# Time-optimal multi-qubit gates: Complexity, efficient heuristic and gate-
time bounds

P. Baßler<EMAIL_ADDRESS>Institute for Theoretical Physics, Heinrich Heine
University Düsseldorf, Germany  M. Heinrich  Institute for Theoretical
Physics, Heinrich Heine University Düsseldorf, Germany  M. Kliesch
<EMAIL_ADDRESS>Institute for Quantum Inspired and Quantum
Optimization, Hamburg University of Technology, Germany

###### Abstract

Multi-qubit entangling interactions arise naturally in several quantum
computing platforms and promise advantages over traditional two-qubit gates.
In particular, a fixed multi-qubit Ising-type interaction together with
single-qubit $\operatorname{\mathsf{X}}$-gates can be used to synthesize
global $\operatorname{\mathsf{ZZ}}$-gates ($\operatorname{\mathsf{GZZ}}$
gates). In this work, we first show that the synthesis of such quantum gates
that are time-optimal is $\mathsf{NP}$-hard. Second, we provide explicit
constructions of special time-optimal multi-qubit gates. They have constant
gate times and can be implemented with linearly many
$\operatorname{\mathsf{X}}$-gate layers. Third, we develop a heuristic
algorithm with polynomial runtime for synthesizing fast multi-qubit gates.
Fourth, we derive lower and upper bounds on the optimal
$\operatorname{\mathsf{GZZ}}$ gate-time. Based on explicit constructions of
$\operatorname{\mathsf{GZZ}}$ gates and numerical studies, we conjecture that
any $\operatorname{\mathsf{GZZ}}$ gate can be executed in a time
$\operatorname{\mathrm{O}}(n)$ for $n$ qubits. Our heuristic synthesis
algorithm leads to $\operatorname{\mathsf{GZZ}}$ gate-times with a similar
scaling, which is optimal in this sense. We expect that our efficient
synthesis of fast multi-qubit gates allows for faster and, hence, also more
error-robust execution of quantum algorithms.

##  1 Introduction

Any quantum computation requires to decompose its logical operations into the
platform’s native instruction set. The performance of the computation depends
heavily on the available instructions and their implementation in the quantum
hardware. In particular for early quantum devices a major challenge is posed
by their short decoherence times, which limits the runtime of a quantum
computation significantly. Therefore, it is not only necessary to optimize the
number of native instructions but also their execution time.

Ising-type interactions give rise to an important and rich class of
Hamiltonians ubiquitous in several quantum computing platforms [1, 2, 3, 4,
5]. Previously, we have utilized these Ising-type interactions in a new
synthesis method [6]. In particular, we have considered the problem of
synthesizing time-optimal multi-qubit gates on a quantum computing platform
that supports the following basic operations:

  1. (I)

single-qubit rotations can be executed in parallel, and

  2. (II)

it offers a fixed Ising-type interaction.

The corresponding synthesis of global $\operatorname{\mathsf{ZZ}}$-gates
($\operatorname{\mathsf{GZZ}}$ gates) is given by the minimization of the
overall gate time, which can be written as a linear program (LP) [6]. This LP
is exponentially large in the number of qubits. It has been unclear whether
such multi-qubit gates can be synthesized in a computationally efficient way
while keeping the gate time optimal.

In this work, we prove that this synthesis problem is $\mathsf{NP}$-hard and
provide a close-to-optimal efficient heuristic solution. To establish this
hardness result, we draw the connection between the synthesis problem and
graph theory. The _cut polytope_ is defined as the convex hull of binary
vectors representing the possible cuts of a given graph [7]. We provide a
polynomial time reduction of the membership problem of the cut polytope to the
synthesis of time-optimal multi-qubit gates. Since this membership problem is
$\mathsf{NP}$-complete [8], the synthesis of time-optimal multi-qubit gates is
$\mathsf{NP}$-hard. This is akin to the $\mathsf{NP}$-hardness of finding an
optimal control pulse for multi-qubit gates using the Mølmer-Sørensen
mechanism [4, 5].

We provide several ways to circumvent the hardness of time-optimal
$\operatorname{\mathsf{GZZ}}$-gates. First, we provide explicit constructions
of time-optimal multi-qubit gates realizing nearest-neighbor coupling under
physically motivated assumptions. Such constructed nearest-neighbor multi-
qubit gates exhibit a constant gate time and can be implemented with only
linearly many single-qubit gate layers. We then use these ideas to define
relaxations of the underlying linear program, leading to a hierarchy of
polynomial-time algorithms for the synthesis of fast multi-qubit gates. By
increasing the level in the hierarchy, this heuristic approach can be adapted
to provide substantially better approximations to the optimal solution at the
cost of higher polynomial runtime. For a small number of qubits, numerical
experiments show that the so-obtained gate times are close to the optimal
solution and come with significant runtime savings.

Among others, we prove bounds on the minimal multi-qubit gate time, and
conjecture that it scales at most linear with the number of qubits. This claim
is supported by a class of explicit constructions of time-optimal multi-qubit
gates achieving the linear upper bound. Moreover, we provide numerical
evidence that these explicit solutions in fact yield the longest gate time for
a small number of qubits.

We expect our results to be useful for the implementation of time-optimal
multi-qubit gates in noisy and intermediate scale quantum (NISQ) devices and
beyond. The polynomial-time heuristic algorithm makes it possible to
efficiently synthesize fast multi-qubit gates for a growing number of qubits.
The here considered multi-qubit gates have been useful in a number of
different applications, some of which we investigated in previous work [6].
Furthermore, they are the main building blocks for a class of Instantaneous
Quantum Polynomial time (IQP) circuits which might be classically hard to
simulate [9]. More recently, it was shown that quantum memory circuits and
boolean functions can be implemented with a constant number of
$\operatorname{\mathsf{GZZ}}$ gates and additional ancilla qubits [10].
Moreover, there is an ancilla-free construction of multi-qubit Clifford
circuits using at most $26$ $\operatorname{\mathsf{GZZ}}$ gates [11]. As noted
in Ref. [11], there is also a shorter implementation for $n\leq 2^{12}$,
requiring only $2(\log_{2}(n)+1)$ $\operatorname{\mathsf{GZZ}}$ gates. This
implementation is based on a decomposition in Ref. [12] and the $\log$-depth
implementation of a $\operatorname{\mathsf{CX}}$ circuit in formula (4) of
Ref. [13].

This paper is structured as follows: We first give a brief introduction to the
synthesis of time-optimal multi-qubit gates [6]. In Section 4 we proof that
the time-optimal multi-qubit synthesis problem is $\mathsf{NP}$-hard. However,
in Section 5 we explicitly construct a certain class of time-optimal multi-
qubit gates with constant gate time. The heuristic algorithm based on the
ideas of the previous section is introduced in Section 6. Section 7 provides
gate time bounds for time-optimal multi-qubit gates. Finally, Section 8
presents numerical evidence for the conjectured linear gate-time scaling and
numerical benchmarks for the heuristic algorithm.

##  2 Synthesizing multi-qubit gates with Ising-type interactions

In this section, we give a short introduction to the synthesis of time-optimal
multi-qubit gates. For more details we refer to the first two sections of Ref.
[6].

On the abstract quantum computing platform with $n$ qubits specified by the
requirements (I) and (II) above, interactions between the qubits are generated
by an Ising-type Hamiltonian

$H_{ZZ}\coloneqq-\sum_{i<j}^{n}J_{ij}\sigma_{z}^{i}\sigma_{z}^{j}\,,$ (1)

where $\sigma_{z}^{i}$ is the Pauli-Z operator acting on the $i$-th qubit.
Note, that diagonal terms, where $i=j$, are excluded since they only change
the Hamiltonian by an energy offset. By $J$ we denote the symmetric matrix
with entries $J_{ij}$ in the upper triangular part and vanishing diagonal. We
call $J$ the _physical coupling matrix_.

Conjugating the Hamiltonian $H_{ZZ}$ with $\operatorname{\mathsf{X}}$ gates on
the qubits indicated by the binary vector
$\boldsymbol{b}\in\mathbb{F}_{2}^{n}$ yields

$\displaystyle\sigma_{x}^{\boldsymbol{b}}H_{ZZ}\sigma_{x}^{\boldsymbol{b}}$
$\displaystyle=-\sum_{i<j}^{n}J_{ij}\sigma_{x}^{b_{i}}\sigma_{x}^{b_{j}}\sigma_{z}^{i}\sigma_{z}^{j}\sigma_{x}^{b_{i}}\sigma_{x}^{b_{j}}$
(2)
$\displaystyle=-\sum_{i<j}^{n}J_{ij}(-1)^{b_{i}}(-1)^{b_{j}}\sigma_{z}^{i}\sigma_{z}^{j}$
$\displaystyle=-\sum_{i<j}^{n}J_{ij}m_{i}m_{j}\sigma_{z}^{i}\sigma_{z}^{j}\eqqcolon
H(\boldsymbol{m})\,,$

where we define the _encoding_ $\boldsymbol{m}\coloneqq(-1)^{\boldsymbol{b}}$
entry wise. The sign of the interaction between qubit $i$ and $j$ is given by
$m_{i}m_{j}\in\\{-1,+1\\}$. We call $H(\boldsymbol{m})$ the _encoded
Hamiltonian_.

Given time steps $\lambda_{\boldsymbol{m}}\geq 0$ during which the encoding
$\boldsymbol{m}$ is used, we consider unitaries of the form

$\prod_{\boldsymbol{m}}\mathrm{e}^{-\mathrm{i}\lambda_{\boldsymbol{m}}H(\boldsymbol{m})}=\mathrm{e}^{-\mathrm{i}\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}H(\boldsymbol{m})}\eqqcolon\mathrm{e}^{-\mathrm{i}H}\,,$
(3)

where we used that the diagonal Hamiltonians $H(\boldsymbol{m})$ mutually
commute. For all possible encodings $\boldsymbol{m}\in\\{-1,+1\\}^{n}$ we
collect the time steps $\lambda_{\boldsymbol{m}}$ in a vector
$\boldsymbol{\lambda}\in\mathbb{R}^{2^{n}}$ and interpret
$t=\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}$ as the total gate time of
the unitary $\mathrm{e}^{-\mathrm{i}H}$, implemented by the sequence of
unitaries (3). Moreover, we use the symmetry

$(-\boldsymbol{m})(-\boldsymbol{m})^{T}=\boldsymbol{m}\boldsymbol{m}^{T}$ (4)

to reduce the degrees of freedom in $\boldsymbol{\lambda}$ from $2^{n}$ to
$2^{n-1}$ by adding up $\lambda_{\boldsymbol{m}}+\lambda_{-\boldsymbol{m}}$ to
a single time step.

The so generated unitary is the time evolution operator under the _total
Hamiltonian_

$H=-\sum_{i<j}^{n}A_{ij}\sigma_{z}^{i}\sigma_{z}^{j}\,,$ (5)

where we have defined the _total coupling matrix_

$A\coloneqq
J\circ\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}\boldsymbol{m}\boldsymbol{m}^{T}\,,$
(6)

and used the linearity of the Hadamard (entry-wise) product $\circ$. By
construction, $A$ is a symmetric matrix with vanishing diagonal. Let us define
the $\binom{n}{2}$-dimensional subspace of symmetric matrices with vanishing
diagonal by

$\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})\coloneqq\left\\{M\in\mathrm{Sym}(\mathbb{R}^{n})\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}M_{ii}=0\
\forall i\in[n]\right\\}\,.$ (7)

For $A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$, we define an
associated multi-qubit gate $\operatorname{\mathsf{GZZ}}(A)$, where
$\operatorname{\mathsf{GZZ}}$ stands for “global $ZZ$ interactions”,

$\displaystyle\operatorname{\mathsf{GZZ}}(A)$
$\displaystyle\coloneqq\exp\left(\mathrm{i}\sum_{i<j}^{n}A_{ij}\sigma_{z}^{i}\sigma_{z}^{j}\right)\,.$
(8)

To determine which matrices can be decomposed as in Eq. 6, we denote the non-
zero index set of a symmetric matrix as
$\operatorname{nz}(M)\coloneqq\\{(i,j)\,|\,M_{ij}\neq 0,i<j\\}$. Then, the
subspace of matrices $A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$
that can be decomposed as in Eq. 6 is exactly given by the condition
$\operatorname{nz}(A)\subseteq\operatorname{nz}(J)$, which we assume to hold
from now on. Thus, all-to-all connectivity enables to decompose any coupling
matrix $A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ but is not
strictly required by our approach. We call the number of encodings
$\boldsymbol{m}$ needed for the decomposition the _encoding cost_ of
$\operatorname{\mathsf{GZZ}}(A)$, and
$\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}$ the _total
$\operatorname{\mathsf{GZZ}}$ time_. Note, that both quantities depend on the
chosen decomposition.

It is convenient to abstract the following analysis from the physical details
given by $J$. For a matrix
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ with
$\operatorname{nz}(A)\subseteq\operatorname{nz}(J)$, its possible
decompositions are in one-to-one correspondence with the decompositions of the
matrix $M\coloneqq A\oslash
J\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ where

$(A\oslash J)_{ij}\coloneqq\begin{cases}A_{ij}/J_{ij},&i\neq j\text{ and
}J_{ij}\neq 0\,,\\\ 0,&\text{otherwise}\,.\end{cases}$ (9)

We further define the linear operator

$\displaystyle\mathcal{V}:\mathbb{R}_{\geq 0}^{2^{n-1}}$
$\displaystyle\rightarrow\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n}),$
(10) $\displaystyle\boldsymbol{\lambda}$
$\displaystyle\mapsto\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}\boldsymbol{m}\boldsymbol{m}^{T},$

represented in the standard basis by a matrix

$V\in\\{-1,+1\\}^{\binom{n}{2}\times(2^{n-1})}\,.$ (11)

Let
$\boldsymbol{v}:\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})\rightarrow\mathbb{R}^{\binom{n}{2}}$
be the (row-wise) vectorization of the upper triangular part of the matrix
input such that the columns of $V$ are given by
$\boldsymbol{v}(\boldsymbol{m}\boldsymbol{m}^{T})$. Our objective is to
minimize the total gate time and the amount of $\boldsymbol{m}$’s needed to
express the matrix $M\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$. To
this end we formulate the following linear program (LP): 1^T λ V λ =v(M) λ
∈R^2^n-1_≥0 , where $\boldsymbol{1}=(1,1,\dots,1)$ is the all-ones vector such
that
$\boldsymbol{1}^{T}\boldsymbol{\lambda}=\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}$.
A _feasible solution_ is an assignment of the variables that fulfills all
constraints of the optimization problem and an _optimal solution_ is a
feasible solution which also minimizes/maximizes the objective function.
Throughout this paper, we indicate optimal solutions by an asterisk, e.g.
$\boldsymbol{\lambda}^{*}$. Note, that the LP (11) has a feasible solution for
any symmetric matrix $M$ with vanishing diagonal [6, Theorem II.2]. The theory
of linear programming then guarantees the existence of an optimal solution
with at most $\binom{n}{2}$ non-zero entries [6, Proposition II.3].

A standard tool in convex optimization is _duality_ [14] which will be used in
Section 7. The dual LP to the LP (11) reads as follows: ⟨M , y⟩ V^T y ≤1 y
∈R^(n2) , with the inner product
$\langle\,\cdot\,,\,\cdot\,\rangle:\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})\times\mathbb{R}^{\binom{n}{2}}\rightarrow\mathbb{R}$,
$\langle M,\boldsymbol{y}\rangle\mapsto\boldsymbol{v}(M)^{T}\boldsymbol{y}$.
Here, inequalities between vectors are to be understood entry-wise. A simple,
but important fact is the following: If $\boldsymbol{\lambda}^{*}$ is an
optimal solution to the primal LP (11), then any feasible solution
$\boldsymbol{y}$ to the dual LP (2) provides a _lower bound_ as $\langle
M,\boldsymbol{y}\rangle\leq\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}$.
Moreover, the feasibility of the LP (11) implies that we have _strong duality_
: if $\boldsymbol{y}^{*}$ is a dual optimal solution, then we have equality
between the optimal values, $\langle
M,\boldsymbol{y}^{*}\rangle=\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}$.

##  3 Main results

We want to highlight our main contributions. First, we provide the hardness
results for the synthesis of time-optimal $\operatorname{\mathsf{GZZ}}$ gates.

######  Theorem (Theorem 6).

The decision version of the LP (11), is $\mathsf{NP}$-complete.

We say that the synthesis of time-optimal multi-qubit gates (solving LP (11))
is $\mathsf{NP}$-hard in the sense of the function problem extension of the
decision problem class $\mathsf{NP}$ [15]. We circumvent the hardness of the
time-optimal multi-qubit gate synthesis by providing an explicit construction
of $\operatorname{\mathsf{GZZ}}$ gates realizing next neighbor coupling with
constant gate time and linear encoding cost. The assumption of a constant
subdiagonal of $J$ is physically motivated and can be realized in an ion trap
[16].

######  Theorem (Theorem 10, informal).

Let the subdiagonal of $J$ and $A$ be a constant with values $c$ and $\varphi$
respectively. Then $\operatorname{\mathsf{GZZ}}(A)$ on $n$ qubits has the
encoding cost of $d\leq 2n$ and constant total $\operatorname{\mathsf{GZZ}}$
time $2\varphi/c$. This total gate time is optimal.

In Section 6 we define Algorithm 3 and introduce LP (6), a polynomial time
heuristic to synthesize $\operatorname{\mathsf{GZZ}}$ gates with small, but
not necessarily minimal, gate time. This heuristic does not rely on any
further assumptions and is applicable for arbitrary
$J,A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$. In Section 8 we show
numerically that this heuristic leads to $\operatorname{\mathsf{GZZ}}$ gate
times close to the optimum while solving the synthesis problem much faster.

Furthermore, we proof lower and upper bounds on the
$\operatorname{\mathsf{GZZ}}$ gate time.

######  Theorem (Theorem 15).

The optimal total gate time of $\operatorname{\mathsf{GZZ}}(A)$ with
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ is lower and upper
bounded by

$\lVert A\oslash
J\rVert_{\ell_{\infty}}\leq\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}\leq\lVert
A\oslash J\rVert_{\ell_{1}}\,.$ (12)

Here, the upper bound scales quadratic in $n$ for a constant matrix
$M=A\oslash J$. This upper bound is loose in the sense that it also holds for
the $\operatorname{\mathsf{GZZ}}$ gates constructed by the heuristic in
Section 6. Therefore, we tighten the upper bound to a linear scaling in $n$.

######  Conjecture (16).

The optimal gate time of $\operatorname{\mathsf{GZZ}}(A)$ with
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ is tightly upper bounded
by

$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}\leq\lVert A\oslash
J\rVert_{\ell_{\infty}}\cdot\begin{cases}n\,,&\text{for odd }n\,,\\\
n-1\,,&\text{for even }n\,.\end{cases}$ (13)

We provide evidence for this conjecture using an explicit construction
(Theorem 21) that realizes the conjectured upper bound for any $n$, as well as
numerical evidence for $n\leq 8$ (Fig. 1). Unfortunately, we were not able to
prove this result and state the challenges in Appendix A.

##  4 Synthesizing time-optimal GZZ gates is $\mathsf{NP}$-hard

In this section, we investigate the complexity of solving the gate synthesis
problem stated as LP (11). We observe that LP (11) is an optimization over the
convex cone generated by

$\mathcal{E}_{n}\coloneqq\left\\{\boldsymbol{m}\boldsymbol{m}^{T}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\boldsymbol{m}\in\\{-1,+1\\}^{n},m_{n}=+1\right\\}\,,$
(14)

which is the set of outer products generated by all possible encodings
$\boldsymbol{m}$. Due to the symmetry Eq. 4 we can uniformly fix the value of
one entry of $\boldsymbol{m}$. We chose the convention $m_{n}=+1$. In the
literature $\mathcal{E}_{n}$ is also known as the _elliptope_ of rank one
matrices [17]. In the following, we consider the polytope

$\operatorname{conv}(\mathcal{E}_{n})\coloneqq\left\\{\sum_{i}\lambda_{i}\boldsymbol{r}_{i}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\sum_{i}\lambda_{i}=1,\lambda_{i}\geq
0,\boldsymbol{r}_{i}\in\mathcal{E}_{n}\ \right\\}\,,$ (15)

and show the connection to graph theory, in particular to the cuts of graphs.

######  Definition 1 (cut polytope [7]).

Let $K_{n}=(V_{n},E_{n})$ be a complete graph with $n$ vertices. Denote
$\delta(X)$ the set of all edges with one endpoint in $X\subset V_{n}$ and the
other endpoint in its complement $\bar{X}$, i.e., $\delta(X)$ defines the cut
between $X$ and $\bar{X}$. Let $\chi^{\delta(X)}\in\\{0,1\\}^{|E_{n}|}$ denote
the characteristic vector of a cut, with $\chi_{e}^{\delta(X)}=1$ if
$e\in\delta(X)$ and $\chi_{e}^{\delta(X)}=0$ otherwise. We define the _cut
polytope_ as the convex hull of the characteristic vectors

$\operatorname{CUT}_{n}\coloneqq\operatorname{conv}\left\\{\chi^{\delta(X)}\in\\{0,1\\}^{|E_{n}|}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}X\subseteq
V_{n}\right\\}\,.$ (16)

######  Lemma 2.

For all $n$, $\operatorname{CUT}_{n}$ is isomorphic to
$\operatorname{conv}(\mathcal{E}_{n})$.

###### Proof.

For each $X\subset V_{n}$ we set $m_{i}=+1$ if $i\in X$ and $m_{i}=-1$ if
$i\in\bar{X}$. Note, that there are $2^{n-1}$ different pairs of $X$ and
$\bar{X}$. We then have $m_{i}m_{j}=-1$ if $i\in X$ and $j\in\bar{X}$ (or the
other way around), and $m_{i}m_{j}=+1$ if $i,j\in X$ or $i,j\in\bar{X}$. So
the characteristic vector can be written as
$\chi_{e}^{\delta(X)}=(1-m_{i}m_{j})/2$ for each edge $e\in E_{n}$ connecting
vertices $i$ and $j$. This is clearly a bijective affine map between the
vertices and thus $\operatorname{CUT}_{n}$ is isomorphic to
$\operatorname{conv}(\mathcal{E}_{n})$. ∎

The following decision problems are membership problems, where the task is to
decide if a given element $\boldsymbol{x}$ belongs to a set or not. In our
case $\boldsymbol{x}$ is a vector and the set is a polytope.

######  Problem 3 ($\operatorname{CUT}_{n}$ membership).

Instance

The adjacency matrix $M\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{Q}^{n}_{\geq
0})$ of a weighted undirected graph with non-negative weights.

Question

Is $M\in\operatorname{CUT}_{n}$?

######  Problem 4 ($\operatorname{conv}(\mathcal{E}_{n})$ membership).

Instance

The matrix $M\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{Q}^{n})$.

Question

Is $M\in\operatorname{conv}(\mathcal{E}_{n})$? That is, does there exist a
decomposition $M=\sum_{i}\lambda_{i}\boldsymbol{r}_{i}$ with
$\boldsymbol{r}_{i}\in\mathcal{E}_{n}$ such that $\lambda_{i}\geq 0$ and
$\sum_{i}\lambda_{i}=1$?

It is well known that the membership problem of the cut polytope, 3, and 4 are
$\mathsf{NP}$-complete [18, 8, 19]. Next, we state the decision version of our
gate synthesis optimization.

######  Problem 5 (time-optimal multi-qubit gate synthesis).

Instance

The matrix $M=A\oslash J\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{Q}^{n})$
and a constant $K\in\mathbb{Q}_{\geq 0}$.

Question

Is there a decomposition $M=\sum_{i}\lambda_{i}\boldsymbol{r}_{i}$ with
$\boldsymbol{r}_{i}\in\mathcal{E}_{n}$ such that $\lambda_{i}\geq 0$ and
$\sum_{i}\lambda_{i}\leq K$?

######  Theorem 6.

5, which is the decision version of the LP (11), is $\mathsf{NP}$-complete.

###### Proof.

A solution of 5 can be verified in polynomial time since there always exists a
decomposition $\sum_{i}\lambda_{i}\boldsymbol{r}_{i}$ of $M$ with minimal
$\sum_{i}\lambda_{i}$ which has at most $\binom{n}{2}=n/2(n-1)$ non-zero terms
[6, Proposition II.3]. Therefore, 5 is $\mathsf{NP}$.

To show that 5 is $\mathsf{NP}$-hard, we construct a polynomial-time mapping
reduction from 4 to 5. Given the matrix
$M\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{Q}^{n})$ and a constant
$K\in\mathbb{Q}_{\geq 0}$ as an instance for 5, let $\lambda_{i}\geq 0$ be the
positive coefficients of the decomposition. If we find
$\sum_{i}\lambda_{i}<K$, then we can always add additional $\lambda$’s such
that equality holds. We choose the additional $\lambda$’s as the coefficients
of the decomposition for the all-zero matrix, see e.g. Lemma 8 below with
$\boldsymbol{k}=\boldsymbol{1}_{n}$ for an explicit construction. We define
the matrix $M^{\prime}\coloneqq M/K$ and the positive coefficients
$\lambda^{\prime}_{i}\coloneqq\lambda_{i}/K$. Then
$M^{\prime}=\sum_{i}\lambda^{\prime}_{i}\boldsymbol{r}_{i}$ with
$\boldsymbol{r}_{i}\in\mathcal{E}_{n}$ and $\sum_{i}\lambda^{\prime}_{i}=1$. ∎

##  5 Time-optimal GZZ synthesis for special instances

Although, solving the LP (11) is $\mathsf{NP}$-hard we present explicit
optimal solutions for certain families of instances which is equivalent to
constructing time-optimal $\operatorname{\mathsf{GZZ}}$ gates. The
constructions of this section yield a qubit-independent total
$\operatorname{\mathsf{GZZ}}$ time which satisfies the optimal lower bound of
Lemma 7 below. Moreover, we show that some $\operatorname{\mathsf{GZZ}}$ gates
can be synthesized with an encoding cost independent of the number of qubits.
However, most of these constructions assume constant values for the elements
of the band-diagonal of the physical coupling matrix $J$. These assumptions
are relaxed throughout this section providing explicit optimal solutions for
physically relevant cases. These results build the foundation for the
heuristic algorithm for fast $\operatorname{\mathsf{GZZ}}$ gate synthesis in
the next section.

By $\lVert M\rVert_{\ell_{p}}$ we denote the $\ell_{p}$-norm of a symmetric
matrix $M$, which is given as the $\ell_{p}$-norm of a vector
$\boldsymbol{v}(M)$ containing all lower/upper triangular matrix elements in
some order. First, we proof a lower bound on the optimal total
$\operatorname{\mathsf{GZZ}}$ time which can be used to verify time-
optimality.

######  Lemma 7.

For any $M\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ the optimal
objective function value of the LP (11) is lower bounded by

$\lVert
M\rVert_{\ell_{\infty}}\leq\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}\,.$ (17)

###### Proof.

The lower bound can be verified by the fact that the matrix representation $V$
of the linear operator in Eq. 10 only has entries $\pm 1$ and that
$\boldsymbol{\lambda}^{*}$ is non-negative. Thus, it holds that $\lVert
M\rVert_{\ell_{\infty}}=\lVert
V\boldsymbol{\lambda}^{*}\rVert_{\ell_{\infty}}\leq\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}$.
∎

Next, we provide calculation rules for the coupling matrix
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ of the
$\operatorname{\mathsf{GZZ}}(A)$ gate. These rules are inherited from matrix
exponentials. Let
$A_{1},A_{2},A_{3}\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ then

  1. (i)

$\operatorname{\mathsf{GZZ}}(A_{1}+A_{2})=\operatorname{\mathsf{GZZ}}(A_{1})\operatorname{\mathsf{GZZ}}(A_{2})$,

  2. (ii)

$\operatorname{\mathsf{GZZ}}(A_{1}\oplus
A_{2})=\operatorname{\mathsf{GZZ}}(A_{1})\otimes\operatorname{\mathsf{GZZ}}(A_{2})$
and

  3. (iii)

$\operatorname{\mathsf{GZZ}}(A_{1}\circ
A_{2})=\operatorname{\mathsf{GZZ}}(A_{3})$, where the coupling matrices can be
decomposed as

$A_{1}\coloneqq\sum_{k=1}^{d_{1}}\lambda_{\boldsymbol{m}_{k}}\boldsymbol{m}_{k}\boldsymbol{m}_{k}^{T},\quad
A_{2}\coloneqq\sum_{k=1}^{d_{2}}\beta_{\boldsymbol{v}_{k}}\boldsymbol{v}_{k}\boldsymbol{v}_{k}^{T},\quad
A_{3}\coloneqq\sum_{k=1}^{d_{3}}\gamma_{\boldsymbol{w}_{k}}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{T}\,,$

with $d_{3}=d_{1}d_{2}$,
$\boldsymbol{w}_{k}=(\boldsymbol{m}_{i}\circ\boldsymbol{v}_{j})_{k=d_{2}i+j}$
and
$\gamma_{\boldsymbol{w}_{k}}=(\lambda_{\boldsymbol{m}_{i}}\beta_{\boldsymbol{v}_{j}})_{k=d_{2}i+j}$
for $j=1,\dots,d_{2}$ and $i=0,\dots,d_{1}-1$.

In Item (iii) we can also express
$\boldsymbol{\gamma}=\boldsymbol{\lambda}\otimes\boldsymbol{\beta}$ with the
Kronecker product. Note, that these generated $\operatorname{\mathsf{GZZ}}$
gates are not necessarily time-optimal.

We denote the $k\times k$ identity matrix by $\mathds{1}_{k}$, the $k$
dimensional all-ones vector by $\boldsymbol{1}_{k}$ and the $k\times k$ matrix
of ones with vanishing diagonal by
$E_{k}\coloneqq\boldsymbol{1}_{k}\boldsymbol{1}_{k}^{T}-\mathds{1}_{k}$. With
the next two lemmas, we bound the encoding cost and total gate time for
special cases of Item (i) and Item (ii), where the total gate time is
constant. These results will provide a basis for all other constructions. The
following result constructs total coupling matrices
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ on arbitrary subsets of
qubits.

######  Lemma 8.

Let the coupling matrix $J$ be constant, i.e. $J_{ij}=c\in\mathbb{R}_{\geq 0}$
for all $i,j\in[n]$ with $i\neq j$. For any $\varphi\in\mathbb{R}$ and
$s\in\mathbb{Z}_{\geq 0}$ the $\operatorname{\mathsf{GZZ}}(A)$ with the matrix

$A=\varphi\bigoplus_{i=1}^{s}E_{k_{i}}\,,$ (18)

with $k_{i}\geq 1$ for $i=1,\dots,s$ has the encoding cost of
$d=2^{\lceil\log_{2}(s)\rceil}\leq 2s$ and constant total
$\operatorname{\mathsf{GZZ}}$ time $\varphi/c$. This total gate time is
optimal.

###### Proof.

W.l.o.g. we set $c=\varphi=1$. We have $n\coloneqq\sum_{i=1}^{s}k_{i}$ qubits.
Note, that $E_{k_{i}}=E_{1}=0$ only contributes to an entry in the diagonal of
$A$, and hence this qubit does not participate in the
$\operatorname{\mathsf{GZZ}}(A)$ gate. We denote the $d\times d$ Hadamard
matrix by $H^{d\times d}$ and the matrix consisting of its first $s$ columns
by $H^{d\times s}$ where $s\leq n$. The orthogonality of the columns of any
Hadamard matrix yields $(H^{d\times s})^{T}H^{d\times s}=d\,\mathds{1}_{s}$.
Replacing each $i$-th column by $k_{i}$ copies of it we obtain the $d\times
n$-matrix $H^{d\times n}_{\boldsymbol{k}}$ from $H^{d\times s}$. Then, we
attain the diagonal block matrix structure

$(H^{d\times n}_{\boldsymbol{k}})^{T}H^{d\times
n}_{\boldsymbol{k}}-d\mathds{1}_{n}=d\bigoplus_{i=1}^{s}E_{k_{i}}\,.$ (19)

We set the diagonal elements on the left-hand side to zero, since they only
contribute to an energy offset in the Hamiltonian. Take each row of
$H^{d\times n}_{\boldsymbol{k}}$ as a possible vector
$\boldsymbol{m}\in\\{-1,1\\}^{n}$ to construct the total coupling matrix

$A=\frac{1}{d}\sum_{\boldsymbol{m}\in\operatorname{rows}(H_{\boldsymbol{k}}^{d\times
n})}\boldsymbol{m}\boldsymbol{m}^{T}\,,$ (20)

i.e., the time steps have been chosen as
$\lambda_{\boldsymbol{m}}=\frac{1}{d}\,[\boldsymbol{m}\in\operatorname{rows}(H^{d\times
n}_{\boldsymbol{k}})]$, with the Iverson bracket $[\,\cdot\,]$. Clearly, we
have $\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}=1$. If we take the
constants $c$ and $\varphi$ into account, we just multiply Eq. 20 by
$\varphi/c$ which gives a total $\operatorname{\mathsf{GZZ}}$ time of
$\varphi/c$. Furthermore, this total $\operatorname{\mathsf{GZZ}}$ time is
optimal since $\lVert M\rVert_{\ell_{\infty}}=\lVert A\oslash
J\rVert_{\ell_{\infty}}=\varphi/c$ satisfies the lower bound of Lemma 7. ∎

The encoding cost of $\operatorname{\mathsf{GZZ}}(A)$ considered in this lemma
can be reduced if redundant encodings
$\boldsymbol{m}=\boldsymbol{m}^{\prime}\in\\{-1,1\\}^{n}$ are present by
adding the corresponding time steps
$\lambda_{\boldsymbol{m}}+\lambda_{\boldsymbol{m}^{\prime}}$. It can be
further reduced by using the Hadamard conjecture [20, 21]. If the Hadamard
conjecture holds, then there exists Hadamard matrices of any dimension
divisible by $4$. It is known that the Hadamard conjecture holds for
dimensions $s\leq 668$ [22, 23], thus the encoding cost can be reduced to
$d=4\lceil s/4\rceil\leq s-4$ in this regime. The encoding cost of the
following Lemma 9, Theorems 10 and 11 and the efficient heuristic in Section 6
can be reduced in the same fashion.

The assumption in Lemma 8 of a constant coupling matrix, $J_{i\neq j}=c$, is
physically unreasonable. Therefore, we relax this assumption for block sizes
$k_{i}=2$ which corresponds to pairwise next neighbor couplings. We use Lemma
8 to construct time-optimal $\operatorname{\mathsf{GZZ}}$ gates for a certain
family of coupling matrices.

######  Lemma 9.

Let $n$ be even and the coupling matrix $J$ be constant on the first
subdiagonal (the other elements are arbitrary), i.e. $J_{ij}=c$ for
$i\in[n-1]$ and $j=i\pm 1$. Then

$\operatorname{\mathsf{GZZ}}(A)=\bigotimes_{i=1}^{n/2}\operatorname{\mathsf{GZZ}}(\varphi
E_{2})\,,$ (21)

has the encoding cost of $d=2^{\lceil\log_{2}(n)-1\rceil}<n$ and constant
total $\operatorname{\mathsf{GZZ}}$ time $\varphi/c$. This total gate time is
optimal.

###### Proof.

Since the matrix $E_{2}$ only has non-zero entries in the first subdiagonal,
the coupling matrix $J$ only needs to be constant there. The claim follows
immediately from Lemma 8 by setting $k_{1}=\dots=k_{n/2}=2$. Clearly, the
total $\operatorname{\mathsf{GZZ}}$ time is optimal since it saturates the
lower bound of Lemma 7. ∎

Lemma 9 guarantees an encoding cost of $d<n$, which is a quadratic saving
compared to the general LP solution with encoding cost $\binom{n}{2}$. We
note, that

$A=\varphi\bigoplus_{i=1}^{n/2}\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}$ (22)

corresponds to parallel $\operatorname{\mathsf{ZZ}}(\varphi)$ gates, which
find applications in simulating molecular dynamics [24]. The assumption of a
constant subdiagonal of $J$ can be realized in an ion trap platform by
applying an anharmonic trapping potential [16].

By combining (i), Lemma 8 and Lemma 9 we obtain the
$\operatorname{\mathsf{GZZ}}$ gate for next neighbor coupling.

######  Theorem 10.

Let the subdiagonal of $J$ be constant, i.e. $J_{ij}=c$ for $i\in[n-1]$ and
$j=i\pm 1$. Then $\operatorname{\mathsf{GZZ}}(A)$ on $n$ qubits with

$A=\varphi\begin{pmatrix}0&1&&&&&\\\ 1&0&1&&&&\\\ &1&0&1&&&\\\ &&1&0&1&&\\\
&&&&\ddots&&\\\ \end{pmatrix}\,,$ (23)

has the encoding cost of $d\leq 2n$ (for $n>4$) and constant total
$\operatorname{\mathsf{GZZ}}$ time $2\varphi/c$. This total gate time is
optimal.

###### Proof.

We set $c=\varphi=1$ w.l.o.g. For now, we assume that $n$ is even. Then

$\displaystyle A$
$\displaystyle=\bigoplus_{i=1}^{n/2}E_{2}+E_{1}\bigoplus\left(\bigoplus_{i=1}^{n/2-1}E_{2}\right)\bigoplus
E_{1}$ (24) $\displaystyle=\begin{pmatrix}0&1&&&&&\\\ 1&0&&&&&\\\ &&0&1&&&\\\
&&1&0&&&\\\ &&&&\ddots&&\\\ \end{pmatrix}+\begin{pmatrix}0&&&&&&\\\
&0&1&&&&\\\ &1&0&&&&\\\ &&&0&1&&\\\ &&&1&0&&\\\ &&&&&\ddots&\\\ \end{pmatrix}$
$\displaystyle=\begin{pmatrix}0&1&&&&&\\\ 1&0&1&&&&\\\ &1&0&1&&&\\\
&&1&0&1&&\\\ &&&&\ddots&&\\\ \end{pmatrix}$

Again, the diagonal entries do not contribute to the interactions. The first
term can be implemented, using Lemma 9, with the encoding cost
$d_{1}=2^{\lceil\log_{2}(n)-1\rceil}$ and the $\operatorname{\mathsf{GZZ}}$
time of $1$. The second term can be implemented, using Lemma 8, as

$\bigotimes_{i=1}^{s}\operatorname{\mathsf{GZZ}}(\varphi E_{k_{i}})\,,$ (25)

where $s=n/2+1$, $k_{1}=k_{s}=1$ and $k_{i}=2$ for $i=2,\dots,n/2$, with the
encoding cost $d_{2}=2^{\lceil\log_{2}(n+2)-1\rceil}$ and the
$\operatorname{\mathsf{GZZ}}$ time of $1$. Adding the encoding costs
$d=d_{1}+d_{2}$ and $\operatorname{\mathsf{GZZ}}$ times of both terms yields
the desired result. If $n$ is not even, then repeat the previous steps for
$n+1$ but in the end reduce the dimension of all the resulting
$\boldsymbol{m}$ by discarding the last entry.

This construction corresponds to a feasible solution of the primal LP (11)
with the objective function value $2$. To show optimality it suffices to
construct a feasible solution for the dual LP (2) with the objective function
value of $2$. First, consider the case $n=3$ with the total coupling matrix

$A=\begin{pmatrix}0&1&0\\\ 1&0&1\\\ 0&1&0\\\
\end{pmatrix}\Rightarrow\boldsymbol{v}(A)^{T}=(1,0,1)\,,$ (26)

and the matrix representation of the linear operator

$V^{T}=\begin{pmatrix}1&1&1\\\ 1&-1&-1\\\ -1&1&-1\\\ -1&-1&1\\\
\end{pmatrix}\,,$ (27)

as in Eq. 10. A feasible dual solution satisfying
$V^{T}\boldsymbol{y}\leq\boldsymbol{1}$ is $\boldsymbol{y}=(1,-1,1)^{T}$.
Thus, we can verify optimality for $n=3$ since the objective function value is
$\langle A,\boldsymbol{y}\rangle=\boldsymbol{v}(A)^{T}\boldsymbol{y}=2$. Now,
we consider arbitrary $n>3$. Extending the dual solution for the case $n=3$
with zeros $\boldsymbol{y}=(1,-1,1,0,\dots,0)^{T}$ does not change the
objective function value $\langle
A,\boldsymbol{y}\rangle=\boldsymbol{v}(A)^{T}\boldsymbol{y}=2$. Such an
extended dual solution is still feasible since $V^{T}$ restricted to the first
three columns only has rows which are already contained in Eq. 27 due to the
symmetry of Eq. 4. ∎

Algorithm 1 is the pseudocode implementation of Theorem 10. It takes the
number of qubits $n$, the constant value $c$ of the subdiagonal of $J$ and the
factor $\varphi$ of $A$ as input and returns the sparse vector
$\boldsymbol{\lambda}$, containing the time steps. The encodings
$\boldsymbol{m}$ are given by the indices of the non-zero elements
$\lambda_{\boldsymbol{m}}\neq 0$.

Algorithm 1 Synthesize $\operatorname{\mathsf{GZZ}}(A)$ as in Theorem 10.

Input: $n,c,\varphi$

Initialize $\boldsymbol{\lambda}=\boldsymbol{0}\in\mathbb{R}^{2^{n}}$

is_odd $\leftarrow false$

if $n$ odd then

$n\leftarrow n+1$

is_odd $\leftarrow true$

Let $H_{1}^{d_{1}\times d_{1}}$ be a Hadamard matrix $\triangleright$ With
$d_{1}=2^{\lceil\log_{2}(n)-1\rceil}$ as in Lemma 9

$H_{1}^{d_{1}\times\frac{n}{2}}\leftarrow\frac{n}{2}\text{ columns of
}H_{1}^{d_{1}\times d_{1}}$ $\triangleright$ It does not matter which columns

$H_{1}^{d_{1}\times n}\leftarrow\text{duplicate columns
}i=1,\dots,\frac{n}{2}\text{ of }H_{1}^{d_{1}\times\frac{n}{2}}$

Let $H_{2}^{d_{2}\times d_{2}}$ be a Hadamard matrix $\triangleright$ With
$d_{2}=2^{\lceil\log_{2}(n+2)-1\rceil}$

$H_{2}^{d_{2}\times(\frac{n}{2}+1)}\leftarrow\frac{n}{2}+1\text{ columns of
}H_{2}^{d_{2}\times d_{2}}$ $\triangleright$ It does not matter which columns

$H_{2}^{d_{2}\times n}\leftarrow\text{duplicate columns
}i=2,\dots,\frac{n}{2}\text{ of }H_{2}^{d_{2}\times(\frac{n}{2}+1)}$

for $j\in\\{1,2\\}$ do

if is_odd then

Delete one column of $H_{j}^{d_{j}\times n}$ $\triangleright$ It does not
matter which column

for $\boldsymbol{m}\in\operatorname{rows}(H_{j}^{d_{j}\times n})$ do

$\lambda_{\boldsymbol{m}}\leftarrow\frac{\varphi}{d_{j}c}$

Output: $\boldsymbol{\lambda}$ $\triangleright$ Can efficiently be saved in a
sparse data format

The following theorem does not require any additional assumptions on $J$. It
shows, how the LP (11) can be supplemented with the explicit solution to
exclude certain qubits.

######  Theorem 11 (Excluding qubits).

Let $N$ be the total number of qubits on the quantum hardware and $n=N-s$ be
the participating qubits in the $\operatorname{\mathsf{GZZ}}$ gate. Synthesize
the $\operatorname{\mathsf{GZZ}}(A)$ gate with
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$, using the LP (11).
Then, the total encoding cost (on $N$ qubits) is at most
$\binom{n}{2}2^{\lceil\log_{2}(s+1)\rceil}$ and the total
$\operatorname{\mathsf{GZZ}}$ time,
$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}$, (on $N$ qubits) is the same as
for the LP run on $n$ qubits.

###### Proof.

Assume w.l.o.g. that all qubits to be excluded are at the end of the qubit
array. Let $k_{1}=n$ and $k_{2},\dots,k_{s+1}=1$. Using Lemma 8 we obtain an
encoding cost of $d_{1}=2^{\lceil\log_{2}(s+1)\rceil}$ and total
$\operatorname{\mathsf{GZZ}}$ time $1$ to generate the matrix

$A_{1}=\begin{pmatrix}E_{n}&0\\\ 0&0\\\ \end{pmatrix}\,.$ (28)

Solving the LP (11) for a matrix
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ yields the encoding cost
$d_{2}=\binom{n}{2}$ and the total $\operatorname{\mathsf{GZZ}}$ time
$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}$. We define the extension of
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ by
$A_{2}\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{N})$. The extension can
be done by appending $s$ arbitrary elements in $\\{-1,+1\\}$ to all vectors
$\boldsymbol{m}\in\\{-1,+1\\}^{n}$ given by the LP (11). Clearly,

$A_{1}\circ A_{2}=\begin{pmatrix}E_{n}&0\\\ 0&0\\\
\end{pmatrix}\circ\begin{pmatrix}A&*\\\ *&*\\\
\end{pmatrix}=\begin{pmatrix}A&0\\\ 0&0\\\ \end{pmatrix}\,.$ (29)

By (iii), the total encoding cost is $d_{1}d_{2}$ and the total
$\operatorname{\mathsf{GZZ}}$ time is

$\sum_{i=1}^{d_{1}}\sum_{j=1}^{d_{2}}1/d_{1}\lambda_{j}^{*}=\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}\,.$
(30)

Consider now an arbitrary coupling matrix $J$. Then

$J\circ A_{1}\circ A_{2}=J\circ\begin{pmatrix}E_{n}&0\\\ 0&0\\\
\end{pmatrix}\circ\begin{pmatrix}\tilde{A}&*\\\ *&*\\\
\end{pmatrix}=\begin{pmatrix}A&0\\\ 0&0\\\ \end{pmatrix}\,,$ (31)

where $\tilde{A}=A\oslash J$ is decomposed by the LP (11). ∎

Algorithm 2 takes the total number of qubits $N$, the total coupling matrix
$A$ (on $n$ qubits) and the set
$\mathcal{Z}\coloneqq\left\\{i\in[N]\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\text{exclude
qubit }i\right\\}$ as input and returns the sparse vector
$\boldsymbol{\gamma}$, containing the time steps. The encodings
$\boldsymbol{w}$ are given by the indices of the non-zero elements
$\gamma_{\boldsymbol{w}}\neq 0$.

Algorithm 2 Excluding qubits.

Input: $N,A,\mathcal{Z}$ (set of qubit indices to be excluded)

Initialize $\boldsymbol{\gamma}=\boldsymbol{0}\in\mathbb{R}^{2^{N}}$

$s\leftarrow|\mathcal{Z}|$

$\overline{\mathcal{Z}}\leftarrow\left\\{i\in[N]\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}i\notin\mathcal{Z}\right\\}$

Let $H_{1}^{d_{1}\times d_{1}}$ be a Hadamard matrix $\triangleright$ With
$d_{1}=2^{\lceil\log_{2}(s+1)\rceil}$

$H_{1}^{d_{1}\times(s+1)}\leftarrow s+1\text{ columns of }H_{1}^{d_{1}\times
d_{1}}$ $\triangleright$ It does not matter which columns

$H_{1}^{d_{1}\times N}\leftarrow$ duplicate one column of
$H_{1}^{d_{1}\times(s+1)}$ $n-1$ times and place them at indices
$i\in\overline{\mathcal{Z}}$

$\boldsymbol{\lambda}^{*}\leftarrow\text{Solve
}A=\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}\boldsymbol{m}\boldsymbol{m}^{T}$
$\triangleright$ using LP (11)

for $\boldsymbol{v}\in\operatorname{rows}(H_{1}^{d_{1}\times N})$ do

for
$\boldsymbol{m}\in\left\\{\boldsymbol{m}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\lambda_{\boldsymbol{m}}^{*}\neq
0\right\\}$ do $\triangleright$ There are at most $\binom{n}{2}$ such
$\boldsymbol{m}$

$\tilde{\boldsymbol{m}}\leftarrow$ extend $\boldsymbol{m}$ with arbitrary
elements from $\\{-1,+1\\}$ at indices $i\in\mathcal{Z}$.

$\boldsymbol{w}\leftarrow\boldsymbol{v}\circ\tilde{\boldsymbol{m}}$

$\gamma_{\boldsymbol{w}}\leftarrow\frac{\lambda_{\boldsymbol{m}}^{*}}{d_{1}}$

Output: $\boldsymbol{\gamma}$ $\triangleright$ Can efficiently be saved in a
sparse data format

We showed explicit constructions of time-optimal $\operatorname{\mathsf{GZZ}}$
gates for total coupling matrices
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ with diagonal block
structure and next neighbor couplings. The resulting
$\operatorname{\mathsf{GZZ}}$ gates have a constant gate time and require only
linear many encodings to be implemented.

##  6 Efficient heuristic for fast GZZ gates

In this section, we build on the results of Lemma 8 to derive a heuristic
algorithm for synthesizing $\operatorname{\mathsf{GZZ}}(A)$ gates with low
total gate time for any $A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$.
This algorithm runs in polynomial time as opposed to the general LP (11),
which we have shown in Theorem 6 to be $\mathsf{NP}$-hard. The runtime saving
is due to the restriction of the elliptope $\mathcal{E}_{n}$ in Eq. 14, with
exponential many elements, to a set with polynomial many elements. This
restriction yields a polynomial sized LP which can be solved in polynomial
time. In practice, the simplex algorithm has a runtime that scales polynomial
in the problem size [25].

Recall the modified Hadamard matrix $H^{d\times n}_{\boldsymbol{k}}$ defined
in the proof of Lemma 8, where we used the rows of $H^{d\times
n}_{\boldsymbol{k}}$ as encodings $\boldsymbol{m}$ to generate block diagonal
target coupling matrices under some assumptions. Here, $s$ is the number of
block matrices on the diagonal of the target matrix,
$\boldsymbol{k}\in\mathbb{N}^{s}$ contains the dimensions for each block and
$d=2^{\lceil\log_{2}(s)\rceil}$ is the required number of encodings to
construct such a block diagonal matrix. From now on, we only consider
$\boldsymbol{k}=(j,1\dots 1)\in\mathbb{N}^{s}$ such that

$(H^{d\times n}_{\boldsymbol{k}})^{T}H^{d\times
n}_{\boldsymbol{k}}-d\mathds{1}_{n}=d\left(E_{j}\oplus E_{1}\oplus\cdots\oplus
E_{1}\right)\,,$ (32)

see Eq. 19. The requirement that $\sum_{i}k_{i}=n$ implies that such a vector
$\boldsymbol{k}$ has $s=n-j+1$ entries. Permuting the columns of $H^{d\times
n}_{\boldsymbol{k}}$ results in the same permutation of the rows and columns
of the right-hand side of Eq. 32. We denote the set of all column-permuted
$H^{d\times n}_{\boldsymbol{k}}$ by $\mathcal{C}^{(j)}$. A specific element of
$\mathcal{C}^{(j)}$ is denoted by $H^{d\times n}_{\boldsymbol{r}}$, where
$\boldsymbol{r}\in\mathbb{N}^{j}$ is an ordered multi-index
$r_{1}<\dots<r_{j}$ indicating which columns of $H^{d\times
n}_{\boldsymbol{r}}$ are identical. For example, $\boldsymbol{r}=(2,5,6)$
indicates that the columns of $H^{d\times n}_{\boldsymbol{r}}$ with indices
$2$, $5$ and $6$ are identical, i.e. replacing column $5$ and $6$ with column
$2$. This notation will be useful later.

######  Definition 12.

For any $j\in\\{2,3,\dots,n\\}$, we define the restricted elliptope

$\mathcal{E}_{n}^{(j)}\coloneqq\left\\{\boldsymbol{m}\boldsymbol{m}^{T}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\boldsymbol{m}\text{
is a row of }H^{d\times n}_{\boldsymbol{r}}\in\mathcal{C}^{(j)}\right\\}.$
(33)

Further we define

$\mathcal{E}_{n}^{[j]}\coloneqq\bigcup_{i=2}^{j}\mathcal{E}_{n}^{(i)}\,.$ (34)

We choose the definition in Eq. 33 similar as in Eq. 14. Next, we show the
size scaling of the restricted elliptopes. This directly translates to the
size and runtime of the heuristic synthesis optimization.

######  Proposition 13.

For any $j\in\\{2,3,\dots,n\\}$, the number of different encodings
$\boldsymbol{m}$ generating the restricted elliptope $\mathcal{E}_{n}^{[j]}$
scales as $\operatorname{\mathrm{O}}(n^{j+1})$.

###### Proof.

Note, that $|\mathcal{C}^{(j)}|=\binom{n}{j}$ since there are $j$ duplicate
columns in $H^{d\times n}_{\boldsymbol{r}}$. The binomial coefficient can be
bounded by $\binom{n}{j}\leq n^{j}/j!$. Since there are
$d=2^{\lceil\log_{2}(n-j+1)\rceil}<2(n-j+1)$ rows of $H^{d\times
n}_{\boldsymbol{r}}$ we have a rough upper bound of the number different
encodings generating the restricted elliptope,
$\big{\lvert}\mathcal{E}_{n}^{(j)}\big{\rvert}\leq
d\binom{n}{j}<2(n-j+1)\binom{n}{j}<2n^{j+1}/j!$. The first inequality is due
to possible redundant encodings in the definition of $\mathcal{E}_{n}^{(j)}$.
Similarly, we can upper bound
$\big{\lvert}\mathcal{E}_{n}^{[j]}\big{\rvert}\leq\sum_{i=2}^{j}\big{\lvert}\mathcal{E}_{n}^{(i)}\big{\rvert}<2\sum_{i=2}^{j}n^{i+1}/i!$
which is a polynomial of order $j+1$. ∎

We denote the convex cone generated by a set $V$ by

$\operatorname{cone}(V)\coloneqq\left\\{\sum_{i}\lambda_{i}v_{i}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\lambda_{i}\geq
0\text{, }v_{i}\in V\right\\}\,.$ (35)

With that, we are ready to present the main result of this section.

######  Theorem 14.

$\operatorname{cone}(\mathcal{E}_{n}^{(2)})=\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$.

###### Proof.

W.l.o.g. we can assume $d=n$ and denote $H^{d\times
d}_{\boldsymbol{r}}\in\mathcal{C}^{(2)}$ by $H^{d}_{(r_{1},r_{2})}$ with the
property

$(H^{d}_{(r_{1},r_{2})})^{T}H^{d}_{(r_{1},r_{2})}-d\mathds{1}_{d}=d\,e_{(r_{1},r_{2})}\,,$
(36)

where $e_{(r_{1},r_{2})}$ is an element of the standard basis for symmetric
matrices with vanishing diagonal. By Eq. 36 it holds
$\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}_{\geq
0}^{n})\subseteq\operatorname{cone}(\mathcal{E}_{n}^{(2)})$, i.e., symmetric
matrices with non-negative entries are in the convex cone.

It is left to show that
$\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}_{<0}^{n})\subseteq\operatorname{cone}(\mathcal{E}_{n}^{(2)})$,
i.e., that also symmetric matrices with negative entries are in the convex
cone. To show this inclusion we define $H^{d}_{(r_{1},-r_{2})}$ similar as
$H^{d}_{(r_{1},r_{2})}$ except the duplicate column at $r_{2}$ is multiplied
by $-1$ such that

$(H^{d}_{(r_{1},-r_{2})})^{T}H^{d}_{(r_{1},-r_{2})}-d\mathds{1}_{d}=-d\,e_{(r_{1},r_{2})}\,.$
(37)

We have to show that for each row
$\boldsymbol{m}\in\operatorname{rows}(H^{d}_{(r_{1},-r_{2})})$ there exist
$\tilde{r}_{1}$ and $\tilde{r}_{2}$ such that
$\boldsymbol{m}\in\operatorname{rows}(H^{d}_{(\tilde{r}_{1},\tilde{r}_{2})})$.
This can be verified straightforwardly for $d=4$ by checking all rows.
W.l.o.g. we show that the hypothesis holds for any $H^{2d}_{(r_{1},-r_{2})}$
by assuming it holds for $H^{d}_{(r_{1},-r_{2})}$. The Sylvester-Hadamard
matrix is constructed inductively according to

$H^{2d}=\begin{pmatrix}H^{d}&H^{d}\\\ H^{d}&-H^{d}\\\ \end{pmatrix}.$ (38)

We consider three cases for $H^{2d}_{(r_{1},-r_{2})}$.

#### Case 1.

For a $H^{2d}_{(r_{1},-r_{2})}$ with identical columns, up to minus sign, at
indices $r_{1},r_{2}\in[d]$ or $r_{1},r_{2}\in[d+1,2d]$ the hypothesis holds
by our assumption by choosing $\tilde{r}_{1},\tilde{r}_{2}\in[d]$ or
$\tilde{r}_{1},\tilde{r}_{2}\in[d+1,2d]$ respectively.

#### Case 2.

Considering the first $d$ rows of $H^{2d}_{(r_{1},-r_{2})}$ with identical
columns, up to minus sign, at indices $r_{1}\in[d]$ and $r_{2}\in[d+1,2d]$.
This case is equivalent to Case 1. with $r_{1},r_{2}\in[d+1,2d]$ since only
the column at $r_{2}$ is negated.

#### Case 3.

Considering the last $d$ rows of $H^{2d}_{(r_{1},-r_{2})}$ with identical
columns, up to minus sign, at indices $r_{1}\in[d]$ and $r_{2}\in[d+1,2d]$.
These rows are included in the last $d$ rows of
$H^{2d}_{(\tilde{r}_{1},\tilde{r}_{2})}$ with
$\tilde{r}_{1},\tilde{r}_{2}\in[d+1,2d]$ and $\tilde{r}_{2}=r_{2}$ since the
column $r_{2}$ is negated which is equivalent to just duplicating a column of
$-H^{d}$.

We have shown that for each row
$\boldsymbol{m}\in\operatorname{rows}(H^{d}_{(r_{1},-r_{2})})$ there exist
$\tilde{r}_{1}$ and $\tilde{r}_{2}$ such that
$\boldsymbol{m}\in\operatorname{rows}(H^{d}_{(\tilde{r}_{1},\tilde{r}_{2})})$.
Therefore,
$\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}_{<0}^{n})\cup\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}_{\geq
0}^{n})=\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})=\operatorname{cone}(\mathcal{E}_{n}^{(2)})$.
The last equality follows from the definition of $\operatorname{cone}$ and
$\mathcal{E}_{n}^{(2)}$. ∎ Theorem 14 shows that the constraint of LP (11) can
always be fulfilled only considering
$\boldsymbol{m}\boldsymbol{m}^{T}\in\mathcal{E}_{n}^{(2)}$. Similar to Eq. 10
we define the restricted linear operator $\mathcal{V}^{[j]}:\mathbb{R}_{\geq
0}^{s}\rightarrow\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n}):\boldsymbol{\lambda}\mapsto\sum_{\boldsymbol{m}}\lambda_{\boldsymbol{m}}\boldsymbol{m}\boldsymbol{m}^{T}$
for all $\boldsymbol{m}\boldsymbol{m}^{T}\in\mathcal{E}_{n}^{[j]}$ with
$h=|\mathcal{E}_{n}^{[j]}|$, represented by a matrix
$V^{[j]}\in\\{-1,+1\\}^{(\binom{n}{2})\times h}$. We define the _restricted
LPj_ to be 1^T λ V^[j] λ =v(M) λ ∈R^h_≥0 . Algorithm 3 summarizes the steps to
construct $\mathcal{E}_{n}^{[j]}$ and therefore the matrix representation of
the restricted linear operator $V^{[j]}$. In practice, Algorithm 3 has to be
executed only once per number of qubits $n$ since the constraints of LP (6)
can be fulfilled for all
$M\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$. This is due to the fact
that $\mathcal{E}_{n}^{(2)}\subseteq\mathcal{E}_{n}^{[j]}$ for any $2\leq
j\leq n$ and
$\operatorname{cone}(\mathcal{E}_{n}^{(2)})=\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$
(Theorem 14). The time and space complexity of Algorithm 3 scales polynomially
in $n$ as shown in Proposition 13. Therefore, the restricted LPj is also of
polynomial size for a fixed $j$. Increasing $j$ leads to better approximations
due to the enlarged search space for the optimal solution. Note, that the
runtime of the mixed integer program (MIP) defined in [6, Section 2.2.2] also
benefits from using $\mathcal{E}_{n}^{[j]}$.

As mentioned in Section 5 the dimension of the Hadamard matrices in Eq. 32 can
be reduced to $d=4\lceil(n-1)/4\rceil\leq n-5$ if we assume that the Hadamard
conjecture holds. Therefore, the runtime of the restricted LPj is reduced as
well if such Hadamard matrices are used. In Section 8 we numerically benchmark
the heuristic algorithm with and without the reduced runtime of the restricted
LPj.

Algorithm 3 Constructing $\mathcal{E}_{n}^{[j]}$.

Input: $n,j$

Initialize $\mathcal{E}_{n}^{[j]}=\emptyset$

for $i=2,\dots,j$ do

Initialize $\mathcal{E}_{n}^{(i)}=\emptyset$

Let $H^{d\times d}$ be a Hadamard matrix $\triangleright$ With
$d=2^{\lceil\log_{2}(n-i+1)\rceil}$

$H^{d\times(n-i+1)}\leftarrow n-i+1\text{ columns of }H^{d\times d}$
$\triangleright$ It does not matter which columns

for all $\boldsymbol{r}$ such that $1\leq r_{1}<\dots<r_{i}\leq n$ do
$\triangleright$ There are $\binom{n}{i}$ such $\boldsymbol{r}$

$H^{d\times n}_{\boldsymbol{r}}\leftarrow$ duplicate one column of
$H^{d\times(n-i+1)}$ $i-1$ times to indices $r_{1},\dots,r_{i}$

$\mathcal{E}_{n}^{(i)}\leftarrow\mathcal{E}_{n}^{(i)}\cup\left\\{\boldsymbol{m}\boldsymbol{m}^{T}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\boldsymbol{m}\text{
is a row of }H^{d\times n}_{\boldsymbol{r}}\right\\}$

$\mathcal{E}_{n}^{[j]}\leftarrow\mathcal{E}_{n}^{[j]}\cup\mathcal{E}_{n}^{(i)}$

Output: $\mathcal{E}_{n}^{[j]}$

##  7 Bounds on the total GZZ time

Our main analytic result (Theorem 15) is that the optimal total
$\operatorname{\mathsf{GZZ}}$ time
$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}$ is lower and upper bounded by the
norms $\lVert M\rVert_{\ell_{\infty}}$ and $\lVert M\rVert_{\ell_{1}}$,
respectively. Note, that for a dense matrix $M$, its norm $\lVert
M\rVert_{\ell_{1}}$ scales quadratic with the number of qubits $n$. We
conjecture an improved upper bound on the total $\operatorname{\mathsf{GZZ}}$
time for dense $M$ that scales at most linear with the number of qubits $n$.
We support this conjecture with explicit solutions for the LP (11) reaching
this bound for any $n$ and numerical results validating the conjecture for
$n\leq 8$.

######  Theorem 15.

The optimal total gate time of $\operatorname{\mathsf{GZZ}}(A)$ with
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ is lower and upper
bounded by

$\lVert
M\rVert_{\ell_{\infty}}\leq\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}\leq\lVert
M\rVert_{\ell_{1}}\,,$ (39)

where $M\coloneqq A\oslash J$. Equality holds for the lower bound for all
matrices $M=C\boldsymbol{m}\boldsymbol{m}^{T}$ for any
$\boldsymbol{m}\in\\{-1,+1\\}^{n}$ and $C\geq 0$.

###### Proof.

The lower bound has been shown in Lemma 7. Equality in the lower bound holds
for $M=C\boldsymbol{m}\boldsymbol{m}^{T}$ by setting
$\lambda_{\boldsymbol{m}}=C=\lVert M\rVert_{\ell_{\infty}}$ and
$\lambda_{\boldsymbol{m}^{\prime}}=0$ for all
$\boldsymbol{m}^{\prime}\neq\boldsymbol{m}$. We use the explicit construction
of the standard basis elements for symmetric matrices from the proof of
Theorem 14 to show the upper bound. To be precise, we have

$\frac{|M_{ij}|}{d}(H^{d\times n}_{\boldsymbol{r}})^{T}H^{d\times
n}_{\boldsymbol{r}}-d\mathds{1}_{n}=M_{ij}e_{(i,j)}\,,$ (40)

where $\boldsymbol{r}=(i,j)$ if $M_{ij}\geq 0$ or
$\boldsymbol{r}=(\tilde{i},\tilde{j})$ if $M_{ij}<0$ as in Eqs. 36 and 37,
respectively. We define $\boldsymbol{\lambda}^{(i,j)}$ with the entries
$\lambda^{(i,j)}_{\boldsymbol{m}}\coloneqq\frac{|M_{ij}|}{d}\,[\boldsymbol{m}\in\operatorname{rows}(H^{d}_{\boldsymbol{r}})]$.
According to Lemma 8 we have
$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{(i,j)}=|M_{ij}|$. Adding
$\boldsymbol{\lambda}^{(i,j)}$ for all $i<j$ yields the upper bound $\lVert
M\rVert_{\ell_{1}}$. ∎

These bounds get tighter the sparser $M$ is. If $M$ has only one non-zero
value, then clearly $\lVert
M\rVert_{\ell_{\infty}}=\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}=\lVert
M\rVert_{\ell_{1}}$. Furthermore, these bounds also hold for the heuristic,
which we presented in Section 6.

Next, we state our conjecture that the optimal gate time of
$\operatorname{\mathsf{GZZ}}(A)$ scales at most linear with the number of
qubits.

######  Conjecture 16.

The optimal gate time of $\operatorname{\mathsf{GZZ}}(A)$ with
$A\in\operatorname{\mathrm{Sym}_{0}}(\mathbb{R}^{n})$ is tightly upper bounded
by

$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}\leq\lVert
M\rVert_{\ell_{\infty}}\cdot\begin{cases}n\,,&\text{for odd }n\,,\\\
n-1\,,&\text{for even }n\,.\end{cases}$ (41)

Hence, it provides a tighter bound for dense $M$ than Theorem 15. To support
the claim of 16 we first construct explicit dual and primal feasible solutions
for the case $M=-E_{n}$ for the LP’s (11) and (2), respectively. Then
optimality is given by showing equality of the objective function values of
the primal and dual problem. We further show that the case $M=-E_{n}$ leads to
the same objective function value as $M=-\boldsymbol{m}\boldsymbol{m}^{T}$ for
any $\boldsymbol{m}\in\\{-1,1\\}^{n}$. Finally, we provide numerical evidence
that the conjecture holds for $n\leq 8$.

For practical purposes it is important to keep in mind that the platform given
$J$ matrix might also scale with the number of qubits resulting in a qubit
dependent constant $\lVert A\oslash J_{n}\rVert_{\ell_{\infty}}=\lVert
M_{n}\rVert_{\ell_{\infty}}$.

###  7.1 Explicit solutions for $M=-\boldsymbol{m}\boldsymbol{m}^{T}$

The following lemma will be used in the proof of the explicit feasible
solution of the dual problem for $M$ being of the form
$M=-\boldsymbol{m}\boldsymbol{m}^{T}$. We can identify
$\boldsymbol{m}\in\\{-1,1\\}^{n}$ with $\boldsymbol{b}\in\mathbb{F}_{2}^{n}$
via $\boldsymbol{m}=(-1)^{\boldsymbol{b}}$ as explained in Section 2.

######  Lemma 17.

It holds that

$\displaystyle P_{\boldsymbol{b}}\coloneqq\sum_{i<j}^{n}(-1)^{b_{i}\oplus
b_{j}}=\binom{n}{2}-2|\boldsymbol{b}|(n-|\boldsymbol{b}|)\,,$ (42)

for any binary vector $\boldsymbol{b}\in\mathbb{F}_{2}^{n}$. We denote the
Hamming weight by $|\boldsymbol{b}|$.

###### Proof.

Let $\boldsymbol{m}=(-1)^{\boldsymbol{b}}$. If the Hamming weight
$|\boldsymbol{b}|$ vanishes, then $\boldsymbol{m}=(+1,\dots,+1)$ and
$P_{\boldsymbol{b}}=\binom{n}{2}$ which is the maximal value. If
$|\boldsymbol{b}|\neq 0$, $\boldsymbol{m}$ contains $|\boldsymbol{b}|$ entries
$-1$, such that the upper triangular part of
$\boldsymbol{m}\boldsymbol{m}^{T}$ contains a rectangle of $-1$’s with length
$|\boldsymbol{b}|$ and width $n-|\boldsymbol{b}|$ so the total amount of
$-1$’s is $|\boldsymbol{b}|(n-|\boldsymbol{b}|)$. Therefore,

$P_{\boldsymbol{b}}=\binom{n}{2}-2|\boldsymbol{b}|(n-|\boldsymbol{b}|)\,.$
(43)

∎

######  Lemma 18 (explicit dual feasible solution).

Let $M=-E_{n}$, then there is an explicit feasible solution $\boldsymbol{y}$
to the dual LP (2) with

$\langle-E_{n},\boldsymbol{y}\rangle=\begin{cases}n\,,&\text{for odd }n\,,\\\
n-1\,,&\text{for even }n\,.\end{cases}$ (44)

###### Proof.

We assume that $y=y_{1}=y_{2}=\dots$. Therefore, it suffices to show that

$y\sum_{i<j}^{n}(-1)^{b_{i}\oplus b_{j}}\leq 1\,.$ (45)

From Lemma 17 we know that
$\boldsymbol{y}=1/\operatorname{min}(P_{\boldsymbol{b}})\boldsymbol{1}$
satisfies the constraint in Eq. 45. The minimum
$\operatorname{min}(P_{\boldsymbol{b}})=-\lfloor n/2\rfloor$ is reached for
$|\boldsymbol{b}|=\lceil n/2\rceil$ or $|\boldsymbol{b}|=\lfloor n/2\rfloor$
which can be verified from the expression of $P_{\boldsymbol{b}}$ in Lemma 17.
Thus, we obtain $\boldsymbol{y}=-1/\lfloor n/2\rfloor\boldsymbol{1}$. The
objective function evaluates to

$\langle-
E_{n},\boldsymbol{y}\rangle=\boldsymbol{-}\boldsymbol{1}^{T}\boldsymbol{y}=\frac{\binom{n}{2}}{\lfloor
n/2\rfloor}=\begin{cases}n\,,&\text{for odd }n\,,\\\ n-1\,,&\text{for even
}n\,.\end{cases}$ (46)

∎

For the construction of a feasible solution to the primal problem we first
require the following result.

######  Lemma 19.

Let $k<n$ be natural numbers. Then

$\sum_{\boldsymbol{b}\in\mathbb{F}_{2}^{n}:|\boldsymbol{b}|=k}|b_{i}\oplus
b_{j}|=2\binom{n-2}{k-1}$ (47)

for all $i,j\in[n]$ with $i\neq j$.

###### Proof.

Consider the case $n=4$ and $k=2$. We get $\binom{4}{2}=6$ binary vectors with
$|\boldsymbol{b}|=2$. It can be easily verified that
$\sum_{|\boldsymbol{b}|=2}|b_{i}\oplus b_{j}|=4=2\binom{2}{1}$. Now, we assume
that for a given $n$ and $k<n$ the Eq. 47 holds. It suffices to verify Eq. 47
for $i\leq n$ and $j=n+1$. We fix $k$, define $n^{\prime}\coloneqq n+1$ and
take a $\boldsymbol{b}\in\mathbb{F}_{2}^{n^{\prime}}$. We have
$\binom{n+1}{k}$ binary vectors with $|\boldsymbol{b}|=k$.

For the case $b_{n+1}=0$ we have $\binom{n}{k}$ such vectors and

$\sum_{\begin{subarray}{c}|\boldsymbol{b}|=k\\\
b_{n+1}=0\end{subarray}}|b_{i}\oplus
b_{n+1}|=\sum_{|\boldsymbol{b}|=k}|b_{i}|=\binom{n-1}{k-1}\,,$ (48)

for all $i\leq n$. There are $\binom{n-1}{k-1}$ different ways in placing
$k-1$ $1$’s.

For the case $b_{n+1}=1$ we have $\binom{n}{k-1}$ such vectors and

$\sum_{\begin{subarray}{c}|\boldsymbol{b}|=k\\\
b_{n+1}=1\end{subarray}}|b_{i}\oplus
b_{n+1}|=\sum_{|\boldsymbol{b}|=k}|b_{i}\oplus 1|=\binom{n-1}{k-1}\,,$ (49)

for all $i\leq n$. There are $\binom{n-1}{k-1}$ different ways in placing
$k-1$ $0$’s.

Combining the two cases $b_{n+1}=0$ and $b_{n+1}=1$ we obtain

$\sum_{|\boldsymbol{b}|=k}|b_{i}\oplus 0|+|b_{i}\oplus
1|=2\binom{n-1}{k-1}=2\binom{n^{\prime}-2}{k-1}\,,$ (50)

for all $i\leq n$. ∎

We motivate the next lemma with the result of the explicit dual feasible
solution for $M=-E_{n}$ from Lemma 18 and the complementary slackness
condition. The complementary slackness condition for a linear program states
that if the $i$-th inequality of the dual problem is a strict inequality for a
feasible solution $\boldsymbol{y}$, then the $i$-th component of a feasible
solution of the primal problem $\boldsymbol{\lambda}$ is zero:

$(V^{T}\boldsymbol{y})_{i}<1\Rightarrow\lambda_{i}=0\,.$ (51)

We use this in the following lemma to construct a feasible solution for the
primal LP (11).

######  Lemma 20 (explicit primal feasible solution).

Let $M=-E_{n}$, then there is an explicit feasible solution
$\boldsymbol{\lambda}$ to the primal LP (11) with

$\boldsymbol{1}^{T}\boldsymbol{\lambda}=\begin{cases}n\,,&\text{for odd
}n\,,\\\ n-1\,,&\text{for even }n\,.\end{cases}$ (52)

###### Proof.

For this proof we define

$k\coloneqq\begin{cases}n\,,&\text{for odd }n\,,\\\ n-1\,,&\text{for even
}n\,.\end{cases}$ (53)

We only consider binary vectors of the set
$S\coloneqq\left\\{\boldsymbol{b}\in\mathbb{F}_{2}^{n}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}|\boldsymbol{b}|=\lceil
n/2\rceil,\lfloor n/2\rfloor\text{ and }b_{n}=0\right\\}$. This set is
motivated by the complementary slackness condition and Lemma 18. It can be
calculated that $|S|=\binom{k}{\lfloor k/2\rfloor}$ using the recurrence
relation

$\binom{n}{k}=\binom{n-1}{k}+\binom{n-1}{k-1}$ (54)

for the binomial coefficients. We show that

$\sum_{\boldsymbol{b}\in S}(-1)^{b_{i}\oplus b_{j}}=-D_{n}\,,$ (55)

for a constant $D_{n}>0$, which we calculate later. If not explicitly stated,
all equations in this proof containing $i$, $j$ hold for all $i,j\in[n]$,
$i\neq j$. We denote $\boldsymbol{\lambda}_{S}$ by all
$\lambda_{\boldsymbol{m}}=\lambda_{\boldsymbol{b}}$ corresponding to the
encoding $\boldsymbol{m}=(-1)^{\boldsymbol{b}}$ with $\boldsymbol{b}\in S$. If
Eq. 55 holds, we can choose a $\boldsymbol{\lambda}_{S}=1/D_{n}\boldsymbol{1}$
resulting in

$\sum_{\boldsymbol{b}\in
S}\boldsymbol{\lambda}_{\boldsymbol{b}}(-1)^{b_{i}\oplus b_{j}}=-1\,,$ (56)

which implies feasibility for $M=-E_{n}$. It is left to show Eq. 55 and
determine $D_{n}$.

First, we consider odd $n$. By definition of $S$ we have that $b_{n}=0$ for
all $\boldsymbol{b}$. Therefore, we obtain $\binom{n-1}{\lfloor
n/2\rfloor}+\binom{n-1}{\lceil n/2\rceil}$ binary vectors with
$|\boldsymbol{b}|=\lceil n/2\rceil$ or $|\boldsymbol{b}|=\lfloor n/2\rfloor$
respectively. Counting the occurrences of “$-1$” in the sum of Eq. 55 is
equivalent to counting the occurrences of “$1$” in the sum

$\displaystyle\sum_{\boldsymbol{b}\in S}|b_{i}\oplus b_{j}|$
$\displaystyle=2\left(\binom{n-3}{\lfloor\frac{n}{2}\rfloor-1}+\binom{n-3}{\lceil\frac{n}{2}\rceil-1}\right)$
(57) $\displaystyle=2\binom{n-2}{\lfloor\frac{n}{2}\rfloor}\,,$

where we used Lemma 19, the recurrence relation for the binomial coefficients
and $\lceil n/2\rceil-1=\lfloor n/2\rfloor$. Counting the occurrences of
“$+1$” in the sum of Eq. 55 yields

$\binom{n-1}{\lfloor\frac{n}{2}\rfloor}+\binom{n-1}{\lceil\frac{n}{2}\rceil}-2\binom{n-2}{\lfloor\frac{n}{2}\rfloor}=\binom{n}{\lfloor\frac{n}{2}\rfloor}-2\binom{n-2}{\lfloor\frac{n}{2}\rfloor}\,.$
(58)

where we used $\binom{n}{\lfloor n/2\rfloor+1}=\binom{n}{\lfloor n/2\rfloor}$
for odd $n$. We now evaluate Eq. 55 for odd $n$

$\displaystyle-D_{n}=\sum_{\boldsymbol{b}\in S}(-1)^{b_{i}\oplus b_{j}}$
$\displaystyle=\binom{n}{\lfloor\frac{n}{2}\rfloor}-4\binom{n-2}{\lfloor\frac{n}{2}\rfloor}$
(59)
$\displaystyle=\binom{n}{\lfloor\frac{n}{2}\rfloor}-\frac{n+1}{n}\binom{n}{\lfloor\frac{n}{2}\rfloor}$
$\displaystyle=\frac{-1}{n}\binom{n}{\lfloor\frac{n}{2}\rfloor}\,.$

The case for even $n$ follows the same steps as for odd $n$, resulting in

$\displaystyle-D_{n}=\sum_{|\boldsymbol{b}|=\frac{n}{2}}(-1)^{b_{i}\oplus
b_{j}}$
$\displaystyle=\binom{n-1}{\frac{n}{2}}-\frac{n}{n-1}\binom{n-1}{\frac{n}{2}}$
(60) $\displaystyle=\frac{-1}{n-1}\binom{n-1}{\frac{n}{2}}\,.$

Equations 59 and 60 show that $D_{n}=1/k\binom{k}{\lfloor k/2\rfloor}$ and
that Eq. 55 holds. Since $|S|=\binom{k}{\lfloor k/2\rfloor}$ the objective
function value is

$\boldsymbol{1}^{T}\boldsymbol{\lambda}_{S}=\binom{k}{\lfloor
k/2\rfloor}\frac{1}{\frac{1}{k}\binom{k}{\lfloor k/2\rfloor}}=k\,.$ (61)

∎

From the equality of the objective functions for the primal and dual problem
from Lemma 20 and Lemma 18 respectively we know that the proposed dual/primal
feasible solutions for $M=-E_{n}$ are in fact optimal solutions. Now, we show
that the $\operatorname{\mathsf{GZZ}}$ gate time with
$M=-\boldsymbol{m}\boldsymbol{m}^{T}$ for any
$\boldsymbol{m}\in\\{-1,+1\\}^{n}$ is the same as for the case $M=-E_{n}$.

######  Theorem 21.

If $A\oslash J\eqqcolon M=-C(\boldsymbol{m}\boldsymbol{m}^{T})$ for any
$\boldsymbol{m}\in\\{-1,+1\\}^{n}$ and $C\geq 0$, then the optimal gate time
of $\operatorname{\mathsf{GZZ}}(A)$ is

$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}=\lVert
M\rVert_{\ell_{\infty}}\begin{cases}n\,,&\text{for odd }n\,,\\\
n-1\,,&\text{for even }n\,.\end{cases}$ (62)

###### Proof.

The statement has been shown for the case $M=-E_{n}$ by constructing an
explicit solution. It is left to show that the cases
$M=-C(\boldsymbol{m}\boldsymbol{m}^{T})$ for a constant $C\geq 0$ yield the
same objective function value. Since the objective function and the
constraints are linear we can w.l.o.g. assume $C=1$. It is clear, that
$\operatorname{sign}(M_{ij})=\operatorname{sign}(-m_{i}m_{j})=-(-1)^{c_{i}\oplus
c_{j}}$ for all $i<j$, with $\boldsymbol{m}=(-1)^{\boldsymbol{c}}$. Let
$\boldsymbol{y}=-y\operatorname{sign}(\boldsymbol{v}(M))$ for a
$y\in\mathbb{R}$, then each constraint of
$V^{T}\boldsymbol{y}\leq\boldsymbol{1}$ of the dual LP (2) reads as

$y\sum_{i<j}(-1)^{b_{i}\oplus b_{j}\oplus c_{i}\oplus
c_{j}}=y\sum_{i<j}(-1)^{\tilde{b}_{i}\oplus\tilde{b}_{j}}\leq 1\,,$ (63)

with $\boldsymbol{\tilde{b}}\coloneqq\boldsymbol{b}\oplus\boldsymbol{c}$
element wise. Consider the ordered set of all
$\boldsymbol{b}\in\mathbb{F}_{2}^{n}$ with $b_{n}=0$, then
$\boldsymbol{\tilde{b}}$ is just a permutation of that set. Due to the
permutation symmetry of the qubits the optimal value of the LP (11) for any
$M=-(\boldsymbol{m}\boldsymbol{m}^{T})$ is the same as for the case
$M=-E_{n}$. Setting $\boldsymbol{y}=-y\operatorname{sign}(\boldsymbol{v}(M))$
with $y=1/\lfloor n/2\rfloor$ as in the proof of Lemma 18 yields an optimal
solution to the dual LP (2) for the case
$M=-(\boldsymbol{m}\boldsymbol{m}^{T})$. ∎

Note that, trivially, the lower bound
$\boldsymbol{1}^{T}\boldsymbol{\lambda}^{*}=\lVert M\rVert_{\ell_{\infty}}$ is
reached if $M=C(\boldsymbol{m}\boldsymbol{m}^{T})$ for any
$\boldsymbol{m}\in\\{-1,+1\\}^{n}$ and $C\geq 0$.

One possibility to prove 16 is to show, that the matrix $M=-E_{n}$ maximizes
the value of the LP (11) among all matrices
$M\in\operatorname{\mathrm{Sym}_{0}}([-1,+1]^{n})$. To this end, consider the
LP ∥ y ∥_ℓ_1 V^T y ≤1 y∈R^(n2) , which is independent of
$M\in\operatorname{\mathrm{Sym}_{0}}([-1,+1]^{n})$. It holds that

$\max_{V^{T}\boldsymbol{y}\leq\boldsymbol{1}}\left(\max_{\lVert
M\rVert_{\ell_{\infty}}\leq 1}\langle
M,\boldsymbol{y}\rangle\right)=\max_{V^{T}\boldsymbol{y}\leq\boldsymbol{1}}\lVert\boldsymbol{y}\rVert_{\ell_{1}}$
(64)

according to
$\lVert\boldsymbol{x}\rVert_{\ell_{1}}=\max_{\lVert\boldsymbol{p}\rVert_{\ell_{\infty}}\leq
1}\boldsymbol{p}^{T}\boldsymbol{x}$. Therefore, the optimal objective value of
LP (21) is an upper bound on all optimal objective values of LP (11). Clearly,
the constructed solution in Lemma 18 is feasible for LP (21). Unfortunately,
proving that this constructed solution is optimal is quite challenging, as we
discuss in Appendix A.

Figure 1:  The optimal objective function value of the dual LP (2) over the
number of qubits. Blue: The range of optimal values for all binary
$M\in\operatorname{\mathrm{Sym}_{0}}(\\{-1,+1\\}^{n})$. Orange: The range of
optimal values for all binary
$M\in\operatorname{\mathrm{Sym}_{0}}(\\{-1,+1\\}^{n})$ without
$M=\pm(\boldsymbol{m}\boldsymbol{m}^{T})$ for any
$\boldsymbol{m}\in\\{-1,+1\\}^{n}$.

##  8 Numerical results

###  8.1 Numeric validation of Conjecture 16 for small $n$

For the numeric validation of the conjecture for small $n$ we solve the dual
LP (2) for all binary $M\in\operatorname{\mathrm{Sym}_{0}}(\\{-1,+1\\}^{n})$
of which there are $2^{\binom{n}{2}}$. For $n=3$ there are only binary
matrices of the form $M=\pm(\boldsymbol{m}\boldsymbol{m}^{T})$ and by Theorem
21 the conjecture holds. Figure 1 shows that 16 holds for $n\leq 8$. For odd
$n\leq 8$ the cases $M=-(\boldsymbol{m}\boldsymbol{m}^{T})$ are in fact the
only cases reaching the upper bound. This can be seen in Fig. 1 by the blue
area exceeding the orange area, which only consists of the optimal values for
all binary matrices without $M=\pm(\boldsymbol{m}\boldsymbol{m}^{T})$ for any
$\boldsymbol{m}\in\\{-1,+1\\}^{n}$.

###  8.2 Numerical benchmark for the heuristic

Figure 2:  Comparing the performance of the original (optimal) LP (11) and the
restricted LPj (6) for $j=2,3,4$. For each line we let the LP’s run for a
fixed time. The black dashed line is the upper bound for the original LP (11)
from 16. The reddish lines show the total $\operatorname{\mathsf{GZZ}}$ times
for the restricted LPj (6) for $j=2,3,4$. The blueish lines show the total
$\operatorname{\mathsf{GZZ}}$ times for the restricted LPj (6) using Hadamard
matrices of dimension $d=4\,k$ to construct the restricted linear operator
$V^{[j]}$. Left: Average case scaling of the total
$\operatorname{\mathsf{GZZ}}$ times. The data points and error bars show the
mean and the standard deviation over $100$ uniformly sampled matrices
$A_{ij}=A_{ji}\in[-1,1]$ for $i<j$. Right: Conjectured worst-case scaling of
the total $\operatorname{\mathsf{GZZ}}$ times. The data points and error bars
show the mean and the standard deviation over $100$ binary matrices
$A=-\boldsymbol{m}\boldsymbol{m}^{T}$ with uniformly sampled encodings
$\boldsymbol{m}\in\\{-1,1\\}^{n}$.

We compare the performance of the restricted LPj (6) to the original LP (11).
To this end we provide numerical results on the average-case and the
conjectured worst-case scaling of the total $\operatorname{\mathsf{GZZ}}$
time. The left plot in Fig. 2 shows the average-case scaling of the total
$\operatorname{\mathsf{GZZ}}$ time. $\operatorname{\mathsf{GZZ}}(A)$ gates
with uniformly sampled matrix elements $A_{ij}=A_{ji}\in[-1,1]$ for all $i<j$
and $i,j\in[n]$ are synthesized. For the worst-case scaling of the right plot
in Fig. 2, $\operatorname{\mathsf{GZZ}}(A)$ gates with binary matrices
$A=-\boldsymbol{m}\boldsymbol{m}^{T}$ with uniformly sampled encodings
$\boldsymbol{m}\in\\{-1,1\\}^{n}$ are synthesized. We note that these are the
matrices with the conjectured worst-case scaling for the LP and that the LPj
might have different worst-case matrices $A$. For convenience, we have set
$M=A$, i.e. setting $J=E_{n}$ which omits the hardware specific time units
given by the quantum platform. For realistic $J$’s and
$\operatorname{\mathsf{GZZ}}(A)$ gate times we refer to [6]. The Python
package CVXPY [26, 27] with the GNU linear program kit simplex solver [28] is
used to solve the LP (11) and the restricted LPj (6).

Figure 2 shows the total $\operatorname{\mathsf{GZZ}}$ time over the number of
qubits $n$ for the restricted LPj (6) and LP (11). We assigned a fixed runtime
($20$ minutes) for each LP to synthesize $\operatorname{\mathsf{GZZ}}(A)$
gates for all sample matrices $A$ and as many $n$ as possible. Clearly, the
runtime of the heuristic algorithm is much shorter than the runtime of the
optimal LP (11). Increasing the hierarchy of the heuristic $j>2$ reduces the
total $\operatorname{\mathsf{GZZ}}$ time significantly while still maintaining
a short runtime. The total $\operatorname{\mathsf{GZZ}}$ time obtained from
the heuristic also seems to scale linear with the number of qubits although
with a different scaling constant. As mentioned in Sections 5 and 6 the size
of the restricted LPj (6), and therefore the runtime, can be further reduced.
This reduction is achieved by using Hadamard matrices of dimension $d=4\,k$
with $k\in\mathbb{N}$ instead of Sylvester-Hadamard matrices of dimension
$d=2^{k}$ in the construction of the restricted linear operator $V^{[j]}$.
This is based on the Hadamard conjecture, which is known to be true for $d\leq
668$ [22, 23]. Surprisingly, using these Hadamard matrices with $d=4\,k$ not
only yields shorter runtime of the heuristic algorithm but also a significant
reduction of the total $\operatorname{\mathsf{GZZ}}$ time compared to the
original restricted LPj (6).

Our numerical results show that the heuristic algorithm approximates well the
optimal total $\operatorname{\mathsf{GZZ}}$ time, while maintaining a short
runtime. This holds true for both, the average-case and the conjectured worst-
case scaling. Therefore, we hope that this heuristic will prove to be an
important tool to implement fast $\operatorname{\mathsf{GZZ}}$ gates in
practice.

##  9 Conclusion

We investigated the time-optimal multi-qubit gate synthesis introduced in Ref.
[6]. We show that synthesizing time-optimal multi-qubit gates in our setting
is $\mathsf{NP}$-hard. However, we also provide explicit solutions for certain
cases with constant gate time and a polynomial-time heuristics to synthesize
fast multi-qubit gates. Our numerical simulations suggest that these
heuristics provide good approximations to the optimal
$\operatorname{\mathsf{GZZ}}$ gate time. Furthermore, tight bounds on the
scaling of the optimal multi-qubit gate times were shown. More precisely, we
showed that the optimal multi-qubit gate time scales at most as $\lVert
A\oslash J\rVert_{\ell_{1}}$, the $\ell_{1}$-norm of the element-wise division
of the total and physical coupling matrices $A$ and $J$, respectively. We also
conjectured that the optimal $\operatorname{\mathsf{GZZ}}$ gate time scales at
most linear with the number of qubits. Our results are practical to estimate
the execution time of a given circuit, where the entangling gates are
implemented as $\operatorname{\mathsf{GZZ}}$ gates. The execution time is a
crucial parameter, in particular, in the NISQ era since it is limiting the
length of a gate sequence due to finite coherence time.

It is our hope to proof the conjectured linear scaling of the optimal
$\operatorname{\mathsf{GZZ}}$ gate time in the near future. Moreover, we would
like to test and verify our proposed time-optimal multi-qubit gate synthesis
methods in an experiment. Depending on the quantum platform we would like to
develop adapted error mitigation schemes for the $\operatorname{\mathsf{GZZ}}$
gates and investigate their robustness against errors.

## Acknowledgements

We are grateful to Lennart Bittel and Arne Heimendahl for fruitful discussions
on complexity theory and convex optimization, respectively. We also want to
thank Frank Vallentin for valuable comments on our conjecture and proof ideas.

This work has been funded by the German Federal Ministry of Education and
Research (BMBF) within the funding program “quantum technologies – from basic
research to market” via the joint project MIQRO (grant number 13N15522) and
the Fujitsu Services GmbH as part of an endowed professorship “Quantum
Inspired and Quantum Optimization”.

## Appendices

##  Appendix A Challenges in proving Conjecture 16

In this section, we want to discuss some obstacles encountered trying to proof
16. 16 holds, if we show that our constructed solutions from Section 7.1 are
optimal solutions for LP (21). Meaning that there is no other feasible
solution resulting in a larger objective function value compared to our
constructed solution.

We tried an inductive proof which turned out to be intricate due to the
additional degrees of freedom in each induction step. Furthermore, we utilized
the connection to graph theory from Lemma 2 to transform the LP (21) to a LP
over the cut polytope. There, the challenge is the affine mapping between the
elliptope and the cut polytope in Eq. 14 and Definition 1, which alters the
optimization problem crucially. In the following, we discuss another approach
based on showing the sufficiency of the Karush–Kuhn–Tucker (KKT) conditions in
more detail.

###  A.1 Concave program

A convex linear program in standard form minimizes a convex objective function
over a convex set. But LP (21) _maximizes_ a convex objective function over a
convex set. Such optimizations are called _concave programs_. It is known that
the maximum is attained at the extreme points of the polytope
$V^{T}\boldsymbol{y}\leq\boldsymbol{1}$ and therefore might have many local
optima [29]. There are several equivalent sufficient conditions for optimality
[30]. Here, we investigate one in detail. First, we define the _conjugate
function_ for a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$

$f^{*}(\boldsymbol{x})\coloneqq\sup\left\\{\boldsymbol{x}^{T}\boldsymbol{y}-f(\boldsymbol{y})\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\boldsymbol{y}\in\mathcal{D}(f)\right\\}\,,$
(65)

with $\mathcal{D}(f)$ the domain of $f$. Furthermore, we need the _support
function_ $h_{A}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ for a closed convex set
$A$

$h_{A}(\boldsymbol{x})\coloneqq\sup\left\\{\boldsymbol{x}^{T}\boldsymbol{y}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\boldsymbol{y}\in
A\right\\}\,.$ (66)

Then the sufficient optimality condition in our case is

$\lVert\boldsymbol{y}^{*}\rVert_{\ell_{1}}=\sup\left\\{h_{A}(\boldsymbol{x})-(\lVert\boldsymbol{x}\rVert_{\ell_{1}})^{*}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\boldsymbol{y}\in\mathbb{R}^{\binom{n}{2}}\right\\}\,,$
(67)

with
$A=\left\\{\boldsymbol{x}\in\mathbb{R}^{\binom{n}{2}}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}V^{T}\boldsymbol{x}\leq\boldsymbol{1}\right\\}$.
The conjugate function of $\lVert\boldsymbol{x}\rVert_{\ell_{1}}$ is

$(\lVert\boldsymbol{x}\rVert_{\ell_{1}})^{*}=\begin{cases}0\,,&\text{if
}\lVert\boldsymbol{x}\rVert_{\ell_{\infty}}\leq 1\,,\\\
\infty\,,&\text{otherwise}\,.\end{cases}$ (68)

Then we have

$\displaystyle\lVert\boldsymbol{y}^{*}\rVert_{\ell_{1}}$
$\displaystyle=\sup\left\\{h_{A}(\boldsymbol{y})\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\lVert\boldsymbol{y}\rVert_{\ell_{\infty}}\leq
1\right\\}$ (69)
$\displaystyle=\sup\left\\{\boldsymbol{x}^{T}\boldsymbol{y}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\lVert\boldsymbol{y}\rVert_{\ell_{\infty}}\leq
1,V^{T}\boldsymbol{x}\leq\boldsymbol{1}\right\\}$
$\displaystyle=\sup\left\\{\lVert\boldsymbol{x}\rVert_{\ell_{1}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}V^{T}\boldsymbol{x}\leq\boldsymbol{1}\right\\}\,,$

which is the same formulation as the original LP (21).

###  A.2 Dualization

If we can formulate the dual LP to the primal LP (21) and find a feasible
solution, then the dual objective function value upper bounds the primal
objective function value by weak duality [31]. The standard form of an
optimization problem with only linear inequality constraints is f(y) A y ≤b.
Then, the Lagrange dual function is given by

$g(\boldsymbol{\lambda})=-\boldsymbol{b}^{T}\boldsymbol{\lambda}-f^{*}(-A^{T}\boldsymbol{\lambda})\,,$
(70)

where $f^{*}$ denotes the conjugate function [14]. For LP (21) we have
$f(\boldsymbol{y})=-\lVert\boldsymbol{y}\rVert_{\ell_{1}}$ (minus sign due to
the minimization in the standard optimization form).

$\displaystyle(-\lVert\boldsymbol{x}\rVert_{\ell_{1}})^{*}$
$\displaystyle=\sup\left\\{\boldsymbol{x}^{T}\boldsymbol{y}+\lVert\boldsymbol{x}\rVert_{\ell_{1}}\nonscript\>|\allowbreak\nonscript\>\mathopen{}\boldsymbol{y}\in\mathbb{R}^{\binom{n}{2}}\right\\}$
(71)
$\displaystyle=\sup\left\\{\boldsymbol{x}^{T}\boldsymbol{y}+\boldsymbol{x}^{T}\boldsymbol{z}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\boldsymbol{y}\in\mathbb{R}^{\binom{n}{2}},\lVert\boldsymbol{z}\rVert_{\ell_{\infty}}\leq
1\right\\}$
$\displaystyle=\sup\left\\{\boldsymbol{x}^{T}\boldsymbol{y}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}\boldsymbol{y}\in\mathbb{R}^{\binom{n}{2}}\right\\}\,,$

which clearly is unbounded if $\boldsymbol{x}\neq\boldsymbol{0}$. Therefore,
we cannot formulate the dual to LP (21).

###  A.3 Invexity

The KKT conditions are optimality conditions for non-linear optimization
problems. Invexity is a generalization of convexity in the sense that the KKT
conditions are necessary and sufficient for optimality [32]. By invexity we
mean that the objective and constraint functions of the optimization problem
are _Type 1 invex_ functions.

######  Definition 22.

Consider the standard form of an optimization problem f(y) g( y ) ≤0 y ∈S,
where $\mathcal{S}\subseteq\mathbb{R}^{m}$ is defined by
$g(\boldsymbol{y})\leq\boldsymbol{0}$. Then $f$ and $g$ are called _Type 1
invex_ functions at point $\boldsymbol{y}^{*}\in\mathcal{S}$ w.r.t. a common
function $\eta(\boldsymbol{y},\boldsymbol{y}^{*})\in\mathbb{R}^{m}$, if for
all $\boldsymbol{y}\in\mathcal{S}$,

$\displaystyle f(\boldsymbol{y})-f(\boldsymbol{y}^{*})$
$\displaystyle\geq\eta(\boldsymbol{y},\boldsymbol{y}^{*})^{T}\nabla
f(\boldsymbol{y}^{*}),$ (72) $\displaystyle-g(\boldsymbol{y}^{*})$
$\displaystyle\geq\eta(\boldsymbol{y},\boldsymbol{y}^{*})^{T}\nabla
g(\boldsymbol{y}^{*})$

hold. It suffices to consider only the active constraints, i.e. the
constraints, where equality holds $g(\boldsymbol{y}^{*})=0$. [32]

Let $K$ be the scaling factor of the conjectured upper bound, i.e.

$K\coloneqq\begin{cases}n\,,&\text{for odd }n\,,\\\ n-1\,,&\text{for even
}n\,.\end{cases}$ (73)

In the case of LP (21) we have
$\mathcal{S}=\left\\{\boldsymbol{y}\in\mathbb{R}^{\binom{n}{2}}\nonscript\>\middle|\allowbreak\nonscript\>\mathopen{}V^{T}\boldsymbol{y}\leq\boldsymbol{1}\right\\}$
and want to show that $\boldsymbol{y}^{*}=-\boldsymbol{1}/\lfloor n/2\rfloor$
is a global optimum. Furthermore, we have

$\displaystyle f(\boldsymbol{y})$
$\displaystyle=-\lVert\boldsymbol{y}\rVert_{\ell_{1}}\,,\hskip
9.49997ptf(\boldsymbol{y}^{*})=-K\,,\,\,\nabla
f(\boldsymbol{y}^{*})=\boldsymbol{1}\,,$ (74) $\displaystyle
g(\boldsymbol{y})$
$\displaystyle=V^{T}\boldsymbol{y}-\boldsymbol{1}\,,\,\,g(\boldsymbol{y}^{*})=0\,,\hskip
16.49995pt\nabla g(\boldsymbol{y}^{*})=V|_{a}\,,$

where $V|_{a}$ are the columns of $V$ such that
$(V|_{a})^{T}\boldsymbol{y}^{*}=\boldsymbol{1}$, i.e. the active constraints.
To show invexity we have to find a common
$\eta(\boldsymbol{y},\boldsymbol{y}^{*})\in\mathbb{R}^{m}$ such that

$\displaystyle\lVert\boldsymbol{y}\rVert_{\ell_{1}}$ $\displaystyle\leq
K-\eta^{T}(\boldsymbol{y},\boldsymbol{y}^{*})\boldsymbol{1}\,\text{ and}$ (75)
$\displaystyle(V^{T}|_{a})\eta(\boldsymbol{y},\boldsymbol{y}^{*})$
$\displaystyle\leq\boldsymbol{0}\,,$

for all $\boldsymbol{y}\in\mathcal{S}$. It is quite challenging to find a
$\eta(\boldsymbol{y},\boldsymbol{y}^{*})\in\mathbb{R}^{m}$ satisfying both
inequalities. In particular, ansätze motivated from the geometry in small
dimensions eventually fail for larger $n$.

##  Appendix B Acronyms

AGF

    average gate fidelity
AQFT

    approximate quantum fourier transform
BOG

    binned outcome generation
CP

    completely positive
CPT

    completely positive and trace preserving
CS

    compressed sensing
DAQC

    digital-analog quantum computing
DFE

    direct fidelity estimation
DM

    dark matter
DD

    dynamical decoupling
EF

    extended formulation
GST

    gate set tomography
GUE

    Gaussian unitary ensemble
HOG

    heavy outcome generation
IQP

    Instantaneous Quantum Polynomial time
KKT

    Karush–Kuhn–Tucker
LDP

    Lamb-Dicke parameter
LP

    linear program
MAGIC

    magnetic gradient induced coupling
MBL

    many-body localization
MIP

    mixed integer program
ML

    machine learning
MLE

    maximum likelihood estimation
MPO

    matrix product operator
MPS

    matrix product state
MS

    Mølmer-Sørensen
MUBs

    mutually unbiased bases
mw

    micro wave
NISQ

    noisy and intermediate scale quantum
POVM

    positive operator valued measure
PVM

    projector-valued measure
QAOA

    quantum approximate optimization algorithm
QFT

    quantum Fourier transform
QML

    quantum machine learning
QMT

    measurement tomography
QPT

    quantum process tomography
QPU

    quantum processing unit
RDM

    reduced density matrix
SFE

    shadow fidelity estimation
SIC

    symmetric, informationally complete
SPAM

    state preparation and measurement
RB

    randomized benchmarking
rf

    radio frequency
RIC

    restricted isometry constant
RIP

    restricted isometry property
TT

    tensor train
TV

    total variation
VQA

    variational quantum algorithm
VQE

    variational quantum eigensolver
XEB

    cross-entropy benchmarking

## References

  * Wang  _et al._ [2001] X. Wang, A. Sørensen, and K. Mølmer,  _Multibit gates for quantum computing_ , Phys. Rev. Lett. 86, 3907 (2001), arXiv:quant-ph/0012055.
  * Monz  _et al._ [2011] T. Monz, P. Schindler, J. T. Barreiro, M. Chwalla, D. Nigg, W. A. Coish, M. Harlander, W. Hänsel, M. Hennrich, and R. Blatt,  _14-qubit entanglement: Creation and coherence_ , Phys. Rev. Lett. 106, 130506 (2011), arXiv:1009.6126.
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver,  _Superconducting qubits: Current state of play_ , Annual Review of Condensed Matter Physics 11, 369 (2020), arXiv:1905.13641.
  * Figgatt  _et al._ [2019] C. Figgatt, A. Ostrander, N. M. Linke, K. A. Landsman, D. Zhu, D. Maslov, and C. Monroe,  _Parallel entangling operations on a universal ion-trap quantum computer_ , Nature 572, 368 (2019), arXiv:1810.11948.
  * Lu  _et al._ [2019] Y. Lu, S. Zhang, K. Zhang, W. Chen, Y. Shen, J. Zhang, J.-N. Zhang, and K. Kim,  _Scalable global entangling gates on arbitrary ion qubits_ , Nature 572, 363 (2019), arXiv:1901.03508.
  * Baßler  _et al._ [2023] P. Baßler, M. Zipper, C. Cedzich, M. Heinrich, P. H. Huber, M. Johanning, and M. Kliesch,  _Synthesis of and compilation with time-optimal multi-qubit gates_ , Quantum 7, 984 (2023), arXiv:2206.06387.
  * Barahona and Mahjoub [1986] F. Barahona and A. R. Mahjoub,  _On the cut polytope_ , Mathematical Programming 36, 157 (1986).
  * Garey and Johnson [2002] M. R. Garey and D. S. Johnson,  _Computers and intractability_ , Vol. 29 (W. H. Freeman and Company, New York, 2002).
  * Bremner  _et al._ [2016] M. J. Bremner, A. Montanaro, and D. J. Shepherd,  _Average-case complexity versus approximate simulation of commuting quantum computations_ , Phys. Rev. Lett. 117, 080501 (2016), arXiv:1504.07999.
  * Allcock  _et al._ [2023] J. Allcock, J. Bao, J. F. Doriguello, A. Luongo, and M. Santha,  _Constant-depth circuits for Uniformly Controlled Gates and Boolean functions with application to quantum memory circuits_ , arXiv:2308.08539 (2023).
  * Bravyi  _et al._ [2022] S. Bravyi, D. Maslov, and Y. Nam,  _Constant-cost implementations of Clifford operations and multiply controlled gates using global interactions_ , Phys. Rev. Lett. 129, 230501 (2022), arXiv:2207.08691.
  * Bravyi and Maslov [2021] S. Bravyi and D. Maslov,  _Hadamard-free circuits expose the structure of the Clifford group_ , IEEE Trans. Inf. Theory 67, 4546 (2021), arXiv:2003.09412.
  * Maslov and Zindorf [2022] D. Maslov and B. Zindorf,  _Depth optimization of CZ, CNOT, and Clifford circuits_ , IEEE Transactions on Quantum Engineering 3, 1 (2022), arxiv:2201.05215.
  * Boyd and Vandenberghe [2009] S. Boyd and L. Vandenberghe,  _Convex Optimization_ (Cambridge University Press, 2009).
  * Rich [2007] E. Rich,  _The problem classes FP and FNP_ , in  _Automata, Computability and Complexity: Theory and Applications_ (Pearson Education, 2007) pp. 510–511.
  * Johanning [2016] M. Johanning,  _Isospaced linear ion strings_ , Appl. Phys. B 122, 71 (2016).
  * Laurent and Poljak [1995] M. Laurent and S. Poljak,  _On a positive semidefinite relaxation of the cut polytope_ , Linear Algebra and its Applications 223-224, 439 (1995).
  * Deza and Laurent [2009] M. M. Deza and M. Laurent,  _Geometry of Cuts and Metrics_, 1st ed., Algorithms and Combinatorics (Springer Berlin Heidelberg, 2009).
  * E.-Nagy  _et al._ [2013] M. E.-Nagy, M. Laurent, and A. Varvitsiotis,  _Complexity of the positive semidefinite matrix completion problem with a rank constraint_ , Springer International Publishing , 105 (2013), arXiv:1203.6602.
  * Paley [1933] R. E. A. C. Paley,  _On orthogonal matrices_ , Journal of Mathematics and Physics 12, 311 (1933).
  * Hedayat and Wallis [1978] A. Hedayat and W. D. Wallis,  _Hadamard matrices and their applications_ , The Annals of Statistics 6, 1184 (1978).
  * Kharaghani and Tayfeh-Rezaie [2005] H. Kharaghani and B. Tayfeh-Rezaie,  _A Hadamard matrix of order 428_ , Journal of Combinatorial Designs 13, 435 (2005).
  * Đoković  _et al._ [2014] D. Ž. Đoković, O. Golubitsky, and I. S. Kotsireas,  _Some new orders of Hadamard and Skew-Hadamard matrices_ , Journal of Combinatorial Designs 22, 270 (2014), arXiv:1301.3671.
  * Cohn  _et al._ [2021] J. Cohn, M. Motta, and R. M. Parrish,  _Quantum filter diagonalization with compressed double-factorized Hamiltonians_ , PRX Quantum 2, 040352 (2021), arXiv:2104.08957.
  * Spielman and Teng [2004] D. A. Spielman and S.-H. Teng,  _Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time_ , Journal of the ACM 51, 385 (2004), arXiv:cs/0111050.
  * Diamond and Boyd [2016] S. Diamond and S. Boyd,  _CVXPY: A Python-embedded modeling language for convex optimization_ , J. Mach. Learn. Res. 17, 1 (2016), arXiv:1603.00943.
  * Agrawal  _et al._ [2018] A. Agrawal, R. Verschueren, S. Diamond, and S. Boyd,  _A rewriting system for convex optimization problems_ , J. Control Decis. 5, 42 (2018), arXiv:1709.04494.
  * Free Software Foundation [2012] Free Software Foundation,  _GLPK (GNU Linear Programming Kit)_ (2012), version: 0.4.6.
  * Phillips and Rosen [1988] A. T. Phillips and J. B. Rosen,  _A parallel algorithm for constrained concave quadratic global minimization_ , Mathematical Programming 42, 421 (1988).
  * Dür  _et al._ [1998] M. Dür, R. Horst, and M. Locatelli,  _Necessary and sufficient global optimality conditions for convex maximization revisited_ , Journal of Mathematical Analysis and Applications 217, 637 (1998).
  * Bazaraa  _et al._ [2013] M. S. Bazaraa, H. D. Sherali, and C. M. Shetty,  _Nonlinear programming: theory and algorithms, 3rd edition_ (John wiley & sons, 2013).
  * Hanson [1999] M. A. Hanson,  _Invexity and the Kuhn–Tucker Theorem_ , Journal of Mathematical Analysis and Applications 236, 594 (1999).

  *[LP]: linear program
  *[NISQ]: noisy and intermediate scale quantum
  *[IQP]: Instantaneous Quantum Polynomial
time

  *[ LP]: linear program
  *[MIP]: mixed integer program
  *[KKT]: Karush–Kuhn–Tucker