# Learning Cross-Lingual IR from an English Retriever

Yulong Li† , Martin Franz‡∗, Md Arafat Sultan‡∗,
Bhavani Iyer‡, Young-Suk Lee‡ Equal contribution. Avirup Sil‡
†IBM Research, ‡IBM Research AI
{yulongl, franzm, bsiyer, ysuklee<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

We present a new cross-lingual information retrieval (CLIR) model trained
using multi-stage knowledge distillation (KD). The teacher and the student are
heterogeneous systems—the former is a pipeline that relies on machine
translation and monolingual IR, while the latter executes a single CLIR
operation. We show that the student can learn both multilingual
representations and CLIR by optimizing two corresponding KD objectives.
Learning multilingual representations from an English-only retriever is
accomplished using a novel cross-lingual alignment algorithm that greedily re-
positions the teacher tokens for alignment. Evaluation on the XOR-TyDi
benchmark shows that the proposed model is far more effective than the
existing approach of fine-tuning with cross-lingual labeled IR data, with a
gain in accuracy of 25.4 Recall@5kt.

##  1 Introduction

Multilingual models are a topic of growing interest in NLP given their
criticality to the universal adoption of AI. Cross-lingual information
retrieval (CLIR) Braschler et al. (1999); Shakery and Zhai (2013); Saleh and
Pecina (2019); Jiang et al. (2020); Asai et al. (2021a); Shi et al. (2021),
for example, can find relevant text in a high-resource language such as
English even when the query is posed in a different, possibly low-resource,
language. In this work, we aim to develop useful CLIR models for this
constrained, yet important, setting where a retrieval corpus is available only
in a single high-resource language (English in our experiments).

A solution to this problem can take one of two general forms. First, machine
translation (MT) of the query into English followed by monolingual (English)
IR Asai et al. (2021a). While such a pipeline approach can provide accurate
predictions, an alternative solution that can tackle the problem purely cross-
lingually, i.e., without involving MT, can be more efficient and cost-
effective. Pre-trained multilingual masked language models (PLMs) such as
multilingual BERT Devlin et al. (2019) and XLM-RoBERTa (XLM-R henceforth)
Conneau et al. (2020) can provide the foundation for such a cross-lingual
solution, as one can simply fine-tune such a model with labeled CLIR data Asai
et al. (2021b).

Here we first run an empirical evaluation of these two approaches on a public
CLIR benchmark Asai et al. (2021a). We use ColBERT111Due to its state-of-the-
art (SOTA) performance outperforming DPR Karpukhin et al. (2020). Khattab and
Zaharia (2020); Khattab et al. (2021) as our IR architecture and XLM-R as the
underlying PLM for both methods (§2). Results indicate that the MT-based
solution can be vastly more effective than CLIR fine-tuning, where we observe
a difference in Recall@5kt of 28.6 (§3). Crucially, the modular design of the
former allows it to leverage additional English-only training data to train
its IR component, which contributes significantly to its performance.

The above results lead naturally to our second research question: Can a more
accurate CLIR model be trained that can operate without having to rely on MT?
To answer the question, instead of viewing the MT-based approach as a
competing one, we propose to leverage its strength via knowledge distillation
(KD) into a single-step CLIR model. KD Hinton et al. (2014) is a powerful
supervision technique typically used to distill the knowledge of a large
teacher model about some task into a smaller student model Sanh et al. (2019).
Here we propose to use KD in a slightly different context, where the teacher
and the student IR models are identical in size, but the former has superior
performance simply due to having access to an MT module and consequently
operating in a high-resource and low-difficulty monolingual environment.

Figure 1: Relevance score distillation. The teacher is shown English triples
but the student’s query input is non-English. Training minimizes the student’s
KL divergence from the teacher’s output softmax distribution over
$S_{q,d^{+}}$ and $S_{q,d^{-}}$ ($\tau$ is the temperature).

We perform two different KD operations (§2.2). The first one directly
optimizes an IR objective, where we use a labeled CLIR dataset consisting of
parallel (English and non-English) questions and their corresponding positive
and negative passages for training. The teacher and the student are shown the
English and non-English versions of a question, respectively, and the training
objective is for the student to match the soft query-passage relevance
predictions of the teacher. Our second KD task is representation learning from
parallel text, where the student learns to encode a non-English text that
matches the teacher’s encoding of the aligned English text at the token level.
The cross-lingual token alignment needed to create the training data for this
task is generated using a novel greedy alignment process which, despite its
noisy nature, proves to be highly effective.

Experimental results (§3) show that our proposed methods can regain most of
the performance loss from the two-stage solution to the purely cross-lingual
model. On an XOR-TyDi test set, the student outperforms the cross-lingual
ColBERT baseline by 25.4 points in terms of Recall@5kt, trailing the teacher
processing English queries by just 3.2 points. Ablation studies also show that
each of our two KD processes contribute significantly towards the final
performance of the student.

Our contributions can be summarized as follows:

  * •

We present an empirical study of the effectiveness of a SOTA IR model
(ColBERT) on cross-lingual IR with and without MT.

  * •

We propose a novel, purely cross-lingual solution that uses knowledge
distillation to learn both improved text representation and IR.

  * •

We demonstrate with a novel cross-lingual alignment algorithm that
distillation using parallel text can strongly augment cross-lingual IR
training.

  * •

We achieve new SOTA results on the XOR-TyDi cross-lingual IR task.

##  2 Method

Here we describe our base IR model (ColBERT) and the proposed KD-based cross-
lingual training algorithms.

###  2.1 The ColBERT Model

The ColBERT Khattab and Zaharia (2020) architecture consists of a shared
transformer-based encoder that separately encodes the input query and
document, followed by a linear compression layer. Each training instance is a
<$q$, $d^{+}$, $d^{-}$> triple, where $q$ is a query, $d^{+}$ is a positive
(relevant) document and $d^{-}$ is a negative (irrelevant) document. ColBERT
first computes a relevance score for the pair $(q,d)$ using Equation 1, where
$d\in\\{d^{+},d^{-}\\}$ and $E_{q_{i}}$ and $E_{d_{j}}$ are the output
embeddings of query token $q_{i}$ and document token $d_{j}$, respectively.
For a given training triple, a pairwise softmax cross-entropy loss is
minimized over the computed scores $S_{q,d^{+}}$ and $S_{q,d^{-}}$.

$S_{q},_{d}:=\sum_{i\in[|q|]}max_{j\in[|d|]}E_{q_{i}}\cdot E^{T}_{d_{j}}$ (1)

For inference, the embeddings of all documents are calculated a priori, while
the query embeddings and the relevance score are computed at runtime.

###  2.2 Knowledge Distillation

Our teacher and student are both ColBERT models that fine-tune the same
underlying multilingual PLM for IR. The teacher is first trained with all-
English triples using the above ColBERT objective. The goal of the subsequent
KD training is to teach the student how to reproduce the behavior of the
teacher when the former is asked a non-English question and the latter its
English translation.

We apply KD at two different stages of the ColBERT workflow: (1) relevance
score computation ($S_{q,d}$ in Equation 1), and (2) encoding (e.g.,
$E_{q_{i}}$). Figure 1 depicts the former in detail, where training minimizes
the KL divergence between the student’s and teacher’s output softmax
distribution (with temperature) over $S_{q,d^{+}}$ and $S_{q,d^{-}}$.

Figure 2: Distillation for representation learning. The student learns to
encode both English and non-English tokens in context that matches the
teacher’s output embeddings for corresponding English tokens.

There is limited availability of labeled training data for CLIR. MT, on the
other hand, is a more established area of research that has produced a large
amount of parallel text. We seek to exploit parallel corpora in our second KD
training stage, where we train the student to compute a representation for a
non-English text that closely matches the teacher’s representation of the
aligned English text. Crucially, since ColBERT computes a single vector for
each individual input token (i.e., a PLM vocabulary item) and not for the
entire input text, our algorithm must support distillation at the token level.

To achieve this, we apply an iterative cross-lingual alignment algorithm.
Assuming $(ne_{1},...,ne_{S})$ to be the ordered tuple of tokens in a non-
English text and $(e_{1},...,e_{T})$ the corresponding tuple from the aligned
English text, each iteration of this algorithm greedily aligns the next
$(ne_{i},e_{j})$ pair with the highest cosine similarity of their output
embeddings. Algorithm 1 implements this idea by repositioning the teacher’s
tokens so that they are position-wise aligned with the corresponding student
tokens. Note that the design choice of using a common multilingual PLM in the
teacher and the student, even though the former is tasked only with handling
English content, is key for the operation of this algorithm as it relies on
the pre-trained PLM’s multilingual representations.

Input:

$v_{T}$: Teacher’s representation (embeddings) of tokenized English (EN) text.

$v_{S}$: Student’s representation of parallel non-EN text.

Output:

$v^{(a)}_{T}$: Reordered teacher output embeddings to reflect position-wise
alignment with $v_{S}$.

Procedure:

$DM\leftarrow cosine\\_distance(v_{T},v_{S})$ //matrix

//get index pairs to swap in $v_{T}$

$swaps\leftarrow[\ ]$

for  _row in rows(DM)_ do

//loop runs $|v_{T}|$ times

$minValue\leftarrow min(DM)$

$i,j\leftarrow index\\_of(minValue)$

//swap rows i and j

$DM[[i,j],:]=DM[[j,i],:]$

//set row and column j to a large value

$DM[j,:]\leftarrow+\infty$

$DM[:,j]\leftarrow+\infty$

$swaps.append((i,j))$

end for

//swap teacher’s output tokens

$v^{(a)}_{T}\leftarrow v_{T}$

for  _s in swaps_ do

$v^{(a)}_{T}[s[0],s[1]]\leftarrow v_{T}[s[1],s[0]]$

end for

Algorithm 1 Cross-lingual alignment.

In addition to cross-lingual alignment, we also perform a similar KD procedure
in which both the teacher and the student are shown the same English text.
This step is useful because ColBERT uses a shared encoder for the query and
the document, necessitating a student that is able to effectively encode text
of both English documents and non-English queries.

Using the alignment information, we train the student by minimizing the
Euclidean distance between its representation of a token (English or non-
English) and the teacher’s representation of the corresponding English token.
Figure 2 shows the KD process for representation learning.

##  3 Experiments and Results

###  3.1 Datasets

Our primary CLIR dataset is XOR-TyDi Asai et al. (2021a), which consists of
examples in seven typologically diverse languages: Arabic, Bengali, Finnish,
Japanese, Korean, Russian and Telugu. The official training set contains
15,221 natural language queries, their short answers, and examples of relevant
(positive) and non-relevant (negative) Wikipedia text snippets with respect to
the queries. For most queries, there are one positive and three negative
examples. We remove the 1,699 (11%) questions that have no answers in the
dataset. In our experiments, we use a random selection of 90% of the remaining
examples for training and the rest for validation.

Following the original XOR-TyDi process, we also obtain additional training
examples by running BM25-based retrieval against a corpus of Wikipedia text
snippets and using answer string match as the relevance criterion. We add
these examples to the original set to obtain three positive and 100 negative
examples per query. For blind testing, we use the 2,113 examples in the
official XOR-TyDi dev set.

System | R@5kt
---|---
$\textnormal{ColBERT}_{cl}$: $\textnormal{FT}_{NQ}$ $\rightarrow$ $\textnormal{FT}_{XOR}$ | 47.7
Teacher: MT + $\textnormal{ColBERT}_{en}$ | 76.3
$\textnormal{ColBERT}_{cl}$ $\rightarrow$ $\textnormal{KD}_{PC}$ $\rightarrow$ $\textnormal{KD}_{XOR}$ | 73.1
Table 1:  Performance on XOR-TyDi. FT: fine-tuning; NQ: the OpenNQ dataset;
PC: parallel corpus; KD: knowledge distillation. Direct fine-tuning of ColBERT
for CLIR ($\textnormal{ColBERT}_{cl}$) lags behind MT + English IR by a large
margin (28.6 points). Our proposed KD methods close the gap by 25.4 points, an
absolute improvement of 89$\%$.

Our monolingual (English) training data are derived from the third fine-tuning
round (ColBERT-QA3) of ColBERT relevance-guided supervision Khattab et al.
(2021) based on OpenNQ data Kwiatkowski et al. (2019). There are about 17.5M
(query, positive, negative) triples in this set. We use an in-house IBM neural
MT system for question translation. Finally, the parallel corpus used in our
KD experiments for representation learning are collected from three sources:
(1) an in-house crawl of Korean, (2) LDC releases (Arabic), and (3)
OPUS222https://opus.nlpl.eu. The corpus has a total of 6.9M passage pairs
which include .9M pairs in Telugu and 1M pairs in each of the other six
languages.

###  3.2 Metric

Our evaluation metric is Recall@5000-tokens (R@5kt) Asai et al. (2021a), which
computes the fraction of questions for which the ground truth short answer is
contained within the top 5,000 tokens of the retrieved passages.

###  3.3 Evaluation on XOR-TyDi

The CLIR baseline for our experiments is a ColBERT model with an underlying
XLM-R PLM, which we first fine-tune with 17.5M NQ examples for one epoch and
then 2.9M XOR-TyDi triples for five epochs. Our student model is initialized
with the parameter weights of the baseline, and is further fine-tuned using
the two KD objectives. The monolingual teacher model—also a ColBERT model on
top of the pre-trained XLM-R—is trained with only the 17.5M NQ triples for one
epoch.

Table 1 compares the performance of our models measured by their R@5kt scores.
The baseline underperforms the MT $\rightarrow$ English IR pipeline by 28.6
points. KD with parallel corpus on the output representations followed by XOR
triples on the query-passage relevance scores achieves an overall improvement
of 25.4 points for the baseline model, which, quite impressively, is only 3.2
points behind the teacher’s score.

###  3.4 Ablation Study

To study the effect of the two different KD operations on our student model’s
performance, we train two additional students. Each of these students goes
through only one of the two KD training steps. Table 2 summarizes the results:
KD with only CLIR examples and only the parallel corpus improves the system’s
score by 15.9 and 20.9 points, respectively. Interestingly, although the
parallel corpus does not contain any IR signal, it contributes more to the
final performance of our model. These results also suggest that our cross-
lingual alignment algorithm does indeed produce useful alignments.

System | R@5kt
---|---
$\textnormal{ColBERT}_{cl}$ $\rightarrow$ $\textnormal{KD}_{PC}$ $\rightarrow$ $\textnormal{KD}_{XOR}$ | 73.1
$\textnormal{ColBERT}_{cl}$ $\rightarrow$ $\textnormal{KD}_{PC}$ | 68.6
$\textnormal{ColBERT}_{cl}$ $\rightarrow$ $\textnormal{KD}_{XOR}$ | 63.6
$\textnormal{ColBERT}_{cl}$ | 47.7
Table 2:  Ablation study results. KD with parallel corpus (PC) and XOR-TyDi IR
triples both play key roles in our student model. Interestingly, the former
has a greater impact on its performance.

##  4 Conclusion

We show that without the help of machine translation (MT) at inference time,
the accuracy of a state-of-the-art IR framework on cross-lingual IR drops
considerably. As a solution, we propose new algorithms to distill the
knowledge of a teacher model that performs monolingual IR on MT output into a
cross-lingual student model capable of operating without MT. We utilize
knowledge distillation (KD) for both IR and text representation, and present
(for the latter) a novel cross-lingual alignment algorithm that only relies on
the underlying masked language model’s multilingual representation
capabilities. In empirical evaluation, our student model recovers most of the
performance drop due to operating in a single-pass cross-lingual mode. Future
work will explore, among other ideas, zero-shot application of our models to
new datasets and utilization of our approach for end-to-end QA.

## References

  * Asai et al. (2021a) Akari Asai, Jungo Kasai, Jonathan Clark, Kenton Lee, Eunsol Choi, and Hannaneh Hajishirzi. 2021a.  XOR QA: Cross-lingual Open-Retrieval Question Answering.  In _NAACL_.
  * Asai et al. (2021b) Akari Asai, Xinyan Yu, Jungo Kasai, and Hannaneh Hajishirzi. 2021b.  One question answering model for many languages with cross-lingual dense passage retrieval.  In _NeurIPS_.
  * Braschler et al. (1999) Martin Braschler, Jürgen Krause, Carol Peters, and Peter Schäuble. 1999\.  Cross-language information retrieval (clir) track overview.  In _TREC_. Citeseer.
  * Conneau et al. (2020) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2020.  Unsupervised Cross-lingual Representation Learning at Scale.  In _ACL_.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  In _NAACL_.
  * Hinton et al. (2014) Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. 2014.  Distilling the knowledge in a neural network.  In _NeurIPS Deep Learning Workshop_.
  * Jiang et al. (2020) Zhuolin Jiang, Amro El-Jaroudi, William Hartmann, Damianos Karakos, and Lingjun Zhao. 2020.  Cross-lingual information retrieval with bert.  _arXiv preprint arXiv:2004.13005_.
  * Karpukhin et al. (2020) Vladimir Karpukhin, Barlas Oguz, Sewon Min, Patrick Lewis, Ledell Wu, Sergey Edunov, Danqi Chen, and Wen-tau Yih. 2020.  Dense passage retrieval for open-domain question answering.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_.
  * Khattab et al. (2021) Omar Khattab, Christopher Potts, and Matei Zaharia. 2021.  Relevance-guided supervision for openqa with colbert.  _Transactions of the ACL_.
  * Khattab and Zaharia (2020) Omar Khattab and Matei Zaharia. 2020.  Colbert: Efficient and effective passage search via contextualized late interaction over bert.  In _SIGIR_.
  * Kwiatkowski et al. (2019) Tom Kwiatkowski, Jennimaria Palomaki, Olivia Redfield, Michael Collins, Ankur Parikh, Chris Alberti, Danielle Epstein, Illia Polosukhin, Matthew Kelcey, Jacob Devlin, Kenton Lee, Kristina N. Toutanova, Llion Jones, Ming-Wei Chang, Andrew Dai, Jakob Uszkoreit, Quoc Le, and Slav Petrov. 2019.  Natural questions: a benchmark for question answering research.  _Transactions of the Association of Computational Linguistics_.
  * Saleh and Pecina (2019) Shadi Saleh and Pavel Pecina. 2019.  An Extended CLEF eHealth Test Collection for Cross-Lingual Information Retrieval in the Medical Domain.  In _European Conference on Information Retrieval_.
  * Sanh et al. (2019) Victor Sanh, Lysandre Debut, Julien Chaumond, and Thomas Wolf. 2019.  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter.  _arXiv preprint arXiv:1910.01108_.
  * Shakery and Zhai (2013) Azadeh Shakery and ChengXiang Zhai. 2013.  Leveraging comparable corpora for cross-lingual information retrieval in resource-lean language pairs.  _Information retrieval_ , 16(1):1–29.
  * Shi et al. (2021) Peng Shi, Rui Zhang, He Bai, and Jimmy Lin. 2021.  Cross-lingual training of dense retrievers for document retrieval.  In _Proceedings of the 1st Workshop on Multilingual Representation Learning_ , pages 251–253.

##  Appendix A Appendix

Following are the hyperparameter combinations used in our different models.
They were selected based on performance on the validation set.

Parameter | Value
---|---
Parameters in ColBERT model |
batch size | 192
accumsteps | 6
linear transfer dim | 128
query_maxlen | 32
doc_maxlen | 180
Parameters in regular training |
lr(NQ) | 1.5e-6
lr(XOR) | 6e-6
Epochs(NQ) | 1
Epochs(XOR) | 5
Parameters in distillation |
loss(XOR/Synthetic data) | KLDiv
loss(Parallel corpus) | MSE
temperature(XOR) | 2
temperature(Synthetic data) | 4
lr(XOR) | 6e-6
lr(Synthetic data) | 1.5e-6
lr(Parallel corpus) | 4.8e-5
Epochs(XOR) | 5
Epochs(Synthetic data) | 1
Epoches(Parallel corpus) | 2
Table 3:  Hyper-parameters used in our test set runs.


We present a framework for the synthesis of phase polynomials that addresses both cases of full connectivity and partial connectivity for NISQ architectures.
In most cases, our algorithms generate circuits with lower CNOT count and CNOT depth than the state of the art or have a significantly smaller running time for similar performances.
We also provide methods that can be applied to our algorithms in order to trade an increase in the CNOT count for a decrease in execution time, thereby filling the gap between our algorithms and faster ones.

§ INTRODUCTION

Quantum circuits optimization is essential to foster the practicability and efficiency of quantum computation.
In particular, to cope with the much-needed compactness of quantum circuits, the synthesis of reversible circuits is being studied thoroughly.
Because the $T$ gate has a high fault-tolerant implementation cost <cit.>, much work has been put into the minimization of the $T$-count <cit.> and the $T$-depth <cit.>.
In contrast, the CNOT gate has a low implementation cost as it is part of the Clifford group <cit.>.
Nonetheless, the usage of metrics based on the $T$ gate have limitations, it turns out that the number of CNOT gates in a circuit is a metric that should not be overlooked as it can have a significant impact on the implementation cost of a circuit <cit.>.

On top of that, quantum computers in the Noisy Intermediate Scale Quantum (NISQ) era <cit.> have architectural constraints.
Concretely, the qubits within these computers are not connected in an all-to-all manner.
It implies that logical gates having an arity of $2$, such as CNOT gates, can only be applied between certain pairs of qubits.
Thus, making a circuit compliant with a given architecture inevitably causes an increase in the CNOT count <cit.>.

A common way of dealing with architectural constraints is to insert SWAP gates to route logical qubits <cit.>.
An alternative is to perform architecture-aware synthesis <cit.>, a method which often produces circuits with a much lower CNOT count while satisfying the architectural constraints.
This approach is typically applied on subsets of circuits that can be represented by high-level constructs such as linear reversible functions.
These circuits can then be put together to form a complete architecture compliant quantum circuit <cit.>.
An important building block in this compilation scheme is the synthesis of circuits composed exclusively of CNOT and $R_Z$ gates.
These circuits can be represented by a high-level construct called phase polynomials.
In this work we tackle the phase polynomials synthesis problem and propose efficient algorithms for both cases of restricted and complete connectivity.

State of the art.
In <cit.>, a SAT-based algorithm that optimally solves the phase polynomials synthesis problem is proposed.
Although this method offers good results regarding the CNOT count, it has an exponential complexity since the SAT problem is NP-complete and is therefore only practical for the synthesis of small phase polynomials.
An efficient heuristic algorithm for phase polynomials synthesis is provided by Amy et al. in <cit.>.
This algorithm, named Gray-Synth, is inspired by Gray code <cit.> and is considered as the current state of the art.
There exists numerous other algorithms for phase polynomials synthesis.
Some of them don't have CNOT minimization as primary objective, as it is the case of the $T$par algorithm <cit.> that aims to parallelize the phase gates of a phase polynomial.
This is also the case, with a lower degree, in <cit.> where automated methods for the optimization of large quantum circuit are given.
As their algorithm scales similarly to the Gray-Synth algorithm and is purposely designed for specific circuits, we will prefer to compare our algorithm with the Gray-Synth algorithm.

Regarding the complexity of the problem, the results presented in <cit.> lead us to think that it is intractable.
The same authors of the Gray-Synth algorithm corroborate this idea by proving the NP-completeness of the problem in some restricted cases <cit.>.

Qubit routing could be used to make the Gray-Synth algorithm compliant with constrained architectures.
This idea was developed and greatly refined by Nash et al. <cit.>, and a modified version of their algorithm has been implemented in the Staq toolkit <cit.>.
An altered version of this algorithm was also recently incorporated by Gheorghiu et al. in a slice-and-build algorithm that optimizes a given quantum circuit while taking into account the connectivity constraints imposed by the physical hardware architecture <cit.>.
In a recent work and with a similar goal, a framework composed of greedy architecture-aware synthesis routines for the compilation of quantum circuits was presented in <cit.>.

A phase polynomial is partially composed of a set of parities which can be stored in a parity table.
As explained in <cit.>, a circuit in which each parity occurs at least once is called a parity network and can be easily modified in order to implement the corresponding phase polynomials.
In all parity network synthesis algorithms the parities are synthesized in an established order, we refer to this order as the parity ordering.
In the Gray-Synth algorithm <cit.> this ordering is inspired by Gray code.
Most of the parity network synthesis algorithms for arbitrary connectivity follow this idea and are based on the Gray-Synth algorithm.
Yet, while the parity ordering defined by the Gray-Synth algorithm is efficient for all-to-all connectivity, it may be unfitted for arbitrary connectivity.
Indeed, most Gray-Synth based algorithms for arbitrary connectivity are not taking the architecture into account when establishing the parity ordering <cit.>.
In other words, the choice of the next parity to synthesize is solely based on the parity table, without taking into account the underlying graph of the architecture.
An algorithm proposed by Arianne Meijer-van de Griend and Ross Duncan <cit.> aims to solve this shortcoming by recursively considering only non-cutting vertices of the underlying graph.
However this algorithm is still based on the parity ordering of the Gray-Synth algorithm which is foremostly designed for all-to-all connectivity.

Our approach.
In this paper we present an efficient alternative to the Gray code inspired parity ordering that helds better results in both cases of all-to-all connectivity and constrained architectures.
In our approach, the parity ordering is defined by a two steps iterative process.
The first step consists in choosing a parity and the second one corresponds to the synthesis of the chosen parity, we iterate until all parities have been synthesized.
Here the parity choice is not bound to a parity ordering uniquely defined upon the parity table as it is the case for Gray code inspired methods, but can also take into account an arbitrary connectivity.
In fact, we will see that this method can be easily adapted to constrained architectures by relying on the commonly used notion of Steiner tree.
This extension to constrained architectures induces an important time cost, nevertheless we will present some techniques to significantly reduce the running time of our algorithm while preserving an important CNOT count reduction when compared to the state of the art.

This paper is organized as follows.
Section <ref> introduces the circuit-polynomial correspondence for quantum circuits over the $\{CNOT, R_z\}$ gate set.
In Section <ref> we present our heuristic algorithm for the synthesis of phase polynomials in the case of full connectivity.
In Section <ref> we extend our algorithm for partial connectivity and we give methods to lower the complexity of our approach.
Benchmarks are given at the end of Sections <ref> and <ref>.

§ PHASE POLYNOMIALS SYNTHESIS

Let $C$ be a quantum circuit operating over $n$ qubits and composed of CNOT and $R_z$ gates.
Such circuit can be best described by exploiting the circuit-polynomial correspondence <cit.>, which associates a phase polynomial and a linear reversible function to $C$.
The action of $C$ on a basis state has the form

$$\lvert \bs x\rangle \mapsto e^{2\pi ip(\bs x)} \lvert g(\bf \bs x)\rangle$$

where $g: \mathbb{F}_2^n \rightarrow \mathbb{F}_2^n$ is a linear reversible function and

$$ p(\bs x) = \sum_{i = 1}^{2^n} \theta_i\  f_{i}(\bs x)$$

is a linear combination of linear Boolean functions $f_{i}: \mathbb{F}_2^n \rightarrow \mathbb{F}_2$.
Any linear Boolean function $f_{i}$ can be written as

$$f_{i}(\bs x) = \bs y^i \cdot \bs x= y^i_1x_1  \oplus y^i_2 x_2 \oplus \ldots y^i_n x_n$$

where $\bs y^i \in \mathbb{F}_2^n$ and $\oplus$ stands for the XOR operation.
The function $p(\bs x)$ is the phase polynomial associated with $C$, and we will refer to the Boolean vectors $\bs y^i$ as the parities of the phase polynomial $p(\bs x)$.
For instance, the circuit represented Figure <ref> performs the mapping

$$ \lvert x_1, x_2, x_3, x_4 \rangle \mapsto e^{ip(x_1, x_2, x_3, x_4)} \lvert x_1 \oplus x_2 \oplus x_3, \; x_1 \oplus x_3 \oplus x_4,\; x_3,\; x_4 \rangle$$

where $p(x_1, x_2, x_3, x_4) = \theta_1(x_1 \oplus x_2) + \theta_2(x_1 \oplus x_2 \oplus x_3) + \theta_3(x_1 \oplus x_3 \oplus x_4)$.
The parities of a phase polynomial can be described by a matrix where each line represents a qubit and each column represents a parity having an associated angle not equal to $0$, we call this matrix the parity table of the phase polynomial and we denote it $P$.
In our example, the parity table of the phase polynomial is

$$P =
\label{matrix}
\begin{pmatrix}
1 & 1 & 1\\
1 & 1 & 0\\
0 & 1 & 1\\
0 & 0 & 1\\
\end{pmatrix}.

Performing the synthesis of the phase polynomial $p(\bs x)$ and the linear reversible function $g(\bs x)$ amounts to constructing a circuit equivalent to $C$.
The synthesis of linear reversible functions is a well studied problem as there exists asymptotically optimal methods <cit.>, as well as efficient heuristic algorithms in both cases of partial and full connectivity <cit.>.
For that reason we will put aside the problem of synthesizing the linear reversible function $g(\bs x)$, and we will focus on the phase polynomials synthesis problem.

@C=1.6em @R=1.4em
x_1          x_1⊕x_2     1.0cm$R_z(\theta_1)$        x_1⊕x_2 ⊕x_3        1.0cm$R_z(\theta_2)$     1                             x_1 ⊕x_2 ⊕x_3

x_2     -1                                 x_1⊕x_3        x_1⊕x_3 ⊕x_4        1.0cm$R_z(\theta_3)$         x_1 ⊕x_3 ⊕x_4

x_3                 -2                                             x_3

x_4                                         -2                     x_4

An annotated circuit composed of CNOT and $R_z$ gates.

Parity networks.
Following Amy et al. <cit.>, we tackle the phase polynomials synthesis problem via the parity network formalism.

A parity network for a parity table $P$ is a CNOT circuit in which each parity $\bs y \in P$ appears at least once.

The $R_z$ gate only modifies the phase and doesn't affect the logical value of the qubit.
Therefore, a parity network can be easily modified by exclusively inserting $R_z$ gates to implement any phase polynomial associated with the parity table $P$.
It implies that the phase polynomials synthesis problem can be reduced to the parity network synthesis problem.
For the remaining of this paper we will consider the parity network synthesis problem and ignore the rotation gates in order to focus on the core of the problem.

Recall that a CNOT gate performs the mapping $\lvert x_i \rangle\lvert x_j \rangle \mapsto \lvert x_i \rangle\lvert x_i \oplus x_j \rangle$.
When applying a CNOT$_{x_i, x_j}$ gate where $x_i$ is the control qubit and $x_j$ is the target qubit, the parity table can be expressed in the new basis by performing the row addition $P_i = P_i \oplus P_j$.
It follows that a parity $\bs y \in P$ is being carried out by a qubit if it satisfies $\sum_{i=1}^n y_i = 1$.

§ PARITY NETWORK SYNTHESIS FOR ALL-TO-ALL CONNECTIVITY

§.§ An efficient heuristic algorithm

In this section we formalize a heuristic algorithm for phase polynomials synthesis in all-to-all connectivity.
Our algorithm is presented in pseudo-code in Algorithm <ref> and an example is provided in Figure <ref>.
The term CNOT$_{i,j}$ refers to a CNOT gate with control $i$ and target $j$ and we define the function $h$ as the Hamming weight of a binary vector or binary matrix, i.e. $h(\bs y) = \sum_{i=1}^n y_i$ where $\bs y \in \mathbb{F}_2^n$ and $h(P) = \sum_{\bs y \in P} h(\bs y)$ where $P$ is a parity table.
Let $S$ be a sequence of row additions, we denote by $\bs y^S$ (resp. $P^S$) the state of $\bs y$ (resp. $P$) after applying the sequence of additions in $S$ onto $\bs y$ (resp. $P$).
Our algorithm follows a two steps iterative process:

* Choose a parity $\bs y \in P$.

* Perform the synthesis of $\bs y$. Remove $\bs y$ from $P$ and go to step $1$.

Step 1.
To know which parity $\bs y \in P$ it would be judicious to choose in step $1$ of our algorithm we first reflect on the minimum CNOT cost induced by the synthesis of $\bs y$.
Let $S_{\bs y}$ be the set of minimum length sequences of additions such that for all $S \in S_{\bs y}$ we have $h(\bs y^S) = 1$.
In other words, $S_{\bs y}$ contains all minimum length sequences of additions that effectively synthesize $\bs y$.
Note that one addition can reduce the value of $h(\bs y)$ by at most $1$.
That being so, it is clear that the length of all $S \in S_{\bs y}$ is equal to $h(\bs y) - 1$, implying a minimum CNOT cost of $h(\bs y) - 1$ for the synthesis of $\bs y$.
Based on this fact it is rather intuitive to choose the parity $\bs y$ for which $h(\bs y)$ is minimum for the step $1$ of our algorithm.

Step 2.
We are left with the second step of the algorithm that raises the following question: which sequence $S \in S_{\bs y}$ of additions should we choose to perform the synthesis of $\bs y$?
A natural choice would be to select the sequence $S \in S_{\bs y}$ that minimizes the value $h(P^S)$.
However we are faced with a challenge as the size of $S_{\bs y}$ is exponential with respect to $h(\bs y)$.
In fact, we will see that the size of $S_{\bs y}$ is greater than the number of spanning arborescences in a complete directed graph composed of $h(\bs y)$ vertices.
lWe define this graph as follows.

Let $P$ be a parity table and $\bs y$ be a parity of $P$, the parity graph associated with $\bs y$ is the complete directed graph $G_{\bs y}=(V, A)$ where $V=\{i \mid y_i = 1\}$ and where each arc $(i, j) \in A$ going from $i$ to $j$ is weighted by $w_{i,j} = h(P_i \oplus P_j) - h(P_j)$.

Let $X$ be a spanning arborescence of $G_{\bs y}$, note that $X$ is composed of $h(\bs y) - 1$ arcs.
For each arc $(i, j) \in X$ going from $i$ to $j$ we associate the addition $P_j = P_j \oplus P_i$.
Now consider a successors-first traversal of $X$ as defined below.

A traversal of an arborescence $X$ is a successors-first traversal if and only if for every vertex $i$ in $X$ the successors of $i$ in $X$ are visited before $i$.

For example, the traversal resulting from a depth-first search postordering is a successors-first traversal.
We can construct a sequence $S$ of additions by following the order of this traversal: if $j$ is the currently visited vertex then we append to $S$ the addition associated with the unique arc $(i, j) \in X$ going from $i$ to $j$.
The length of $S$ is equal to $h(\bs y) - 1$ and $h(\bs y^S) = 1$, thus $S \in S_{\bs y}$ and we say that $X$ is the spanning arborescence associated with $S$.
A different successors-first traversal would give us a different sequence of additions for the same spanning arborescence $X$, and for every sequence $S \in S_{\bs y}$ there is a corresponding spanning arborescence in $G_{\bs y}$.
Hence there is a surjection between $S_{\bs y}$ and the set of spanning arborescences in $G_{\bs y}$, in particular we have $|S_{\bs y}| \geq |\{X \mid X\text{ is a spanning arborescence of }G_{\bs y}\}| = n^{n-1}$ where $n=h(\bs y)$.
As a result we have to choose $S\in S_{\bs y}$ among an exponential number of possibilities.

To cope with this problem we can first notice that if two sequences $S_i, S_j$ are associated with the same arborescence, then they are equivalent in the sense that $P^{S_i} = P^{S_j}$.
As our only metric is the number of CNOT we don't make any distinction between these equivalent sequences, and we can equivalently refer to $S \in S_{\bs y}$ or its associated spanning arborescence in $G_{\bs y}$.
Recall that each arc $(i, j) \in A$ going from $i$ to $j$ is weighted by $w_{i,j} = h(P_i \oplus P_j) - h(P_j)$, then we have
$$h(P^S) = h(P) + \sum_{(i, j) \in X} h(P_i \oplus P_j) - h(P_j) = h(P) + \sum_{(i, j) \in X} w_{i,j} $$
where $S \in S_{\bs y}$ and $X$ is the spanning arborescence associated with $S$ in $G_{\bs y}$.
In consequence, minimizing $h(P^S)$ amounts to minimizing $\sum_{(i, j) \in X} w_{i, j}$, which is optimally satisfied when $X$ is the minimum weight spanning arborescence of $G_{\bs y}$.
The minimum weight spanning arborescence problem is a well-known problem, an algorithm proposed by Robert Endre Tarjan <cit.> solves it with a complexity of $\mathcal{O}(|V|^2)$ for complete graphs where $|V|$ is the number of vertices.

Correctness and complexity.
Let $P$ be a parity table of size $n \times m$.
Our algorithm terminates when $P$ is empty and a parity $\bs y \in P$ is removed from $P$ at each iteration.
Therefore our algorithm performs $m$ iterations and finishes.
The algorithm starts with an empty circuit $C$ and at each iteration $C$ is extended to synthesize a parity $\bs y \in P$ not yet occurring in $C$.
Hence the constructed circuit $C$ is a parity network for $P$ and our algorithm is correct.

Choosing the parity in the step $1$ of our algorithm has a cost of $\mathcal{O}(mn)$.
For the step $2$, constructing the graph $G_{\bs y}$ has a complexity of $\mathcal{O}(mn^2)$, whereas computing its minimum weight spanning arborescence and performing the for loop over the depth-first search both have a smaller complexity of $\mathcal{O}(n^2)$ and $\mathcal{O}(mn)$ respectively.
Both steps are performed $m$ times so the overall complexity of our algorithm is $\mathcal{O}(m^2n^2)$.

Parity table $P$
Circuit synthesizing a parity network associated with $P$
$C \leftarrow$ new empty circuit

$P$ non-empty
$\displaystyle \bs y \leftarrow \min \operatorname*{argmin}_{\bs y} \{h(\bs y) \mid \bs y \in P\}$

$P \leftarrow P \setminus \bs y$

$G_{\bs y} \leftarrow (\{i \mid y_i = 1\}, \{(i, j, h(P_i \oplus P_j) - h(P_j)) \mid y_i=1, y_j=1, i \neq j\})$

$X \leftarrow$ MinimumWeightSpanningArborescence$(G_{\bs y})$

$i \in\ $DepthFirstSearchPostordering$(X)$
$j \leftarrow $ direct predecessor of $i$ in $X$

$C \leftarrow C :: CNOT_{i, j}$

$P_i \leftarrow P_i \oplus P_j$

Parity network synthesis

\[
\begin{tikzpicture} \matrix (m)[ matrix of math nodes, left delimiter  = (, right delimiter = ), ] {
\bf 0 & 0 & 1 & 1 & 1 & 1 & 1\\
\bf 1 & 1 & 0 & 0 & 1 & 1 & 1\\
\bf 1 & 1 & 0 & 1 & 0 & 0 & 1\\
\bf 0 & 1 & 1 & 1 & 0 & 1 & 1\\
\bf 1 & 0 & 1 & 1 & 1 & 1 & 1\\
} ;
\end{tikzpicture}
\]

\[
\begin{tikzpicture}[>=Stealth,->,line width=.7pt] \matrix [matrix of math nodes, column sep={2.0cm,between origins},
row sep={2.5cm,between origins}, nodes={circle, draw, minimum size=1cm}] {
& |(2)| 2 &   & |(3)| 3 & \\
&  & |(5)| 5 & \\
\draw[transform canvas={xshift=2pt,yshift=-2pt},shorten >= 2pt,shorten <= -2pt] (2) to node[transform canvas={xshift=-2pt}, below]{$-1$} (3);
\draw[transform canvas={xshift=-2pt,yshift=2pt},shorten >= 2pt,shorten <= -2pt, line width=0.55mm] (3) to node[transform canvas={xshift=2pt}, above]{$-2$} (2);
\draw[transform canvas={xshift=2pt,yshift=-2pt},shorten >= -1pt] (5) to node[transform canvas={xshift=-2pt, yshift=2pt}, below right]{$0$}  (3);
\draw[transform canvas={xshift=-2pt,yshift=2pt},shorten >= -1pt] (3) to node[transform canvas={xshift=2pt, yshift=-2pt}, above left]{$-2$} (5);
\draw[transform canvas={xshift=2pt,yshift=2pt},shorten <= -1pt] (5) to node[transform canvas={xshift=-2pt, yshift=-2pt}, above right]{$-2$} (2);
\draw[transform canvas={xshift=-2pt,yshift=-2pt},shorten <= -1pt, line width=0.55mm] (2) to node[transform canvas={xshift=2pt, yshift=2pt}, below left]{$-3$} (5);
\end{tikzpicture}
\]

Parity table and graph $G_{\bs y}$. The minimum weight spanning arborescence of $G_{\bs y}$ is represented by bold arcs.

\[
\begin{tikzpicture} \matrix (m)[ matrix of math nodes, left delimiter  = (, right delimiter = ), ] {
\bf 0 & 0 & 1 & 1 & 1 & 1 & 1\\
\bf 0 & 0 & 0 & 1 & 1 & 1 & 0\\
\bf 1 & 1 & 0 & 1 & 0 & 0 & 1\\
\bf 0 & 1 & 1 & 1 & 0 & 1 & 1\\
\bf 0 & 1 & 1 & 1 & 0 & 0 & 0\\
} ;
\end{tikzpicture}
\]

\[
\begin{array}{c}
\Qcircuit @C=1.8em @R=1.4em {
\lstick{x_1} & \qw & \qw & \qw & \rstick{\!\!\!\!\!\!\!\!\!\!x_1}\\
\lstick{x_2} & \targ  & \ctrl{1} & \qw & \rstick{\!\!\!\!\!\!\!\!\!\!x_2 \oplus x_5}\\
\lstick{x_3} & \qw & \targ & \qw  & \rstick{\!\!\!\!\!\!\!\!\!\!x_2 \oplus x_3 \oplus x_4}\\
\lstick{x_4} & \qw & \qw & \qw & \rstick{\!\!\!\!\!\!\!\!\!\!x_4}\\
\lstick{x_5} & \ctrl{-3} & \qw & \qw  & \rstick{\!\!\!\!\!\!\!\!\!\!x_5}\\
& & & &\\
\end{array}
\]

State of the parity table after the synthesis of $\bs y$, the additions performed are $P_5 = P_5 \oplus P_2$ and $P_2 = P_2 \oplus P_3$. And circuit corresponding to the synthesis of $\bs y$ with respect to the minimum weight spanning arborescence of $G_{\bs y}$.

Example for 1 iteration of Algorithm <ref>. The chosen parity $\bs y$ is represented in bold.

Further optimizations.
The parity $\bs y$ selected in the step $1$ of our algorithm satisfies $\bs y = \operatorname*{argmin} \{h(\bs y) \mid \bs y \in P\}$.
Yet it often arises that multiple parities satisfy this property.
A second selection criteria is needed in order to choose among these parities.
After experimenting with different methods, we found that what gives the best results is to simply choose the parity representing the smallest integer.
This way the parities will be ordered according to their most significant bit as it is done in <cit.>, and the parities will be processed in an order analogous to the Gray code.
As a result and in the same manner as the Gray-Synth algorithm, our algorithm is optimal when the given parity table contains all possible parities.

As previously mentioned, if two sequences $S_i, S_j$ are associated with the same spanning arborescence $X$, then they are equivalent in the sense that $P^{S_i} = P^{S_j}$.
Nonetheless, the circuits produced by $S_i$ and $S_j$ can have different depths.
Thus, by choosing wisely among these equivalent sequences we could enhance the depth performances of our algorithm without affecting the CNOT count of the outputed circuit.
This opportunity is interesting as depth is often considered as a second metric in CNOT circuits synthesis.
Despite that, we decided not to implement such optimization as finding the depth-optimal traversal of $X$ induces a computational overhead.
As a side note, we want to mention that one could also prioritize depth minimization over CNOT count minimization by finding the sequence of additions $S$ that minimizes $h(P^S)$ among all the sequences in $S_{\bs y}$ with a corresponding CNOT circuit of depth $\lceil \log_2 h(\bs y) \rceil$.

§.§ Benchmarks

For a parity table $P$ of size $n \times m$ we call density the value $100\times\frac{m}{2^n - 1}$, which is the percentage representation of the ratio between the number of parities in $P$ and the number of possible parities for $n$ qubits.
To evaluate the perfomances of our algorithm we generate random parity tables for 7, 10, 13 and 16 qubits with a density varying from $1\%$ to $100\%$.
We compare our results to the state of the art, namely the Gray-Synth algorithm <cit.> and the Lazy-Synth framework <cit.>.
Throughout this paper, we have always set the Lazy-Synth depth parameter to $3$ for its recursive search.
Our algorithm as well as the Gray-Synth algorithm have been implemented in Python whereas the Lazy-Synth algorithm has been implemented in C++.
The standard deviation $\sigma$ is not represented in our benchmarks as it is particularly low ($\sigma < 10^{-2}$).

We first discuss the CNOT count and depth performances of our algorithm in comparison to the Gray-Synth algorithm as presented in Figure <ref>.
As expected, both algorithms are optimal for full density parity tables.
That being said, our algorithm has a better CNOT depth than the Gray-Synth algorithm for all other densities and it also has a better CNOT count for all parity tables with a density between $1\%$ and $95\%$.
Moreover, we observe that the percentage of CNOT gained over the Gray-Synth algorithm increases as the number of qubits gets larger.
It is therefore reasonable to project that our algorithm outperforms the Gray-Synth algorithm for all phase polynomials acting on larger number of qubits and whose parity table has a density between $1\%$ and $95\%$.

As it can be seen in Figure <ref>, for $10$ and $13$ qubits the Lazy-Synth algorithm is solving the parity network synthesis problem with less CNOT than our algorithm for parity tables having a density under $50\%$.
However, this performance gap never exceeds $6\%$ and comes with a tremendous time cost.
The runtimes for the $3$ algorithms are shown in Table <ref>, we weren't able to execute the Lazy-Synth algorithm for $16$ qubits due to its important computational time.
As expected by the complexity analysis, our algorithm is slower than the Gray-Synth algorithm but still has a decent enough running time to be applied on large parity tables.
We provide methods to further reduce the running time of our algorithm while preserving good performances in Section <ref>.
Also, the complexity of our algorithm is majored by the creation of the graph $G_{\bs y}$, yet this task can be easily parallelized to significantly reduce the computational time of the algorithm.

We also performed benchmarks on much lower densities in order to test our algorithm on a higher number of qubits, the results are depicted in Figure <ref>.
We can see that our algorithm still offers a significant CNOT reduction when compared to the Gray-Synth algorithm for higher number of qubits.
Interestingly, the CNOT count ratio is increasing when the number of qubits increases while the CNOT depth ratio does the opposite.
The computational times for this setting are shown in Table <ref>.
In a similar way as in Table <ref>, we can see that the computational time of our algorithm grows faster than the computational time of the Gray-Synth algorithm as the size of the parity table increases.

Overall, with regards to the CNOT metrics, our algorithm offers much better results than the Gray-Synth algorithm and its performances are similar to the Lazy-Synth algorithm but with a viable execution time for large phase polynomials, thus achieving an efficient performance over time ratio.

CNOT count and depth of the circuits generated by our algorithm divided by the CNOT count and depth of the circuits outputed by the Gray-Synth algorithm.
Each point is averaged over $1000, 100, 10$ and $10$ randomly generated parity tables for $7, 10, 13$ and $16$ qubits respectively.

CNOT count and depth of the circuits generated by our algorithm divided by the CNOT count and depth of the circuits outputed by the Lazy-Synth algorithm.
Each point is averaged over $1000, 100, 10$ and $10$ randomly generated parity tables for $7, 10, 13$ and $16$ qubits respectively.

Average computational time in seconds of the Gray-Synth algorithm, the proposed algorithm and the Lazy-Synth algorithm.

2*Density     3c|7 qubits     3c|10 qubits     3c|13 qubits     2c|16 qubits

Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed

1%     -     -     -    0.006     0.011     0.045    0.058     0.094     0.798    0.583     0.839

20%    0.007     0.018     0.056    0.079     0.156     1.393    0.775     1.606     94.539    8.909     37.710

40%    0.011     0.032     0.113    0.111     0.303     4.283    1.167     3.708     329.136    19.042     119.998

60%    0.015     0.048     0.185    0.146     0.447     8.729    1.622     6.541     688.243    30.944     253.264

80%    0.017     0.063     0.279    0.176     0.591     14.985    2.030     10.121     1195.277    43.895     431.766

100%    0.020     0.078     0.380    0.201     0.751     22.193    2.402     14.143     1748.428    56.923     652.193

CNOT count and depth of the circuits generated by our algorithm divided by the CNOT count and depth of the circuits outputed by the Gray-Synth algorithm. The number of parities ranges from $1\%$ of $n^2$ to $n^2$ where $n$ is the number of qubits.
Each point is averaged over $100$ randomly generated parity tables.

Average computational time in seconds of the Gray-Synth algorithm and the proposed algorithm for $n^2$ parities where $n$ is the number of qubits.

2*Algorithm     5c|Number of qubits

20     30     40     50     60

Gray-Synth     0.749     3.696     13.406     30.7     71.3

Proposed     0.815     4.924     22.598     81.2     243

§ PARITY NETWORK SYNTHESIS FOR PARTIAL CONNECTIVITY

Many synthesis algorithms are first designed for all-to-all connectivity before being adapted to the case of partial connectivity.
We follow the same process in this work and now demonstrate how our algorithm can be extended to perform architecture-aware synthesis.

§.§ Extending the heuristic function

As we have seen in Section <ref>, in the case of all-to-all connectivity the cost function $\mathcal{C}: \mathbb{F}_2^n \rightarrow \mathbb{N}$ that gives the minimum CNOT cost for the synthesis of a parity $\bs y$ is equal to $h(\bs y)-1$.
This is utterly incorrect in the case of partial connectivity as the function $h$ doesn't even take into account the connectivity graph of the architecture.
Hence, to extend our algorithm for constrained architectures we first need to determine the value of our heuristic function $\mathcal{C}$ for the case of partial connectivity.
For this purpose we will rely on Steiner trees, which are commonly used when it comes to designing architecture-aware synthesis algorithms <cit.>.

Steiner trees.
Given a graph $G=(V, E)$ and a set of vertices $S \subseteq V$, the Steiner tree problem consists of finding the minimum tree $T=(V_T, E_T)$, called Steiner tree, such that $T$ is a subgraph of $G$ and $S \subseteq V_T$.
The vertices in $S$ are called terminals and the vertices in $V_T \setminus S$ are called Steiner nodes.
To perform the synthesis of a parity $\bs y$, a Steiner tree is required since we must interconnect the vertices in $S=\{i \mid y_i = 1\}$.
Once the Steiner tree $T$ is computed, two steps must be carried out:

* Fill-in all Steiner nodes such that for all $u \in V_T$ we have $y_u = 1$.

* Perform the sequence of additions associated with a successors-first traversal of any spanning arborescence of $T$.

An example of this process is given Figure <ref>, and the pseudo-code to perform the fill-in step is provided in Algorithm <ref>.
We use the notation $G[X]$ to refer to the induced subgraph of $G$ formed by the subset of vertices $X$.
For the second step it is exactly the same procedure as the for loop in Algorithm <ref> for an arbitrary spanning arborescence of $T$.
The fill-in step requires $|V_T \setminus S|$ additions and the second step takes $|V_T| - 1$ additions, hence the total CNOT cost for the synthesis of $\bs y$ is $\mathcal{C}(\bs y) = |V_T \setminus S| + |V_T| - 1 = 2|V_T| - |S| - 1$.
Solving optimally the Steiner tree problem would minimize $C(\bs y)$, unfortunately finding the minimal Steiner tree is NP-hard <cit.>.
We will therefore rely on an approximation algorithm proposed by Takashashi et al. <cit.> that has an approximation ratio of $2 - 2/|S|$.
The algorithm is presented in pseudo-code in Algorithm <ref>, it starts by adding an arbitrary vertex of $S$ to the Steiner tree $T$ and then constructs $T$ by iteratively adding the shortest path between $T$ and one vertex of $S$ not yet in $T$.
Its runtime is $\mathcal{O}(|S||V|^2)$, although it can be lowered to $\mathcal{O}(|V|^2)$ if all the shortest paths are provided.

$T$, $S$
$C \leftarrow$ new empty circuit

$F \leftarrow T[V_T \setminus S]$

$F$ non-empty
$u \leftarrow$ a leaf of $F$

$v \leftarrow $ any vertex in $S \cap T.neighbors(u)$

$C \leftarrow C :: CNOT_{u, v}$

$S \leftarrow S \cup \{u\}$

$F \leftarrow F[V_{F} \setminus \{u\}]$

Arbitrary fill-in

$S$, $paths$
$u \leftarrow $ any vertex in $S$

$S \leftarrow S \setminus \{u\}$

$T \leftarrow $ Graph($\{u\}, \emptyset$)

$S$ non-empty
$\displaystyle u, v \leftarrow \operatorname*{argmin}_{u,v} \{|paths_{u,v}| \mid u \in T, v \in S\}$

$T \leftarrow T \cup paths_{u, v}$

$S \leftarrow S \setminus \{v\}$

Steiner tree

\[
\begin{tikzpicture}[scale=1.22, anchor=west]
\tikzstyle{w}=[circle,draw=black!80, thick, font=\normalsize, minimum size=0.7cm, inner sep=0pt]
\tikzstyle{b}=[circle,fill=gray,draw=black!80, thick, font=\normalsize, minimum size=0.7cm, inner sep=0pt]

\node[b] (v1) at (0,0) {1};
\node[w] (v2) at (1.25,0) {2};
\node[w] (v3) at (2.5,0) {3};
\node[w] (v4) at (0,-1.25) {4};
\node[w] (v5) at (1.25,-1.25) {5};
\node[b] (v6) at (2.5,-1.25) {6};
\node[w] (v7) at (0,-2.5) {7};
\node[b] (v8) at (1.25,-2.5) {8};
\node[w] (v9) at (2.5,-2.5) {9};
\node[b] (v10) at (0,-3.75) {10};
\node[w] (v11) at (1.25,-3.75) {11};
\node[w] (v12) at (2.5,-3.75) {12};

\draw[very thick] (v1) -- (v4);
\draw[very thick] (v4) -- (v5);
\draw[very thick] (v5) -- (v6);
\draw[very thick] (v5) -- (v8);
\draw[very thick] (v8) -- (v11);
\draw[very thick] (v10) -- (v11);
\foreach \from/\to in {v1/v2, v2/v3, v4/v5, v5/v6, v7/v8, v8/v9, v2/v5, v3/v6, v4/v7, v5/v8, v6/v9, v7/v10, v11/v12, v12/v9}
\draw (\from) -- (\to);
\end{tikzpicture}
\]

\[
\begin{array}{c}
\Qcircuit @C=1.0em @R=1.18em {
\lstick{x_1}    & \targ & \qw & \qw & \qw & \qw & \qw & \qw & \qw & \targ & \qw \\%& \rstick{x_1 \oplus x_6 \oplus x_8 \oplus x_{10}}\\
\lstick{x_4}    & \ctrl{-1} & \targ & \qw & \qw & \qw & \qw & \qw & \targ & \ctrl{-1} & \qw \\%& \rstick{x_4 \oplus x_6 \oplus x_8 \oplus x_{10}}\\
\lstick{x_5}    & \qw & \ctrl{-1} & \qw & \qw & \qw & \targ & \targ & \ctrl{-1} & \qw & \qw \\%& \rstick{x_5 \oplus x_{6} \oplus x_8 \oplus x_{10}}\\
\lstick{x_6}    & \qw & \qw & \qw & \qw & \qw & \qw & \ctrl{-1} & \qw & \qw & \qw \\%& \rstick{x_6}\\
\lstick{x_8}    & \qw & \qw & \targ & \qw & \targ & \ctrl{-2} & \qw & \qw & \qw & \qw \\%& \rstick{x_8 \oplus x_{10}}\\
\lstick{x_{10}} & \qw & \qw & \qw & \ctrl{1} & \qw & \qw & \qw & \qw & \qw & \qw \\%& \rstick{x_{10}}\\
\lstick{x_{11}} & \qw & \qw & \ctrl{-2} & \targ & \ctrl{-2} & \qw & \qw & \qw & \qw & \qw \\%& \rstick{x_{10} \oplus x_{11}}\\
\end{array}
\]

Synthesis example for the parity $x_1 \oplus x_6 \oplus x_8 \oplus x_{10}$ on a grid architecture. The gray nodes are the terminals and the Steiner tree is represented by bold edges. The chosen root for the spanning arborescence is the vertex $1$.

§.§ Architecture-aware algorithm

In this section we present an architecture-aware version of Algorithm <ref> by relying on the heuristic function $\mathcal{C}$ defined in Section <ref>.
The pseudo-code of the algorithm is provided in Algorithm <ref>.
We reuse some notations of Section <ref>: $P$ is a parity table, $\bs y$ is a parity of $P$ and $G_{\bs y}$ is the parity graph associated with $\bs y$.
We define $P^X$ as being the state of the parity table $P$ after performing the sequence of additions associated with any successors-first traversal of the arborescence $X$.
We denote by $T_{\bs y}$ the Steiner tree of $\bs y$ where the terminals are the vertices in the set $S=\{i \mid y_i = 1\}$.

We consider again the two steps process described in Section <ref> and modify it to take into account the architecture's connectivity.
Namely, we choose the parity $\bs y$ that minimizes $\mathcal{C}(\bs y)$ and we perform the synthesis of $\bs y$ in a architecture-aware manner such as described in Section <ref>.
In most cases, the number of minimum size fill-in of $T_{\bs y}$ is exponential with respect to the number of terminals.
As we want our algorithm to be scalable we get over this step by doing an arbitrary fill-in as presented in Algorithm <ref>.
Then, we must perform the sequence of additions associated with a successors-first traversal of any spanning arborescence $X$ of $T_{\bs y}$.
Our heuristic function $\mathcal{C}$ isn't merely based on the Hamming weight anymore, it implies that we cannot rely again on the minimum weight spanning arborescence to choose among the spanning arborescences of $T_{\bs y}$ as we have done in Section <ref>.
However, since $T_{\bs y}$ is a tree it only has $|V_{T_{\bs y}}|$ different spanning arborescences (a tree can have exactly one spanning arborescence rooted at each of its vertices), we can compute all the possibilities to choose the one that gives the best results and still have a polynomial time algorithm.

Let $\bs c^X = \text{sort}\{\mathcal{C}(\bs y) \mid \bs y \in P^X\}$ be the cost vector sorted in ascending order with respect to a spanning arborescence $X$ of $T_{\bs y}$.
As we don't rely on the minimum weight spanning arborescence anymore, we are not compelled to choose the spanning arborescence $X$ of $T_{\bs y}$ that minimizes $\sum_{\bs y \in P^X} \mathcal{C}(\bs y)$.
Our experiments have shown that choosing the spanning arborescence $X$ of $T_{\bs y}$ such that $\bs c^X_1$ is minimal leads to better results.
If several arborescences satisfy this property then we choose among them by taking the spanning arborescence $X$ such that $\bs c^X_2$ is minimal.
We repeat the process for $\bs c^X_i$ where $i \in [3, \ldots, m]$ until we only have one spanning arborescence left or until the end of the cost vectors is reached.
The chosen spanning arborescence $X$ of $T_{\bs y}$ then satisfies

\begin{align*}
X &= \operatorname*{argmin}\{\bs c^X \mid X \text{ is a spanning arborescence of }T_{\bs y}\}\\
&= \operatorname*{argmin}\{\text{sort}\{\mathcal{C}(\bs y) \mid \bs y \in P^X\} \mid X \text{ is a spanning arborescence of }T_{\bs y}\}\}.
\end{align*}

A parity table $P$ and a connectivity graph $G$
Circuit synthesizing a parity network associated with $P$ and compliant with the architectural constraints described by $G$
$C \leftarrow$ new empty circuit

$paths \leftarrow $ShortestPaths$(G)$

$P$ non-empty
$\displaystyle \bs y \leftarrow \operatorname*{argmin}_{\bs y} \{\mathcal{C}(\bs y) \mid \bs y \in P\}$

$P \leftarrow P \setminus \bs y$

$S \leftarrow \{i \mid y_i=1\}$

$T \leftarrow$ SteinerTree($S$, $paths$)

$CNOT_{i,j} \in\ $FillIn($T$, $S$)
$P_i \leftarrow P_i \oplus P_j$

$C \leftarrow C :: CNOT_{i,j}$

$\displaystyle X \leftarrow \operatorname*{argmin}_X \{$sort$\{\mathcal{C}(\bs y) \mid \bs y \in P^X\} \mid X$ is a spanning arborescence of $T\}$

$i \in\ $DepthFirstSearchPostordering$(X)$
$j \leftarrow $ direct predecessor of $i$ in $X$

$C \leftarrow C :: CNOT_{i, j}$

$P_i \leftarrow P_i \oplus P_j$

Architecture-aware parity network synthesis

Complexity analysis.
Let $n$ be the number of qubits and $m$ the number of parities.
Algorithm <ref> compute the Steiner tree of a parity with a complexity of $\mathcal{O}(n^2)$.
The complexity of Algorithm <ref> is majored by the task of finding the optimal spanning arborescence $X$, which takes $\mathcal{O}(mn^3)$ operations as it requires to compute $\mathcal{O}(mn)$ Steiner trees.
The algorithm performs $m$ iterations so the overall complexity is $\mathcal{O}(m^2n^3)$.
Some methods to further reduce the complexity of our algorithm are described in Section <ref>.

Further optimizations.
We refer to the implementation cost of a Steiner tree $T_{\bs y}$ as the value $\mathcal{C}(\bs y)$.
In most cases, in order to find the spanning arborescence $X$ satisfying $X = \operatorname*{argmin}\{\text{sort}\{\mathcal{C}(\bs y) \mid \bs y \in P^X\}\} = \operatorname*{argmin}\{\bs c^X\}$, we don't actually need to know all the values of the vector $\bs c^X$.
Our experiments show that the CNOT performances of our algorithm doesn't change when we only consider the $K=10$ first values of $\bs c^X$.
Indeed, it is rather rare to have $2$ arborescences $X_1$, $X_2$ such that $\bs c^{X_1}_i = \bs c^{X_2}_i \  \forall i \in [1, \ldots, K=10]$.
Also, when it happens, then choosing $X_1$ or $X_2$ doesn't make a significant difference in the algorithm performances.
Consequently, to determine the values $\bs c^{X}_i \ \forall i \in [1, \ldots, K]$ for an arborescence $X$, we only have to compute the $K$ Steiner trees of $\{T_{\bs y} \mid \bs y \in P^X\}$ having a minimum implementation cost $\mathcal{C}(\bs y)$.

In order to find the $K$ minimum cost Steiner trees we propose an algorithm that simultaneously constructs all the Steiner trees and stops when the $K$ minimum cost Steiner trees are found.
In this algorithm, only $K$ Steiner trees will be completely computed and all the others will be partially constructed, thus saving a considerable amount of time.
The pseudo-code of the algorithm is given in Algorithm <ref>, it takes as input the shortest paths of the graph, an integer $K$ and a set $S$ of sets $S_1, \ldots, S_m$ of terminals from which we want to find the $K$ minimum cost Steiner trees.

$S$, $paths$, $K$
$R \leftarrow$ new empty list

$L \leftarrow$ new list of empty stacks

$cost \leftarrow 0$

$i \in \{1, \ldots, |S|\}$
$u \leftarrow $ any vertex in $S_i$

$S_i \leftarrow S_i \setminus \{u\}$

$T_i \leftarrow $ Graph($\{u\}, \emptyset$)

$L$ non-empty
$Q \leftarrow L$.popleft

$Q$ non-empty
$i \leftarrow Q.$pop

$S_i$ is empty
$R \leftarrow R :: (cost,\; T_i)$

$|R| = K$

$\displaystyle u, v \leftarrow \operatorname*{argmin}_{u,v} \{|paths_{u,v}| \mid u \in T_i, v \in S_i\}$

$T_i \leftarrow T_i \cup paths_{u, v}$

$S_i \leftarrow S_i \setminus \{v\}$

$j \leftarrow 1 + 2 \times (|paths_{u, v}|-2)$

$cost \leftarrow cost + 1$

$K$ minimum cost Steiner trees

§.§ Lowering the complexity

In this section we describe some methods that can be applied to our algorithms in order to reduce their complexity and achieve the best trade-off possible between performances and running time constraints.
We demonstrate these methods in the case of constrained architectures, although they can also be applied in the same way on Algorithm <ref> for all-to-all connectivity.

Greedy method.
For this method we consider a greedier version of Algorithm <ref> by choosing an arbitrary spanning arborescence instead of the one satisfying $X = \operatorname*{argmin} \{\text{sort}\{\mathcal{C}(\bs y) \mid \bs y \in P^X\}$.
The algorithm then simply consists in choosing the parity $\bs y$ minimizing $\mathcal{C}(\bs y)$, performing the synthesis of $\bs y$ optimally and reiterating until all parities have been synthesized.
The complexity of this algorithm is $\mathcal{O}(m^2n^2)$.

Sliding window.
Another way of lowering the complexity of our algorithm is to only consider the parities that are in the range of a sliding window of size $\alpha$ that runs over the parity table.
At each iteration of our algorithm the parity which has just been synthesized is removed from the window and another parity not yet considered is added to the window.
With this method the complexity of our algorithm is lowered to $\mathcal{O}(\alpha mn^3)$.
By combining the greedy and sliding window methods, the complexity of the algorithm is further lowered to $\mathcal{O}(\alpha mn^2)$.

§.§ Benchmarks

We compare our algorithms with the algorithm proposed by Meijer-van de Griend et al. <cit.>, the phase polynomials synthesis algorithm implemented in Staq <cit.> and the Lazy-Synth algorithm <cit.>.
Beside Algorithm <ref>, we benchmark 3 modified versions of it stemming from Section <ref>: with a sliding window of size $50$, with the greedy method and with both the greedy method and a sliding window of size $50$.
All algorithms have been implemented in Python except for the Lazy-Synth algorithm which has been implemented in C++.
The standard deviation $\sigma$ is not represented in our benchmarks as it is particularly low ($\sigma < 10^{-2}$).

Figure <ref> compares the algorithms on a $3 \times 3$ grid with a density ranging from $1\%$ to $100\%$.
In this setting, our algorithm and the Lazy-Synth algorithm both offer the best CNOT performances, our algorithm being faster.
In Figure <ref>, for the Melbourne architecture, the Lazy-Synth algorithm provides the best performances with regards to the CNOT metrics but has an important running time.
All the variants of Algorithm <ref> have a better CNOT count and CNOT depth than the algorithms proposed by Meijer-van de Griend et al. and Staq.
Morever, the greedy version with a sliding window of size $50$ has a similar running time.
We tested the algorithms on various architectures and didn't notice any significant changes.

We also did benchmarks for a constant number of parities with an increasing number of qubits, the results are shown in Figure <ref>.
We didn't run the sliding window variants for this benchmark because their usefulness is limited by the small constant number of parities.
We can notice that the algorithm proposed by Meijer-van de Griend et al. has a smaller CNOT depth than the other algorithms when the number of qubits is greater than $36$, however it also has the worst CNOT count.
Algorithm <ref> and the Lazy-Synth algorithm are producing again the smallest circuits.
Nevertheless, as it can be seen in Table <ref>, they both have a rapidly growing computational time.
The greedy method seems to offer the best compromise.
Indeed, it produces circuits with significantly better CNOT count than the algorithms proposed by Meijer-van de Griend et al. and Staq, and has a much more tolerable running time than Algorithm <ref> and the Lazy-Synth algorithm.

Grid $3 \times 3$

14 qubits IBM Melbourne

Scaling of the circuit size, circuit depth and computational time on a $3 \times 3$ grid architecture and the 14 qubits IBM's Melbourne architecture.
Each point is averaged over $1000$ randomly generated parity tables.

Scaling of the circuit size and circuit depth on square grids for $100$ parities.
Each point is averaged over $100$ randomly generated parity tables.

Average computational time in seconds on square grids for $100$ parities.

2*Algorithm     8c|Grid size

$3\times 3$     $4\times 4$     $5\times 5$     $6\times 6$     $7\times 7$     $8\times 8$     $9\times 9$     $10 \times 10$

Proposed     0.303     1.439     5.763     18.650     50.745     137.250     339.595     768.073

Greedy     0.088     0.209     0.498     1.070     2.095     4.270     8.254     15.098

Staq     0.038     0.078     0.157     0.268     0.421     0.648     0.990     1.616

Meijer-van de     2*0.066     2*0.217     2*0.517     2*0.941     2*1.577     2*2.584     2*3.960     2*5.977

Griend et al.

Lazy-Synth     0.560     5.524     80.920     523.639     2381.802     -     -     -

§ CONCLUSION

We presented heuristic algorithms for the synthesis of circuits over the $\{CNOT, R_Z\}$ gate set.
We covered both cases of full connectivity and partial connectivity for NISQ architectures.
When compared to the state of the art, the benchmarks have shown that our algorithms are producing circuits of smaller or comparable size.
State of the art algorithms yielding analogous CNOT count and CNOT depth performances are outperformed when it comes to the execution time.

Our framework could be further expanded by adding depth and width to our algorithms, which would result in an even smaller CNOT count at the cost of a runtime increase.
This methods was utilized in the Lazy-Synth framework, instead of going in the same direction we presented methods to lower the complexity of our algorithms while preserving as much performances as possible, with the aim of providing a more comprehensive toolkit for the synthesis of phase polynomials.

Finally, the modularity of our approach makes it easily adaptable for the synthesis of other circuits such as sequences of Pauli rotations, a generalization of phase polynomials.
This is not the case for the Gray-Synth algorithm and all its architecture-aware derivatives.
Indeed, the Gray-Synth algorithm relies on the fact that the row $j$ of the parity table isn't affected when applying the gate $CNOT_{i, j}$.
This invariant isn't true when considering sequences of Pauli rotations.
Therefore the Gray-Synth algorithm needs to be coupled with phase polynomials extraction routines in order to compile a generic circuit.

§ ACKNOWLEDGMENTS

This work was supported in part by the French National Research Agency (ANR) under the research project SoftQPRO ANR-17-CE25-0009-02, and by
the DGE of the French Ministry of Industry under the research project PIAGDN/QuantEx P163746-484124.

§ BENCHMARKS DATA

§.§ All-to-all connectivity

CNOT count averaged over 1000, 100, 10 and 10 randomly generated parity tables for 7, 10, 13 and 16 qubits respectively.

2*Density     3c|7 qubits     3c|10 qubits     3c|13 qubits     2c|16 qubits

Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed

1%     -     -     -     35     24     26     308     192     196     2540     1501

5%     14     10     11     135     91    91     1123     721    689     9096     5556

10%     25     18     19     230     157    156     1896     1243    1177     15190     9551

15%     34     25     26     306     217    212     2524     1703    1619     20245     12980

20%     42     32    32     376     271    266     3078     2116    2033     24740     16221

25%     50     38     39     437     324    318     3573     2529    2437     28687     19421

30%     56     44     45     494     372    368     4024     2915    2852     32324     22581

35%     63     50     51     548     421    418     4443     3309    3263     35760     25768

40%     69     56     57     596     468     470     4839     3706    3669     38842     28903

45%     74     61     62     640     517     520     5218     4080    4074     41766     32005

50%     79     67     68     685     565     570     5548     4474     4488     44512     35104

55%     84     72     74     725     611     622     5873     4875     4888     47139     38284

60%     89     78     80     762     659     673     6197     5251     5311     49544     41473

65%     94     84     86     800     707     722     6487     5645     5712     51929     44668

70%     98     89     92     836     756     774     6757     6038     6118     54157     47838

75%     102     95     98     866     804     824     7017     6436     6523     56276     51041

80%     106     100     104     900     850     874     7283     6835     6942     58293     54265

85%     110     106     110     929     899     927     7515     7231     7345     60249     57477

90%     113     111     116     958     944     975     7742     7619     7754     62068     60712

95%    116     116     122    986     987     1027    7969     7981     8162    63836     63895

100%    120     120     128    1013     1013     1077    8178     8178     8572    65519     65519

CNOT depth averaged over 1000, 100, 10 and 10 randomly generated parity tables for 7, 10, 13 and 16 qubits respectively.

2*Density     3c|7 qubits     3c|10 qubits     3c|13 qubits     2c|16 qubits

Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed

1%     -     -     -     32     14     21     298     114     151     2523     857

5%     13     8     9     130     60     71     1112     479     532     9079     3678

10%     23     13     15     223     110     122     1882     856     896     15173     6669

15%     32     19     21     299     154     167     2510     1213     1264     20225     9344

20%     40     24     26     369     200     214     3065     1563     1628     24718     12114

25%     47     29     31     431     242     257     3560     1918     1996     28668     15119

30%     54     34     36     488     286     302     4006     2275     2364     32305     18016

35%     60     40     41     541     331     343     4431     2668     2709     35742     21167

40%     66     45     47     590     378     390     4828     3051    3045     38814     24299

45%     72     51    51     634     426     431     5205     3436    3403     41744     27466

50%     77     56    56     678     476    474     5538     3841    3767     44490     30661

55%     82     62    61     719     528    515     5860     4269    4101     47110     33863

60%     87     68    66     755     578    561     6183     4622    4433     49513     37056

65%     92     75    71     794     629    600     6475     5049    4808     51904     40347

70%     95     81    76     829     684    645     6740     5471    5144     54134     43575

75%     100     88    82     860     736    686     7006     5869    5465     56244     46892

80%     104     94    86     895     792    727     7273     6303    5814     58269     50189

85%     108     101    92     923     847    774     7505     6738    6204     60223     53503

90%     112     107    97     953     903    814     7730     7212    6509     62043     56948

95%     115     114    102     982     961    856     7959     7714    6828     63818     60866

100%     120     120    107     1013     1013    901     8178     8178    7218    65519     65519

Runtime in seconds averaged over 1000, 100, 10 and 10 randomly generated parity tables for 7, 10, 13 and 16 qubits respectively.

2*Density     3c|7 qubits     3c|10 qubits     3c|13 qubits     2c|16 qubits

Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed     Lazy-Synth     Gray-Synth     Proposed

1%     -     -     -    0.006     0.011     0.045    0.058     0.094     0.798    0.583     0.839

5%    0.002     0.005     0.024    0.025     0.048     0.217    0.270     0.412     9.207    2.387     5.376

10%    0.004     0.010     0.034    0.047     0.085     0.522    0.464     0.785     28.997    4.495     13.714

15%    0.005     0.014     0.045    0.064     0.122     0.894    0.644     1.197     53.802    6.662     24.442

20%    0.007     0.018     0.056    0.079     0.156     1.393    0.775     1.606     94.539    8.909     37.710

25%    0.008     0.021     0.069    0.092     0.194     1.950    0.798     2.056     138.939    11.294     53.946

30%    0.009     0.025     0.082    0.091     0.221     2.568    0.943     2.561     187.390    13.733     72.766

35%    0.010     0.029     0.097    0.101     0.261     3.385    1.046     3.127     258.234    16.311     94.937

40%    0.011     0.032     0.113    0.111     0.303     4.283    1.167     3.708     329.136    19.042     119.998

45%    0.012     0.036     0.129    0.120     0.327     5.207    1.285     4.361     405.330    22.268     150.989

50%    0.013     0.040     0.145    0.129     0.368     6.215    1.398     5.077     489.180    25.036     181.316

55%    0.013     0.044     0.165    0.138     0.409     7.333    1.508     5.829     569.902    28.043     215.597

60%    0.015     0.048     0.185    0.146     0.447     8.729    1.622     6.541     688.243    30.944     253.264

65%    0.015     0.051     0.207    0.153     0.471     10.186    1.731     7.374     819.486    34.068     293.270

70%    0.016     0.055     0.228    0.162     0.515     11.751    1.837     8.243     943.779    37.294     336.222

75%    0.017     0.059     0.253    0.168     0.547     13.318    1.935     9.213     1072.623    40.457     382.445

80%    0.017     0.063     0.279    0.176     0.591     14.985    2.030     10.121     1195.277    43.895     431.766

85%    0.018     0.067     0.301    0.182     0.640     16.741    2.149     11.117     1312.032    47.108     483.807

90%    0.018     0.070     0.328    0.189     0.695     18.431    2.215     12.084     1452.226    50.472     539.097

95%    0.019     0.075     0.355    0.196     0.710     20.374    2.313     13.160     1596.560    53.741     596.925

100%    0.020     0.078     0.380    0.201     0.751     22.193    2.402     14.143     1748.428    56.923     652.193