Engineering the Cost Function of a Variational Quantum Algorithm for Implementation on Near-Term Devices

The Hamiltonian of the TFIM for a one-dimensional ring of $N$ qubits is given by

where the transverse field direction was chosen for convenience of implementation in superconducting systems (compare to [20] with $\textsf{X}\leftrightarrow\textsf{Z}$). We use natural units (i.e. $\hbar=k_{b}=1$) throughout the manuscript. Consider the special case of generating Thermofield Double (TFD) states in a four-qubit system which was recently demonstrated experimentally, using superconducting qubits by [17]. In this case, the intra-system Hamiltonian reduces to

which describes interactions within each of the two subsystems ($A$ and $B$) with $g$ being the transverse field strength. The ultimate objective is to have the full system undergo unitary evolution such that it will yield a thermal state (or Gibbs state) on subsystem $A$, if it is considered in isolation. In practice, this can be studied by performing a partial trace over subsystem $B$ [11]. Conversely, this technique can be viewed as a purification of the Gibbs state, resulting in a TFD in the full system [19]. The TFD state $\ket{\xi}$ at an inverse temperature $\beta=T^{-1}$ is thus defined as

where $\mathcal{N}$ is a normalization factor, and

$\ket{\xi(0)}$ is the TFD state at $\beta=0$ or $T\rightarrow\infty$, which should be a pairwise maximally entangled Bell state (since tracing subsystem $B$ out of this full state yields a maximally-mixed state for subsystem $A$). For simulation purposes we thus set the initial state to be

where the qubits are labeled as $\ket{A_{1},A_{2},B_{1},B_{2}}$ and qubits $\left\{A_{i},B_{i}\right\}$ are pairwise maximally entangled through applied unitaries.

The protocol for generation of TFD states is described in [20], and we follow a similar path to invoke the variational ansatz motivated by the quantum alternating operator ansatz [5]. This involves alternating between the subsystem Hamiltonian $\mathcal{H}_{A}+\mathcal{H}_{B}$ and the entangling Hamiltonian $\mathcal{H}_{AB}$. Here the subsystem $B$ Hamiltonian is defined as

Near-term quantum computing systems have stringent coherent operations limits and it is desirable to minimize the number of steps required for a given workload [21]. Hence, here we focus on the case of increasing efficiency of single-step TFD generation under the given Hamiltonians (eqs. 2 and 7). With our restricted model, the quantum circuit for TFD generation is described by eq. 8 with $\left\{\alpha_{1},\alpha_{2},\gamma_{1},\gamma_{2}\right\}$ as variational parameters.

For numerical convenience it is typical to utilize the error in fidelity $\mathcal{E}$ with respect to the target as the cost function during optimization [19],

Evaluation of the latter expression requires access to the actual wavefunction of the generated state. Experimental determination of the wavefunction is typically limited to either estimation based on extensive tomographic reconstruction [1], or an array of sophisticated weak measurements [8]. Thus it is undesirable in its present form for evaluation in an actual hybrid quantum-classical experimental system. In principle, it is possible to modify eq. 9 to use the density matrix of the generated state, but this also is cumbersome and reasons are discussed further in the next section.

An alternative to eq. 9 is using the free energy of the system as the cost function during optimization [19]. In this case, the free energy $F_{A}$ is calculated on the reduced density matrix for subsystem $A$ as

where $E_{A}$, $S_{A}$, $\varsigma_{A}=\Tr_{B}\outerproduct{\psi_{\vec{\alpha},\vec{\gamma}}}{\psi_{\vec{\alpha},\vec{\gamma}}}$ are the energy, von Neumann entropy, and the reduced density operator for subsystem $A$, respectively, and $T$ is the system temperature. In this case quantum state tomography (QST) [1] of subsystem $A$ is necessary to calculate $F_{A}$. Recently, methods to approximate the von Neumann entropy have also been proposed [2].

Typically, QST of a system requires a complete set of measurements related to the number of unknowns of the system size [6]. Given a system of $n_{\textrm{q}}$ qubits, the number of unique density matrix elements is given by $\left(2^{2n_{\textrm{q}}}-1\right)$. Thus, this is the total number of unique measurements required to fully characterize the qubit system (e.g. for a two-qubit subsystem, this gives 15 measurements). As system size increases, QST becomes prohibitively expensive and renders free energy $F_{A}$ a poor choice for large-scale optimization problems.

To formulate a more experiment-friendly expression, it is possible to hypothesize a cost function based on the form of $F_{A}$ as described in section III-B. This would involve substituting alternative expressions for $E_{A}$ and $S_{A}$ that require fewer measurements than QST. It is straightforward to infer that $F_{A}$ will be dominated by $S_{A}$ at higher temperatures, and by $E_{A}$ at lower temperatures. Given the simplicity of the TFIM Hamiltonian, the expression for the low temperature regime is easily obtained as follows,

where the $\expectationvalue{\cdot}$ notation indicates the expectation value of the relevant correlator with respect to the evolved wavefunction $\ket{\psi_{\vec{\alpha},\vec{\gamma}}}$ in eq. 8. The reduction in the number of terms is not strictly due to the use of correlators, but it is rather the underlying symmetry of the problem that allows us to determine $E_{A}$ with only three system measurements for subsystem $A$. However, this method of constructing the cost function allows an explicit represention using tangible measurements for the hybrid quantum-classical optimization algorithm.

There is no straightforward method to derive an expression for $S_{A}$ based on correlator expectation values as before. However, it is trivial to verify that $\ket{\xi(0)}$ is the ground state of the negative of the inter-system Hamiltonian, $-\mathcal{H}_{AB}$. Given that we expect $S_{A}$ to dominate at $T\rightarrow\infty$ and $\ket{\xi(0)}$ is the infinite temperature TFD state, it is natural to hypothesize an approximate form of $S_{A}$ to be generalized to $\mathcal{H}_{AB}$ when considering the full system (i.e. both $A$ and $B$) [14]. This results in the following expression for a cost function $\mathcal{C}_{0}(T)$, which is more amenable for practical implementation in a quantum-classical optimization algorithm.

where the system temperature $T$ is a parameter, and $\left\{\alpha_{1},\alpha_{2},\gamma_{1},\gamma_{2}\right\}$ are optimization variables. Here, we have generalized the cost function to the full system and included $\mathcal{H}_{B}$ in the energy of the system at low temperatures.

When the strength of the transverse field is set to $g=1$, a critical point between antiferromagnetic and paramagnetic quantum phases is expected [3, 13]. Hence, for cost function performance comparison purposes, we will primarily explore the case of $g=1$. The case of $g\neq 1$ will be considered for completeness in appendix C. We simulate TFD state generation using Differential Evolution, which is a global optimization algorithm [18] supported by Wolfram Mathematica for non-linear optimization. The optimization is performed over a wide inverse temperature range of six orders of magnitude to ensure complete coverage.

During optimization we minimize the cost functions corresponding to $F_{A}$ and $\mathcal{C}_{0}$ separately, and see that there is excellent agreement between them for extreme low and high temperatures (see fig. 1). We have chosen trace distance $\mathcal{T}$ and fidelity $\mathcal{F}$ as defined below [11], as proximity measures comparing the ideal TFD state and the optimally generated state with the different cost functions.

Here, $\rho=\outerproduct{\xi(\beta)}{\xi(\beta)}$ is the density matrix corresponding to the ideal TFD state, $\sigma=\outerproduct{\Psi(\beta)}{\Psi(\beta)}$ represents the density matrix for the circuit-generated state $\ket{\Psi(\beta)}$ utilizing the relevant cost function, and $\rho_{A}/\sigma_{A}$ are corresponding subsystem states after tracing out $B$, respectively.

For perfect TFD state preparation, we would expect $\mathcal{T}=0$ and $\mathcal{F}=1$. However, in the range $10^{-2}<\beta<10^{2}$, the performance of $\mathcal{C}_{0}$ as a cost function is poor, and demonstrates difficulty in reaching the target TFD state. This is somewhat expected, given that we hypothesized the form of the simple cost function based on extreme high/low temperature behavior. Also note that for $\beta\rightarrow\infty$, $\mathcal{F}\neq 1$ and $\mathcal{T}\neq 0$ indicating that even while using free energy as the cost function we do not construct the ideal TFD state. This is also expected since $T\rightarrow 0$ is the most difficult regime to generate the TFD state based on the given protocol [19]. This indicates that the depth of the circuit is most likely insufficient to construct a better state approaching the TFD state. However depending on the measure used, we find that better approximations to the TFD states are possible using different engineered cost functions.

Given the shortcomings of $\mathcal{C}_{0}$ at intermediate $\beta$ values, we began our new approach to construct a better cost function by generalizing eq. 12 for $g=1$ by adding coefficients to the expression containing correlators as follows,

where

Here $\zeta$ and $\tau$ are parameters which we optimize to find a cost function $\mathcal{C}_{1}$ that can yield better approximations to TFD states across the full inverse temperature range. For simplicity and clarity, instead of performing nested optimizations over $\left\{\zeta,\tau\right\}$ and $\left\{\vec{\alpha},\vec{\gamma}\right\}$, we execute TFD generation for a range of $\zeta$ and $\tau$ values (see appendix A for details). By varying $1.4<\zeta<1.9$, and $1.2<\tau<1.7$, and evaluating the minimization quantity of interest $\Xi$, we see that the best agreement between $\ket{\xi(\beta)}$ and $\ket{\Psi(\beta)}$ is obtained for $\zeta=1.6$ and $\tau=1.48$ (see fig. 2). The improvements from using $\mathcal{C}_{1}$ are discussed in section IV-C.

To further the work towards engineering a better cost function, we used an alternative approach based on density matrix elements’ closeness to generate a cost function that shows significant improvement over both the hypothesized cost function $\mathcal{C}_{0}$, and its enhanced version $\mathcal{C}_{1}$. In this case, we studied the characteristics of the density matrices of the target TFD state, $\rho$, and the single-step circuit-generated state, $\sigma$. Note that we are not tracing out subsystem $B$ of these matrices to obtain the thermal state. Instead we directly compare the density matrices of the target and generated states. Also we relax the assumption of $g=1$, and keep $g$ as a parameter throughout the analysis.

First, we observe the redundancies present in the ideal TFD state and obtain 15 unique real elements for $\rho$:

where the density matrix elements are labeled consistent with the notation in eq. 5 (see appendix B for details). Similarly, we observe the redundancies and Hermiticity in $\sigma$ to obtain 10 unique complex elements for off-diagonal elements, and 5 unique real elements for the diagonal elements:

Inspired by trial optimizations, we find that $\sigma$ is explicitly symmetrized by choosing particular values for the inter-system variational parameters $\vec{\alpha}$,

resulting in 14 unique real elements in $\sigma$. With this choice for $\vec{\alpha}$, it is found that $\sigma_{35}=\sigma_{06}$, indicating that this constrained $\sigma$ will be limited in fully matching $\rho$.

Following symmetrization, a cost function is constructed explicitly by calculating the sum of the square of differences between the density matrix elements for $\rho$ and $\sigma$:

where $r_{i}\in\mathcal{R}$, $s_{i}\in\mathcal{S}$, and $a_{i}\in\{0,1\}$ is a weight used to prune elements. We find that choosing

generates a simple cost function which yields good results across the full temperature range and for different $g$ values. For this system of four qubits, we find that it is beneficial to substitute variational angles in the engineered cost function $\mathcal{C}_{2}$ with the operator expectation values as was the case in $\mathcal{C}_{0}$ and $\mathcal{C}_{1}$. This will enable the experiments to be performed using identical quantum processor measurements, while modifying the classical processor evaluation to improve efficiency. Thus, the non-zero elements for the engineered cost function $\mathcal{C}_{2}$ are explicitly given by eq. 61. Note that given the symmetric nature of the evolution of the subsystems, it is sufficient to include only the intra-system expectation values for one subsystem.

We now evaluate the performance of the various cost functions when generating the TFD states for a wide temperature range. Fidelity and trace distances are calculated for the traced out subsystem as described in section III-D. In fig. 3, we observe that $\mathcal{C}_{1}$ yields vastly superior results compared to the original cost function $\mathcal{C}_{0}$, especially at intermediate temperatures. We also find that $\mathcal{C}_{2}$ performs significantly better than $\mathcal{C}_{1}$ at high temperatures.

For $\beta>1$, we note that the two measures ($\mathcal{F}$ and $\mathcal{T}$) offer somewhat inconsistent results. When considering trace distance $\mathcal{T}$ as the proximity measure, $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ seems to give better results compared to both $\mathcal{C}_{0}$ and $F_{A}$. In appendix C we study a few cases of $g\neq 1$ and find that $\mathcal{C}_{2}$ outperforms all other cost functions (including $F_{A}$) irrespective of the proximity measure used. We attribute this anomaly to the definition of the measures themselves, and further study of this aspect is beyond the scope of this work.