Corresponding author<EMAIL_ADDRESS>
# Optimization of the Variational Quantum Eigensolver for Quantum Chemistry
Applications

R.J.P.T. de Keijzer  Eindhoven University of Technology, P. O. Box 513, 5600
MB Eindhoven, The Netherlands  [  V.E. Colussi  Eindhoven University of
Technology, P. O. Box 513, 5600 MB Eindhoven, The Netherlands  INO-CNR BEC
Center and Dipartimento di Fisica, Università di Trento, 38123 Povo, Italy  B.
S̆korić  Eindhoven University of Technology, P. O. Box 513, 5600 MB Eindhoven,
The Netherlands  S.J.J.M.F. Kokkelmans  Eindhoven University of Technology, P.
O. Box 513, 5600 MB Eindhoven, The Netherlands

(August 27, 2024)

###### Abstract

This work studies the variational quantum eigensolver algorithm, which is
designed to determine the ground state of a quantum mechanical system by
combining classical and quantum hardware. Two methods of reducing the number
of required qubit manipulations, prone to induce errors, for the variational
quantum eigensolver are studied. First, we formally justify the multiple
$\mathbb{Z}_{2}$ symmetry qubit reduction scheme first sketched by Bravyi et
al. [arXiv:1701.08213 (2017)]. Second, we show that even in small, but non-
trivial systems such as H2, LiH, and H2O, the choice of entangling methods
(gate based or native) gives rise to varying rates of convergence to the
ground state of the system. Through both the removal of qubits and the choice
of entangler, the demands on the quantum hardware can be reduced. We find that
in general, analyzing the VQE problem is complex, where the number of qubits,
the method of entangling, and the depth of the search space all interact. In
specific cases however, concrete results can be shown, and an entangling
method can be recommended over others as it outperforms in terms of difference
from the ground state energy.

##  I Introduction

Presently, quantum computing is in the noisy intermediate-scale quantum (NISQ)
era Preskill (2018), where the available universal quantum computers cannot
outperform their classical counterparts except for a few, specific, well-
designed cases Arute  _et al._ (2019). However, even in the NISQ era quantum
computers can be used in practice. One such application is the variational
quantum eigensolver (VQE) algorithm, first proposed in a paper by Peruzzo and
McClean in 2014 Peruzzo  _et al._ (2014). The algorithm is hybrid, e.g. at
certain points a quantum processor unit (QPU) is addressed. The VQE algorithm
has shown proof of concept for small molecules with several different designs
of qubits Kandala  _et al._ (2017); Hempel  _et al._ (2018); Aspuru-Guzik and
Walther (2012); O’Brien  _et al._ (2018); Arute  _et al._ (2020). The aim of
the algorithm is to determine the lowest eigenvalue of a quantum Hamiltonian,
which can be equated to finding the ground state energy of a molecule. In
recent years VQE has been implemented on many different qubit systems and has
become a highly active area of research (c.f. Refs. Cerezo  _et al._ (2020);
Bharti  _et al._ (2021) for recent comprehensive overviews).

To illustrate the VQE algorithm, a schematic of the relevant circuit diagram
is shown in Fig. 1. To determine the ground state energy of a Hamiltonian the
VQE algorithm maps the fermionic degrees of freedom of molecules onto a set of
qubits. These qubits can then be manipulated into a trial state with a certain
energy by applying a sequence of entanglement operators $U_{ent}$ and
individual rotations $U_{i,j}$. The length of this sequence is the depth $d$
which plays an important role in determining which part of the state space can
be addressed by the VQE. The entanglement operator $U_{ent}$ is a multi-qubit
gate that ensures entangled states can be reached and thus its choice is
critical in exploiting the quantum nature of the method. This operator is
necessary for utilizing the computational powers of the multi-qubit structure.
By applying a problem dependent Hamiltonian to the trial state, its energy can
be measured using a QPU. Based on these measurements the classical part of the
computer proposes a new trial state (new rotations) by means of a classical
optimization algorithm Kuroiwa and Nakagawa (2021); IBM-Research-Editorial-
Staff (2017). Because of the limited hardware of quantum computers, each qubit
manipulation has a non-negligible probability of error. Therefore, it is
important to minimize the number of necessary manipulations.

Figure 1: Schematic circuit diagram of the VQE algorithm. The qubits are
manipulated by a sequence of entanglement operators $U_{ent}$ and individual
rotations $U_{i,j}$ into a prepared trial state.

This paper examines two ways of achieving this minimization: (i) by reducing
the number of qubits used to describe the problem as proposed by Bravyi et al.
Bravyi  _et al._ (2017) and (ii) by reducing the number of qubit manipulations
using entanglement methods native to the physical system. Both of these
methods are directed to optimizing the quantum hardware rather than the
classical side.

The layout of this paper is as follows. In Secs. II and III the VQE algorithm
and the method of simulation are described. Sec. IV gives a constructive proof
of the multiple $\mathbb{Z}_{2}$ qubit reduction symmetries scheme by Bravyi
et al. Bravyi  _et al._ (2017) Sec. V describes multi-qubit entangling methods
used in the VQE analysis. Finally, in Sec. VI we compare and contrast the
entanglement methods and their results on small molecule VQE test case
problems.

##  II Quantum measurements and simulation setup

The VQE algorithm exploits the variational principle and aims to find a set of
parameters minimizing the energy of a quantum system. The states of a system
are mapped to qubits and the classical optimization algorithm optimizes the
variational calculation. To do so for Hamiltonians of small molecules requires
a characterization of their electron spin-orbitals. For this the Hartree-Fock
method is applied, which reformulates the Hamiltonian in first-quantization
form McArdle  _et al._ (2018). The Born-Oppenheimer approximation is applied
which ensures that the nuclear contribution to the Hamiltonian is constant
McArdle  _et al._ (2018). The Hamiltonian $H_{1}$ under the Born-Oppenheimer
approximation becomes

$\displaystyle H_{1}$
$\displaystyle=-\sum_{i}\frac{\nabla^{2}_{\vec{r}_{i}}}{2}-\sum_{i,j}\frac{Z_{i}}{|\vec{R}_{i}-\vec{r}_{j}|}$
(1)
$\displaystyle+\sum_{i,j>i}\frac{1}{|\vec{r}_{i}-\vec{r}_{j}|}+\sum_{i,j>i}\frac{Z_{i}Z_{j}}{|\vec{R}_{i}-\vec{R}_{j}|},$

where, $\vec{R}_{i}$ and $Z_{i}$ are the position and charge of nucleus $i$.
$\vec{r}_{i}$ is the position of electron $i$. In VQE, the Hamiltonian is
described in Hartree units 111 Hartree units are also known as atomic units
where the reduced constant of Plank $\hbar$, the electron mass $m_{e}$, the
elementary charge $e$ and Coulomb’s constant $k_{e}=1/4\pi\epsilon_{0}$ are
equal to unity. In this system the Bohr radius
$a_{0}=4\pi\epsilon_{0}\hbar^{2}/m_{e}e^{2}$ is also equal to unity. Hartree
units are often used in molecular level calculations McQuarrie (2008).. The
Hamiltonian is reformulated in second-quantization form McClean  _et al._
(2016) $H_{2}$, by projecting it onto a finite set of orthogonal spin-orbital
modes $\\{\phi_{1}\\}_{i=1}^{n}$, as

$H_{2}=V_{nn}+\sum_{p,q}h_{pq}a^{\dagger}_{p}a_{q}+\frac{1}{2}\sum_{p,q,r,s}h_{pqrs}a^{\dagger}_{p}a^{\dagger}_{q}a_{r}a_{s},$
(2)

where the nuclei-nuclei interaction potential
$V_{nn}=\sum_{i,j>i}Z_{i}Z_{j}/|\vec{R}_{i}-\vec{R}_{j}|$ is constant and the
annihilation and creation operators $a_{i}^{\dagger}$ and $a_{i}$ work on a
set of orthogonal spin-orbital modes $\\{\phi_{1}\\}_{i=1}^{n}$ determined by
the linear combination of atomic orbitals (LCAO) method and which induce the
fermionic algebra McQuarrie (2008). The number of considered orbitals
determines the complexity of the problem and with that the number of required
qubits. In many practical cases only a select few low energy orbitals are
considered while all higher energy atomic orbitals are ignored. The molecular
integrals $h_{pq}$ and $h_{pqrs}$ are referred to as the one-electron and two-
electron integrals respectively

$\displaystyle h_{pq}$ $\displaystyle=\int
S^{*}_{p}(\omega)S_{q}(\omega)d\omega$ (3)
$\displaystyle\enskip\times\int_{\mathbb{R}^{3}}\phi^{*}_{p}(\vec{r})\left(\frac{\nabla^{2}_{\vec{r}}}{2}-\sum_{i}\frac{Z_{i}}{|\vec{R_{i}}-\vec{r}|}\right)\phi_{q}(\vec{r})d\vec{r},$

$\displaystyle h_{pqrs}=$ $\displaystyle\int\int
S^{*}_{p}(\omega_{1})S_{s}(\omega_{1})S^{*}_{q}(\omega_{2})S_{r}(\omega_{2})d\omega_{1}d\omega_{2}$
(4)
$\displaystyle\times\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\phi^{*}_{p}(\vec{r_{1}})\phi^{*}_{q}(\vec{r_{2}})\phi_{r}(\vec{r_{2}})\phi_{s}(\vec{r_{1}})}{|\vec{r_{1}}-\vec{r_{2}}|}d\vec{r_{1}}d\vec{r_{2}},$

where $S_{i}(\omega)$ is the spin part of the spin orbitals 222In this work
the integrals in Eq. 3 and Eq. 4 have been determined using quantum
computation libraries OpenFermion McClean  _et al._ (2017) and Psi4 Smith  _et
al._ (2020) in the STO-3G basis..

Some important concepts in quantum computing are introduced below. These are
required to understand the VQE method and are used in proofs that follow.

######  Definition 1.

Pauli Matrices

$I:=\begin{pmatrix}1&0\\\ 0&1\\\ \end{pmatrix},\enskip
X:=\begin{pmatrix}0&1\\\ 1&0\\\ \end{pmatrix},\enskip
Y:=\begin{pmatrix}0&-i\\\ i&0\\\ \end{pmatrix},\enskip
Z:=\begin{pmatrix}1&0\\\ 0&-1\\\ \end{pmatrix}.$ (5)

######  Definition 2.

$m$-fold Single Qubit Operator

$\sigma^{\rho}_{j}:=I^{\otimes(j-1)}\otimes\rho\otimes I^{\otimes(m-j)},\quad
j\in{1,...,m},\quad\rho\in\\{X,Y,Z\\}.$ (6)

These $m$-fold single qubit operator work on a $m$-qubit state but only change
the state of one of the qubits.

######  Definition 3.

$m$-fold Pauli operator space $P_{m}$

$P_{m}:=\\{\pm 1,\pm i\\}\times\\{I,X,Y,Z\\}^{\otimes m}.$ (7)

$P_{m}$ forms a group under operator multiplication since the Pauli operators
$P_{1}$ form a group. Note that since the Pauli matrices are a complete set
for the Hilbert space of complex $2\times 2$ matrices, $P_{m}$ is a complete
set for the Hilbert space of complex $2^{m}\times 2^{m}$ matrices

The Jordan-Wigner transformation can be used to map the creation and
annihilation operators of Eq. (2) to combinations of Pauli matrices Coleman
(2015). The transformation is given by

$a_{j}=I^{\otimes j-1}\otimes\sigma^{+}\otimes Z^{\otimes m-j},$ (8)
$a^{\dagger}_{j}=I^{\otimes j-1}\otimes\sigma^{-}\otimes Z^{\otimes m-j},$ (9)

where $m$ is the number of qubits. The $\sigma^{\pm}$ operators are defined as

$\sigma^{\pm}:=\frac{X\pm iY}{2}.$ (10)

It can be shown that the transformed operators $a_{j}$ and $a^{\dagger}_{j}$
obey the canonical commutation relations for fermions. Since this
transformation is local the electron configuration can be determined from the
qubit state. Since the number of required qubits is a key factor determining
the required effort for calculation on both physical systems and in
simulations it is desirable to reduce the number of qubits where possible. The
Hamiltonian as constructed in the Jordan-Wigner transformation commutes with
the number spin up and number spin down operators $N_{\uparrow}$ and
$N_{\downarrow}$ so we have

$[H,N_{\uparrow}]=[H,N_{\downarrow}]=0.$ (11)

This symmetry can be used to taper off 2 qubits based respectively on particle
and spin number as described by Bravyi et al. Bravyi  _et al._ (2017). This
effectively projects the Hamiltonian on those states describing the right
number of electrons and spin. This method can be extended to taper off more
qubits if additional symmetries are present in the Hamiltonian, which is the
subject of Sec. IV.

Because $P_{m}$ is a complete set for the Hilbert space of complex
$2^{m}\times 2^{m}$ matrices, any $m$-qubit Hamiltonian $H$ can be written as
a linear combination of elements of the $m$-fold Pauli operators in $P_{m}$ as

$\displaystyle H$
$\displaystyle=\sum_{k}h_{k}\boldsymbol{\sigma}^{k},\quad\boldsymbol{\sigma}^{k}=\sigma^{k}_{1}\otimes\sigma^{k}_{2}\otimes...\otimes\sigma^{k}_{m}\in
P_{m},\quad$ (12) $\displaystyle h_{k}$ $\displaystyle=\langle
H,\boldsymbol{\sigma}^{k}\rangle_{\mathcal{F}}=\sum_{i,j}\overline{H^{T}_{ij}}\boldsymbol{\sigma}^{k}_{ij}=tr(H^{\dagger}\boldsymbol{\sigma}^{k}).$

This is a finite dimensional Fourier decomposition with the Frobenius inner
product $\langle\cdot,\cdot\rangle_{\mathcal{F}}$ on matrix representations of
the operators. Note that in Eq. (12) $\sigma_{ij}^{k}$ is a matrix element,
not a Pauli matrix working on a qubit. Any rotation on a single qubit can be
written as

$R(\alpha,\beta,\gamma,\delta)=e^{i\alpha}Z_{\beta}X_{\gamma}Z_{\delta},$ (13)

where $\alpha,\beta,\gamma,\delta\in[0,2\pi]$. The operators $X_{\phi}$ and
$Z_{\phi}$ denote a rotation of an angle $\phi$ around the given axis and are
given by

$\displaystyle X_{\phi}=\cos(\phi/2)I+\sin(\phi/2)X,$ (14) $\displaystyle
Y_{\phi}=\cos(\phi/2)I-i\sin(\phi/2)Y,$ $\displaystyle
Z_{\phi}=\cos(\phi/2)I+\sin(\phi/2)Z.$

##  III Trial State Initialization and Measurement

In order to initialize a trial state the individual qubits have to be put in a
requested state. Following Kandala et al. Kandala  _et al._ (2017), we assume
control over the rotations of the individual qubits and let the qubits
entangle by a passive always-on interaction described by some drift
Hamiltonian $H_{0}$. As long as the rotations can be performed in a relatively
small time frame compared to the time between two rotations $t_{0}$ we can
model the evolution of the system as a sequence of rotations and entanglement
operations $U_{ent}=e^{iH_{0}t_{0}}$. This method is called the hardware
efficient ansatz Gogolin  _et al._ (2021) and is especially relevant in NISQ
machines where there will likely be a passive interaction Preskill (2018). If
such an interaction is not present and the qubit system possesses some form of
control in entanglement operations then there is more freedom in choosing what
$U_{ent}$ should look like. The depth $d$ of a state preparation is defined as
the length of the sequence of rotations and entanglement operations, or
equivalently the number of entanglement operations applied. Each state on the
Bloch sphere of a single qubit can be reached with a $ZXZ$-rotation. Such a
rotation on a qubit $q$ at a depth $i$ can be written as

$U_{q,i}(\vec{\theta})=Z_{\theta_{1}^{q,i}}X_{\theta_{2}^{q,i}}Z_{\theta_{3}^{q,i}},$
(15)

where $\vec{\theta}$ has three elements for every qubit and depth pair. In
total $3d+3$ rotations would be performed on every qubit. However, in this
paper the initial state will always be the vacuum state. The first
$Z$-rotation can therefore be omitted. Thus the parameter vector
$\vec{\theta}$ is an element of the search space $[0,2\pi]^{\otimes D}$, where
$D=(3d+2)m$. To reach the trial state described by the parameter vector
$\vec{\theta}$ the initial state $|\psi_{init}\rangle$ is rotated and
entangled. The entanglement operator is $U_{ent}$. The prepared trial state
after all rotations and entanglements will be

$\displaystyle|\Psi(\vec{\theta})\rangle=$
$\displaystyle\left(\prod^{m}_{q=1}U_{q,d}(\vec{\theta})\times
U_{ent}\right)\times\left(\prod^{m}_{q=1}U_{q,d-1}(\vec{\theta})\times
U_{ent}\right)\times...$ (16)
$\displaystyle\times\left(\prod^{m}_{q=1}U_{q,0}(\vec{\theta})\right)|\psi_{init}\rangle.$

With the trial state prepared, the expectation value of $H$ can be measured on
a quantum computer as

$\displaystyle\langle H\rangle_{\vec{\theta}}=$
$\displaystyle\langle\Psi(\vec{\theta})|H|\Psi(\vec{\theta})\rangle=\sum_{k}h_{k}\langle\Psi(\vec{\theta})|\boldsymbol{\sigma}^{k}|\Psi(\vec{\theta})\rangle$
(17)
$\displaystyle=\sum_{k}h_{k}\langle\Psi(\vec{\theta})|\sigma_{1}^{k}\otimes\sigma_{2}^{k}\otimes...\otimes\sigma_{m}^{k}|\Psi(\vec{\theta})\rangle.$

In quantum computers these expectation values require sufficient measurements.
In simulation the expectation values require an inner product of matrices and
vectors.

##  IV $\mathbb{Z}_{2}$ symmetry qubit reduction scheme

This section formally justifies the qubit removal process for multiple
$\mathbb{Z}_{2}$ symmetries as described by Bravyi et al. Bravyi  _et al._
(2017). In this process the $\mathbb{Z}_{2}$ symmetries in Hamiltonians Bravyi
_et al._ (2017), often originating from geometric symmetries in the molecule,
are exploited to reduce the number of qubits necessary to describe the
Hamiltonian. These $\mathbb{Z}_{2}$ symmetries include the spin and electron
number symmetries and have recently been linked to more physically intuitive
point group symmetries Setia  _et al._ (2020). The reductions will
substantially lower computation power and time. The qubit reductions based on
$\mathbb{Z}_{2}$ symmetries related to spin and electron number have found
substantial use in recent works Cerezo  _et al._ (2020); Cao  _et al._ (2019);
Ganzhorn  _et al._ (2019) and have been included in quantum computing
libraries, e.g. OpenFermion McClean  _et al._ (2017). Implementations of the
additional $\mathbb{Z}_{2}$ symmetries discussed in this work have are limited
in number, most likely because there are not many or even no such symmetries
present in the small systems that can analyzed in classical simulations.
However, in larger systems they can lead to substantial reduction in hardware
requirement Bravyi  _et al._ (2017). The main idea of the qubit reduction
scheme is to transform a $m$-qubit Hamiltonian in such a way that the Pauli
operators in the decomposition will all have either $I$ or $X$ as the last
$r\in\mathbb{N}$ factors (e.g. YXXZIX and ZXYZXX). If the $i$-th term of every
Pauli operator in the decomposition is $I$ or $X$ then $[H,\sigma^{x}_{i}]=0$,
so the Hamiltonian $H$ and $\sigma^{x}_{i}$ commute and thus share common
eigenvectors. One can therefore replace the $i$-th factors of the operators
with the eigenvalues $\pm 1$ and thus taper off a qubit Bravyi  _et al._
(2017). We ignore the phases of the $m$-fold Pauli operators because these do
not influence commutation. Doing so every $m$-fold Pauli operator becomes its
own inverse.

In the following proof, we place the qubit reduction scheme sketched in
Sec.VII of the paper by Bravyi et al. Bravyi  _et al._ (2017) on a
mathematically rigorous level. We begin first with some definitions and lemmas
required to provide a constructive proof of the scheme.

######  Definition 4.

Symmetry Group

A symmetry group $S$ of a group $G$ is an abelian subgroup of $G$ such that
$-I\not\in S$.

######  Definition 5.

Center and Centralizer

Let $G$ be a group. The center $Z(G)$ of $G$ is the set of elements commuting
with every element of $G$.

$Z(G)=\\{z\in G\quad|\quad zg=gz\quad\forall g\in G\\}.$ (18)

The centralizer $C_{G}(A)$ of $A\subseteq G$ is the set of all elements in $G$
which commute with every element of $A$.

$C_{G}(A)=\\{g\in G\quad|\quad ag=ga\quad\forall a\in A\\}.$ (19)

######  Lemma 1.

Symmetry group commuting with $2^{m}\times 2^{m}$ Hamiltonian
Let $S\subseteq P_{m}$ be a symmetry group. $S$ commutes with a $2^{m}\times
2^{m}$ Hamiltonian $H=\sum_{k}h_{k}\boldsymbol{\sigma}^{k}$ if $S$ commutes
with every Pauli operator in the decomposition of H. Thus, let
$H=\sum_{k}h_{k}\boldsymbol{\sigma}^{k}$. If $\forall k$ $S$ commutes with
$\boldsymbol{\sigma}^{k}$ then $S$ commutes with $H$.

###### Proof.

$sH=s\sum_{k}h_{k}\boldsymbol{\sigma}^{k}=\sum_{k}h_{k}s\boldsymbol{\sigma}^{k}=\sum_{k}h_{k}\boldsymbol{\sigma}^{k}s=Hs$
∎

######  Lemma 2.

Redefining of generators
Let the generators $\tau_{1},\tau_{2},...,\tau_{r}$ generate the abelian group
$S=\langle\tau_{1},\tau_{2},...,\tau_{r}\rangle\subset P_{m}$ and let
$i,j\in\\{1,....r\\}$ with $i\neq j$. Then
$S=\langle\tau_{1},\tau_{2},...,\tau_{i},...,\tau_{j},...,\tau_{r}\rangle=\langle\tau_{1},\tau_{2},...,\tau_{i},...,\tau_{j-1},\tau_{i}\tau_{j},\tau_{j+1},...,\tau_{r}\rangle$

###### Proof.

Since $\tau_{i},\tau_{j}\in S$ it holds that $\tau_{i}\tau_{j}\in S$.
Therefore,
$\langle\tau_{1},\tau_{2},...,\tau_{i},...,\tau_{j-1},\tau_{i}\tau_{j},\tau_{j+1},...,\tau_{r}\rangle\subseteq\langle\tau_{1},\tau_{2},...,\tau_{i},...,\tau_{j},...,\tau_{r}\rangle$.
Now if $\tau_{i}\tau_{j},\tau_{i}\in S$ then
$\tau_{i}\tau_{i}\tau_{j}=\tau_{j}\in S$. Therefore
$\langle\tau_{1},\tau_{2},...,\tau_{i},...,\tau_{j},...,\tau_{r}\rangle\subseteq\langle\tau_{1},\tau_{2},...,\tau_{i},...,\tau_{j-1},\tau_{i}\tau_{j},\tau_{j+1},...,\tau_{r}\rangle$.
∎

######  Lemma 3.

Commutation of product
Let $A,B_{1},B_{2}$ and $B_{3}$ be operators such that $A$ commutes with
$B_{1}$ and anti-commutes with $B_{2}$ and $B_{3}$. Then $B_{1}B_{2}$ and
$B_{2}B_{1}$ anti-commute with A. $B_{2}B_{3}$ and $B_{3}B_{2}$ commute with
A.

###### Proof.

$AB_{1}B_{2}=B_{1}AB_{2}=-B_{1}B_{2}A$ and
$AB_{2}B_{1}=-B_{2}AB_{1}=-B_{2}B_{1}A$. So indeed $B_{1}B_{2}$ and
$B_{2}B_{1}$ anti-commute with A. Furthermore,
$AB_{2}B_{3}=-B_{2}AB_{3}=B_{2}B_{3}A$ and
$AB_{3}B_{2}=-B_{3}AB_{2}=B_{3}B_{2}A$. So indeed $B_{2}B_{3}$ and
$B_{3}B_{2}$ anti-commute with A. ∎

An important result of stabilizer theory Brun (2019) is that a symmetry group
commuting with a $m$-qubit Hamiltonian will at most require $m$ generators.
This is a key fact in the qubit reduction scheme (the exact proof of this is
given by Fujii Fujii (2015)). A proof sketch would take the following form.
Each symmetry group commuting with a $m$-qubit Hamiltonian is a symmetry group
of $P_{m}$. Each symmetry group generator projects a $m$-qubit basis state
onto one of its eigenvalues $\pm 1$. Therefore, each of the generators divides
the Hilbert space of $m$-qubits into two subspaces. The generators share a
common set of eigenvectors as they all commute. The generators thus partition
the $m$-qubit space into $2^{\mathrm{\\#generators}}$ subspaces. Since there
are only $2^{m}$ basis states, $\mathrm{\\#generators}\leq m$.

Now the algorithm for the qubit reduction is given below. For notation
purposes the single qubit operators $X,Y$ and $Z$ on qubit $i$ are sometimes
denoted as $\sigma_{i}^{X},\sigma_{i}^{Y}$ and $\sigma_{i}^{Z}$ respectively.
The first goal is to find a symmetry group $S$, commuting with $H$, which has
a maximal number of generators. First off, every Pauli operator
$\boldsymbol{\sigma}^{k}\in P_{m}$ will be encoded as a binary row-vector
$(a_{x}|a_{z})$ of length $2m$ as follows

$\sigma(a_{x}|a_{z})=\sigma(a_{x1},...,a_{xm}|a_{z1},...,a_{zm})=\prod_{i=1}^{m}(\sigma_{i}^{x})^{a_{xi}}(\sigma_{i}^{z})^{a_{zi}},$
(20)

where $a_{xi},a_{zi}\in\\{0,1\\}$. As an example, $IZXY=(0011|0101)$. Using
the standard inner product on $2m$-dimensional vectors in $\mathbb{Z}_{2}$,
resulting in a value 0 or 1, the following commutation relation holds

$\sigma(a_{x}|a_{z})\sigma(b_{x}|b_{z})=(-1)^{a_{x}\cdot b_{z}+a_{z}\cdot
b_{x}}\sigma(b_{x}|b_{z})\sigma(a_{x}|a_{z}).$ (21)

Following this, a scalar product $\times$ is defined on the vector space of
binary vectors of length $2m$ as

$a\times b:=a_{x}\cdot b_{z}+a_{z}\cdot b_{x}\in\\{0,1\\}.$ (22)

This vector product is symmetrical and linear in both $a$ and $b$. By
evaluating this vector product it can easily be determined if the
corresponding operators will commute or anti-commute, see Eq. (21).

The set of $k$ Pauli operators
$\boldsymbol{\sigma}^{1},...,\boldsymbol{\sigma}^{k}$ appearing in the
decomposition of $H$ can be represented by vectors
$(a_{x}^{1}|a_{z}^{1}),...,(a_{x}^{k}|a_{z}^{k})$ which can in turn be
represented in a $k\times 2m$ matrix $G$ as

$G=\left[\begin{array}[]{ccc|ccc}\cdots&a_{x}^{1}&\cdots&\cdots&a_{z}^{1}&\cdots\\\
\cdots&a_{x}^{2}&\cdots&\cdots&a_{z}^{2}&\cdots\\\ &\vdots&&&\vdots&\\\
\cdots&a_{x}^{r}&\cdots&\cdots&a_{z}^{r}&\cdots\\\ \end{array}\right]$ (23)

Now if a Pauli operator $\sigma(b_{x}|b_{z})$ were to commute with every term
in $H$ it must be that $G\times(b_{x}|b_{z})=0$, where the vector product
$\times$, defined in Eq. (22), is taken for every row in $G$. From the theory
of stabilizer codes the concept of a parity check matrix $E$ is borrowed
Gottesman (1997). This matrix $E$ is defined as

$E=\left[\begin{array}[]{ccc|ccc}\cdots&a_{z}^{1}&\cdots&\cdots&a_{x}^{1}&\cdots\\\
\cdots&a_{z}^{2}&\cdots&\cdots&a_{x}^{2}&\cdots\\\ &\vdots&&&\vdots&\\\
\cdots&a_{z}^{r}&\cdots&\cdots&a_{x}^{r}&\cdots\\\ \end{array}\right]$

Note that $EG^{T}=0$. Now if the operator $\sigma(b_{x}|b_{z})$ were to be
part of the symmetry group $S$ then $(b_{x}|b_{z})\in\text{ker}(E)$. Next a
basis for $\text{ker}(E)$ can be constructed. Let
$d:=\text{dim}(\text{ker}(E))$. These basis vectors $b^{1},...,b^{d}$ give
rise to Pauli operators
$\boldsymbol{\sigma}^{1},...,\boldsymbol{\sigma}^{d}\in P_{m}$. However, it is
not guaranteed that these Pauli operators mutually commute and thus will form
a set of generators for a symmetry group $S$. Commutation can be enforced by
applying the algorithm below to obtain a set of new vectors
$g^{1},...,g^{r}\in\text{span}(b^{1},...,b^{d})$ with $r\leq d$. This is the
part not formally justified in the paper by Bravyi et al Bravyi  _et al._
(2017). The operators corresponding to these vectors will be linearly
independent and commuting. We are thus looking for generators of a maximal
Abelian subgroup of $\text{ker}(E)$. Hence

$\forall i,j\in\\{1,...,r\\}:g^{i}\times g^{j}=0.$ (24)

The algorithm for finding maximal abelian subgroups Kurosh (2014) is described
as follows

  1. 1.

Define $G=\langle\sigma(b^{1}),...\sigma(b^{d})\rangle$.

  2. 2.

Now compute the center $C$ of $G$ and set $A:=C$. Since $A$ is abelian
$A\subseteq C_{G}(A)$.

  3. 3.

If $C_{G}(A)=A$ then $A$ is a maximal abelian subgroup of $G$. Else, take an
arbitrary $a\in C_{G}(A)\backslash A$. Define $A:=\langle A\cup\\{a\\}\rangle$
and repeat step 3.

This procedure must halt as $G$ is finitely generated and $A$ keeps growing
while still $A\subseteq G$ holds. As mentioned previously, it is known that
there are maximally $m$ generators in $S$ Fujii (2015). Thus
$S=\langle\tau_{1},\tau_{2},...\tau_{r}\rangle$ where $r\leq m$.

At this point in the proof, Ref. Bravyi  _et al._ (2017) claims that in all
considered examples the symmetry generators $\tau_{1},...,\tau_{r}$ are
$Z$-type Pauli operators and in addition that a subset of qubits
$q(1),...,q(k)$ can be chosen such that

$\sigma_{q(i)}^{x}\tau_{j}=(-1)^{\delta_{i,j}}\tau_{j}\sigma_{q(i)}^{x}.$ (25)

However, in the following, we find that this claim is nontrivial, although we
will show that the assumption that the generators are $Z$-type is unnecessary.

######  Definition 6.

Generator transform
Let $i\in{1,...,r}$ and let $\tau_{i}$ be a generator of a symmetry group
$S=\langle\tau_{1},...,\tau_{r}\rangle$. Let $\tau_{i}$ anti-commute with a
single qubit operator $\sigma_{j}^{\rho}$, where $j\in{1,...m}$ and
$\rho\in\\{X,Y,Z\\}$. Define new operators $\tilde{\tau}_{k}$ as

$\tilde{\tau_{k}}=\begin{cases}\tau_{k}\quad\quad\quad\mathrm{if}\quad
k=i\vee\tau_{k}\sigma_{j}^{\rho}=\sigma_{j}^{\rho}\tau_{k}\\\
\tau_{i}\tau_{k}\quad\quad\mathrm{else}\end{cases}.$ (26)

These new operators are linearly independent, mutually commuting and generate
$S$ according to lemma 2. Furthermore, these new generators were transformed
in such a way that $\tilde{\tau}_{i}$ anti-commutes with $\sigma_{j}^{\rho}$
while all other $\tilde{\tau}_{k}$ commute with $\sigma_{j}^{\rho}$ according
to lemma 3.

For each $n\leq r$ the qubit reduction scheme requires a triple
$\\{A_{n},B_{n},V_{n}\\}$ with $A_{n}=\\{q(i)|i=1,...,n\quad
i\not=j\rightarrow q(i)\not=q(j)\\}$ a set of qubit-indices,
$B_{n}=\\{\rho(i)|i=1,...,n\quad\rho(i)\in\\{X,Z\\}\\}$ a set of labels and
$V_{n}=\\{\nu_{i}|i=1,...,r,\quad\nu_{i}=\tau_{1}^{\alpha_{1}}...\tau_{r}^{\alpha_{r}},\alpha_{j}\in\\{0,1\\}\\}$
a set of generators of $S$. This triple $\\{A_{n},B_{n},V_{n}\\}$ should
satisfy the predicate

$\displaystyle\Pi(\\{A_{n},B_{n},V_{n}\\}):$ $\displaystyle\quad\forall
i\in\\{1,...,n\\}\hskip 5.0pt\forall j\in\\{1,...,r\\}$ (27)
$\displaystyle\sigma^{\rho(i)}_{q(i)}\nu_{j}=(-1)^{\delta_{i,j}}\nu_{j}\sigma^{\rho(i)}_{q(i)}$

This requirement states that each qubit $i$ should anti-commute only with
generator $i$.

######  Lemma 4.

Assume that $\Pi(\\{A_{n},B_{n},V_{n}\\})$ is satisfied. Let
$\tau_{n+1}=P_{1}P_{2}...P_{m}$ where $P_{i}\in\\{I,X,Y,Z\\}$. There exists
$j\in\\{1,...,m\\}$ such that $j\not\in A_{n}$ and $P_{j}\not=I$.

###### Proof.

Proof by contradiction. Assume that $\forall i\not\in A_{n}:P_{i}=I$. Define
the subgroup $S_{l}^{n}<P_{m}$ for $l\leq n$ as

$\displaystyle S_{l}^{n}\left(\\{A_{n},B_{n},V_{n}\\}\right)=$ (28)
$\displaystyle\\{\sigma=P_{1}P_{2}...P_{m}|\hskip 5.0ptj\not\in
A_{n}\Rightarrow P_{j}=I,\hskip 5.0pt\forall i\in\\{1,...,l\\}$
$\displaystyle\sigma\tau_{i}=\tau_{i}\sigma,\hskip 5.0pt\forall
w\in\\{1,...,n\\}\hskip
5.0pt\sigma\sigma^{\rho(w)}_{q(w)}=\sigma^{\rho(w)}_{q(w)}\sigma\\}$

$S_{l}^{n}$ is a subgroup of $P_{m}$ as for every $l\leq n$ it is true that
$I^{\otimes m}\in S_{l}$ and if $\sigma_{a},\sigma_{b}\in S_{l}^{n}$ then
$\sigma_{a}\sigma_{b}\in S_{l}^{n}$ (the inverses are of course trivially
satisfied as every element is its own inverse). One can see that
$\tau_{n+1}\in S^{n}_{n}$ by construction. From the condition $\forall
w\in\\{1,...,n\\}\quad\sigma\sigma^{p_{q(w)}}_{q(w)}=\sigma^{p_{q(w)}}_{q(w)}\sigma$
it follows that $|S^{n}_{0}|=2^{n}$, as in every position in $A_{n}$ only two
operators out of $\\{I,X,Y,Z\\}$ are allowed. It is easy to see that
$S^{n}_{l+1}$ is a subgroup of $S^{n}_{l}$. Lagrange’s theorem gives

$|S^{n}_{l}|/|S^{n}_{l+1}|=|S^{n}_{l}/S^{n}_{l+1}|\in\mathbb{N}.$ (29)

Note now that $\sigma^{\rho(l+1)}_{q(l+1)}$ is in $S^{n}_{l}$ but not in
$S^{n}_{l+1}$, since it must anti commute with $\tau_{l+1}$. So,
$|S^{n}_{l}|/|S^{n}_{l+1}|\geq 2$. Then $|S^{n}_{0}|/|S^{n}_{n}|\geq 2^{n}$.
Therefore, $|S_{n}|=1$ and thus $S_{n}={I}$, the only group with one element.
But by construction it was known that $\tau_{n+1}\in S^{n}_{n}$. Since $I$ can
not be a generator a contradiction is reached. Therefore, there must indeed
exist a $j\in{1,...,m}$ such that $j\not\in A_{n}$ and $P_{j}\not=I$. ∎

######  Theorem 5.

Let $S=\langle\tau_{1},\tau_{2},...\tau_{r}\rangle$ be a symmetry group. For
every $n\leq r$ there exists a triple $\\{A_{n},B_{n},V_{n}\\}$ such that
$\Pi(\\{A_{n},B_{n},V_{n}\\})$ is satisfied.

###### Proof.

The proof is given by induction over $n$

Base case $n=1$: consider $\tau_{1}=P_{1}P_{2}...P_{m}$ where
$P_{i}\in\\{I,X,Y,Z\\}$. Since $\tau_{1}$ is a generator there must be a
$P_{j}\not=I$. The first qubit index is chosen as $q(1)=j$ thus
$A_{1}=\\{j\\}$. The first label $\rho(1)$ is chosen as

$\rho(1)=\begin{cases}X\quad\quad\mathrm{if}\quad P_{j}\in\\{Y,Z\\}\\\
Z\quad\quad\mathrm{if}\quad P_{j}=X\end{cases}.$ (30)

Thus, $B_{1}=\\{\rho(1)\\}.$ The generator transform from definition 6 can now
be applied to the original generators of $S$ to yield new generators
$\tilde{\tau}_{1},...,\tilde{\tau}_{r}$ which form $V_{1}$. For ease of
argumentation the new generators $\tilde{\tau}_{1},...,\tilde{\tau}_{r}$ are
renamed to $\tau_{1},...,\tau_{r}$. Lemmas 2 and 3 now show that
$\Pi(\\{A_{1},B_{1},V_{1}\\})$ is satisfied.

Induction cases $2\leq n+1\leq r$.

From lemma 4 there exists a $j\in{1,...m}$ such that $j\not\in A_{n}$ and
$P_{j}\not=I$. One can then choose $q(n+1)=j$ and

$\rho(n+1)=\begin{cases}X\quad\quad\mathrm{if}\quad P_{j}\in\\{Y,Z\\}\\\
Z\quad\quad\mathrm{if}\quad P_{j}=X\end{cases}.$ (31)

The generator transform of definition 6 can be applied to the generators in
$V_{n}$ to yield new generators $\tilde{\tau}_{1},...,\tilde{\tau}_{r}$. These
generators form $V_{n+1}$. For ease of argumentation the new generators
$\tilde{\tau}_{1},...,\tilde{\tau}_{r}$ are renamed to
$\tau_{1},...,\tau_{r}$. $\Pi(\\{A_{n+1},B_{n+1},V_{n+1}\\})$ is now
satisfied. ∎

The existence and construction of a triple $\\{A_{r},B_{r},V_{r}\\}$ such that
$\Pi(\\{A_{r},B_{r},V_{r}\\})$ is satisfied has been proven.

Using the found triple $\\{A_{r},B_{r},V_{r}\\}$ one can define the operators
$U_{i}$ for $i=1,...r,$ according to

$U_{i}=\begin{cases}\frac{1}{2}(\sigma^{z}_{q(i)}+\tau_{i})(\sigma^{z}_{q(i)}+\sigma^{x}_{q(i)})\quad\mathrm{if}\quad\rho(i)=Z\\\
\frac{1}{\sqrt{2}}(\sigma^{x}_{q(i)}+\tau_{i})\quad\quad\quad\quad\quad\quad\hskip
1.99997pt\mathrm{if}\quad\rho(i)=X\end{cases}.$ (32)

Below certain properties of the defined operators $U_{i}$ are shown. Since
Pauli matrices are Hermitian and unitary
$\sigma^{\rho(i)\dagger}_{q(i)}=\sigma^{\rho(i)}_{q(i)}$,
$\tau_{i}^{\dagger}=\tau_{i}$ and
$\sigma^{\rho(i)\dagger}_{q(i)}\sigma^{\rho(i)}_{q(i)}=\tau_{i}^{\dagger}\tau_{i}=1$.
Using the implied commutation relations the following holds

$\displaystyle U_{i}^{\dagger}U_{i}=1,\quad
U_{i}\sigma^{x}_{q(i)}U_{i}^{\dagger}=\tau_{i},$ (33) $\displaystyle
U_{i}\sigma^{x}_{q(j)}U_{i}^{\dagger}=\sigma^{x}_{q(j)}\quad\mathrm{for}\quad
j\not=i.$

This allows for the definition of the unitary operator
$U=U_{1}U_{2}...U_{r}W$, where $W$ is a permutation operator assigning qubit
$i$ to qubit $q(i)$. As a combination of unitary operators $U$ is also
unitary. The following relation follows from Eq. (33)
$U\sigma^{x}_{i}U^{\dagger}=\tau_{i}.$

The transformed Hamiltonian $H^{\prime}$ can be defined as

$H^{\prime}=U^{\dagger}HU=\sum_{k}h_{k}U^{\dagger}\boldsymbol{\sigma}^{k}U.$
(34)

Define $\boldsymbol{\eta}^{k}=U^{\dagger}\boldsymbol{\sigma}^{k}U$. Since
$\forall i,j\in\\{1,...,r\\}$ $[\boldsymbol{\sigma^{i}},\tau_{j}]=0$ it must
be that $\forall i,j\in\\{1,...,r\\}$
$[\boldsymbol{\eta^{i}},\sigma^{x}_{j}]=0$. Therefore, $H^{\prime}$ commutes
with the last $r$ qubits and thus shares eigenvectors with
$\sigma^{x}_{m-(r-1)},...\sigma^{x}_{m}$. The last $r$ qubits can thus be
replaced by the $X$ operator eigenvalues $\pm 1$.

Since the proof presented in this section is constructive, it should be
feasible to implement this scheme in practice. For instance, in Ref. Bravyi
_et al._ (2017) this scheme has been applied successfully to Hamiltonians of
small molecules such as BeH2 and HCl to reduce the number of qubits necessary
to describe the problem by four qubits, which is two more than just those from
electron and spin number. In the remainder of the paper, we focus on the
importance of the role that entangling methods play in the convergence of the
VQE. We analyze only relatively simple systems that posses no additional
$\mathbb{Z}_{2}$ symmetries other than the electronic and spin symmetries
present.

##  V Multi-qubit entangling methods

In the preparation of the trial state the entanglement operator $U_{ent}$
ensures that the full Hilbert space of trial states can be reached. In this
section, a representative group of multi-qubit entangling methods will be
described that will be used in the VQE simulations of Sec. VIII. Some of the
entanglers are gate based while others are native to system interactions a_o .
With regards to the latter, we study in particular the Rydberg quantum
computing system and its accompanying Rydberg interaction Saffman (2016);
Adams  _et al._ (2019); Jaksch  _et al._ (2000); Morgado and Whitlock (2020a).
The goal here is to show that entanglement gates based on these interactions
can reach similar energy differences with the ground state as standard multi
qubit entanglement gates in the form of combinations of CNOT gates. The
advantage of these interaction based entanglement operations is that they are
native to a system and thus do not require extensive constructions or highly
specific schemes. To illustrate the comparable convergence between these
methods we look at a small ensemble of interaction based/native entanglement
operations and compare this to standard combinations of CNOT gates.

Figure 2: Quantum circuit diagram of the CNOT chain and pairs methods for
$m=4$.

The standard entanglement operator on two qubits is the CNOT gate. A way to
entangle more than two qubits is by applying the CNOT gate in a chain or
pairwise as is shown in Fig. 2.

Since qubit operations are prone to errors it is desirable to use as few as
possible to entangle all qubits. The chain and pair methods respectively use
$m-1$ and $m(m-1)/2$ CNOT operations for $m$ qubits. A way of reducing the
number of operations is using entanglement interactions native to the qubit
system. These interactions can entangle all qubits using only one operation
and thus when considered as a single gate, are a more efficient way of
entangling all qubits. Note that this efficiency only holds if the interaction
is native to the qubit system, on other systems it might have to be
constructed out of multiple other interactions.

A method proposed by Jaksch et al. Jaksch  _et al._ (2000) uses the Rydberg
interaction. The Rydberg interaction is a long range interaction which, when
one Rydberg atom (qubit) in the excited state, blocks others from that state.
Thus, it can function much like a CNOT gate. However, because of the long
range interaction it also allows for a Cm-1NOT gate where one qubit controls
all $m-1$ other qubits, thus multi-qubit entanglement is achieved Walker and
Saffman (2012); Morgado and Whitlock (2020b). Another physical multi-qubit
entangling method is the Krawtchouk chain Groenland and Schoutens (2018). The
Krawtchouk chain is an entangling method based on the physical Hamiltonian of
a Rydberg system which can be described as an interacting spin chain. Such a
spin chain could consist of trapped ions or ultra-cold Rydberg atoms Clark
_et al._ (2005). The Krawtchouk chain allows for the construction of the PSTm
gate which mirrors the left and right side of a chain of qubits and multiplies
it by a factor $i$ (for instance $|01101\rangle\rightarrow i|10110\rangle$).
This effect is often called perfect state transfer (PST).

The Krawtchouk chain is an important method to construct other multi-qubit
gates. Several studies have shown that the Krawtchouk Hamiltonian gate can be
transformed, using only few one and two qubit rotations Groenland and
Schoutens (2018); Clark  _et al._ (2005). Two examples of these
transformations are the iSWAP2 gate and the Cm-1NOT gate. The iSWAP2 gate
swaps the last two qubits and multiplies the state by a phase $i$. This
swapping is controlled by the first $m-2$ qubits.

In the quantum chemistry problems considered, some combinations of orbitals
interact more than others. Similarly, in some of these entanglement schemes,
such as CNOT pairs and Cm-1NOT, the qubits play different roles, e.g. control
versus target. This means that different assignments of orbitals to qubits may
result in different convergence. Based on prior information on the problem
from quantum chemistry, it might thus be beneficial to assign certain qubits
to certain orbitals. Since this paper does not focus on initial chemistry
problem analysis, we simply assign orbital 1 to qubit 1, etc.

##  VI Results

In this section, we investigate the behaviour of the discussed multi-qubit
entangling methods and classical optimization algorithms on VQE problems of
H2, LiH and H2O. The simultaneous perturbation stochastic approximation (SPSA)
algorithm, which is described in App. A, is used as the classical optimization
algorithm throughout. For problems with a small number of qubits, such as
considered here, the dimension of the state space is limited. In other more
complex problems it may be advantageous instead to start off with a specific
initial state to check if the VQE algorithm is able to improve upon this. For
consistency, we choose the initial state to be the vacuum state
$|\psi_{init}\rangle=|0...0\rangle$ in this work. When the found energies
$E_{VQE}$ of two methods are compared to the ground state energy $E_{g}$, the
method with the lower difference $E_{VQE}-E_{g}$ is said to outperform the
other.

Figure 3: A typical energy difference from the ground state energy versus
iteration number for the entangling methods PSTm, Cm-1NOT and iSWAP2 compared
to the standard CNOT chain. The VQE problem describes a LiH molecule (6
qubits) at depth 6, with classical optimization SPSA.

Figure 4: Differences from the exact ground energy vs. depth for $10^{3}$
iteration SPSA simulations comparing the entangling methods for H2 (a) and LiH
(b) at depths 1 to 6. Energy differences are calculated using the geometric
average over several internal configurations of the molecule. The spread of
each data point is given by the geometric standard deviation. Depth values are
slightly offset for readability.

Specifically, H2 is modelled with a 1s orbital for each atom resulting in a 4
qubit system which is reduced to a 2 qubit system by the $\mathbb{Z}_{2}$
symmetries. We run the simulations over interatomic distances from
0.392$a_{0}$ to $0.931a_{0}$ in steps of 0.049$a_{0}$, with $a_{0}$ the Bohr-
radius. The LiH molecule is aligned on the $x$-axis and we consider the $1s$,
$2s$ and $2p_{x}$ orbitals on the Li atom and a $1s$ on the H atom resulting
in 8 qubits which is again reduced to 6. Simulations are run over interatomic
distances from 0.72$a_{0}$ to $1.71a_{0}$, in steps of 0.09$a_{0}$. For both
molecules the resulting errors with the ground state energy are geometrically
averaged over the chosen bond length.

To investigate the dependence of the VQE results on entangling methods, we
first analyze the ground state of LiH for d=6 using the PSTm, iSWAP2, Cm-1NOT
and CNOT chain gates (discussed in Sec. V) as shown in Fig. 3. Here it is
found that methods based on the Rydberg interaction are able to perform
comparably to the CNOT chain method (see Sec. V.). Repeating these simulations
for LiH and H2 as a function of depth $d$ we get the resulting energy
differences given in Fig. 4. These results seem to further support that
Rydberg interaction entangling methods can compare to the CNOT chain method in
terms of performance. The optimal entangling method however appears to be
problem dependent. The entanglement depth $d$ can open up the search space for
the trial state allowing lower energy differences to be reached as can be seen
in many of the LiH cases of Fig. 4. However, it can also overcomplicate the
problem by introducing too many rotation parameters and thus take more
iterations to converge as we see in H2 cases. Another interesting fact is seen
for PSTm and iSWAP2 at H2, where the entangling method is not fit for the
problem regardless of the depth, and an energy difference as low as the CNOT
chain and CkNOT differences is not reached.

To see that the preference of one entangling method over the other is not the
result of the geometric averaging of the energy difference we solve the VQE
problem for H2O (8 qubits) at bond length $0.99a_{0}$ with angles between the
H atoms ranging from $\varphi=0.2\pi$ to $\varphi=0.7\pi$ in steps of $0.1\pi$
(the angle with the lowest ground state energy for H2O is $\varphi\approx
0.58\pi$). Doing so allows us to determine the influence of a slightly changed
internal configuration of the molecule on the VQE convergence and entangler
performance. The H2O molecule lies in the $xy$-plane, the 1s orbital is filled
for the O atom, the $2s$, $2p_{x}$ and $2p_{y}$ are active, and a 1s orbital
is considered active on each of the H atoms. After reduction, this results in
8 qubits. Fig. 5 shows the results per angle while Fig. 6 shows the energy
differences geometrically averaged over all angles together with the geometric
standard deviation. The results in Fig. 5 further support the claim that
Rydberg based entanglers can perform as well as the CNOT chain. Furthermore,
in Fig. 5 it is seen that for all angles the iSWAP2 and PSTm gates tend to
converge faster than the CNOT chain and Cm-1NOT gates. We also see that the
Cm-1NOT method tends to converge the slowest. At the $\varphi=0.5\pi$ angle we
see the CNOT chain getting stuck in a local minimum and not able to decrease
the energy any further. The fact that iSWAP2 and PSTm perform consistently
well at slightly different internal configurations shows that some
entanglement methods are more suited to a certain molecular problem than
others.

Figure 5: Iteration number versus energy difference with the ground state
energy for H2O with depth $d=6$ at different angles $\varphi$ between the H
atoms, optimized with the SPSA algorithm: (a) $\varphi=0.2\pi$, (b)
$\varphi=0.3\pi$, (c) $\varphi=0.4\pi$, (d) $\varphi=0.5\pi$, (e)
$\varphi=0.6\pi$, and (f) $\varphi=0.7\pi$.

Figure 6: Iteration number versus energy difference with the ground state
energy for H2O with depth $d=6$ geometrically averaged over the different
angles $\varphi$ between the H atoms, optimized with the SPSA algorithm. In
(a) results are shown for the CNOT chain and PSTm entanglement methods, in (b)
results are shown for the ISWAP2 and Cm-1NOT entanglement methods. Spread
obtained via geometric variance.

##  VII Conclusion

In this paper, we have focused on optimizing the VQE algorithm for quantum
chemistry applications. First, a constructive proof was shown for a qubit
reduction scheme using the $\mathbb{Z}_{2}$ symmetries of Hamiltonians, based
on the paper by Bravyi et al. Bravyi  _et al._ (2017), which effectively
reduces the total number of necessary qubits. The multiple $\mathbb{Z}_{2}$
qubit reduction can thus be an important ingredient for reducing the total
number of qubit manipulations for VQE problems.

After this the performance of the VQE algorithm as a function of the
entangling method choice and depth was examined. Entangling methods based on
both physical interactions, such as PSTm and Cm-1NOT gave a convergence that
was comparable to that of the CNOT chain method in the tested small molecule
VQE problems. Looking at individual problems there seems to be preference for
certain entangling operations over others: in H2 the CNOT chain and CkNOT
performed well while in LiH and H2O the CkNOT and PSTm gates were preferred.
We thus see that even in VQE problems with low number of qubits the choice of
entanglement operation plays an important role in the speed of convergence.

To further investigate the link between entangling method and VQE problems,
VQE calculations should be performed for more complicated molecular systems to
see if a certain preference is maintained at a higher number of qubits. In
this regime the barren plateau problem McClean  _et al._ (2018) will have to
accounted for. Another interesting point would be the prior analysis of the
quantum chemistry problem to generate more optimal assignments from orbitals
to qubits. Furthermore, our findings motivate future studies on the link
between a molecular VQE problem and the entangling method. A measure of which
entangling operation would do well prior to calculation would be valuable
knowledge concerning the demands on hardware. Such a concept would be
foundational towards the design of efficient VQE circuitry.

## Acknowledgements

We thank Jasper Postema and Gijs Groeneveld for discussions. This research is
financially supported by the Dutch Ministry of Economic Affairs and Climate
Policy (EZK), as part of the Quantum Delta NL programme, and by the
Netherlands Organisation for Scientific Research (NWO) under Grants No.
680-47-623 and 680-92-18-05.

## Conflict of interest

The authors have no conflicts to disclose.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, _et al._ , Nature 574, 505 (2019).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 4213 (2014), arXiv:1304.3061 [quant-ph] .
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature (London) 549, 242 (2017), arXiv:1704.05018 [quant-ph] .
  * Hempel  _et al._ (2018) C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Phys. Rev. X 8, 031022 (2018).
  * Aspuru-Guzik and Walther (2012) A. Aspuru-Guzik and P. Walther, Nature Physics 8 (2012), 10.1038/nphys2253.
  * O’Brien  _et al._ (2018) T. O’Brien, P. Rożek, and A. Akhmerov, Physical review letters 120, 220504 (2018).
  * Arute  _et al._ (2020) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, _et al._ , Science 369, 1084 (2020).
  * Cerezo  _et al._ (2020) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, “Variational quantum algorithms,” (2020), arXiv:2012.09265 [quant-ph] .
  * Bharti  _et al._ (2021) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, “Noisy intermediate-scale quantum (nisq) algorithms,” (2021), arXiv:2101.08448 [quant-ph] .
  * Kuroiwa and Nakagawa (2021) K. Kuroiwa and Y. O. Nakagawa, Phys. Rev. Research 3, 013197 (2021).
  * IBM-Research-Editorial-Staff (2017) IBM-Research-Editorial-Staff, “How to measure a molecule’s energy using a quantum computer,”  (2017).
  * Bravyi  _et al._ (2017) S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, arXiv e-prints , arXiv:1701.08213 (2017), arXiv:1701.08213 [quant-ph] .
  * McArdle  _et al._ (2018) S. McArdle, S. Endo, A. Aspuru-Guzik, S. Benjamin, and X. Yuan, arXiv e-prints , arXiv:1808.10402 (2018), arXiv:1808.10402 [quant-ph] .
  * Note (1) Hartree units are also known as atomic units where the reduced constant of Plank $\hbar$, the electron mass $m_{e}$, the elementary charge $e$ and Coulomb’s constant $k_{e}=1/4\pi\epsilon_{0}$ are equal to unity. In this system the Bohr radius $a_{0}=4\pi\epsilon_{0}\hbar^{2}/m_{e}e^{2}$ is also equal to unity. Hartree units are often used in molecular level calculations McQuarrie (2008).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016), arXiv:1509.04279 [quant-ph] .
  * McQuarrie (2008) D. McQuarrie,  _Quantum Chemistry_, v. 1 (University Science Books, 2008) Chap. 10, pp. 266–290.
  * Note (2) In this work the integrals in Eq. 3 and Eq. 4 have been determined using quantum computation libraries OpenFermion McClean  _et al._ (2017) and Psi4 Smith  _et al._ (2020) in the STO-3G basis.
  * Coleman (2015) P. Coleman, “Simple examples of second quantization,” in  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015) p. 71–94.
  * Gogolin  _et al._ (2021) C. Gogolin, G.-L. Anselmetti, D. Wierichs, and R. M. Parrish, New Journal of Physics (2021).
  * Setia  _et al._ (2020) K. Setia, R. Chen, J. E. Rice, A. Mezzacapo, M. Pistoia, and J. Whitfield, “Reducing qubit requirements for quantum simulation using molecular point group symmetries,” (2020), arXiv:1910.14644 [quant-ph] .
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Chemical reviews 119, 10856 (2019).
  * Ganzhorn  _et al._ (2019) M. Ganzhorn, D. J. Egger, P. Barkoutsos, P. Ollitrault, G. Salis, N. Moll, M. Roth, A. Fuhrer, P. Mueller, S. Woerner, _et al._ , Physical Review Applied 11, 044092 (2019).
  * McClean  _et al._ (2017) J. R. McClean, K. J. Sung, I. D. Kivlichan, Y. Cao, C. Dai, E. Schuyler Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner, T. Hardikar, V. Havlíček, O. Higgott, C. Huang, _et al._ , arXiv e-prints , arXiv:1710.07629 (2017), arXiv:1710.07629 [quant-ph] .
  * Brun (2019) T. A. Brun, “Quantum error correction,” (2019), arXiv:1910.03672 [quant-ph] .
  * Fujii (2015) K. Fujii,  _Quantum Computation with Topological Codes_, 1st ed., 8, Vol. 1 (Springer Singapore, 2015).
  * Gottesman (1997) D. Gottesman,  _Stabilizer codes and quantum error correction_ , Ph.D. thesis, California Institute of Technology (1997).
  * Kurosh (2014) A. G. Kurosh,  _Theory of groups_ (Amer Mathematical Society, 2014).
  * (29) Another widely used method for VQE that can create entanglement is the universal coupled cluster method (UCC). Here a different Ansatz and parameters are used and the state is evolved via $e^{T(\boldsymbol{\theta})-T(\boldsymbol{\theta})^{\dagger}}\left|\psi_{0}\right\rangle$, where $T=\Sigma_{i,j}\theta_{i}^{j}a_{i}^{\dagger}a_{j}+\Sigma_{i,j,k,l}\theta_{i,j}^{k,l}a_{i}^{\dagger}a_{j}^{\dagger}a_{k}a_{l}$ is the cluster operator, the parametrized sum of all possible one and two excitation operators. Since in this method the entanglement operations are parametrized and we do not consider that type of control, we choose not to include this in our research Fedorov  _et al._ (2021).
  * Saffman (2016) M. Saffman, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 202001 (2016).
  * Adams  _et al._ (2019) C. S. Adams, J. D. Pritchard, and J. P. Shaffer, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 012002 (2019).
  * Jaksch  _et al._ (2000) D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Phys. Rev. Lett. 85, 2208 (2000), arXiv:quant-ph/0004038 [quant-ph] .
  * Morgado and Whitlock (2020a) M. Morgado and S. Whitlock, “Quantum simulation and computing with rydberg-interacting qubits,” (2020a), arXiv:2011.03031 [quant-ph] .
  * Walker and Saffman (2012) T. G. Walker and M. Saffman, in  _Advances in Atomic, Molecular, and Optical Physics_, Advances In Atomic, Molecular, and Optical Physics, Vol. 61, edited by P. Berman, E. Arimondo, and C. Lin (Academic Press, 2012) pp. 81 – 115.
  * Morgado and Whitlock (2020b) M. Morgado and S. Whitlock, “Quantum simulation and computing with rydberg-interacting qubits,” (2020b), arXiv:2011.03031 [quant-ph] .
  * Groenland and Schoutens (2018) K. Groenland and K. Schoutens, Phys. Rev. A 97, 042321 (2018), arXiv:1707.05144 [quant-ph] .
  * Clark  _et al._ (2005) S. R. Clark, C. M. Alves, and D. Jaksch, New Journal of Physics 7, 124 (2005).
  * Clark  _et al._ (2005) S. R. Clark, C. Moura Alves, and D. Jaksch, New Journal of Physics 7, 124 (2005), arXiv:quant-ph/0406150 [quant-ph] .
  * McClean  _et al._ (2018) J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nature Communications 9, 4812 (2018).
  * Smith  _et al._ (2020) D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, _et al._ , The Journal of Chemical Physics 152, 184108 (2020), https://doi.org/10.1063/5.0006002 .
  * Fedorov  _et al._ (2021) D. A. Fedorov, B. Peng, N. Govind, and Y. Alexeev, “Vqe method: A short survey and recent developments,” (2021), arXiv:2103.08505 [quant-ph] .
  * Spall (1998) J. C. Spall, IEEE Transactions on Aerospace and Electronic Systems 34, 817 (1998).
  * Nicholas (2015) P. J. Nicholas, in  _A Dividing Rectangles Algorithm for Stochastic Simulation Optimization_ (2015).
  * Piao Tan and A. Floudas (2009) M. Piao Tan and C. A. Floudas, “Determining the optimal number of clustersdetermining the optimal number of clusters,” in  _Encyclopedia of Optimization_, edited by C. A. Floudas and P. M. Pardalos (Springer US, Boston, MA, 2009) pp. 687–694.
  * Schuld  _et al._ (2019) M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Phys. Rev. A 99, 032331 (2019).
  * Uvarov  _et al._ (2020) A. Uvarov, J. D. Biamonte, and D. Yudin, Phys. Rev. B 102, 075104 (2020).

##  Appendix A Classical Optimization Method SPSA

Alongside the multi-qubit entangling method the classical optimization method
plays an important role in the VQE algorithm, see Fig. 1. This paper uses the
SPSA method which is widely used in VQE research Cerezo  _et al._ (2020);
Kandala  _et al._ (2017); Bharti  _et al._ (2021)

SPSA starts with an input initial trial state $|\Psi(\vec{\theta_{0}})\rangle$
with parameter vector $\vec{\theta_{0}}\in[0,2\pi]^{\otimes D}$. At every
iteration step a small perturbation is done from the parameter vector towards
a random direction to obtain two new parameter vectors
$\vec{\theta}_{k,\pm}$=$\vec{\theta}_{k}+c_{k}\vec{\Delta_{k}}$. Here
$\vec{\Delta_{k}}$ is a $(3d+2)m$-dimensional random vector that is generated
by a sufficiently random generator and $c_{k}$ is a yet to be defined
constant. In this paper $\vec{\Delta_{k}}$ is determined by the Rademacher
distribution where every vector entry has a probability of $1/2$ of being
either 1 or -1. The gradient $\vec{g_{k}}$ is estimated by

$\vec{g_{k}}=\frac{\langle\Psi(\vec{\theta_{k,+}})|H|\Psi(\vec{\theta_{k,+}})\rangle-\langle\Psi(\vec{\theta_{k,-}})|H|\Psi(\vec{\theta_{k,-}})\rangle}{2c_{k}}\vec{\Delta_{k}}.$
(35)

To decrease the energy of the trial state at iteration $k$ a step against the
gradient direction is taken such that

$\vec{\theta}_{k+1}=\vec{\theta}_{k}-a_{k}\vec{g}_{k},$ (36)

where $a_{k}$ is a weight factor dependent on the iteration step. The factors
$a_{k}$ and $c_{k}$ are given by

$a_{k}=\frac{a}{k^{A}},\quad c_{k}=\frac{c}{k^{\Gamma}}.$ (37)

The constants $A,c$ and $\Gamma$ are best chosen to be 0.602, 0.01 and 0.101
respectively to guarantee the smoothest gradient descent Spall (1998). The
constant $a$ is dependent on the gradient around the initial position and is
defined according to

$a=\frac{2\pi}{5}\frac{c}{\left\langle|\langle\Psi(\vec{\theta}_{1,+})|H|\Psi(\vec{\theta}_{1,+})\rangle-\langle\Psi(\vec{\theta}_{1,-})|H|\Psi(\vec{\theta}_{1,-})\rangle|\right\rangle_{\vec{\Delta}_{1}}},$
(38)

where $\left\langle...\right\rangle_{\vec{\Delta}_{1}}$ is the expectation
value over the distribution of $\vec{\Delta}_{1}$. This is the way the average
slope around the initialization vector is calculated in order to calibrate the
optimization process. In practice this expectation value is approximated using
Monte-Carlo methods. A step is taken in the direction opposite of the gradient
in order to, conceivably, reach a lower energy value. The SPSA algorithm is a
local optimization method.

It should be noted that many other classical optimization techniques have been
developed and applied to VQE in recent years such as DIRECTNicholas (2015);
Piao Tan and A. Floudas (2009) and SGD (c.f. Ref. Cerezo  _et al._ (2020) for
a comprehensive overview). In particular the parameter shift rule for
gradients allows for gradients as in Eq. (35) to be determined precisely by
the QPU Schuld  _et al._ (2019). However, since the focus of this work is on
the convergence w.r.t the entanglement methods, we do not discuss further the
details of these other methods.

Recent developments in VQE research have shown the issue of barren plateaus
arising when the number of qubits (and thus the complexity of the problem) is
increased McClean  _et al._ (2018); Uvarov  _et al._ (2020). These plateaus
are large flat parts of the parameter-energy landscape arising in VQE problems
which can greatly decrease the effectiveness of a gradient descend method.
Since our goal is to show how different entanglement methods and optimizers
can influence the convergence of VQE problems even for problems with a small
number of qubits this problem is not encountered in the present work.