Improved iterative quantum algorithm for ground-state preparation

The inputs in our algorithm are the Hamiltonian $\hat{\mathcal{H}}=\sum_{k=0}^{K-1}h_{k}\hat{\mathcal{H}}_{k}$, a tensor product state $|0^{\otimes\tilde{k}}\rangle$ and the initial state $|\phi(0)\rangle$. A well-decided initial state $|\phi(0)\rangle$ should have a sufficient overlap with the ground state $|u_{0}\rangle$ such that $c_{0}^{(0)}=|\langle u_{0}|\phi(0)\rangle|=\mathcal{O}[\textrm{poly}(1/n)]$ He2021Quantum ; Note1 . Denote $S$ the largest iterative step. Our algorithm iteratively run the following procedures for $s=0,1,\cdots,S-1$.

(a) Performing the unitary operator $U(\boldsymbol{\theta}_{\textrm{opt}})\otimes I_{N}$ on the state $|\Psi(s)\rangle=|0^{\otimes\tilde{k}}\rangle|\phi(s)\rangle$, we obtain the state $|\Psi_{a}(s)\rangle=|\boldsymbol{y}\rangle|\phi(s)\rangle$.

(b) Apply the controlled unitary $\Lambda_{\tilde{k}}G=\sum_{k=0}^{K}|k\rangle\langle k|\otimes G_{k}$ on the whole system. The state $|\Psi_{a}(s)\rangle$ is transformed into the state

(c) Implementing the unitary $U^{{\dagger}}(\boldsymbol{\theta}_{\textrm{opt}})\otimes I_{N}$ on the state $|\Psi_{b}(s)\rangle$, we obtain

where the ancillary state of $|\Psi_{\perp}(s)\rangle$ is orthogonal to $|0\cdots 0\rangle$.

(d) Measure the ancillary system with the projector $M=|0^{\otimes\tilde{k}}\rangle\langle 0^{\otimes\tilde{k}}|\otimes I_{N}$. The state of the whole system after the measurement is

The success probability of preparing the state $|\phi(s+1)\rangle$ is

Furthermore, the following Lemma gives a critical property of the success probability that $P(s+1)$ would not decay exponentially to zero with the iteration (see the proof in Appendix A).

Although QAA provides a quadratic speedup on the measurement sampling complexity, $\Theta[P^{-1/2}(s+1)]$, this step is expensive since the Hamiltonian simulation has a sophisticated circuit structure and requires many copies of pure states $|\Psi_{c}(s)\rangle\langle\Psi_{c}(s)|$ Llyod2014Quantum . Instead, we perform the classical Monte Carlo sampling Montanaro2015Quantum and as a result it needs roughly $\Theta[P^{-1}(s+1)]$ copies of $|\Psi_{c}(s)\rangle$ to generate the state $|\phi(s+1)\rangle$.

We consider two common metrics for the solution quality, the energy and fidelity. The energy is given as an expectation value of the Hamiltonian $\hat{\mathcal{H}}$ on the evolved state $|\phi(s)\rangle$ such that

The largest step $S$ is determined by the convergence criterion such as the energy difference between two consecutive steps, $|E(s)-E(s-1)|\leq\varepsilon$, where $\varepsilon$ is an error threshold Gomes2021Adaptive . The fidelity between two the pure states $|\phi(s)\rangle$ and $|u_{0}\rangle$ is defined as $F(|\phi(s)\rangle,|u_{0}\rangle)=\sqrt{|\langle\phi(s)|u_{0}\rangle|^{2}}$ Nielsen2000Quantum .

In order to enact the non-unitary operator $G$ (Eq. 4) on an arbitrary state $|\phi(s)\rangle$, we consider the implementation of the unitary operation $\Lambda_{\tilde{k}}G=\sum_{k=0}^{K}|k\rangle\langle k|\otimes G_{k}$ in a similar way adopted in Barenco1995Elementary ; Wei2020A ; Liang2022Quantum . Since $\Lambda_{\tilde{k}}G$ is a block diagonal unitary, it can be implemented by applying $K+1$ $\tilde{k}$-fold controlled $n$-qubit gates, $\Lambda_{\tilde{k}}G=\prod_{k=0}^{K}\mathcal{C}_{\tilde{k},|k\rangle}G_{k}$, where the unitary

The index $\tilde{k}$ in $\mathcal{C}_{\tilde{k},|k\rangle}G_{k}$ denotes the number of control qubits, see Fig. (1.c). The index $|k\rangle$ in $\mathcal{C}_{\tilde{k},|k\rangle}G_{k}$ denotes the controlled state. Since the Hamiltonian $\hat{\mathcal{H}}$ is $l$-local, each operator $\mathcal{C}_{\tilde{k},|k\rangle}G_{k}$ can then be decomposed into at most $l$ $\tilde{k}$-fold controlled single-qubit gates $\mathcal{C}_{\tilde{k},|k\rangle}\sigma_{a}^{(b)}$, where $a\in\{x,y,z\}$, $\sigma_{a}^{(b)}$ denotes the Pauli operator $\sigma_{a}$ acting on $b$th qubit of the working system, and the factor $s_{k}$ has been absorbed in one of the operator $\sigma_{a}^{(b)}$. Notice that if $s_{k}=-1$, the operator $s_{k}\sigma_{a}^{(b)}$ can be achieved by the unitary $R_{x}(2\pi)\sigma_{a}^{(b)}$ associated with a rotation operator with respect to the $x$ axe.

For an arbitrary unitary $\mathcal{P}\in\textrm{U}(2)$, the controlled operator

can be decomposed into three controlled unitary operators $\mathcal{C}_{\tilde{k}-1,|(K-1)/2\rangle}Q$, $\mathcal{C}_{1,|1\rangle}Q$ and $\mathcal{C}_{1,|1\rangle}Q^{{\dagger}}$, as well as two Toffolli gates over $\tilde{k}-1$ qubits, where $Q^{2}=\mathcal{P}$ Bergholm2005Quantum . The cost of simulating the two Toffolli gates is $\mathcal{O}(\tilde{k})$ which counts the number of single qubit and CNOT gates Xin2017Quantum . For an arbitrary unitary $\mathcal{C}_{\tilde{k},|k\rangle}\mathcal{P}$, $k=0,\cdots,K$, we utilize the Pauli operator $\sigma_{x}$ to realize the transformation between the operators $\mathcal{C}_{\tilde{k},|k\rangle}\mathcal{P}$ and $\mathcal{C}_{\tilde{k},|K\rangle}\mathcal{P}$. Particularly, given the binary representation $k=k_{1}k_{2}\cdots k_{\tilde{k}}$, we have

where $[\sigma_{x}]=\sigma_{x}^{!k_{1}}\otimes\sigma_{x}^{!k_{2}}\otimes\cdots\sigma_{x}^{!k_{\tilde{k}}}\otimes I_{2}$ with $!k_{i}$ representing the NOT operator that returns $1$ when $k_{i}=0$ and returns $0$ when $k_{i}=1$. Denoting $\mathcal{T}$ the cost of simulating the unitary $\mathcal{C}_{\tilde{k},|k\rangle}\mathcal{P}$, we have $\mathcal{T}=\mathcal{O}(\tilde{k}^{2})$. Thus, the total gate complexity of simulating $\Lambda_{\tilde{k}}G$ is roughly $\mathcal{T}_{\textrm{total}}=l(K+1)\mathcal{T}=\mathcal{O}[l(K+1)\widetilde{k}^{2}]$.

Algorithm 1 presents the detailed process of variational quantum algorithm to prepare the state $|\boldsymbol{y}\rangle$.

In step (1), the parameterized quantum circuit $U(\boldsymbol{\theta})$ is a hardware-efficient ansatz Kandala2017Hardware ; Havlivcek2019Supervised ; Commeau2020Variational . For instance, As shown in Fig. (1b), the PQC consists of single qubit gates $R_{y}(\theta)=e^{-\iota\frac{\theta}{2}\sigma_{y}}$ depending on the tunable parameter $\theta$, the Pauli operator $\sigma_{y}$ and the controlled Pauli gates $\sigma_{z}$ applying to the neighbor qubits. The probability $\boldsymbol{p}_{z}$ is generated by measuring the given trial state multiple times. Let the total number of samples be $N_{\textrm{sample}}$ and the number of samples of state $|z\rangle$ be $N_{z}$. Thus the probability is $\boldsymbol{p}_{z}=N_{z}/N_{\textrm{sample}}$.

In step (2), we define a cost function

where $|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$ and $|\cdot|$ denotes absolute value. The first term

is dubbed as the Kullback-Leibler divergence (KL) which quantifies the amount of information lost when changing from the vector distribution $\boldsymbol{y}^{[2]}=(\boldsymbol{y}_{0}^{2},\cdots,\boldsymbol{y}_{K}^{2})$ to another distribution $\boldsymbol{p}=(\boldsymbol{p}_{0},\cdots,\boldsymbol{p}_{K})$. The second term ensures the obtained states learned by the quantity $F_{\textrm{KL}}(\boldsymbol{\theta})$ have positive local phases. For example, we aim to prepare a single qubit state $|\boldsymbol{y}\rangle=\frac{1}{\sqrt{5}}|0\rangle+\frac{2}{\sqrt{5}}|1\rangle$ corresponding to a classical vector $\boldsymbol{y}=(\frac{1}{\sqrt{5}},\frac{2}{\sqrt{5}})$. Variational optimizing the cost function $F_{\textrm{KL}}$ generates four states,

The desired state $|\boldsymbol{y}_{0}\rangle$ can be obtained by optimizing the second term. Note that the second term can also be

It is clear that $F(\boldsymbol{\theta})$ is positive and is zero if and only if $U(\boldsymbol{\theta}_{\textrm{opt}})|0^{\otimes\tilde{k}}\rangle=|\boldsymbol{y}\rangle$. Much small cost values indicate a high fidelity preparation of the state $|\boldsymbol{y}\rangle$. Thus, the cost function is faithful and operationally meaningful Cerezo2021Variational . Here the inner product $\langle+^{\otimes\tilde{k}}|\psi(\boldsymbol{\theta})\rangle$ is computed via the Hadamard test Aharonov2009A . Updating the circuit parameter by performing a classical optimization procedure such as the Nelder-Mead algorithm Nelder1965A , we produce the optimal parameter $\boldsymbol{\theta}_{\textrm{opt}}$ until the cost function converges. The learning scheme maps the classical vector into a set of parameters $\boldsymbol{\theta}_{\textrm{opt}}$ such that $\boldsymbol{y}\mapsto\{\boldsymbol{\theta}_{\textrm{opt}}\}$. Loading the parameter $\boldsymbol{\theta}_{\textrm{opt}}$ into NISQ devices equipped with a parameterized quantum circuit $U(\boldsymbol{\theta})$, we then produce the state $|\boldsymbol{y}\rangle\approx|\tilde{\boldsymbol{y}}\rangle=U(\boldsymbol{\theta}_{\textrm{opt}})|0^{\otimes\tilde{k}}\rangle$.

In this section, we report strategies to determine the learning rate $\mu$. Theorem 1 provides a theoretical benchmark for analyzing the learning rate $\mu$ (the proof is provided in Appendix B).

We here remark that from Theorem 1 the learning rate $\mu$ can be an arbitrary positive number when the eigenvalues satisfy $E_{N-1}+E_{0}\leq 0$. One may wonder whether one can add diagonal constant shift to the Hamiltonian $\hat{\mathcal{H}}$ to switch the condition $E_{N-1}+E_{0}\leq 0$. Let $\hat{\mathcal{H}}^{{}^{\prime}}=\hat{\mathcal{H}}+\tau I_{N}$ with eigenvalues $E_{i}^{{}^{\prime}}=E_{i}+\tau$ for some constant $\tau\in\mathbb{R}$, such that $E_{N-1}^{{}^{\prime}}+E_{0}^{{}^{\prime}}\leq 0$ but $E_{N-1}+E_{0}>0$. Based on Theorem 1, the corresponding gradient operator is given by $G^{{}^{\prime}}=I_{N}-2\mu^{{}^{\prime}}\hat{\mathcal{H}}^{{}^{\prime}}$ with an arbitrary positive valued learn rate $\mu^{{}^{\prime}}$. However, $E_{N-1}^{{}^{\prime}}+E_{0}^{{}^{\prime}}\leq 0$ means that the shift constant $\tau\leq-\frac{E_{N-1}+E_{0}}{2}$. Therefore, although appropriate shift constant can avoid the selection of the learning rate, the fundamental convergence behavior does not modify.

We also remark that the condition (14) is impractical since the eigenvalues of $\hat{\mathcal{H}}$ are required to be estimated in advance. Thus, we next provide a practical strategy to choose $\mu$ by demonstrating the equivalence between the imaginary time evolution (ITE) and the algorithm proposed in Sec. II.1.

ITE is an iterative computational method to solve the ground-state of many-body quantum systems Aulicino2021State . Consider the time-independent Schrödinger equation in imaginary time ($t\rightarrow-\iota t$),

where $|\phi(t)\rangle$ is a quantum state at time $t$ and $\hat{\mathcal{H}}$ is the Hamiltonian. The formal solution of Eq. (15) can be expressed as $|\phi(t)\rangle=e^{-\hat{\mathcal{H}}t}|\phi(0)\rangle$, where $|\phi(0)\rangle$ is an initial state at time $t=0$. Suppose the quantum state $|\phi(0)\rangle$ is expanded in the eigenbasis of $\hat{\mathcal{H}}$ given in (1),

The (unnormalized) sequence states are given by

in the limitation of large $t$ with non-zero overlap $(c_{0}^{(0)}\neq 0)$ Lehtovaara2006Solution . The trial energy

Namely, the ITE always converges to the ground state after long time iterations.

Consider a small time $\Delta t=t/N_{t}$. The exponential operator $e^{-\hat{\mathcal{H}}t}=e^{-\hat{\mathcal{H}}N_{t}\Delta t}$ and each non-unitary operator $e^{-\hat{\mathcal{H}}\Delta t}$ have good approximations in terms of the following $J$th-order product formula Childs2021Theory ,

where the Taylor series is truncated at order $J$ Berry2015Simulating . In particular, with the first-order truncation of the Taylor series, $e^{-\hat{\mathcal{H}}\Delta t}$ is approximated as a non-unitary operator

with error $\mathcal{O}(\Delta t^{2})$. By representing $\tilde{G}$ as an LCU, the non-unitary operation $\tilde{G}$ can be implemented on a quantum computer Berry2015Simulating ; Long2006General . We repeatedly apply the non-unitary operator $\tilde{G}$ on a random initial state $|\phi(0)\rangle$. When the time $t$ is large enough, the updated state can be seen as an approximate ground state.

We observe that the two operators $\tilde{G}=I_{N}-\Delta t\hat{\mathcal{H}}$ and $G=I_{N}-2\mu\hat{\mathcal{H}}$ are equivalent when the small time step $\Delta t=2\mu$. Thus, in practice we select the learning rate $\mu$ in terms of $\Delta t$ such that $\mu=\Delta t/2$. $\Delta t$ should be sufficient small. Fox example, let the error of the first-order approximation of operator $e^{-\hat{\mathcal{H}}\Delta t}$ be $\epsilon=\Delta t^{2}$. This means that we should select $\mu=\Delta t/2=\sqrt{\epsilon}/2$. In the next section, we find that $\epsilon\leq 10^{-2}$ is an amenable error.

In this subsection, we apply the proposed algorithm to simulate the ground state of two models including the deuteron molecule Dumitrescu2018Cloud ; Shehab2019Toward ; Aydeniz2020Practical and the Heisenberg spin$-1/2$ model with external magnetic field Bespalova2021Hamiltonian .