Filtering variational quantum algorithms for combinatorial optimization

Given an n-qubit Hamiltonian $\mathcal{H}$, we define a filtering operator $F\equiv f(\mathcal{H};\tau )$ via a real-valued function f(E; τ) of the energy E and a parameter τ > 0. We require that the function f²(E; τ) is strictly decreasing on the interval given by the complete spectrum of the Hamiltonian E ∈ [E_min, E_max]. Filtering operators are Hermitian and commute with the Hamiltonian by definition.

When a quantum state |ψ⟩ is sampled in the eigenbasis {|λ_x ⟩: x = 0, 1, ..., 2ⁿ − 1} of a Hamiltonian, then the probability distribution of eigenvectors is given by P_ψ (λ_x ) = |⟨ψ|λ_x ⟩|². The application of the filtering operator on the quantum state produces a new quantum state $\vert F\psi \rangle =F\vert \psi \rangle /\sqrt{{\langle {F}^{2}\rangle }_{\psi }}$ which generates a probability distribution that depends on the energy E_x of each eigenstate |λ_x ⟩:

$$\begin{equation}{P}_{F\psi }({\lambda }_{x})=\frac{{f}^{2}({E}_{x};\tau )}{{\langle {F}^{2}\rangle }_{\psi }}{P}_{\psi }({\lambda }_{x}).\end{equation} \tag{ 1 }$$

For non-eigenstates |ψ⟩ the action of the filtering operator increases the probability of all overlapping eigenstates |λ_x ⟩ ($\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\enspace {P}_{\psi }({\lambda }_{x}) > 0$) for which ${f}^{2}({E}_{x};\tau ) > {\langle {F}^{2}\rangle }_{\psi }$ and decreases the probability of all overlapping eigenstates |λ_y ⟩ for which ${f}^{2}({E}_{y};\tau )< {\langle {F}^{2}\rangle }_{\psi }$. Since f²(E; τ) is strictly decreasing as a function of E, such eigenstates exist for any non-eigenstate |ψ⟩. Hence, under the application of the filtering operator, the probability of sampling eigenstates with low energy increases and the probability of sampling eigenstates with high energy decreases. This naturally leads to a reduction of the average energy. In particular, if the ground state |λ₀⟩ has a finite overlap (P_ψ (λ₀) > 0) with the initial state, the probability of sampling it increases with every application of the filtering operator. After sufficiently many applications the ground state is produced.

Some filtering function definitions are given in figure 2. The Chebyshev filtering operator approximates the Dirac delta operator $\delta (\mathcal{H})$ using the following expansion in terms of Chebyshev polynomials up to order τ [27]:

$$\begin{equation}\delta (\mathcal{H})\approx f(\mathcal{H};\tau )=\sum\limits _{r=0}^{\lfloor \tau /2\rfloor }{(-1)}^{r}\frac{2-{\delta }_{r,0}}{\pi }{g}_{2r}^{(\tau )}{T}_{2r}(\mathcal{H})\end{equation} \tag{ 2 }$$

$$\begin{equation}{g}_{s}^{(\tau )}=\frac{(\tau -s+1)\mathrm{cos}\enspace \frac{\pi s}{\tau +1}+\mathrm{sin}\enspace \frac{\pi s}{\tau +1}\enspace \mathrm{cot}\enspace \frac{\pi }{\tau +1}}{\tau +1}.\end{equation} \tag{ 3 }$$

Here δ_r,s is the Kronecker delta and the Chebyshev polynomials are defined via the recursive formula T_s+1(x) = 2xT_s (x) − T_s−1(x) with T₀(x) = 1, T₁(x) = x.

Zoom In Zoom Out Reset image size

The parameter τ in the filtering operator definition is inspired by the time step parameter of ITE. In fact, for the exponential filtering operator in figure 2 τ is precisely the imaginary time step. This parameter interpolates the action of the filtering operator between two limits. For vanishing values of τ → 0 the filtering operator becomes the identity operator and for sufficiently large values of τ → ∞ the filtering operator becomes a projector onto the ground state.

For our selection of filters we took inspiration from several works. The inverse filter is inspired by the inverse iteration procedure which is a common ingredient in numerical routines that calculate eigenvalues and eigenvectors of matrices [28–30]. The cosine filter was previously used in non-variational algorithms to achieve faster ground state preparation [31] and to analyze finite energy intervals [32]. The Chebyshev filter considered here was also employed in powerful tensor network algorithms for the study of thermalization [33–35].

The QVF algorithm approximates the repeated action of a filtering operator on some initial quantum state by successively optimizing the variational parameters of a parameterized quantum circuit. The algorithm starts at optimization step t = 0 by preparing an initial state |ψ₀⟩ that has finite overlap with the ground state: ${P}_{{\psi }_{0}}({\lambda }_{0}) > 0$. Then the algorithm proceeds iteratively and in each optimization step t ⩾ 1 approximates the state |F_t ψ_t−1⟩—that results from exactly applying the filtering operator F_t to the state |ψ_t−1⟩—by a state |ψ_t ⟩. The subscript t in F_t indicates that the filtering operator can change at each optimization step. The algorithm stops after an initially chosen number of optimization steps.

In order to approximate the application of the filtering operator, we prepare a parameterized quantum circuit ansatz |ψ( θ )⟩ that depends on a vector of m parameters θ = (θ₁, ..., θ_m ). At optimization step t we search for the parameters that minimize the Euclidean distance between the parameterized quantum state and |F_t ψ_t−1⟩:

$$\begin{equation}{\mathcal{C}}_{t}(\boldsymbol{\theta })=\frac{1}{2}{\Vert}\vert \psi (\boldsymbol{\theta })\rangle -\vert {F}_{t}{\psi }_{t-1}\rangle {{\Vert}}^{2}=1-\frac{\mathrm{R}\mathrm{e}\langle {\psi }_{t-1}\vert {F}_{t}\vert \psi (\boldsymbol{\theta })\rangle }{\sqrt{\langle {F}_{t}^{2}{\rangle }_{{\psi }_{t-1}}}}.\end{equation} \tag{ 4 }$$

The final vector of parameters obtained at the end of the minimization of equation (4) defines the quantum state |ψ_t ⟩ ≡ |ψ( θ _t )⟩ at optimization step t. The cost function in equation (4) can be minimized with the help of the Hadamard test, which needs one additional ancilla qubit and several additional controlled operations, as we explain in appendix A.

To avoid the additional quantum resources required by the Hadamard test in QVF, in the following we develop F-VQE. The F-VQE algorithm uses a specific gradient-based procedure that requires essentially the same circuits as VQE.

The partial derivative of the cost function in equation (4) with respect to one parameter θ_j is derived in appendix B.1:

$$\begin{equation}\frac{\partial {\mathcal{C}}_{t}(\boldsymbol{\theta })}{\partial {\theta }_{j}}=-\frac{\mathrm{R}\mathrm{e}\langle {\psi }_{t-1}\vert {F}_{t}\vert \psi (\boldsymbol{\theta }+\pi {\boldsymbol{e}}_{j})\rangle }{2\sqrt{{\langle {F}_{t}^{2}\rangle }_{{\psi }_{t-1}}}}.\end{equation} \tag{ 5 }$$

Here the state |ψ( θ + π e _j )⟩ is produced by the same ansatz circuit except that the vector of angles is shifted by an amount π along the direction e _j of parameter θ_j . If the gradient is evaluated at the current vector of parameters θ _t−1, then the parameter-shift rule [36, 37] yields:

$$\begin{equation}{\left.\frac{\partial {\mathcal{C}}_{t}(\boldsymbol{\theta })}{\partial {\theta }_{j}}\right\vert }_{{\boldsymbol{\theta }}_{t-1}}=-\frac{{\langle {F}_{t}\rangle }_{{\psi }_{t-1}^{j+}}-{\langle {F}_{t}\rangle }_{{\psi }_{t-1}^{j-}}}{4\sqrt{{\langle {F}_{t}^{2}\rangle }_{{\psi }_{t-1}}}}.\end{equation} \tag{ 6 }$$

Here the three circuits |ψ_t−1⟩ and $\vert {\psi }_{t-1}^{j\pm }\rangle \equiv \vert \psi ({\boldsymbol{\theta }}_{t-1}\pm \frac{\pi }{2}{\boldsymbol{e}}_{j})\rangle$ are generated by the ansatz with different parameter vectors. Note that the expectation value in the denominator is the same for all partial derivatives at fixed t.

The F-VQE algorithm takes advantage of this favorable case as follows. At optimization step t, F-VQE performs a single gradient-descent update:

$$\begin{equation}{\boldsymbol{\theta }}_{t}={\boldsymbol{\theta }}_{t-1}-\eta \sum\limits _{j=1}^{m}{\left.\frac{\partial {\mathcal{C}}_{t}(\boldsymbol{\theta })}{\partial {\theta }_{j}}\right\vert }_{{\boldsymbol{\theta }}_{t-1}}{\boldsymbol{e}}_{j},\end{equation} \tag{ 7 }$$

where η > 0 is the learning rate. Then F-VQE moves on to the next cost function ${\mathcal{C}}_{t+1}(\boldsymbol{\theta })$ and proceeds identically. For each optimization step this algorithm requires the evaluation of 2m + 1 circuits.

The expectation value ${\langle {F}_{t}\rangle }_{\psi }$ of filtering operators can be efficiently evaluated by sampling the quantum state in the Hamiltonian eigenbasis. If each eigenstate |λ_x ⟩ is sampled M_x times from a total of M samples, the filtering operator expectation value can be approximated via the Monte Carlo estimator f(E; τ) as:

$$\begin{equation}{\langle {F}_{t}\rangle }_{\psi }\approx \frac{1}{M}\sum\limits _{x}{M}_{x}f({E}_{x};\tau ).\end{equation} \tag{ 8 }$$

In this article, we represent combinatorial optimization problems by diagonal quadratic unconstrained binary optimization (QUBO) Hamiltonians [8–10] for which the eigenbasis is the computational basis and energies can be efficiently computed. Therefore the expectation value of a filtering operator can be approximated by sampling the quantum state in the computational basis.

At each optimization step t the samples used to compute ${\langle {F}_{t}^{2}\rangle }_{{\psi }_{t-1}}$ in equation (6) are also used to compute the average energy ${\langle \mathcal{H}\rangle }_{{\psi }_{t-1}}$. As t increases, the average energy is expected to decrease and the probability of sampling the ground state is expected to increase. Thus F-VQE provides the average energy and a growing chance of sampling a low energy eigenstate or even the ground state at no extra cost during the optimization.

The gradient in F-VQE is equivalent to the one in VQE under certain assumptions. We derive the VQE gradient in appendix B.2. If the Hamiltonian in the VQE gradient of equation (B7) is replaced by −F_t , the new VQE gradient evaluated at the point |ψ( θ )⟩ = |ψ_t−1⟩ coincides with the F-VQE gradient in equation (6) up to a positive multiplicative factor. However, we emphasize that the corresponding VQE cost function −⟨ψ( θ )|F_t |ψ( θ )⟩ is different from the F-VQE cost function in equation (4) where the dependence on the parameters is of the form −Re⟨ψ_t−1|F_t |ψ( θ )⟩. We note that the F-VQE gradient in equation (5), in general, coincides only at the point |ψ( θ )⟩ = |ψ_t−1⟩ with the gradient of the modified VQE cost function. We can also see in equations (B6) and (B8) that the second derivatives do not coincide, not even at that point |ψ( θ )⟩ = |ψ_t−1⟩. Therefore both algorithms explore parameter landscapes with different curvatures.

Both the cost function in equation (4) and its gradient in equation (6) depend on the parameter τ via the expectation value of the filtering operator in equation (8). We dynamically adapt τ to keep the gradient norm as close as possible to some desired large and fixed value at every optimization step. This can prevent the gradient from vanishing and enable us to determine its value more accurately with a fixed number of measurements.

In F-VQE we employ the following heuristic to dynamically adapt τ. At each optimization step the dependence $g(\tau )\equiv {\Vert}\nabla {\left.{\mathcal{C}}_{t}(\tau )\right\vert }_{{\boldsymbol{\theta }}_{t-1}}{{\Vert}}^{2}$ between the gradient norm and τ is used to select the value τ_t that returns a gradient norm as close as possible below a certain threshold g_c > 0. For each optimization step such value τ_t is obtained by solving the implicit equation g(τ_t ) = g_c . Note that g(0) = 0 since for τ = 0 the filtering operator becomes the identity operator and the gradient norm vanishes. For large values of τ the gradient norm saturates at a finite value that is determined by the overlap of the gradient circuits with the ground state. Taking this into account, we select τ_t in the following way. We evaluate the gradient norm for increasing values of τ until either (i) an upper bound τ_u > τ_t is found such that g(τ_u ) > g_c or (ii) g(τ) converges to a constant. In the first case (i) we search for τ_t in the range [0, τ_u ] up to a certain precision. In the second case (ii) we select from the tried values the one that provided a gradient norm closest below the threshold. Note that for the Chebyshev filtering operator only positive integer values of τ are allowed.

We emphasize that our heuristic is different from a simple re-scaling of the gradient by a constant. As shown in figure 3 each partial derivative changes non-trivially as a function of τ. Moreover, in the simulations we observe that the gradient norm has a consistent dependence on τ across different optimization steps and problems.

Zoom In Zoom Out Reset image size

The F-VQE algorithm is based on expectation value computations of filtering operators and it is natural to ask whether this method can benefit from using causal cones. The causal cone of an observable is the quantum circuit composed of only those qubits and gates in the ansatz circuit that have an actual effect on the expectation value. In general, causal cones allow us to simplify the computation of expectation values for local observables. The simplification follows from the fact that outside the causal cone of a local operator unitary gates cancel with their adjoints. Causal cones are a crucial ingredient in various tensor network methods, e.g. based on the multiscale entanglement renormalization ansatz [38]. We have previously used them to make variational quantum algorithms more hardware-efficient for the simulation of the time evolution of quantum many-body systems [23].

Figures 4(b) and (c) show the causal cones for observables with support on two neighboring qubits and two non-neighboring qubits, respectively. Only the qubits and gates inside the causal cone need to be prepared experimentally. Note that for observables with support on distant qubits, the causal cone splits into two separable causal cones that can be independently realized in hardware. Therefore causal cones can reduce the required number of gates and qubits when the observables have small support.

Zoom In Zoom Out Reset image size

Inspecting the filtering operators in figure 2, we see that the exponential, power, and Chebyshev filters can make use of causal cones. The exponential filter $\mathrm{exp}(-\tau \mathcal{H})$ is equivalent to a product of two-local terms that can be processed independently if additional approximations are made [23]. The power filter ${(1-\mathcal{H})}^{\tau }$ for integer values of τ is equivalent to a sum of at most 2τ-local terms so that the entire expectation value can be determined from the sum of simpler expectation values. Similarly the expectation value of the Chebyshev filter can be calculated using the sum of expectation values of at most 2τ-local observables, as the Chebyshev filtering operator is a polynomial in $\mathcal{H}$ of degree τ.

We use 25 random MaxCut instances for each problem size of n ∈ {5, 7, 9, 11, 13, 23} qubits. These problem sizes n correspond to graphs with N = n + 1 vertices. Each instance is defined on a random three-regular weighted simple undirected and connected graph $\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$ where $\mathcal{V}=\left\{1,\enspace 2,\enspace \dots ,\enspace N\right\}$ is the set of vertices, $\mathcal{E}\subset \mathcal{V}\times \mathcal{V}$ is the set of edges between different vertices, and $\mathcal{W}=\left\{{w}_{e}\in [0,1]:\enspace e\in \mathcal{E}\right\}$ is the set of random weights uniformly distributed in the range [0, 1] for all edges. A cut is represented by variables z_v = {+1, −1}, with $v\in \mathcal{V}$, that are +1 for the vertices in one subset of the cut and −1 for the vertices in the other subset of the cut. This formulation has the obvious symmetry of swapping labels +1, −1 for the two subsets. We break this symmetry by assigning +1 to the last vertex v = N and thereby reduce the number of variables to n = N − 1. Then the MaxCut problem consists in solving the optimization problem $({z}_{1}^{\ast },{z}_{2}^{\ast },\dots ,{z}_{n}^{\ast })={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{({z}_{1},{z}_{2},\dots ,{z}_{n})}C({z}_{1},{z}_{2},\dots ,{z}_{n})$, with cost function

$$\begin{equation}C({z}_{1},{z}_{2},\dots ,{z}_{n})=\sum\limits _{e=\left\{u,N\right\}\in \mathcal{E}}{w}_{e}\frac{1-{z}_{u}}{2}+\sum\limits _{e=\left\{u,v\ne N\right\}\in \mathcal{E}}{w}_{e}\frac{1-{z}_{u}{z}_{v}}{2}.\end{equation} \tag{ 9 }$$

The Hamiltonian formulation of the MaxCut problem is obtained via any real coefficients a and b > 0 as

$$\begin{equation}\mathcal{H}=a1-bC({Z}_{1},{Z}_{2},\dots ,{Z}_{n}),\end{equation} \tag{ 10 }$$

where each variable z_u in the MaxCut cost function is replaced by the Pauli operator Z_u acting on qubit u ∈ {1, ..., n}. A ground state of $\mathcal{H}$ in equation (10) is the computational state $\vert {\lambda }_{0}\rangle ={\bigotimes}_{v=1}^{n}\vert (1-{z}_{v}^{\ast })/2\rangle$. The approximation ratio of an n-qubit quantum state |ψ⟩ is defined as

$$\begin{equation}\langle {\alpha \rangle }_{\psi }=\frac{\langle C{({Z}_{1},{Z}_{2},\dots ,{Z}_{n})\rangle }_{\psi }}{{\mathrm{max}}_{({z}_{1},{z}_{2},\dots ,{z}_{n})}C({z}_{1},{z}_{2},\dots ,{z}_{n})}.\end{equation} \tag{ 11 }$$

Before we apply filtering operators to MaxCut Hamiltonians, we re-scale the energy range to [0, 1] using the coefficients a and b in equation (10). To achieve this, we compute lower and upper bounds of the MaxCut cost function in equation (9). We choose the upper bound to be the optimum cost of the semidefinite programming relaxation of the MaxCut problem [41]. We fix the lower bound to the minimum cost 0, which corresponds to the trivial cut z_u = +1 for all u ∈ {1, ..., n}.

F-VQE. This algorithm uses the parameterized quantum circuit shown in figure 4(a). An initial state |ψ₀⟩ = |+⟩^⊗n is prepared by setting to π/2 the parameters in the last rotation on each qubit and setting the remaining parameters to 0. Parameters are iteratively updated using analytical gradient descent as described in section 2.3. At each optimization step, the value of the parameter τ is adapted according to the procedure explained in Section IID. We choose a threshold of g_c = 0.1 for the gradient norm and solve the implicit equation with a precision 0 < g_c − g(τ_t ) < 0.01. For the learning rate η we choose the inverse of the Hessian's diagonal. As shown in appendix B.1, at each optimization step t all diagonal elements of the Hessian have the same value which can be computed by means of the circuit |ψ_t−1⟩. This is a quasi-Newton method [42, 43] that uses only the diagonal of the Hessian matrix and can be realized without additional cost.

HE-ITE. In relation to using causal cones with F-VQE, we have chosen to concentrate our analysis on the exponential filter. In this case the F-VQE method is equivalent to the HE-ITE algorithm called angle update in [23]. We adapt this algorithm to the general QUBO Hamiltonians with long-range interactions considered here. We use the ansatz circuit depicted in figure 4(a) with p = 1 layer. Figures 4(b) and (c) show examples of causal cones in HE-ITE. A total of 70 time steps are performed with a fixed imaginary time step τ = 1.0. For the 23-qubit problems, we choose the number of measurement shots dependent on the number n_cone of qubits in the cone as ${2}^{{n}_{\text{cone}}+2}$.

VQE. For this algorithm we employ the same ansatz as in F-VQE and HE-ITE, shown in figure 4(a). The computation of the analytical gradient for this ansatz requires two quantum circuits per parameter as described in appendix B.2. For the VQE cost function the diagonal of the Hessian can be obtained using the same circuits that are needed for the analytical gradient, similar to F-VQE. Thus it is possible to apply the same quasi-Newton method to the VQE optimization. However, we found that this heuristic performs worse than simply fixing a learning rate for all optimization steps and partial derivatives. More specifically, in the simulations we found that the Hessian diagonal elements frequently vanish so that the parameter update often diverges with this heuristic. The analysis of this phenomenon goes beyond the scope of this work and therefore we choose to fix a learning rate η for all optimization steps and partial derivatives. Comparing the performance of the values η = 1, 0.1 and 0.01 using our five-qubit MaxCut instances, we conclude that η = 1 performs best and hence use this fixed value for η in our simulations and experiments.

QAOA. The parameterized quantum circuit for QAOA is:

$$\begin{equation}\vert \psi (\boldsymbol{\gamma },\boldsymbol{\beta })\rangle =U({\gamma }_{p},{\beta }_{p})\dots U({\gamma }_{1},{\beta }_{1})\vert +{\rangle }^{\otimes n},\end{equation} \tag{ 12 }$$

$$\begin{equation}U({\gamma }_{j},{\beta }_{j})=\mathrm{exp}\left(-i{\gamma }_{j}\sum\limits _{q=1}^{n}{X}_{q}\right)\mathrm{exp}\left(-i{\beta }_{j}\mathcal{H}\right).\end{equation} \tag{ 13 }$$

Here $\vert +\rangle =(1/\sqrt{2})(\vert 0\rangle +\vert 1\rangle )$, X_q denotes the Pauli operator X acting on qubit q, and the ansatz is defined by the m = 2p parameters in the vectors γ = (γ₁, ..., γ_p ) and β = (β₁, ..., β_p ). We initialize the parameters randomly in the range [0, π]. We optimize the parameters using analytical gradient descent. As shown in appendix B.3, the computation of the analytical gradient for the QAOA ansatz requires just the ansatz circuit with various parameter sets. In the QAOA simulations we do not use the quasi-Newton method to determine the learning rate, as this would need additional circuits. Instead we fix a learning rate η for all optimization steps and partial derivatives. We choose η = 1 as it is the best performing learning rate for the five-qubit MaxCut instances where we compared the performance of η = 1, 0.1, and 0.01.

In the following, we present and analyze the numerical and experimental results of the F-VQE algorithms and compare them with those of HE-ITE, VQE and QAOA. For each MaxCut instance we pay special attention to two benchmark quantities: the approximation ratio and the probability of measuring the ground state. These quantities have some dependence. A large probability of sampling the ground state P_ψ (λ₀) ≈ 1 implies a large approximation ratio ⟨α⟩_ψ ≈ 1. However, a quantum state |ψ⟩ given by a superposition of low-energy excited states can exhibit a large approximation ratio ⟨α⟩_ψ ≈ 1 but low ground state probability P_ψ (λ₀) ≈ 0.

We compare the performance of various filtering operators in F-VQE in figure 5. Here filtering operators are sorted from left to right by performance and the inverse filter is the best performing one. Figure 5(a) shows for each algorithm considered in this article the final approximation ratio averaged over all 25 MaxCut instances for each problem size. We observe that the best performing filters achieve the largest approximation ratios and are more reliable in obtaining such approximation ratios, which can be gathered from the small values of the corresponding standard deviations. Figure 5(b) shows the distribution of the minimum number of optimization steps required to achieve an approximation ratio above 0.75. Here all filters show a similar performance: they require five or fewer optimization steps with little deviation from the median. Figure 5(c) shows the fraction of MaxCut instances where the algorithms obtain a probability of sampling the ground state above 0.25. The probability of sampling it at least once with M measurement shots is then 1 − (1 − 0.25)^M . The best filters achieve this probability for a larger fraction of MaxCut instances.

Zoom In Zoom Out Reset image size

Let us now compare the best performing filtering operator, the inverse filter, with VQE and QAOA. F-VQE requires less optimization steps than VQE and QAOA to achieve larger and more consistent approximation ratios. This can be seen in figures 5(a) and (b) as well as in figure 1(a). Additionally, as shown in figure 5(c), F-VQE converges to the ground state more often for all problem sizes.

The HE-ITE algorithm also achieves approximation ratios close to the optimum and frequently converges to the ground state. Moreover, the evolution of the approximation ratio shown in figure 1(a) almost overlaps with that for F-VQE. Similarly figure 1(b) shows that HE-ITE obtains an average approximation ratio close to the optimum for 25 MaxCut instances of 23 qubits. Importantly, the quantum circuits never use more than six qubits. The inset in figure 1(b) shows the average fraction of circuits used for each qubit count. We observe that the majority of circuits require just four qubits.

The performance of HE-ITE depends on the imaginary time step τ and to analyze it we have run HE-ITE for a long total imaginary time of 100. Figure 6 compares for three values of τ the final approximation ratios and the fraction of MaxCut instances where a probability of sampling the ground state above 0.25 is achieved. We conclude that HE-ITE improves systematically by choosing smaller values of τ.

Zoom In Zoom Out Reset image size

To demonstrate the experimental feasibility, we run F-VQE with the inverse filter on the Quantinuum H1 trapped-ion quantum processor and solve a random nine-qubit MaxCut instance. Experimental settings (see table 1(c)) are similar to those used for the numerical simulations, with some differences described in detail in appendix C. The main difference is that we use only p = 1 layer and remove all rotation gates except the first and the last ones for each qubit. Figure 7 shows the approximation ratio and the probability of measuring the ground state at each optimization step. The final approximation ratio is 0.9844 ± 0.0062 and the probability of sampling the ground state after the final optimization step is 0.928 ± 0.024. Here the value after ± indicates a 95% confidence interval.

Zoom In Zoom Out Reset image size

In the quantum ansatz circuits |ψ( θ )⟩ of figure 4(a) parameters θ = (θ₁, θ₂, ..., θ_m ) are present only in rotation gates of the form R_G (θ_j ) = exp(−iθ_j G/2), where G is a single-qubit Pauli operator or a tensor product of Pauli operators. As a function of a single parameter, the circuit can be expressed in terms of two fixed unitaries V_A , V_B and the rotation gate corresponding to the parameter:

$$\begin{equation}\vert \psi ({\theta }_{j})\rangle ={V}_{A}{R}_{G}({\theta }_{j}){V}_{B}\vert 0\rangle ,\end{equation} \tag{ B1 }$$

where |0⟩ = |0⟩^⊗n is the initial register state. Hence the first and second derivatives with respect to a parameter θ_j are:

$$\begin{equation}\frac{\partial \vert \psi ({\theta }_{j})\rangle }{\partial {\theta }_{j}}=\frac{1}{2}{V}_{A}{R}_{G}({\theta }_{j})(-iG){V}_{B}\vert 0\rangle =\frac{1}{2}\vert \psi ({\theta }_{j}+\pi )\rangle ,\end{equation} \tag{ B2 }$$

$$\begin{equation}\frac{{\partial }^{2}\vert \psi ({\theta }_{j})\rangle }{\partial {\theta }_{j}^{2}}=\frac{1}{4}{V}_{A}{R}_{G}({\theta }_{j}){(-iG)}^{2}{V}_{B}\vert 0\rangle =-\frac{1}{4}\vert \psi ({\theta }_{j})\rangle ,\end{equation} \tag{ B3 }$$

where we have used that −iG = R_G (π).

Then the first and second derivatives of the cost function in equation (4) are:

$$\begin{equation}\frac{\partial {\mathcal{C}}_{t}(\boldsymbol{\theta })}{\partial {\theta }_{j}}=-\frac{\mathrm{R}\mathrm{e}\langle {\psi }_{t-1}\vert {F}_{t}\vert \psi (\boldsymbol{\theta }+\pi {\boldsymbol{e}}_{j})\rangle }{2\sqrt{{\langle {F}_{t}^{2}\rangle }_{{\psi }_{t-1}}}},\end{equation} \tag{ B4 }$$

$$\begin{equation}\frac{{\partial }^{2}{\mathcal{C}}_{t}(\boldsymbol{\theta })}{\partial {\theta }_{j}^{2}}=\frac{\mathrm{R}\mathrm{e}\langle {\psi }_{t-1}\vert {F}_{t}\vert \psi (\boldsymbol{\theta })\rangle }{4\sqrt{{\langle {F}_{t}^{2}\rangle }_{{\psi }_{t-1}}}}.\end{equation} \tag{ B5 }$$

The first equation corresponds to equation (5). Since the filtering operator F_t is Hermitian, when the first equation is evaluated for the vector of parameters θ _t−1 that produces the state |ψ_t−1⟩ = |ψ( θ _t−1)⟩ we can use the parameter-shift rule to express the numerator as a sum of two circuits, which leads to equation (6). When the second equation is evaluated for θ _t−1, it results in:

$$\begin{equation}{\left.\frac{{\partial }^{2}{\mathcal{C}}_{t}(\boldsymbol{\theta })}{\partial {\theta }_{j}^{2}}\right\vert }_{{\boldsymbol{\theta }}_{t-1}}=\frac{\langle {F}_{t}{\rangle }_{{\psi }_{t-1}}}{4\sqrt{\langle {F}_{t}^{2}{\rangle }_{{\psi }_{t-1}}}},\end{equation} \tag{ B6 }$$

which requires only the quantum circuit for |ψ_t−1⟩.

The cost function is the average energy ${\langle \mathcal{H}\rangle }_{\psi (\boldsymbol{\theta })}$ and the ansatz circuit |ψ( θ )⟩ is the same as the one employed by F-VQE shown in figure 4(a). Given that each parameter is present in only one rotation gate, the parameter-shift rule can be applied to express the analytical gradient as:

$$\begin{equation}{\left.\frac{\partial {\langle \mathcal{H}\rangle }_{\psi (\boldsymbol{\theta })}}{\partial {\theta }_{j}}\right\vert }_{{\boldsymbol{\theta }}_{t-1}}=\frac{1}{2}\left({\langle \mathcal{H}\rangle }_{{\psi }_{t-1}^{j+}}-{\langle \mathcal{H}\rangle }_{{\psi }_{t-1}^{j-}}\right).\end{equation} \tag{ B7 }$$

As with F-VQE, the circuits $\vert {\psi }_{t-1}^{j\pm }\rangle \equiv \vert \psi \left({\boldsymbol{\theta }}_{t-1}\pm \frac{\pi }{2}{\boldsymbol{e}}_{j}\right)\rangle$ are implemented by shifting the parameter θ_j by an amount ±π/2.

Using the same methods, the second derivative with respect to each parameter can be evaluated as

$$\begin{equation}{\left.\frac{{\partial }^{2}\langle {\mathcal{H}\rangle }_{\psi (\boldsymbol{\theta })}}{\partial {\theta }_{j}^{2}}\right\vert }_{{\boldsymbol{\theta }}_{t-1}}=\frac{1}{2}\langle {\mathcal{H}\rangle }_{{\psi }_{t-1}^{j++}}-\frac{1}{2}\langle {\mathcal{H}\rangle }_{{\psi }_{t-1}},\end{equation} \tag{ B8 }$$

where |ψ_t−1⟩ = |ψ( θ _t−1)⟩ is the state with no shifts, and $\vert {\psi }_{t-1}^{\enspace j++}\rangle =\vert \psi \left({\boldsymbol{\theta }}_{t-1}\pm \pi {\boldsymbol{e}}_{j}\right)\rangle$ is the state with a + π shift.

The QAOA ansatz in equation (12) is not of the previous form: here parameters multiply sums of tensor products of Pauli operators so that the partial derivatives become sums of circuit evaluations. Given a Hamiltonian $\mathcal{H}={\sum }_{k=1}^{K}{h}_{k}{Z}_{{Q}_{k}}$ with real coefficients h_k and ${Z}_{{Q}_{k}}={\bigotimes}_{q\in {Q}_{k}}{Z}_{q}$, the ansatz derivatives are:

$$\begin{equation}\frac{\partial \vert \psi (\boldsymbol{\gamma },\boldsymbol{\beta })\rangle }{\partial {\gamma }_{j}}=\sum\limits _{q=1}^{n}{\tilde{U}}_{(j,p)}(-i{X}_{q}){\tilde{U}}_{(1,j-1)}\vert + {\rangle }^{\otimes n},\end{equation} \tag{ B9 }$$

$$\begin{equation}\frac{\partial \vert \psi (\boldsymbol{\gamma },\boldsymbol{\beta })\rangle }{\partial {\beta }_{j}}=\sum\limits _{k=1}^{K}{h}_{k}{\tilde{U}}_{(j+1,p)}(-i{Z}_{{Q}_{k}}){\tilde{U}}_{(1,j)}\vert + {\rangle }^{\otimes n}.\end{equation} \tag{ B10 }$$

Here ${\tilde{U}}_{a,b}=U({\boldsymbol{\gamma }}_{b},{\boldsymbol{\beta }}_{b})\dots U({\boldsymbol{\gamma }}_{a+1},{\boldsymbol{\beta }}_{a+1})U({\boldsymbol{\gamma }}_{a},{\boldsymbol{\beta }}_{a})$, where U( γ _j , β _j ) is defined in equation (13), contains all QAOA circuit layers from a to b. Therefore the partial derivatives of the QAOA cost function are:

$$\begin{equation}\frac{\partial \langle {\mathcal{H}\rangle }_{\psi (\boldsymbol{\gamma },\boldsymbol{\beta })}}{\partial {\gamma }_{j}}=2\mathrm{R}\mathrm{e}\langle \psi (\boldsymbol{\gamma },\boldsymbol{\beta })\vert \mathcal{H}\frac{\partial \vert \psi (\boldsymbol{\gamma },\boldsymbol{\beta })\rangle }{\partial {\gamma }_{j}},\end{equation} \tag{ B11 }$$

$$\begin{equation}\frac{\partial \langle {\mathcal{H}\rangle }_{\psi (\boldsymbol{\gamma },\boldsymbol{\beta })}}{\partial {\beta }_{j}}=2\enspace \mathrm{R}\mathrm{e}\langle \psi (\boldsymbol{\gamma },\boldsymbol{\beta })\vert \mathcal{H}\frac{\partial \vert \psi (\boldsymbol{\gamma },\boldsymbol{\beta })\rangle }{\partial {\beta }_{j}}.\end{equation} \tag{ B12 }$$

If the parameter-shift rule is applied to every term individually, we obtain the analytical gradient for the QAOA ansatz in terms of ansatz circuits for various parameter sets:

$$\begin{equation}\frac{\partial \langle {\mathcal{H}\rangle }_{\psi (\boldsymbol{\gamma },\boldsymbol{\beta })}}{\partial {\gamma }_{j}}=\sum\limits _{q=1}^{n}\left(\langle {\mathcal{H}\rangle }_{{\psi }_{X}^{(j,q)+}}-\langle {\mathcal{H}\rangle }_{{\psi }_{X}^{(j,q)-}}\right),\end{equation} \tag{ B13 }$$

$$\begin{equation}\frac{\partial \langle {\mathcal{H}\rangle }_{\psi (\boldsymbol{\gamma },\boldsymbol{\beta })}}{\partial {\beta }_{j}}=\sum\limits _{k=1}^{K}{h}_{k}\left(\langle {\mathcal{H}\rangle }_{{\psi }_{\mathcal{H}}^{(j,k)+}}-\langle {\mathcal{H}\rangle }_{{\psi }_{\mathcal{H}}^{(j,k)-}}\right).\end{equation} \tag{ B14 }$$

The evaluation of all gradient components requires the following 2p(n + K) circuits that are defined by inserting a rotation gate in the ansatz circuit:

$$\begin{equation}\vert {\psi }_{X}^{(j,q)\pm }\rangle ={\tilde{U}}_{(j,p)}{R}_{{X}_{q}}\left(\pm \frac{\pi }{2}\right){\tilde{U}}_{(1,j-1)}\vert +{\rangle }^{\otimes n},\end{equation} \tag{ B15 }$$

$$\begin{equation}\vert {\psi }_{\mathcal{H}}^{(j,k)\pm }\rangle ={\tilde{U}}_{(j+1,p)}{R}_{{Z}_{{Q}_{k}}}\left(\pm \frac{\pi }{2}\right){\tilde{U}}_{(1,j)}\vert +{\rangle }^{\otimes n}.\end{equation} \tag{ B16 }$$