§ INTRODUCTION

§ METHODS

§.§ The Adaptive Derivative-Assembled Pseudo-Trotter Variational Quantum Eigensolver

The adaptive derivative-assembled pseudo-Trotter variational quantum eigensolver (ADAPT-VQE)<cit.> is a VQE-inspired algorithm designed to approximate the ground state wave-function and ground state energy of a given Hamiltonian. Unlike many other classical variational quantum eigensolvers such as the various flavours of trotterised unitary coupled cluster <cit.> however, ADAPT-VQE does not specify a fixed ansatz for the sought-after ground state at the beginning of the algorithm. Instead, ADAPT-VQE functions by first fixing a set of admissible Hermitian generators (the so-called operator pool). The ansatz wave-function is then grown iteratively by parametrically exponentiating a carefully selected Hermitian generator, appending this exponentiated generator to the previous ansatz wave-function, and then variationally tuning the new ansatz wave-function. Since the selection procedure is tailored to the specific Hamiltonian system under consideration (see below), one usually hopes to obtain a more compact ansatz than the one generated by non-adaptive VQEs while still retaining the practical advantages of the VQE for near-term quantum hardware.

The general workflow of the ADAPT-VQE algorithm is as follows. Given the qubit representation of an input Hamiltonian $H$, a pool of admissible Hermitian generators $\mathbb{P}$, and a stopping criterion:

* Boot the qubits to an initial state $\ket{\Psi^{(0)}}$.

* At the $m^{\rm{th}}$ iteration, identify the Hermitian generator $B_m \in \mathbb{P}$ such that the action of the parameterised unitary operator $\exp(-\imath \theta_m B_m), ~\theta_m \in [-\pi, \pi)$ on the current ansatz $\ket{\Psi^{(m-1)}}$ is likely to produce a new wave-function with the largest drop in energy. This identification is done by computing the gradient, at $\theta=0$
, of the expectation value of the Hamiltonian, i.e.,

\begin{align}\label{eq:new_1}
B_m = \argmax_{B\in\mathbb P}  \left \vert\frac{\partial}{{\partial \theta}} \braket{\Psi^{(m-1)}|\exp(\imath \theta B)H\exp(-\imath \theta B)|\Psi^{(m-1)}}\big\vert_{\theta=0} \right \vert.
\end{align}

Note that the criterion (<ref>) is simply a heuristic, and there is no guarantee that the Hermitian generator $B_m$ selected through this criterion will indeed lead to the parameterised unitary operator whose action on the current ansatz $\ket{\Psi^{(m-1)}}$ results in the largest drop in energy. This point will be the subject of further discussion in Section <ref>.

* Exit the iterative process if the stopping criterion is met (see below for more explanation). Otherwise, append the resulting parametrised unitary operator to the left of the current ansatz wave-function $\ket{\Psi^{(m-1)}}$, i.e., define

\begin{align*}
\ket{\widetilde{\Psi^{(m)}}}:=& \exp(-\imath \theta_mB_m)\ket{\Psi^{(m-1)}}
=\exp(-\imath \theta_mB_m)\exp(-\imath \theta_{m-1}'B_{m-1})\ldots \exp(-\imath \theta_1'B_1)\ket{\Psi^{(0)}}.
\end{align*}

* Run a classical VQE routine by optimising all parameters $\theta_m, \theta_{m-1}, \ldots, \theta_1$ in the new ansatz wave-function $ \ket{\widetilde{\Psi^{(m)}}}$ so as to minimize the expectation value of the Hamiltonian, i.e., solve the optimisation problem

\begin{align}
\vec{\theta}^{\rm{opt}} :=& (\theta_1', \ldots, \theta_{m-1}', \theta_m')\\ \nonumber
:=& \underset{\theta_1, \ldots, \theta_{m-1}, \theta_{m}}{\operatorname{argmin}}\Big \langle\Psi^{(0)}\Big| \prod_{k=1}^{k=m} \exp(\imath \theta_kB_k)H\prod_{k=m}^{k=1} \exp(-\imath \theta_kB_k)\Big \vert \Psi^{(0)} \Big \rangle,
\end{align}

and define the new ansatz wave-function $\ket{\Psi^{(m)}}$ using the newly optimized parameters $\theta_1', \ldots, \theta_m'$, i.e., define

\begin{align*}
\ket{\Psi^{(m)}}:= \prod_{k=m}^{k=1}\exp(-\imath \theta_k'B_k)\ket{\Psi^{(0)}} .
\end{align*}

Let us emphasize that although we also denote the newly optimized parameters at the current $m^{\rm th}$ iteration by $\theta_1',\ldots \theta_m'$, these optimized values are not necessarily the same as those used to define $\ket{\Psi^{(m-1)}}$ and referenced in Step 3 above.

* Return to Step 2 with the updated ansatz $\ket{\Psi^{(m)}}$.

Let us remark here that a common choice of stopping criterion is to impose a pre-defined threshold tolerance $\epsilon>0$ on the magnitude of the gradients computed in Step 2 above, i.e., exit the ADAPT-VQE algorithm at iteration $m$ if

\begin{align*}
\max_{B\in\mathbb P}  \left \vert\frac{\partial}{{\partial \theta}} \braket{\Psi^{(m-1)}|\exp(\imath \theta B)H\exp(-\imath \theta B)|\Psi^{(m-1)}}\big\vert_{\theta=0} \right \vert < \epsilon.
\end{align*}

An obvious alternative option is to impose a maximal iteration count on the number of ADAPT-VQE steps or a minimal decrease of the expectation value between two iterates

\begin{equation*}
\Big \langle\Psi^{(0)}\Big| \prod_{k=1}^{k=m-1} \exp(\imath \theta'_kB_k)H\prod_{k=m-1}^{k=1} \exp(-\imath \theta'_kB_k)\Big \vert \Psi^{(0)} \Big \rangle - \Big \langle\Psi^{(0)}\Big| \prod_{k=1}^{k=m} \exp(\imath \theta'_kB_k)H\prod_{k=m}^{k=1} \exp(-\imath \theta'_kB_k)\Big \vert \Psi^{(0)} \Big \rangle< \varepsilon
\end{equation*}

Next, let us discuss some commonly used operator pools for ADAPT-VQE.

§.§ Operator Pools for ADAPT-VQE

As one may expect by studying the ADAPT-VQE workflow, the success of the algorithm is strongly impacted by the choice of the Hermitian generator pool $\mathbb P$. As an extreme example, if all generators in the operator pool commute with the Hamiltonian, then the algorithm will terminate at the first iteration thus resulting in no improvement of the initial guess. The goal of this section is to briefly present some popular operator pools. While a great variety of operator pools have been introduced in the literature <cit.>, we will limit ourselves to a `chemically-inspired' pool which is popular for simulating quantum chemical systems <cit.>, a simplified version of this chemically-inspired pool which (empirically) leads to lower quantum gate counts <cit.>, and finally a so-called minimal operator pool which possesses some useful mathematical properties. Before doing so however, let us first clarify some additional details concerning the ADAPT-VQE algorithm.

ADAPT-VQE was primarily developed for quantum chemistry applications, i.e., for application to molecular systems typically described by a second-quantized molecular Hamiltonian mapped to a qubit representation through the Jordan-Wigner transformation. In this setting, an obvious choice for the initial state in the ADAPT-VQE algorithm is the Hartree-Fock wave-function. Indeed, in the standard formalism where each qubit is used to represent a specific spin-orbital, we can write $\ket{0}_p$ and $\ket{1}_p$ to denote states corresponding to an empty and occupied spin-orbital $p$ respectively. With this notation, the reference Hartree-Fock state for a system having $n$ electrons in $N$ spin-orbitals can be expressed as $\ket{\Psi_{\text{\rm HF}}}:= \ket{1_{0}\ldots 1_{n-1}0_{n}\ldots 0_{N-1}}$, which is straightforward to represent on quantum circuits. Note that this is an example of a case where we have access to an initial state that is both simple to represent on quantum architecture and also yields a wave-function having a reasonable overlap with the sought-after ground state wave-function. Of course, for arbitrary Hamiltonians, such an efficient choice for the initial state might not be possible, in which case we might have to rely on random initialisations, for instance.

The Qubit Excitation-based Pool <cit.>

The first commonly used operator pool is the Qubit excitation-based (QEB) pool which is inspired by the popular coupled cluster method from computational quantum chemistry. The QEB pool consists of so-called single-qubit and double-qubit excitation operators which take the form

\begin{align}\label{eq:single}
\end{align}

\begin{equation}\label{eq:double}
\begin{split}
X_rY_sX_pX_q+Y_rX_sX_pX_q+Y_rY_sY_pX_q+Y_rY_sX_pY_q -X_rX_sY_pX_q-X_rX_sX_pY_q-Y_rX_sY_pY_q-X_rY_sY_pY_q\right).
\end{split}
\end{equation}

Here $p, q, r, $ and $s$ denote qubit indices and $X_p$ and $Y_p$ are the usual one-qubit Pauli gates acting on qubit $p$. Thus, the single-qubit generator $A_{pq}$ acts between the single qubits $p$ and $q$ while the double-qubit generator $A_{pqrs}$ acts between the qubit pairs $(p, q)$ and $(r, s)$.

Given an $N$-qubit system, it is easy to see that the QEB operator pool a priori has $\mathcal{O}(N^4)$ elements. In practice however, not all possible excitation operators of the form (<ref>)-(<ref>) are included in the operator pool. Instead, the QEB operator pool is limited to those single-qubit and double-qubit excitation operators which preserve important symmetries in the system such as spin or the number of particles. Additionally, it is readily checked that the parametric exponentiation of a QEB operator is easy to calculate, and the resulting unitary operators have well-known CNOT-optimised circuits (see, e.g., Figure <ref>
in the appendix).

The Qubit Hardware-efficient Pool <cit.>

While the Qubit excitation-based pool provides excellent performance in numerical simulations on quantum simulators, the practical implementation of QEB-based ansatz wave-functions on near-term quantum hardware remains challenging. This is primarily due to the fact that the number of CNOT gates required to construct the associated QEB circuits, while significantly smaller than the CNOT counts for classical “fixed-ansatz” approaches, is still far too high. The so-called Qubit hardware-efficient pool <cit.> addresses this issue by considering instead a pool consisting of modified single and double excitation operators of the form

\begin{align}\label{eq:single_2}
\widetilde{X}_{pq}=\frac{1}{2}X_qY_p
\end{align}

\begin{equation}\label{eq:double_2}
\begin{split}
\widetilde{X}^{(1)}_{pqrs}=&\frac{1}{8} X_rY_sX_pX_q, \quad \widetilde{X}^{(2)}_{pqrs}=\frac{1}{8} Y_rX_sX_pX_q, \quad \widetilde{X}^{(3)}_{pqrs}=\frac{1}{8} Y_rY_sY_pX_q, \quad \widetilde{X}^{(4)}_{pqrs}=\frac{1}{8} Y_rY_sX_pY_q,\\
\widetilde{X}^{(5)}_{pqrs}=&\frac{1}{8} X_rX_sY_pX_q, \quad \widetilde{X}^{(6)}_{pqrs}=\frac{1}{8} X_rX_sX_pY_q, \quad \widetilde{X}^{(7)}_{pqrs}=\frac{1}{8} Y_rX_sY_pY_q, \quad \widetilde{X}^{(8)}_{pqrs}=\frac{1}{8} X_rY_sY_pY_q,
\end{split}
\end{equation}

where $p, q, r, s$ again denote qubit indices and $X_p$ and $Y_p$ are one-qubit Pauli gates acting on qubit $p$. Note that not all modified double qubit-excitation operators of the form (<ref>) are added to the operator pool since, for instance, the operator $\widetilde{X}^{(1)}_{pqrs}=X_rY_sX_pX_q$ and $\widetilde{X}^{(7)}_{pqrs}=Y_rX_sY_pY_q$ are related by a global rotation <cit.>. Let us also emphasize here that the qubit hardware-efficient pool is not particle conserving and thus violates an important symmetry of the system.

Numerical experiments involving the qubit hardware-efficient pool demonstrate that the resulting ansatz wave-functions can indeed be represented on quantum architectures using a lower CNOT count than that required for the QEB-pool based ansätze. The tradeoff, however, is that the qubit hardware-efficient pool is several times larger than than QEB-pool, and this produces a significant computational overhead when evaluating operator gradients for use in the operator selection criterion. The next operator pool attempts to resolve this issue.

Minimal Hardware-efficient Pool <cit.>

It is well-known that when the Hamiltonian under study is real-valued, then the ground state eigenfunction can be expressed as a real linear combination of real basis vectors. For such settings therefore, it is possible to show the existence of a so-called minimal operator pool that allows the transformation of any real-valued wave-function (in particular, the Hartree-Fock reference state) to another real-valued wave-function (in particular, the sought-after ground state eigenfunction). More precisely, given an $N$-qubit system, we may define the operator pool

\begin{align*}
\mathbb P=\{Y_p\}_{p=0}^{N-2} ~\cup ~\{Z_pY_{p+1}\}_{p=0}^{N-2},
\end{align*}

where $Y_p$ and $Z_p$ are the usual one-qubit Pauli gates acting on qubit $p$. Then, for any two real-valued wave-functions $\ket\Phi$ and $\ket\Psi$, there exist $\theta_1, \ldots, \theta_M \in [-\pi, \pi)$ and Hermitian generators $B_1, \ldots, B_M \in \mathbb{P}$ such that

\begin{align*}
\ket\Phi=\prod_{k=M}^{1}\exp(-\imath\theta_k B_k)\ket\Psi.
\end{align*}

In other words, for a precise choice of parameters and generators in the pool $\mathbb{P}$, we can construct the sought-after ground-state eigenfunction by applying the parametrised, exponentiated generators to the Hartree-Fock reference state.

The key advantage of this so-called minimal hardware-efficient pool $\mathbb P$ is that it consists of only $2N-2$ elements, and the operators in this pool can be parametrically exponentiated using very simple circuits that require a minimal number of CNOT gates.

§.§ Resource Saving Enhancements and the GGA-VQE Algorithm

As described in detail in Section <ref>, a core step in the ADAPT-VQE algorithm is the selection of an optimal Hermitian generator $B$ from the operator pool whose addition to the current ansatz can produce a new ansatz wave-function with the largest drop in energy. Current implementations of ADAPT-VQE make this choice through a heuristic criterion based on evaluating certain gradients of the expectation value of the Hamiltonian. More precisely, for a given pool of operators $\mathbb{P}$, at the $m^{\rm th}$ iteration, one computes (c.f., Equation (<ref>))

\begin{align}\label{eq:new_11}
B_m = \underset{B\in\mathbb P}{\text{\rm argmax}}  \left \vert\frac{\partial}{{\partial \theta}} \braket{\Psi^{(m-1)}|\exp(\imath \theta B)H\exp(-\imath \theta B)|\Psi^{(m-1)}}\big\vert_{\theta=0} \right \vert,
\end{align}

As explained above I suggest to get rid of $m$ in $\theta_m$

\begin{equation*}
B_m = \underset{B\in\mathbb P}{\text{\rm argmax}}  \left \vert\frac{\partial}{{\partial \theta}} \braket{\Psi^{(m-1)}|\exp(\imath \theta B)H\exp(-\imath \theta B)|\Psi^{(m-1)}}\big\vert_{\theta=0} \right \vert,
\end{equation*}

where $\ket{\Psi^{(m-1)}}$ denotes the ansatz wave-function at the ${m-1}$ iteration. Of course, Equation (<ref>) is still a heuristic, and it is conceivable that there exists another Hermitian generator $\widetilde{B}_m$ such that

\begin{align*}
\min_{\theta \in [-\pi, \pi)}  \braket{\Psi^{(m-1)}|\exp(\imath \theta \widetilde{B}_m)H\exp(-\imath \theta \widetilde{B}_m)|\Psi^{(m-1)}}&<\min_{\theta \in [-\pi, \pi)}  \braket{\Psi^{(m-1)}|\exp(\imath \theta {B}_m)H\exp(-\imath \theta {B_m})|\Psi^{(m-1)}}\\
\intertext{while}
\left \vert\frac{\partial}{{\partial \theta}} \braket{\Psi^{(m-1)}|\exp(\imath \theta \widetilde{B}_m)H\exp(-\imath \theta \widetilde{B}_m)|\Psi^{(m-1)}}\big\vert_{\theta=0} \right \vert &< \left \vert\frac{\partial}{{\partial \theta}} \braket{\Psi^{(m-1)}|\exp(\imath \theta {B}_m)H\exp(-\imath \theta {B}_m)|\Psi^{(m-1)}}\big\vert_{\theta=0} \right \vert.
\end{align*}

A representative example of this situation is displayed in Figure <ref>.

Illustration of a situation where the ADAPT-VQE selection criterion does not pick the optimal operator leading to the largest energy drop.

Given that at each iteration of ADAPT-VQE, we face the task of optimising a multi-dimensional objective function which is very noisy due to the quality of the current quantum hardware, the selection of the wrong operator to append to the current ansatz wave-function can be a costly mistake <cit.>. In this section, we introduce an energy sorting algorithm that allows the exact selection of the locally, optimal Hermitian generator from any of the three operator pools that we have introduced in Section <ref>. In other words, we show that it is possible, using a few measurements on the quantum device, to exactly solve the optimisation problem

\begin{align}\label{eq:new_opt}
B_m= \underset{B \in \mathbb{P}}{\text{argmin}}  \min_{\theta \in [-\pi, \pi)} \mathcal{L}(B, \theta, \ket{\Psi^{(m-1)}}) :=  \underset{B \in \mathbb{P}}{\text{argmin}} \min_{\theta \in [-\pi, \pi)}\braket{\Psi^{(m-1)}|\exp(\imath \theta B)H\exp(-\imath \theta B)|\Psi^{(m-1)}},
\end{align}

where, for notational convenience, we have introduced the objective function $\mathcal{L}(B, \theta, \ket{\Psi^{(m-1)}})$, which, in some literature, is also referred to as a landscape function.

The main idea of our energy sorting algorithm is to take advantage of a result that is largely known in the VQE literature on `fixed-ansatz' methods, namely that the expectation value of a Hermitian operator with respect to an ansatz wave-function $\ket{\Psi(\theta)}$ parametrised by a single quantum gate can be expressed in terms of elementary trigonometric functions of $\theta$ <cit.>. While this result is typically expressed for ansatz wave-function parametrised by rotation gates, the core arguments can be extended to wave-functions parametrised by Hermitian generators belonging to any of the three operator pools that we have introduced in Section <ref>. Indeed, taking advantage of the precise functional form of the Hermitian generators introduced in Section <ref>, a simple calculation shows that

* For any generator $B$ in the Qubit-Excitation-Based (QEB) operator pool, it holds that $B^3=B$ with $I$ denoting the identity matrix (see Appendix <ref>).

* For any generator $B$ in the Qubit hardware-efficient and minimal hardware-efficient pools, it holds that $B^2=I$ .

The above simple relations now imply (see also Appendix <ref>) that for any generator $B$ in the Qubit-Excitation-Based (QEB) operator pool and any $\theta \in [-\pi, \pi)$, it holds that

\begin{align}\label{eq:new_has_1}
\exp(-\imath \theta B)=I+(\cos(\theta)-1)B^2-\imath\sin(\theta) B,
\end{align}

and for any generator $B$ in the Qubit hardware-efficient and minimal hardware-efficient pools and any $\theta \in [-\pi, \pi)$, it holds that

\begin{align}\label{eq:new_has_2}
\exp(-\imath \theta B)=\cos(\theta)I-\imath\sin(\theta) B.
\end{align}

Using Equations (<ref>) and (<ref>), it is easy to establish that for any wave-function $\ket{\phi}$ and any $\theta \in [-\pi, \pi)$ the objective function $\mathcal{L}(B, \theta, \ket{\phi})$ has the analytical form

\begin{equation}\label{eq:landscape}
\mathcal{L}(B, \theta, \ket{\phi}) =\begin{cases} \braket{\phi|H|\phi}
+\big(\cos(\theta)-1\big)(\braket{\phi|\{H,B^2\}|\phi}-2\braket{\phi|BHB|\phi} ) \hspace{2cm} \text{if } B^3=B,\\[0.5em]
+\big(1-\cos(\theta)\big)^2\big(\braket{\phi|B^2HB^2|\phi}-\braket{\phi|BHB|\phi} \big)\\[0.5em]
+\sin(\theta)(\cos(\theta)-1)\braket{\phi|\imath B[H,B]B|\phi} \\[0.5em]
\cos^2(\theta)\braket{\phi|H|\phi} +\dfrac{\sin(2\theta)}{2}\braket{\phi|\imath [B,H]|\phi} \hspace{4.15cm}\text{if } B^2=I.\\[0.75em]
\end{cases}
\end{equation}

where $\{\cdot, \cdot\}$ and $[\cdot, \cdot]$ denote the anti-commutator and commutator respectively. A demonstration of this result can be found in Appendix <ref>.

Equation (<ref>) implies that for any Hermitian generator $B$ from our operator pools and any arbitrary wave-function $\ket{\phi}$, the objective function $\mathcal{L}(B, \theta, \ket{\phi})$ can be expressed in terms of elementary trigonometric functions of $\theta$. An important consequence of this expression is that, if we now evaluate the objective function $\mathcal{L}(B, \theta, \ket{\phi})$ at certain well-chosen angles, we can obtain a linear system of equations for the unknown operator expectation values. More precisely,

For the QEB pool ($B^3=B$):

Since $\braket{\phi|H|\phi}$ can be measured directly on the quantum device, four measurements of $ℒ(B, θ, |ϕ⟩)$ at well-chosen $θ= θ^(1), θ^(2), θ^(3), θ^(4)$ yields a linear system for the four unknowns $⟨ϕ|{H,B^2}|ϕ|-⟩2⟨ϕ|BHB|ϕ|$⟩, $\braket{\phi|B^2HB^2|\phi}-\braket{\phi|BHB|\phi}$, $⟨ϕ|iB[H,B]B|ϕ|$⟩ and $\braket{\phi|\imath[B,H]|\phi}$.\vspace{2mm}

\item[\rm \textbf{For the hardware efficient pools ($B^2=I$):}]~

Since $⟨ϕ|H|ϕ|$⟩ can be measured directly on the quantum device, two measurements of $\mathcal{L}(B, \theta, \ket{\phi})$ at well-chosen $\theta= \theta^{(1)}, \theta^{(2)}$ yields a linear system for the two unknowns $\braket{\phi|i[B,H]|\phi}$ and $⟨ϕ|BHB|ϕ|$⟩.

In other words, using a minimal number of measurements and by solving a very small linear system, we can compute all terms involving $B, H$, and $\ket{\phi}$ in the expression (<ref>) for the objective function $\mathcal{L}(B, \theta, \ket{\phi})$. Since the dependency of this function on $\theta$ is through elementary trigonometric functions, we can thus express the objective function $\mathcal{L}(B, \theta, \ket{\phi})$ analytically for any generator $B$, any angle $\theta$, and any wave-function $\ket{\phi}$.
This allows us to solve the optimisation problem (<ref>) up to arbitrary precision for any Hermitian generator $B$ from our operator pool, and thereby obtain the locally, optimal generator that should be added to the current ansatz wave-function $\ket{\Psi^{(m-1)}}$. Note that if we assume an operator pool of size $M$, a total of $4M+1$ measurements will be required to screen all Hermitian generators from the qubit-excitation based pool while a total of $2M+1$ measurements will be required to screen all Hermitian generators from the two hardware efficient pools.

Two important remarks are now in order. First, as mentioned previously, analytical expressions such as (<ref>) for landscape functions $\mathcal{L}(B, \theta, \ket{\phi})$ of the form (<ref>) are present in the existing literature on `fixed-ansatz' methods <cit.>, albeit not- to the best of our knowledge- for the specific operator pools that we have considered in the current study. However, these landscape functions are almost exclusively used to perform an operator-by-operator, analytical optimisation of a structurally fixed, parameterised quantum circuit. Indeed, the only numerical method we are aware of that extends the applicability of analytical landscapes functions beyond simple iterative optimisation of a fixed ansatz circuit is the Rotoselect algorithm <cit.>, and Rotoselect still assumes a fixed structure for the parametrised quantum circuit, in which only the choice of $X, Y,$ or $Z$ rotation gate (and not the `location') can be varied according to the associated landscape function.

Second, let us point out that the landscape function (<ref>) assumes the addition of a single operator to the current ansatz wave-function at each iteration of the adaptive algorithm. If the pool of potential unitary operators is commutative, then the specific order in which operators are chosen is unimportant, and it is therefore sufficient to consider a sequential application of the representation (<ref>) of $\mathcal{L}(B, \theta, \ket{\phi})$ to determine, at each iteration, the optimal operator to append to the current ansatz. On the other hand, if the Hermitian generators belonging to the operator pool do not commute (which is often the case), then the ordering of the operators is important, and it is potentially useful to consider landscape functions based on the simultaneous addition of $d>1$ operators to the current ansatz wave-function at each iteration. We now briefly discuss this generalisation.

Multi-dimensional analytical landscape functions

Let us consider an adaptive procedure in which $d$ unitary operators, constructed using $d$ Hermitian generators from a given operator pool $\mathbb{P}$ are to be appended to the current ansatz wave-function $\ket{\Psi^{(m-1)}}$ at iteration $m$. We are now interested in determining the ordered $d$-tuple of Hermitian generators $(B_{m_d}, \ldots, B_{m_1})$ such that

\begin{align}\label{eq:new_opt_general}
(B_{m_d}, \ldots, B_{m_1})=& \underset{B_d, \ldots, B_1 \in \mathbb{P}} {\text{argmin}}  \hspace{2mm}\min_{\substack{\theta_d, \ldots \theta_1 \\[0.2em] \in [-\pi, \pi)}}\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)\\ \nonumber
:=&  \underset{B_d, \ldots, B_1 \in \mathbb{P}}{\text{argmin}} \hspace{2mm} \min_{\substack{\theta_d, \ldots \theta_1\\[0.2em] \in [-\pi, \pi)}} \braket{\Psi^{(m-1)}|\exp(\imath \theta_1 {B}_1)\ldots \exp(-\imath \theta_d B_d) H\exp(\imath \theta_d B_d)\ldots\exp(-\imath \theta_1 B_1)|\Psi^{(m-1)}}.
\end{align}

In order to obtain an analytical representation of the $d$-dimensional objective function $\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)$, the fundamental idea is to appeal once again to Equations (<ref>) and (<ref>) and expand each exponential in $\theta_j, ~ j \in \{1, \ldots, d\}$ as a sum of a sine and cosine function of $\theta_j$. This expansion allows us to conclude that the  $d$-dimensional objective function $\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)$ can be written in the general form

\begin{align*}\label{eq:landscape_2}
&\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)= \\ \nonumber
&\begin{cases}  \Bigg \langle \Psi^{(m-1)}\Big \vert \underset{j=1}{\overset{j=d}{\prod}} \left(I+(\cos(\theta_j)-1)B_j^2+\imath\sin(\theta_j) B_j\right) H\underset{j=d}{\overset{j=1}{\prod}} \left(I+(\cos(\theta_j)-1)B_j^2-\imath\sin(\theta_j) B_j\right)\Big \vert \Psi^{(m-1)}\Bigg \rangle \hspace{0.7cm} \text{if } B^3=B,\\[2em]
\Bigg \langle \Psi^{(m-1)}\Big \vert \underset{j=1}{\overset{j=d}{\prod}} \left(\cos(\theta_j)I+\imath\sin(\theta_j) B_j\right) H\underset{j=d}{\overset{j=1}{\prod}} \left(\cos(\theta_j)I-\imath\sin(\theta_j) B_j\right)\Big \vert \Psi^{(m-1)}\Bigg \rangle\hspace{4.15cm}\text{if } B^2=I.
\end{cases}
\end{align*}

In other words, the $d$-dimensional objective function $\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)$ can be written as a polynomial of the variables $\big\{1, \cos(\theta_j), \sin(\theta_j)\colon j \in \{1, \ldots, d\} \big\}$ with the exact structure of the polynomial depending on the properties of the operator pool $\mathbb{P}$. As a representative example, the landscape function for hardware efficient pools with $d=2$, after some simplifications, is of the form:

\begin{equation}
\label{eq:landscape_2}
\begin{split}
\mathcal{L}\big((B_2, \theta_2), (B_1, \theta_1), \ket{\phi}\big) =&
\braket{\phi|H|\phi}\\[0.25em]
+&\frac{\cos(2\theta_1)}{2}\Big(\braket{\phi| H-B_1HB_1|\phi} +\frac{\cos(2\theta_2)}2\braket{\phi| H-B_1HB_1-B_2HB_2+B_2B_1HB_1B_2|\phi}\\[0.25em]
&\hspace{1cm}+\frac{\sin(2\theta_2)}2\braket{\phi| i[B_2,H-B_1HB_1]|\phi}\Big)\\
\langle\phi| i[B_1,H]|\phi\rangle+\frac{\cos(2\theta_2)}2\langle \phi|
\end{split}
\end{equation}

Consequently, a total of 7 measurements on a quantum device are required to deduce an analytical expression for the two-dimensional landscape function $\mathcal{L}\big((B_2, \theta_2), (B_1, \theta_1), \ket{\phi}\big)$ for any Hermitian generators $B_1, B_2$ belonging to either of the two hardware efficient operator pools. Since the selection of the best two operator to append to the current ansatz wave-function requires comparing all pairs of Hermitian generators, we conclude that for an operator pool of size $M$, at most $6M^2+1$ measurements are required to determine the locally optimal pair of unitary operators that should be appended to the current ansatz wave-function at each iteration in order to achieve the largest drop in expectation value of the underlying Hermitian operator.

In the case of a general $d$-dimensional objective function $\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)$, similar arguments yield that

* for the Qubit-Excitation-Based (QEB) operator pool of size $M$, we require $\mathcal{O}(5^d M^d)$ measurements to determine the locally optimal $d$-tuple of unitary operators that should be appended to the current ansatz wave-function at each iteration in order to achieve the largest drop in expectation value of the underlying Hermitian operator;

* for the Qubit hardware-efficient and minimal hardware-efficient pools of size $M$, we require $\mathcal{O}(3^d M^d)$ measurements to determine the locally optimal $d$-tuple of unitary operators that should be appended to the current ansatz wave-function at each iteration in order to achieve the largest drop in expectation value of the underlying Hermitian operator.

While this procedure can become computationally intractable for moderately large $d$, various simplifications are possible that can lead to more tractable gradient-free adaptive algorithms involving multi-operator selection and optimisation, and it is likely that such methods offer an advantage when formulating a greedy gradient-free adaptive VQE for a complex Hamiltonian using a non-commutative operator pool. As representative examples, given an ansatz wave-function $\ket{\Psi^{(m-1)}}$, at iteration $m$:

* We may use the energy sorting algorithm based on one-dimensional landscape functions to classify, in descending order of importance, the best $d$ Hermitian generators $(B_d, \ldots, B_1)$ whose addition to the current ansatz wave-function can result in the largest drops in the expectation value of the underlying Hermitian operator, i.e.,

\begin{align*}
\min_{\theta_d \in [-\pi, \pi)}\braket{\Psi^{(m-1)}|\exp(\imath \theta B_d)H\exp(-\imath \theta B_d)|\Psi^{(m-1)}} \geq \ldots \geq \min_{\theta_1 \in [-\pi, \pi)}\braket{\Psi^{(m-1)}|\exp(\imath \theta B_1)H\exp(-\imath \theta B_1)|\Psi^{(m-1)}}.
\end{align*}

We can then switch to the analytical expression of $d$-dimensional landscape function $\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)$ in order to compute the optimal parameters $(\theta_d^*, \ldots, \theta^*_1)$ such that

\begin{align*}
(\theta_d^*, \ldots, \theta_1^*)=&\text{argmin}_{\substack{\theta_d, \ldots \theta_1 \\[0.2em] \in [-\pi, \pi)}}\mathcal{L}\Big( (B_d, \theta_d),\ldots, (B_1, \theta_1), \ket{\Psi^{(m-1)}}\Big)\\
=&\text{argmin}_{\substack{\theta_d, \ldots \theta_1 \\[0.2em] \in [-\pi, \pi)}}\braket{\Psi^{(m-1)}|\exp(\imath \theta_1 {B}_1)\ldots \exp(-\imath \theta_d B_d) H\exp(\imath \theta_d B_d)\ldots\exp(-\imath \theta_1 B_1)|\Psi^{(m-1)}}
\end{align*}

In other words, we may use one-dimensional landscape functions to identity the best $d$ Hermitian generators to add to the current ansatz wave-function and we may employ the $d$-dimensional landscape functions to perform the analytical optimisation. The number of quantum measurements required by this procedure scales as $3^d M$ (resp. $5^d M$) for a hardware efficient (resp. QEB) operator pool of size $M$.

* Taking the newly obtained ansatz wave-function $\ket{\Psi^{(m)}}$ after iteration $m$ as structurally fixed, we may perform Rotoselect-style <cit.> backwards and forwards optimisation sweeps over all parameterised unitary operators $\exp(\imath \theta_1^* {B}_d), \ldots \exp(\imath \theta_d^* {B}_1)$. In particular, thanks to the analytical expression (<ref>) for the $d$-dimensional landscape function, each iteration in these optimisation sweeps can involve $d$ parametrised unitary operators simultaneously. The number of quantum measurements required by a single sweep utilising  $d$-dimensional landscape functions scales as $3^d M$ (resp. $5^d M$) for a hardware efficient (resp. QEB) operator pool of size $M$.

Such considerations will be the subject of numerical investigations in forthcoming works. Since the goal of the current study is the successful implementation of adaptive variational quantum algorithms on current generation NISQ devices, we will, for the moment, limit ourselves to the computationally cheap case of one-dimensional landscape functions, which suffices for the relatively simple Hamiltonians that we consider in the sequel. In this one-dimensional setting, the energy sorting algorithm that we have introduced is the basis of the following greedy gradient-free adaptive VQE, which we dub GGA-VQE.

The Greedy Gradient-free Adaptive Variational Quantum Eigensolver (GGA-VQE)

Given the qubit representation of an input Hamiltonian $H$, a pool of admissible Hermitian generators $\mathbb{P}$, and a stopping criterion:

* Boot the qubits to an initial state $\ket{\Psi^{(0)}}$.

* At the $m^{\rm{th}}$ iteration, use the energy sorting algorithm detailed above to identify the Hermitian generator $B_m \in \mathbb{P}$ that solves the optimisation problem (<ref>), i.e.,

\begin{align}\label{eq:energy_sorting_opt}
B_m=  \underset{B \in \mathbb{P}}{\text{argmin}} \min_{\theta \in [-\pi, \pi)}\braket{\Psi^{(m-1)}|\exp(\imath \theta \widetilde{B})H\exp(-\imath \theta \widetilde{B})|\Psi^{(m-1)}}.
\end{align}

* Exit the iterative process if the stopping criterion is met. Otherwise, append the resulting parametrised unitary operator to the left of the current ansatz wave-function $\ket{\Psi^{(m-1)}}$, i.e., define the new ansatz wave-function

\begin{align*}
\ket{\Psi^{(m)}}:=& \exp(-\imath \theta_m'B_m)\ket{\Psi^{(m-1)}}=\exp(-\imath \theta_m'B_m)\exp(-\imath \theta_{m-1}'B_{m-1})\ldots \exp(-\imath \theta_1'B_1)\ket{\Psi^{(0)}},
\end{align*}

where the angle $\theta_m'$ is obtained in the process of solving the optimisation problem (<ref>).

* Return to Step 2 with the updated ansatz $\ket{\Psi^{(m)}}$.

It is important to emphasise that, in contrast to the classical ADAPT-VQE procedure, the GGA-VQE algorithm described above does not involve a global optimisation of all parameters in the current ansatz at each iteration. Instead, at each iteration, we use the energy sorting algorithm to identify the locally optimal Hermitian generator $B$ as well as the optimal angle $\theta$ which should be used to construct the new ansatz wave-function– a process which involves the optimisation of one-dimensional, elementary trigonometric functions. In particular, in contrast to the classical ADAPT-VQE, we avoid entirely the need to optimise a multi-dimensional and extremely noisy cost function involving the system Hamiltonian. The resulting huge savings in quantum resources suggest that the GGA-VQE is particularly suited for implementation on near-term quantum devices. Let us also note that the $d$-dimensional landscape functions introduced in Section <ref> above can be used to develop natural generalisations of the GGA-VQE algorithm, which we dub GGA-VQE(d), that are likely to be particularly suited for the ground state preparation of strongly correlated systems.

It can now readily be seen that the main computational bottle-neck in the GGA-VQE algorithm is the energy sorting procedure which requires $\mathcal{O}(M)$ measurements of the system Hamiltonian for an operator pool of size $M$. It is therefore natural to ask if the number of measurements required to perform the energy sorting can be further reduced, at least for certain types of Hamiltonians. In the next section, we answer this question affirmatively by showing that for a certain class of Ising Hamiltonians, it is possible to perform energy sorting using a number of measurements independent of both the size of the operator pool and the number of qubits.

§.§ GGA-VQE for a Transverse-field Ising Model

While the ADAPT-VQE algorithm is predominantly applied to compute the ground state energies of molecular systems, there is, in principle, no restriction in applying the method to obtain ground state energies for more general Hamiltonians <cit.>. The goal of this section is to describe in detail, the application of the GGA-VQE algorithm that we have introduced in Section <ref> to an open boundary transverse-field Ising Hamiltonian <cit.>. Ising Hamiltonians of this type are of great importance in condensed-matter physics since they are among the simplest models capable of representing different phases of matter, depending on the value of various systems parameters <cit.>. As the Ising Hamiltonian is well-known theoretically, it also presents a good first test for computational experiments prior to tackling more complex molecular Hamiltonians.

Given an $N$-qubit register, we consider the transverse-field Ising Hamiltonian given by

\begin{align}\label{eq:Ising_Hamiltonian}
H=h\sum_{p=0}^{N-1}X_p+J\sum_{p=0}^{N-2}Z_p Z_{p+1},
\end{align}

where $X_p$ and $Z_p$ denote the usual $Y$ and $Z$ Pauli matrices acting on qubit $p$, and $h, J >0$ are system parameters. The physical constant $h$ models the intensity of a magnetic field directed along the $x$-axis, whereas the constant $J$ models the strength of the nearest-neighbour interactions. If $J<0$, neighbouring spins tend to align, and the opposite is true if $J> 0$. Note that in this model, each qubit represents a spin-state.

Since the Ising Hamiltonian is real valued, a natural choice of operator pool is the minimal hardware-efficient pool introduced in Section <ref>, which is given by

\begin{align}\label{eq:Ising_Pool}
\mathbb P=\{Y_p\}_{p=0}^{N-1} ~\cup ~\{Z_pY_{p+1}\}_{p=0}^{N-1}.
\end{align}

Let us now recall from Section <ref> that implementing the GGA-VQE algorithm requires us to solve, at each iteration, a minimisation problem so as to identify the optimal Hermitian generator which should be used to construct the new ansatz wave-function. The objective function associated with this minimisation problem (see Equation (<ref>)) is given by

\begin{align*}
\mathcal{L}(B, \theta, \ket{\Psi^{(m-1)}})=\braket{\Psi^{(m-1})|\exp(\imath \theta B)H\exp(-\imath \theta B)|\Psi^{(m-1)}},
\end{align*}

where $B \in \mathbb{P}$ is any Hermitian generator from our operator pool, the parameter $\theta \in [-\pi,\pi)$, and $\ket{\Psi^{(m-1)}}$ denotes the previous ansatz wave-function.

It can now be shown (see the Appendix for a detailed demonstration) that for the Ising Hamiltonian defined through Equation (<ref>) and the minimal hardware-efficient pool $\mathbb{P}$ given by (<ref>), the objective function $\mathcal{L}(B, \theta, \ket{\Psi^{(m-1)}})$ has the following simple structure:

\begin{equation}
\mathcal{L}(B, \theta, \ket{\Psi^{(m-1)}}) = \begin{cases}
\Big\langle \Psi^{(m-1)} \Big\vert H \Big\vert \Psi^{(m-1)}\Big \rangle \hspace{7cm}\text{ if }~ B= Y_p,\\[0.5em]
\Psi^{(m-1)}\Big\vert hZ_p-J(X_pZ_{p+1}+Z_{p-1}X_p\delta_{p>0})\Big\vert \Psi^{(m-1)}\Big\rangle\\[0.5em]
-2\sin^2(\theta)\Big\langle\Psi^{(m-1)}\Big\vert hX_p+J(Z_pZ_{p+1}+Z_{p-1}Z_p\delta_{p>0})\Big\vert \Psi^{(m-1)}\Big\rangle.\\[2em]
\Big\langle \Psi^{(m-1)}  \Big\vert H \Big\vert \Psi^{(m-1)} \Big\rangle  \hspace{7cm}\text{ if }~ B= Z_pY_{p+1}. \\[0.5em]
+\sin(2\theta)\Big\langle \Psi^{(m-1)}\Big \vert h(Z_pZ_{p+1}-Y_pY_{p+1})\Big\vert\Psi^{(m-1)}\Big \rangle\\[0.5em]
-\sin(2\theta)\Big\langle \Psi^{(m-1)} \Big\vert J(X_{p+1}+Z_pX_{p+1}Z_{p+2}\delta_{p+2<n})\Big\vert \Psi^{(m-1)}  \Big\rangle\\[0.5em]
-2\sin^2(\theta)\Big\langle \Psi^{(m-1)} \Big\vert hX_p+hX_{p+1}+JZ_pZ_{p+1}+JZ_{p+1}Z_{p+2}\delta_{p+2<n}\Big\vert \Psi^{(m-1)}
\Big\rangle.
\end{cases}
\label{eq:ising_landscape}
\end{equation}

A close study of the right-hand side of Equation (<ref>) now indicates that many terms involving the expectation values of the Pauli matrices can be measured directly and simultaneously on the quantum device without the need to run over all possible Hermitian generators in $\mathbb{P}$.

* The terms containing only tensor products of $Z$ operators can readily be measured in the computational basis.

* The terms containing only tensor products of $X$ (resp. $Y$) operators can be measured by applying a Hadamard (resp. $S^\dag\equiv \text{diag}(1,-\imath)$ and a Hadamard) gate on each qubit.

* The remaining terms are of the form $X_pZ_{p+1}$ or $Z_{p-1}X_{p}Z_{p+1}$. Terms of this form can be measured by applying a Hadamard gate on qubit $p$. The terms corresponding to $p$ even commute and can therefore be measured simultaneously. The same holds true for the $p$ odd terms which can thus also be measured simultaneously.

Consequently, it is possible, at each iteration of the Greedy-ADAPT-VQE algorithm to construct exactly five quantum circuits whose measurements allow us to recover an analytical expression for all objective functions $\mathcal{L}(B, \theta, \ket{\Psi^{(m-1)}}), ~B\in \mathbb{P}$ in terms of elementary trigonometric functions of $\theta$. We have thus achieved a radical reduction in the number of required measurements from $4N-3$ (for the minimal hardware-efficient pool of size $2N-2$) to five.

We end this section by noting that a simple choice of initial state for the energy-sorting frozen core ADAPT-VQE procedure is given by the ground-state of the non-interacting Hamiltonian $\sum_{p=0}^{N-1} X_p$, i.e.,

\begin{align*}
\ket{\Psi^{(0)}}=\ket-^{\otimes n}
\end{align*}

Assuming that the system parameters satisfy $|h|>|J|$, it is not unreasonable to expect the ground state of the true interacting Hamiltonian to be a perturbation of $\ket{\Psi^{(0)}}$.

§.§ Overlap-GGA-VQE for Molecular Systems

Unfortunately, despite the various resource-saving enhancements that we have introduced, the application of the GGA-VQE algorithm to molecular systems of chemical interest seems to be out of reach on current quantum hardware. Indeed, in contrast to the Ising model introduced in Section <ref> for which the expectation value of the Hamiltonian can be be computed through measurements of just two quantum circuits, computing the expectation value of an $N$-spin orbital molecular Hamiltonian typically requires $\mathcal{O}(N^4)$ measurements– a process which introduces an overwhelming amount of hardware and measurement noise. Fortunately, we have at our disposal an alternative and considerably simpler adaptive algorithm that can be used to explore the limits of the current quantum hardware for the simulation of molecular systems. This is the so-called Overlap-ADAPT-VQE algorithm introduced in <cit.>.

Overlap-ADAPT-VQE is a hybrid quantum/classical algorithm in the spirit of ADAPT-VQE which aims to construct compact approximations of target wave-functions through an iterative procedure. As discussed in an earlier contribution <cit.>, Overlap-ADAPT-VQE seeks to improve the construction of adaptive ansatz wave-function for a VQE procedure in the following two ways:

* First, it can be used as a compression strategy to generate more compact ansatz wave-functions which can be represented using shallower quantum circuits and fewer CNOT gates. More precisely, given a target wave-function represented on a quantum device, the Overlap-ADAPT-VQE can generate a compact approximation of this target wave-function. This compact approximation can then be used as an initialisation for a subsequent second adaptive VQE procedure with the aim of further improving the quality of the ansatz whilst reducing the overall quantum circuit-depth.

* Second, Overlap-ADAPT-VQE can be used to generate a high-quality initialisation for a subsequent adaptive VQE algorithm on a quantum device. More precisely, by taking a moderately accurate, classically computed wave-function as the target, Overlap-ADAPT-VQE can produce a high-fidelity, compact approximation of this classical wave-function on a quantum device. This compact approximation can then be used as the initialisation for a subsequent adaptive VQE procedure which helps alleviate the issue of initial barren plateaus.

In contrast to ADAPT, the Overlap-ADAPT procedure does not require the measurement of the expectation value of the Hamiltonian. Instead, at each iteration of Overlap-ADAPT, we measure the overlap between the current ansatz wave-function and the target wave-function to be approximated– a measurement that is much simpler to achieve.
To be more precise, the general workflow of the Overlap-ADAPT-VQE algorithm is as follows.

Given a target wave-function $\ket{\Psi_{\rm ref}}$, a pool of admissible Hermitian generators $\mathbb{P}$, and a maximal operator count $p$:

* Boot the qubits to an initial state $\ket{\Psi^{(0)}}$.

* At the $m^{\rm{th}}$ iteration, identify the Hermitian generator $B_m \in \mathbb{P}$ such that the action of the parameterised unitary operator $\exp(\imath \theta_m B_m), ~\theta_m \in [-\pi, \pi)$ on the current ansatz $\ket{\Psi^{(m-1)}}$ is likely to produce a new wave-function having the largest overlap with the target wave-function. This identification is done by maximising a specific gradient involving the current ansatz wave-function at $\theta_m=0$, i.e.,

\begin{align}\label{eq:new_2}
B_m = \argmax_{B\in\mathbb P} \left \vert \frac{\partial}{{\partial \theta}} \left \vert \left\langle \exp(\imath \theta B)\Psi^{(m-1)}\Big\vert \Psi_{\rm ref}\right \rangle \right \vert^2 \Big\vert_{\theta=0} \right\vert.
\end{align}

Note that, as in the classical ADAPT-VQE procedure, the criterion (<ref>) is a heuristic, and there is no guarantee that the Hermitian generator $B_m$ selected through this criterion will indeed lead to the parametrised unitary operator whose action on the current ansatz $\ket{\Psi^{(m-1)}}$ results in the greatest increase in overlap.

* Append the resulting parametrised unitary operator to the left of the current ansatz wave-function $\ket{\Psi^{(m-1)}}$, i.e., define

\begin{align*}
\ket{\widetilde{\Psi^{(m)}}}:=& \exp(-\imath \theta_mB_m)\ket{\Psi^{(m-1)}}\\
=&\exp(-\imath \theta_mB_m)\exp(-\imath \theta_{m-1}'B_{m-1})\ldots \exp(-\imath \theta_1'B_1)\ket{\Psi^{(0)}}.
\end{align*}

* Run a classical VQE routine by optimising all parameters $\theta_m, \theta_{m-1}, \ldots, \theta_1$ in the new ansatz wave-function $ \ket{\widetilde{\Psi^{(m)}}}$ so as to maximise its overlap with the target wave-function i.e., solve the optimisation problem

\begin{align}
\vec{\theta}^{\rm{opt}} :=& (\theta_1', \ldots, \theta_{m-1}', \theta_m')\\ \nonumber
:=& \underset{\theta_1, \ldots, \theta_{m-1}, \theta_{m}}{\operatorname{argmax}}\left\vert \Big \langle \prod_{k=m}^{k=1} \exp(-\imath \theta_kB_k)\Psi^{(0)}\Big \vert  \Psi_{\rm ref}  \Big\rangle\right \vert^2,
\end{align}

and define the new ansatz wave-function $\ket{\Psi^{(m)}}$ using the newly optimised parameters $\theta_1', \ldots, \theta_m'$, i.e., define

\begin{align*}
\ket{\Psi^{(m)}}:= \prod_{k=m}^{k=1}\exp(-\imath \theta_k'B_k)\ket{\Psi^{(0)}} .
\end{align*}

Let us emphasise that although we also denote the newly optimised parameters at the current $m^{\rm th}$ iteration by $\theta_1',\ldots \theta_m'$, these optimised values are not necessarily the same as those used to define $\ket{\Psi^{(m-1)}}$ and referenced in Step 3 above.

* If the total number of operators in the updated ansatz is equal to $p$, exit the iterative process. Otherwise go to Step 2 with the updated ansatz $\ket{\Psi^{(m)}}$.

It is not difficult to see that the Overlap-ADAPT-VQE procedure can also be viewed as finding the maximizer of a Hamiltonian $H$ given by $H=\ket\Psi_{\rm ref}\bra\Psi_{\rm ref}$. Thus, the formalism developed in Sections <ref>-<ref> can readily be adapted to fit the framework of Overlap-ADAPT-VQE, and in particular, we can define an Overlap-GGA-VQE algorithm. In order to be able to take advantage of the energy sorting algorithm in this setting, we will, additionally, describe how to compute the expectation of this type of Hamiltonian, i.e., how to compute the overlap between two arbitrary states.

The Compute-Uncompute Method

One method to compute the overlap between two states represented on $N$-qubit quantum registers is to use the so-called compute-uncompute method. Indeed, assume we have knowledge of two quantum circuits $U_{\Psi}$ and $U_{\Phi}$ such that $U_{\Psi}\ket{0}=\ket{\Psi}$ and $U_{\Phi}\ket{0}=\ket{\Phi}$, where $\ket0$ denotes the initial (usually Hartree-Fock) state. Then, the overlap $|\braket{\Phi|\Psi}|^2$ can be computed as the expectation value of the projector on the zero state $\ket0\bra0=\left(\frac{I+Z}2\right)^{\otimes n}$ with respect to the state $U_{\Phi}^\dag\ket\Psi$. Indeed, we have

\begin{align*}
|\braket{\Phi|\Psi}|^2=|\braket{0|U_{\Phi}^\dag U_{\Psi}|0}|^2=\bra0 U_{\Psi}^\dag U_{\Phi}\ket0\bra0 U_{\Phi}^\dag U_{\Psi}\ket0.
\end{align*}

The compute-uncompute method has the advantage of not requiring any additional qubits beyond those required to represent the circuits $U_{\Psi}$ and $U_{\Phi}$. It does, however, require combining the individual circuits $U_{\Phi}$ and $U_{\Psi}$ into a single quantum circuit which therefore has twice the depth of the initial circuits.

The Hadamard SWAP-Test

The overlap $|\braket{\Phi|\Psi}|^2$ may also be computed through the so-called SWAP test method for which the associated circuit is shown below. The essential idea of this method is to construct a circuit containing an ancillary qubit such that the probability $p(0)$ of measuring $0$ on the ancillary qubit is related to the overlap through the relation

\begin{align*}
\end{align*}

The SWAP test circuit has the advantage of having the same circuit depth as that of the individual circuits $U_{\Psi}$ and $U_{\Phi}$ representing the states $\ket\Psi$ and $\ket\Phi$ respectively. On the other hand, this efficiency in circuit depth comes at the cost of doubling the number of qubits and requiring $N$ SWAP gates.

\[
\begin{array}{c}
\Qcircuit @C=4.5em @R=1em @!R {
\nghost{{q}_{0} :  } & \lstick{\ket0} & \gate{\mathrm{H}} & \ctrl{1} & \gate{\mathrm{H}} & \meter & \qw & \qw\\
\nghost{{q}_{1} :  } & \lstick{\ket\Phi} & \qw & \qswap & \qw & \qw & \qw & \qw\\
\nghost{{q}_{2} :  } & \lstick{\ket\Psi} & \qw & \qswap \qwx[-1] & \qw & \qw & \qw & \qw\\
\nghost{\mathrm{{c} :  }} & \lstick{\mathrm{{c} :  }} & \lstick{/_{_{1}}} \cw & \cw & \cw & \dstick{_{_{\hspace{0.0em}0}}} \cw \ar @{<=} [-3,0] & \cw & \cw\\}
\end{array}\]

Hadamard SWAP-Test circuit to compute the overlap $|\braket{\Phi|\Psi}|^2$.

§.§ Hybrid Measurements of Quantum Simulations

Before proceeding to the actual numerical results, let us briefly describe the specific outcomes that we wish to explore using QPU implementations of these adaptive algorithms. The primary objective of such adaptive procedures is to yield a wave-function ansatz that accurately represents the ground state of the physical system under study. Our goal in executing such adaptive algorithms on the QPU therefore is to obtain an ordered set of operators (together with corresponding optimal parameters) whose application to the initial state, yield a state that exhibits a high fidelity to the true ground state of the physical system under study. Note that the fidelity of two quantum states $\ket{\psi}$ and $\ket{\phi}$ is defined as the overlap squared of these two states, i.e., $F(\ket{\psi}, \ket{\phi}) = |\langle{\psi}{\ket{\phi}}|^2$.

To achieve the sought-after state preparation, the GGA-VQE algorithm minimises the variational energy of the ansatz wave-function, while the Overlap-GGA-VQE algorithm maximises the overlap (fidelity) of the ansatz with an accurate target state. In order to evaluate the success of our QPU implementations, we must therefore measure the fidelities of the ansatz wave-functions generated by these adaptive algorithms and compare the results obtained on the QPU to those obtained from an HPC simulator. Large differences between the results obtained from the QPU implementations and those obtained from the classical HPC simulations would naturally indicate that the measurement and hardware noise on the quantum device is too great for our algorithm to succeed.

It is crucial to note, however, that any divergence between the QPU and HPC simulator results, while ultimately due to device noise, can nevertheless arise from two different sources:

* The divergence could be due to algorithmic failure. In other words, the noise on the quantum device results in the adaptive algorithm either selecting the wrong parametrised unitary operators to append to the current ansatz, or the wrong `optimal' angles, or both.

* The divergence could be due to ansatz evaluation errors. In other words, while the adaptive algorithm has succeeded (at least to a certain degree) in selecting appropriate parametrised unitary operators and corresponding optimal angles, the actual measurement of the fidelity on the quantum device is a failure due to hardware and measurement noise.

In our opinion, the ansatz evaluation error is largely a reflection of the limits of the current QPU hardware and not indicative of the failure of our adaptive algorithms. Consequently, we will primarily focus on evaluating possible algorithmic failures in the QPU implementations of our adaptive algorithms. An obvious strategy to make this evaluation, is to retrieve the ansatz wave-function yielded by the QPU-implemented GGA-VQE or Overlap-GGA-VQE methods, represent this ansatz wave-function using classical methods on a HPC simulator, and measure the sought-after observables. We refer to this approach as 'hybrid' observable evaluation in the sequel. It is important to emphasise that the construction of the ansatz wave-function, i.e., the choice and order of unitary operators as well as the corresponding angles are all determined by calculations on the QPU. It is only the final measurement of the ansatz wave-function that takes place on the simulator. Of course, it is also of interest to determine the degree of ansatz evaluation errors and we will, where possible, evaluate the energy or fidelity of the generated ansatz wave-function directly on the quantum computer.

A schematic overview of the GGA-VQE algorithm with hybrid observable measurement. We emphasise that the HPC simulator is used only to evaluate the final ansatz wave-function obtained by executing GGA-VQE on the QPU.

§ RESULTS

The algorithmic procedures in this research have been performed using an in-house code developed within the Amazon Braket SDK. All quantum computations were performed using Amazon Braket and have been executed on an IonQ Aria 25-qubits trapped-ion quantum computer which incorporates built-in error mitigation techniques, and each observable is evaluated using 2500 shots. Classical simulations were conducted using our internal Hyperion multi-GPU-accelerated quantum simulator <cit.> with all simulator computations being run on a single NVIDIA DGX A100 node.

§.§ The GGA-VQE algorithm applied to the Ising Model

For our first set of numerical experiments, we apply the GGA-VQE algorithm described in Section <ref> to the transverse-field Ising model described in Section <ref>. We set the system parameters of the Ising Hamiltonian to $h=0.5$, $J=0.2$, which ensures that the two-body interactions in this Ising Hamiltonian play an important role.

Energy convergence of the GGA-VQE algorithm with respect to the number of iterations. The blue reference curve denotes the energy of a classically simulated ansatz. The green and orange curves denote the hybrid and QPU energy evaluations of the GGA-VQE ansatz wave-function respectively. Note that the hybrid evaluation is carried out by retrieving the GGA-VQE ansatz wave-function generated by the QPU, re-implementing it on the Hyperion HPC simulator, and then evaluating the variational energy.

Figure <ref> illustrates the convergence of the hybrid energy evaluations of the GGA-VQE ansatz wave-function with respect to the number of algorithm iterations. We remind the reader that these hybrid energy evaluations are obtained by first running the GGA-VQE algorithm on the IonQ Aria QPU, retrieving the resulting ansatz wave-function and re-implementing it on the Hyperion HPC simulator, and then evaluating the variational energy on the HPC simulator. For reference, we also plot the corresponding energy curve obtained by executing the GGA-VQE ansatz directly on the HPC Hyperion simulator using $10^6$ samples per measured circuit. In addition, as an indication of the measurement and hardware noise, we also plot the GGA-VQE energies obtain by direct measurement on the QPU.

Figure <ref> clearly indicates that QPU-implemented GGA-VQE procedure successfully provides an ansatz wave-function that closely matches the ground state. Moreover, the QPU implementation and HPC simulator implementation of the GGA-VQE algorithm seem highly consistent despite significant noise in the quantum evaluation of observables, as noticeable from the QPU energy evaluation curve in Figure <ref>. Indeed, the greedy, gradient-free operator selection procedure that we have introduced in this study, which relies on a function extrapolation using five noisy evaluations on the QPU, is able to build an ansatz with an energy error below $2.50\times 10^{-2}$ eV and a fidelity exceeding 98% with the exact ground state (see Figure <ref>).

A convergence plot of the fidelity of the GGA-VQE ansatz wave-function produced by the QPU and re-implemented in the Hyperion HPC simulator (hybrid evaluation approach) with the exact ground state of this Ising model obtained using a diagonalisation procedure on the Hyperion simulator. The figure on the right is a zoomed-in version of the figure on the left to better appreciate the fidelity of the GGA-VQE ansatz wave-function.

Expected energy drop of each Hermitian generator from the operator pool during the GGA-VQE iterative procedure on the QPU. The minimal pool operators, numbered from 1 to 48, are listed in the same order as defined in Equation (<ref>). All energy differences are expressed in eV.

To better illustrate the outstanding robustness of the GGA-VQE procedure with respect to QPU noise, we depict in Figure <ref> the expected maximal energy drop of each Hermitian generator from the chosen minimal operator pool throughout the iterative procedure. We observe maximal energy drops of approximately $1.5\times 10^{-2}$ eV for the first 24 iterations followed by a great decrease in the potential energy drops from iteration 25 onwards. This is consistent with the energy curve displayed in Figure <ref> which steadily decreases for the first 24 iterations and then reaches a plateau. It is important to note that while Figure <ref> displays a decrease in the hybrid evaluation of energy between the 24th and 25th iterations, in our opinion, this decrease cannot be attributed to the newly introduced operator. Rather, it is simply a fortuitous consequence of measurement noise. Indeed, the reference simulator curve does not demonstrate any such energy decrease after the 24th iteration. Note that Figure <ref> also suggests a possible explanation for the noise-resilience of the GGA-VQE algorithm. Indeed, we see that the pool of potential unitary operators that can be appended to the current ansatz wave-function can broadly be divided into two camps, namely, operators whose addition to the current ansatz will lead to an approximate energy drop of $1.5\times 10^{-2}$ eV, and operators whose addition to the current ansatz will not meaningfully lower the expectation energy. Thus, the presence of hardware noise can result in a less optimal operator from the first camp being appended to the current ansatz but it will likely not yield a useless operator from the second camp.

For further confirmation, we examine the energy landscapes associated with certain Hermitian generators from the operator pool, extrapolated using five noisy measurements on the QPU. We compare these noisy QPU-based landscapes  with the reference landscapes obtained from the Hyperion HPC simulator. Our results, displayed in Figure <ref>, indicate a nearly perfect match for an operator that enables an energy drop ($Z_0Y_1$), as well as for an operator that does not improve the ansatz ($Y_0$) at the first GGA-VQE iteration. This finding explains the remarkable resilience to QPU noise of operator selection in the GGA-VQE procedure, and suggests that the QPU-implemented algorithm can consistently pick the optimal operator and associated parameter, resulting in a gradual reduction of the variational energy of the ansatz wave-function and thus convergence towards the ground state.

One-dimensional energy landscapes of the operators $Y_0$ and $Z_0Z_1$ when applied to the initial state. The orange curve is extrapolated from five noisy circuit evaluations following the method described in <ref>. The black curve is the exact energy landscape obtained from an HPC simulation. The energy differences are expressed in eV, and the angles are given in radians.

§.§ Overlap-GGA-VQE for a Stretched HF molecule

For our next set of numerical experiments, we consider the application of the Overlap-GGA-VQE algorithm for the approximation of the ground state eigenfunction of the HF molecule at a bond distance of $2.5$ Å. We consider an active space of 8 electrons in 10 spin orbitals in the minimal STO-3g basis set, thus freezing the lowest $1s$ orbital as doubly occupied. The Hartree-Fock state can therefore be represented as $\ket{\phi}=\ket{1111111100}$, which requires 10 qubits. The target wave-function for the Overlap-GGA-VQE process is obtained using a QEB-ADAPT-VQE procedure carried out on the Hyperion HPC simulator until convergence at the chemical accuracy level. The resulting QEB-ADAPT-VQE target wave-function, which has an error of about $1.4$ mHa, is constructed using four generators from the Qubit Excitation-based pool (see Section <ref>) leading to a total CNOT circuit count of 32. The purpose of applying the Overlap-GGA-VQE algorithm is to obtain a high-fidelity approximation of this target wave-function using fewer CNOT gates.

We employ a subset $\mathbb{P}$ of the qubit hardware-efficient pool introduced in Section <ref>. More precisely, we define an index set $P$ for pairs of qubits given by

\begin{align*}
P =\{(4,0),(8,0),(5,1),(9,1),(5,0),(7,0),(7,1)\}.
\end{align*}

Corresponding to this index set $P$, we define the sub-pool $\mathbb{P}$ of qubit hardware-efficient operators as

\begin{align*}
\mathbb P=\{X_{\textbf{p}}= \frac{1}{2}X_{p_2}X_{p_1} \colon \textbf{p}=(p_1, p_2) \in {P}  \}.
\end{align*}

In other words, $\mathbb{P}$ consists of a collection of single excitation qubit hardware-efficient operators (recall Equation (<ref>)).
Equipped with the operator pool $\mathbb{P}$, we apply the Overlap-GGA-VQE algorithm to the target QEB-ADAPT-VQE wave-function.  It is important to note that, for the current HF system, the initial Hartree-Fock state exhibits no overlap with the QEB-ADAPT-VQE target.

Overlap difference induced by each operator from the operator pool $\mathbb{P}$ when applied to the ansatz in the QPU execution (Fig. <ref>), and in a classical simulation (Fig. <ref>). The steps correspond to the iterations in the Overlap-GGA-VQE procedure. At step 2, none of the operators in the pool can significantly improve the overlap with the target state, indicating that convergence has been achieved.

Ansatz fidelity of the Overlap-GGA-VQE ansatz with the target state. The blue bars correspond to the fidelity of a classically simulated ansatz with the target state.  The green and orange bars denote the hybrid and QPU fidelity evaluations respectively. Note that the hybrid evaluation is carried out by retrieving the Overlap-GGA-VQE ansatz wave-function generated by the QPU, re-implementing it on the Hyperion HPC simulator, and then evaluating the variational energy.

Figure <ref> shows the maximal overlap increases that can be achieved by the application of different operators from the pool on the initial Hartree-Fock state. These overlap values are obtained through extrapolation of the one-dimensional objective function defined in Equation (<ref>) for each operator from the pool $\mathbb{P}$. The extrapolations require, for each operator, two noisy overlap evaluations in the QPU. It is readily seen however, that despite the possible hardware and measurement noise on the QPU, the QPU extrapolations are nearly identical to the extrapolations obtained using an HPC simulator implementation of the Overlap-GGA-VQE algorithm displayed in Figure <ref>. This numerical result supports the noise robustness of our operator selection approach. Indeed, we see that in both the QPU implementation and the simulator implementation, the Overlap-GGA-VQE algorithm correctly identifies the operators acting on the qubit pair (5, 0) as the ones that should lead to the largest increase in the overlap, and thus the largest increase in fidelity with the target. For this simple example, convergence is reached after a single iteration as, at the second step, no generator from the operator pool can further improve the overlap with the target state.

The final fidelities of the QPU-implemented and simulator-implemented Overlap-GGA-VQE ansatz wave-functions are plotted in Figure <ref>. We remind the reader that the term `hybrid fidelity evaluation' refers to the classically recomputed fidelity of the QPU-generated ansatz (see the discussion at the start of this section). This hybrid fidelity evaluation precisely matches the value of the fidelity obtained through the pure HPC simulator implementation of Overlap-GGA-VQE, regardless of the choice of overlap measurement technique.

For completeness, we have also plotted the fidelities obtained from the QPU implementation of the Overlap-GGA-VQE procedure through direct measurement on the QPU. In this case, we see that the Swap test exhibits a higher noise level than the compute-uncompute method which indicates that– at least for this hardware– a deeper circuit involving only 10 qubits is less affected by device noise than a shorter circuit that requires gates and connectivity across 20 qubits. Note, however, that in both cases, the hybrid evaluation of the Overlap-GGA-VQE ansatz is highly accurate as indicated in the green bars in Figure <ref>. Indeed, the QPU implementation of the Overlap-GGA-VQE procedure manages to provided an ansatz wave-function that achieves a fidelity of over 99% with a chemically accurate target wave-function while using only 2 CNOT gates.

§ DISCUSSION

§ DATA AVAILABILITY

Data generated during the study is available upon request from the authors (E-mail:

§ CODE AVAILABILITY

The code used during the study is available upon request from the authors (E-mail:

§ ACKNOWLEDGEMENTS

This work has been funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant No 810367), project EMC2 (J.-P. P. and Y.M.). Support from the PEPR EPiQ and HQI programs is acknowledged. We thank Amazon Braket (G. Tourpe) for partial funding of the computations on the Aria IonQ machine and the whole team for their help in the setup of our computations on Amazon Braket's SDK (Software Development Kit). GPU computations have been performed at GENCI (IDRIS, Orsay, France) on grant no A0130712052.

§ APPENDIX