# Computational Supremacy of Quantum Eigensolver by Extension of Optimized
Binary Configurations

Hayun Park  Department of Liberal Studies, Kangwon National University,
Samcheok, 25913, Republic of Korea  Hunpyo Lee  Email<EMAIL_ADDRESS>Department of Liberal Studies, Kangwon National University, Samcheok, 25913,
Republic of Korea  Quantum Sub Inc., Samcheok, 25913, Republic of Korea

###### Abstract

We developed a quantum eigensolver (QE) which is based on an extension of
optimized binary configurations measured by quantum annealing (QA) on a D-Wave
Quantum Annealer (D-Wave QA). This approach performs iterative QA measurements
to optimize the eigenstates $|\psi\rangle$ without the derivation of a
classical computer. The computational cost is $\eta ML$ for full eigenvalues
$E$ and $|\psi\rangle$ of the Hamiltonian $\hat{H}$ of size $L\times L$, where
$M$ and $\eta$ are the number of QA measurements required to reach the
converged $|\psi\rangle$ and the total annealing time of many QA shots,
respectively. Unlike the exact diagonalized (ED) algorithm with $L^{3}$
iterations on a classical computer, the computation cost is not significantly
affected by $L$ and $M$ because $\eta$ represents a very short time within
$10^{-2}$ seconds on the D-Wave QA. We selected the tight-binding $\hat{H}$
that contains the exact $E$ values of all energy states in two systems with
metallic and insulating phases. We confirmed that the proposed QE algorithm
provides exact solutions within the errors of $5\times 10^{-3}$. The QE
algorithm will not only show computational supremacy over the ED approach on a
classical computer but will also be widely used for various applications such
as material and drug design.

###### pacs:

71.10.Fd,71.27.+a,71.30.+h

_Introduction_ —Fast computation to find the full eigenvalues $E$ and
eigenstates $|\psi\rangle$ of a Hermitian matrix $\hat{H}$ of size $L\times L$
is a highly desirable technique which has been widely applied across various
applications, including many-body electronic systems, electronic structures,
novel material design and drug design. The full values of $E$ and
$|\psi\rangle$ can be computed using an exact diagonalization (ED) method,
where the computational cost demands iterations of $L^{3}$ on a classic
computer based on bits Weise2008 . Another approach using a classic computer
is to employ an optimization approach. The lowest eigenvalue
$E^{\text{Lowest}}$ and its $|\psi\rangle$ can be determined using the
optimization equation
$E^{\text{Lowest}}=\text{min}_{|\psi\rangle}\langle\psi|\hat{H}|\psi\rangle$,
where $|\psi\rangle$ retains continuous variables. A gradient descent (GD)
method is then employed to solve the continuous-variable optimization problem.
The computational cost for determining $E^{\text{Lowest}}$ and its
$|\psi\rangle$ requires the derivation of $N\times L$, where $N$ is the number
of iterations required to reach the converged $E^{\text{Lowest}}$. $N$
generally increases with $L$. The excited eigenvalue and eigenstate are also
computed using the same optimization method with the constraint of the
orthonormal condition between the eigenstates. The overall computational cost
for determining the full $E$ and $|\psi\rangle$ in all the energy states is
approximately $N^{2}\times L^{2}$. Therefore, it is difficult to calculate the
full $E$ and $|\psi\rangle$ of $\hat{H}$ with massive $L$ in a feasible time
using both ED and GD approaches on classic computer.

The variational quantum eigensolver (VQE) algorithm on gate-type quantum
computers with qubit circuits is another promising tool for determining
$E^{\text{Lowest}}$ and its $|\psi\rangle$ of $\hat{H}$ Preskill2018 ;
Higgott2019 ; Bharti2022 . These are measured by the operation of qubits on a
circuit in combination with derivations on a classical computer. Many degrees
of freedom for expressing continuous variables $|\psi\rangle$ require a deep
circuit, which is accompanied with large quantum errors. Therefore, another
qubit is required to minimize the larger quantum errors. Present gate-type
quantum computers are only composed of approximately hundreds of qubits, which
is insufficient to reduce quantum errors for large $\hat{H}$, even though they
show limitless possibilities with increasing qubits in the future Preskill2018
.

The quantum annealing (QA) approach on adiabatic quantum computers is another
alternative for fast computation of optimization problems Kadowaki1998 ;
Johnson2011 ; King2022 . The QA approach measures the possible
$|\psi^{\text{Binary}}\rangle$ with only the optimized binary configurations
through the quantum adiabatic process in the transverse-field Ising model,
where $|\psi^{\text{Binary}}\rangle$ is the optimized binary spin
configuration. This has already been performed on D-Wave quantum annealer
(D-Wave QA) with 5000 qubits, which can compete with classic computers in
terms of computational speed and accuracy for combinatorial optimization
problems of moderate size Amin2018 ; Inoue2021 . It can be also employed as
quantum simulators for exploring phase transition and dynamic behaviors
King2021(1) ; Kairys2020 ; King2021 ; Park2022 ; Irie2021 ; Albash2018 ;
Kumar2024 ; Teplukhin2020 . On the other hand, the method only works well on
binary combinatorial optimization problems and is limited to continuous-
variable optimization problems. Therefore, an extension of the approach to
continuous-variable optimization problems that can be applied to more
important applications is required in D-Wave QA.

In this Letter, we describe the development of a full quantum eigensolver (QE)
that is based on the extension of optimized binary configurations determined
using QA on the D-Wave QA. This method determines the possible $E^{\text{QE}}$
and $|\psi^{\text{QE}}\rangle$ values with continuous variables in all energy
states through iterative QA measurements. The computational cost is $\eta ML$
for full $E^{\text{QE}}$ and $|\psi^{\text{QE}}\rangle$ of $\hat{H}$ of size
$L\times L$, where $M$ and $\eta$ are the number of QA measurements required
to reach the converged $|\psi\rangle$ and the total annealing time of many QA
shots, respectively. This means that, because $\eta$ is very short, $L$ and
the other parameters in $\hat{H}$ do not significantly affect the
computational time. We selected a one-dimensional (1D) ionic non-interacting
tight-binding $\hat{H}$ whose exact $E$ is known in the entire energy
spectrum. We measured the possible $E^{\text{QE}}$ and
$|\psi^{\text{QE}}\rangle$ in various cases of both metallic and insulating
phases using our QE algorithm for D-Wave QA. We compared them with the exact
$E$ and confirmed that it provides an exact $E$ within the errors of $5\times
10^{-3}$. We believe that the proposed QE approach exhibits computational
superiority over ED with $L^{3}$ iterations and GD with $N^{2}L^{2}$
derivations on a classical computer.

Figure 1:  (Color online) (a) Zephyr graph with $20$ couplers between qubits
on D-Wave QA Advantage2 prototype2.3. (b) Number of qubits as a function of
$L$ for embedding Hamiltonian $\hat{H}$ required for quantum annealing
measurement of the fully connected systems. The original $\hat{H}$ means the
topology without the chains.

_Algorithm of Quantum Eigensolver (QE)_ —$E^{\text{Lowest}}$ and its
$|\psi\rangle$ of $\hat{H}$ are determined by the optimization equation given
as

$E^{\text{Lowest}}=\text{min}_{|\psi\rangle}\langle\psi|\hat{H}|\psi\rangle,$
(1)

where $|\psi\rangle$ retains the continuous variable. The idea of our
algorithm is to separate $|\psi\rangle$ into
$|\psi\rangle=|\psi^{\prime}\rangle+|\phi\rangle$. Both
$|\psi^{\prime}\rangle$ and $|\phi\rangle$ are determined by iterative QA
measurements. The detailed computational procedure is followed as: (i) We
first determine $E^{\prime}$ and $|\psi^{\prime}\rangle$ in
$E^{\prime}=\text{min}_{|\psi^{\prime}\rangle}\langle\psi^{\prime}|\hat{H}|\psi^{\prime}\rangle$
by QA measurement. The initial $|\psi^{\prime}\rangle$ is only expressed as
possible $|\psi^{\text{Binary}}\rangle$ that shows very rough solution of the
optimized binary configuration. (ii) We modify $\hat{H}$ into
$\hat{H}-E^{\prime}I$, where $I$ is the identity matrix. Now we rewrite
$|\phi\rangle$ into $|\phi\rangle=|\psi\rangle-|\psi^{\prime}\rangle$, and
optimize $|\psi\rangle$ in $\hat{H}-E^{\prime}\hat{I}$. The optimization
equation for $|\psi\rangle$ is given as

$E_{|\psi\rangle}=\text{min}_{|\psi\rangle}[(\langle\psi|-\langle\psi^{\prime}|)\hat{H}^{(ii)}(|\psi\rangle-|\psi^{\prime}\rangle)],$
(2)

where $\hat{H}^{(ii)}$ is $\hat{H}^{(ii)}=\hat{H}-E^{\prime}\hat{I}$.
$|\psi^{\prime}\rangle$ were already known from QA measurement in (i)
procedure. We measure $|\psi\rangle$ in Eq. (2) through QA approach again.
Note that $\langle\psi^{\prime}|\hat{H}|\psi\rangle$ and
$\langle\psi^{\prime}|\hat{H}|\psi\rangle^{\text{T}}$ are only diagonal parts
in the matrix of Eq. (2). (iii) Next, $|\psi^{\prime}\rangle$ is replaced into
$|\psi\rangle$ measured in (ii) procedure, and (i) procedure is runing again.
$|\psi^{\prime}\rangle$ is no longer $|\psi^{\text{Binary}}\rangle$. (iv) This
process is repeated until $|\psi^{\prime}\rangle$ converges into
$|\psi\rangle$ without assistance of derivation on classic computer.

Figure 2:  (Color online) (a) The lowest eigenvalues $E^{\text{Lowest}}/t_{1}$
and (b) $|\cos\theta|$ as a function of $t_{2}/t_{1}$. The results are
determined by the optimized binary configurations (Binary) and quantum
eigensolver (QE) approaches. $|\cos\theta|$ is computed by
$|\cos\theta|=\langle\psi^{\text{Exact}}|\psi^{\text{Compare}}\rangle$, where
$|\psi^{\text{Compare}}\rangle$ is $|\psi^{\text{Exact}}\rangle$,
$|\psi^{\text{Binary}}\rangle$ or $|\psi^{\text{QE}}\rangle$. Here, all
eigenstates are normalized.

The excited eigenvalue and its $|\psi\rangle$ are determined using the
modified Hamiltonian $H^{\text{Excited}}$, which is given by
$\hat{H}^{\text{Excited}}=\hat{H}-w|\Psi\rangle\langle\Psi|$, where
$|\Psi\rangle$ is the eigenstate of a one-level-lower eigenvalue Negre2022 .
We impose an orthonormal condition on $|\psi\rangle$ and $|\Psi\rangle$. $w$
is a constant estimated as $w=h-l$, where $h$ and $l$ are the highest and
lowest eigenvalues of $\hat{H}$, respectively. Finally, the excited $E$ and
its $|\psi\rangle$ in $\hat{H}^{\text{Excited}}$ are computed using the same
optimization method as described in Eq. (1). The computational expense is
$\eta ML$ for the full $E$ and $|\psi\rangle$ of $\hat{H}$. Notably, $L$ and
other parameters in $\hat{H}$ do not significantly affect the computational
time, contrary to the ED algorithm with $L^{3}$ iterations and GD optimization
with $N^{2}L^{2}$ algorithm on a classical computer, because $\eta$ is a very
short time within $10^{-2}$ seconds for the D-Wave QA.

_Hamiltonian for inspection of QE_ —We selected the one-dimensional ionic
$t_{1}-t_{2}$ tight-binding $\hat{H}$ that contains the exact $E$ in all
energy states to demonstrate the accuracy and efficiency of our algorithm on
the D-Wave QA. The Hamiltonian $\hat{H}$ is given by

$\displaystyle\hat{H}=$ $\displaystyle-
t_{1}\sum_{<i,j>}(c_{i}^{\dagger}c_{j}+\text{H.C})-t_{2}\sum_{<i,j^{\prime}>}(c_{i}^{\dagger}c_{j^{\prime}}+\text{H.C})$
$\displaystyle+\sum_{i=0}^{L}\left[\Delta(-1)^{(i)}-\mu\right]n_{i},$ (3)

where $c_{i}^{\dagger}$ and $c_{i}$ are the electron creation and annihilation
operators at site $i$, respectively. H.C is a Hermitian conjugate. $t_{1}$ and
$t_{2}$ are the nearest neighbors at site $j$ and next nearest-neighbor
hopping at site $j^{\prime}$, respectively. $L$ is the number of lattices
under periodic boundary condition. $\mu$ and $\Delta$ are the chemical and
ionic potentials, respectively.

Figure 3:  (Color online) Deviations $D$ as a function of number of iteration
for convergence in QE algorithm with different size $L$. $D$ are determined by
$D=|E^{\text{Exact}}-E^{\text{QE}}|$ at the lowest states with
$t_{2}/t_{1}=1.0$. Inset shows the subtle $D$ after convergence of QA
measurements using the QE approach. We confirmed that $D$ were below $5\times
10^{-3}$ in all cases.

_Quantum Annealing Measurements_ —We used D-Wave QA Advantage2 prototype2.3,
where the qubits are designed in a Zephyr structure, for the QA measurements,
as shown in Fig. 1 (a). Under Zephyr topology, each qubit is connected to $20$
different qubits via a coupler. The architecture of the original $\hat{H}$ (or
$\hat{H}^{\text{Excited}}$) topologically matches that of a Zephyr graph by
embeddings with chains of ferromagnetic (FM) order between qubits. Here, the
embedding Hamiltonian $\hat{H}^{\text{Embedding}}$ is given by
$\hat{H}^{\text{Embedding}}=\hat{H}+\hat{H}^{\text{Chain}}$, where
$\hat{H}^{\text{Chain}}$ is the chain part required for embedding. More qubits
are required than in the original ones to form chains. In addition, the
elements in the excited $\hat{H}^{\text{Excited}}$ are completely occupied. We
used the ’dwave.embedding.zephyr.find clique embedding’ library provided by
D-Wave Ocean Package for embedding of the fully connected systems. Fig. 1 (b)
shows the number of qubits as a function of $L$ for the original $\hat{H}$ and
$\hat{H}^{\text{Embedding}}$. The optimal chain coupling for an appropriate FM
order in the chains was determined using the method proposed by our group. The
annealing time was set to $10^{-4}$ seconds per a QA shot. The total annealing
time $\eta$ with $100$ QA shots was set to $10^{-2}$ seconds.

Figure 4:  (Color online) Energy dispersion $\epsilon(\bf{k})/t_{1}$ as a
function of momentum $k/a$ for (a) metal with $\Delta/t_{1}=0.0$ and (b)
insulator with $t_{2}/t_{1}=0.0$. Other parameters are $t_{2}/t_{1}=1.0$ and
$\Delta/t_{1}=0.6$ for metal and insulator, respectively. The lattice constant
$a$ is set to $1$.

_The lowest eigenvalues and eigenstates on metallic phase_ —We first measured
$E^{\text{Lowest}}$ and its $|\psi\rangle$ of $\hat{H}$ with $\Delta=0.0$ as a
function of $t_{2}/t_{1}$ using a QA with optimized binary configurations and
QE approaches on the D-Wave QA. $E^{\text{Binary}}$ is determined using

$E^{\text{Binary}}=\text{min}_{|\psi^{\text{Binary}}\rangle}\langle\psi^{\text{Binary}}|\hat{H}|\psi^{\text{Binary}}\rangle.$
(4)

Fig. 2 (a) shows $E^{\text{Lowest}}/t_{1}$ as a function of $t_{2}/t_{1}$.
From $t_{2}/t_{1}=0.0$ to $t_{2}/t_{1}=0.3$, $E^{\text{Binary}}$ and
$E^{\text{QE}}$ are equal to the exact $E^{\text{Exact}}$. After passing
through this region, $E^{\text{QE}}$ remained equal to $E^{\text{Exact}}$
within errors of $5\times 10^{-3}$, whereas $E^{\text{Binary}}$ was much
higher than $E^{\text{Exact}}$. We also computed $|\cos\theta|$, which is
given by the dot product of $|\psi^{\text{Binary}}\rangle$,
$|\psi^{\text{QE}}\rangle$ and $|\psi^{\text{Exact}}\rangle$ to confirm the
similarity of the eigenstates. Here, $|\psi^{\text{Exact}}\rangle$ was
computed by the NumPy library. The results are present in Fig. 2 (b). We
confirmed that $|\cos\theta|$ between exact and QE eigenstates is nearly $1$.
It means that $|\psi^{\text{QE}}\rangle$ and $|\psi^{\text{Exact}}\rangle$
agree well, as indicated by $E^{\text{Lowest}}/t_{1}$ in Fig. 2 (a).

In Fig. 3 we present the deviations $D$ as a function of the number of
iterations in the QE algorithm for systems with $L=30$, $42$ and $60$ at
$t_{2}/t_{1}=1.0$. Here, the deviations $D$ are computed by
$D=|E^{\text{Exact}}-E^{\text{QE}}|$. Regardless of $L$, the convergence
condition under which $D$ disappears was satisfied after $50$ iterations. It
means that only $50\times 10^{-2}$ and $L\times 50\times 10^{-2}$ seconds were
enough to reach converged results in the QE algorithm of Eq. (1) and to
measure full $E^{\text{QE}}$ in all energy states, respectively. The inset of
Fig. 3 shows the subtle $D$ observed in the QA measurements using the QE
approach, after convergence. We confirmed that $D$ were below $5\times
10^{-3}$ in all cases.

_Energy dispersion including information of the full eigenvalues of all
spectrum_ —Finally, we computed the energy dispersion $\epsilon({\bf k})$ in
two cases with metallic and insulating phases at $\mu=0.0$, because
$\epsilon({\bf k})$ includes $E$ in all energy states in the momentum space.
Here, $\epsilon({\bf k})$ are expressed as $\epsilon({\bf k})=-2t_{1}\cos{\bf
k}-2t_{2}\cos 2{\bf k}$ and $\epsilon({\bf k})=\pm\sqrt{(-2t_{1}\cos{\bf
k})^{2}+\Delta^{2}}$ for metal with $\Delta/t_{1}=0$ and insulater with
$t_{2}/t_{1}=0$. We also employed $t_{2}/t_{1}=1.0$ and $\Delta/t_{1}=1.0$ for
metallic and insulating phases, respectively. Here, $\bf k$ means the momentum
space. The lattice constant $a$ is set to $1$. The results of $\epsilon(\bf
k)$ that contain the full $E$ values for the entire spectra for the metallic
and insulating phases of different $L$ values are shown in Figs. 4(a) and (b),
respectively. We confirmed that the $E$ values determined using the QE method
exactly matche exact $\epsilon({\bf k})$ in both the metallic and insulating
phases.

_Conclusion_ —The QA approach only provides $|\psi^{\text{Binary}}\rangle$ and
its corresponding $E^{\text{Binary}}$. Thus, an extension of the optimized
binary configuration is required to capture $|\psi\rangle$ with continuous
variables. We developed a full QE based on iterative QA measurements on the
D-Wave QA. The approach adjusted $|\psi\rangle$ correctly from the initial
$|\psi^{\text{Binary}}\rangle$ without requiring derivation using a classical
computer. The computational cost is $\eta ML$ for full $E$ and $|\psi\rangle$
for $\hat{H}$ with $L$. $\eta$ was set to be below the maximum value of
$10^{-2}$ seconds. Thus, the computational time was not significantly affected
by $L$ and $M$. This differs to the ED algorithm with $L^{3}$ iterations and
the GD approach with derivations of $N^{2}L^{2}$ on a classical computer. We
considered two cases with metallic and insulating phases in the one-
dimensional non-interacting ionic tight-binding $\hat{H}$ that contains the
exact $E$ in the entire spectra to confirm the efficiency and accuracy of the
QE algorithm on D-Wave QA. We determined that the iterations of the QA
measurements using the QE method converged well. In addition, we confirmed
that the QE method provided an exact $E$ within an error of $5\times 10^{-3}$.
Finally, we believe that the proposed QE will not only demonstrate
computational supremacy over various numerical approaches on classic
computers, but also be widely used for various applications, such as novel
material and drug design.

_Acknowledgements_ —This work was supported by Institute of Information and
communications Technology Planning Evaluation (IITP) grant funded by the
Korean government (MSIT) (No. RS-2023-0022952422282052750001) and by the
quantum computing technology development program of the National Research
Foundation of Korea (NRF) funded by the Korean government (Ministry of Science
and ICT(MSIT)) (No. 2020M3H3A1110365).

## References

  * (1) A. Weiße and H. Fehske, Exact Diagonalization Techniques, Computational Many- Particle Physics Lecture Notes in Physics, 739, (2008).
  * (2) J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * (3) Oscar Higgott, Daochen Wang, and Stephen Brierley, Variational Quantum Computation of Excited States, Quantum 3, 156 (2019).
  * (4) Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik, Noisy intermediate-scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022).
  * (5) Tadashi Kadowaki and Hidetoshi Nishimori, Quantum annealing in the transverse Ising model, Phys. Rev. E 58, 5355 (1998).
  * (6) M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, et al., Quantum annealing with manufactured spins, Nature 473, 194 (2011).
  * (7) A.D. King, S. Suzuki, J. Raymond et al., Coherent quantum annealing in a programmable 2000 qubit Ising chain, Nature Phys. 102, 1745 (2022).
  * (8) M.H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and Roger Melko, Quantum Boltzmann Machine, Phys. Rev. X 8, 021050 (2018).
  * (9) Inoue, D., Okada, A., Matsumori, T. et al., Traffic signal optimization on a square lattice with quantum annealing, Sci Rep 11, 3303 (2021).
  * (10) A. D. King, C. D. Batista, J. Raymond, T. Lanting, I. Ozfidan, G. Poulin-Lamarre, H. Zhang, and M. H. Amin, Quantum Annealing Simulation of Out-of-Equilibrium Magnetization in a Spin- Chain Compound, PRX Quantum 2, 030317 (2021).
  * (11) Paul Kairys, Andrew D. King, Isil Ozfidan, Kelly Boothby, Jack Raymond, Arnab Banerjee, and Travis S. Humble, Simulating the Shastry-Sutherland Ising Model Using Quantum Annealing, PRX Quantum 1, 020320 (2020).
  * (12) King, A.D., Raymond, J., Lanting, T. et al., Scaling advantage over path-integral Monte Carlo in quantum simulation of geometrically frustrated magnets, Nat. Commun. 12, 1113 (2021).
  * (13) H. Park, and H. Lee, Frustrated Ising Model on D-wave Quantum Annealing Machine J. Phys. Soc. Jpn. 91, 074001 (2022).
  * (14) Irie, H., Liang, H., Doi, T. et al. Hybrid quantum annealing via molecular dynamics. Sci Rep 11, 8426 (2021).
  * (15) Tameem Albash and Daniel A. Lidar, Demonstration of a Scaling Advantage for a Quantum Annealer over Simulated Annealing, Phys. Rev. X 8, 031016 (2018).
  * (16) V. Kumar, N. Baskaran, V. S. Prasannaa, K. Sugisaki, D. Mukherjee, K. G. Dyall, and B. P. Das, Computation of relativistic and many-body effects in atomic systems using quantum annealing, Phys. Rev. A 109, 042808 (2024).
  * (17) A. Teplukhin, B. K. Kendrick, S. Tretiak, and P. A. Dub, Electronic structure with direct diagonalization on a D-wave quantum annealer, Scientific Reports 10, 20753 (2020).
  * (18) C. F. A. Negre, A. Lopez-Bezanilla, Y. Zhang, P. D. Akrobotu, S. M. Mniszewski, S. Tretiak, and P. A. Dub, Toward a QUBO-Based Density Matrix Electronic Structure Method, Journal of Chemical Theory and Computation 18 (7), 4177 (2022).