*]Facility for Rare Isotope Beams and Department of Physics and Astronomy, Michigan State University, East Lansing, 48824 MI, USA

# Quantum techniques for eigenvalue problems

Dean Lee<EMAIL_ADDRESS>[

###### Abstract

This article is a brief introduction to quantum algorithms for the eigenvalue
problem in quantum many-body systems. Rather than a broad survey of topics, we
focus on providing a conceptual understanding of several quantum algorithms
that cover the essentials of adiabatic evolution, variational methods, phase
detection algorithms, and several other approaches. For each method, we
discuss the potential advantages and remaining challenges.

###### keywords:

adiabatic evolution, quantum approximate optimization algorithm, unitary
coupled cluster, Thouless theorem, Jordan-Wigner transformation, Trotter-
Suzuki expansion, variational methods, phase estimation, iterative phase
estimation, rodeo algorithm

##  1 Introduction

Quantum computing has the potential to address many of the unsolved problems
of quantum many-body physics. By allowing for arbitrary linear combinations of
tensor products of qubits, one can store exponentially more information than
classical bits. This opens the possibility of calculations of strongly-
interacting systems with many degrees of freedom without the need for Monte
Carlo methods and their accompanying problems associated with sign
oscillations [1, 2, 3, 4, 5, 6]. Furthermore, qubits naturally evolve with
unitary real-time dynamics, providing access to non-equilibrium processes,
which are often well beyond the reach of first-principles calculations using
classical computers. But there are also great challenges to realizing the
promise of quantum computing. One of the main problems is the fact that the
quantum computing devices available today have significant limitations due to
gate errors, qubit decoherence, faulty measurement readout, small numbers of
qubits, and limited qubit connectivity. These problems severely limit the
class of problems that one can address at present. Nevertheless, significant
advances are being made in quantum hardware performance and scale [7, 8, 9],
and it is useful to consider the design and performance of quantum algorithms
as quantum resources grow and become more reliable.

There is an excellent and comprehensive review on quantum computing and
quantum many-body systems in Ref. [10]. Instead of writing another review with
similarly broad scope, in this article we instead focus on several algorithms
of interest for eigenvalue problems. The aim is to provide a readable
introduction for novice readers with enough detail to demonstrate the concepts
and execution of each method. We should note that there are many useful
algorithms of relevance to eigenvalue problems that we do not cover here.
These include cooling algorithms [11, 12, 13, 14, 15], coupled heat bath
approaches [16, 17], dissipative open system methods [18, 19, 20], spectral
combing [21], symmetry projection techniques [22, 23, 24, 25, 26], linear
combinations of unitaries [27, 28, 29, 30, 31], and imaginary time evolution
[32, 33, 34, 35, 36].

In the following, we start with a review of the adiabatic theorem and the
performance of adiabatic evolution for the preparation of eigenstates. After
this, we cover the broad class of variational methods. We discuss gradient
calculation techniques for optimization and several specific variational
algorithms. Thereafter we present several phase detection algorithms. These
include phase estimation, iterative phase estimation, and the rodeo algorithm.
We then conclude with a summary and outlook for the future.

##  2 Adiabatic Evolution

The adiabatic theorem states that if a quantum state is an eigenstate of an
initial Hamiltonian $H(0)=H_{0}$, then the quantum state will remain trapped
in an exact eigenstate of the instantaneous Hamiltonian $H(t)$ in the limit
that the time dependence of $H(t)$ is infinitely slow [37, 38]. If this
evolution has only finite duration, then the error will scale inversely with
the total time evolution, $T$. We can use quantum adiabatic evolution to
prepare the eigenstates of any Hamiltonian $H_{1}$ by preparing an exact
eigenstate of some simple initial Hamiltonian $H_{0}$. For the purpose of
analysis, it is convenient to scale out the dependence on the total duration
of time $T$ and work with the rescaled variable $s=t/T$. We then make a smooth
interpolation $H(s)$ with $s$ ranging from $s=0$ to $s=1$, with $H(0)=H_{0}$
and $H(1)=H_{1}$ [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]. Let us
define the adiabatic evolution operator

$U(s)={\cal T}\exp\left[-iT\int_{0}^{s}H(s^{\prime})ds^{\prime}\right],$ (1)

where ${\cal T}$ indicates time ordering where operators at later times are
placed on the left. In the limit of large time $T$, the unitary transformation
$U(1)$ will map any eigenstate of $H_{0}$ to an eigenstate of $H_{1}$. In Ref.
[51], it is observed that the unitarily-transformed Hamiltonian,

$H^{\prime}(1)=U^{\dagger}(1)H_{1}U(1),$ (2)

is a Hamiltonian whose eigenvalues are equal to $H_{1}$ but whose eigenvectors
are equal to $H_{0}$. For this reason, the term “Hamiltonian translator” was
used to describe the unitary transformation $U(1)$. Suppose we start from the
Hamiltonian $H(0)$ and perform a perturbation theory expansion in the
difference, $H^{\prime}(1)-H(0)$,

$H^{\prime}(1)=H(0)+[H^{\prime}(1)-H(0)].$ (3)

Since $H(0)$ and $H^{\prime}(1)$ share the same eigenvectors, we find that
first-order perturbation for the energy is exact and all other terms in
perturbation theory for the energy or wave function vanish.

Let us now consider the one-parameter eigenvector $\ket{\psi(s)}$, which is an
instantaneous eigenvector of $H(s)$ for $s$ in the interval $[0,1]$. Let
$\Delta(s)$ be the spectral gap between $\ket{\psi(s)}$ and the rest of the
energy spectrum of $H(s)$. In computing the spectral gap, we can ignore
sectors that are orthogonal to $\ket{\psi(s)}$ due to symmetries that are
respected by $H(s)$. We note that $\ket{\psi(0)}$ is an eigenstate of $H_{0}$.
Let us define $\ket{\psi_{U}(s)}$ as $U(s)\ket{\psi(0)}$.

We use the symbol $||\cdot||$ to denote the operator norm. Building upon the
work of Ref. [52, 53, 54], Jansen et al. [55] derived the rigorous bound that

$T\geq\frac{1}{\delta}\left\\{\int^{s}_{0}\left[\frac{||\partial_{s}^{2}H(s^{\prime})||}{\Delta^{2}(s)}+7\frac{||\partial_{s}H(s^{\prime})||^{2}}{\Delta^{3}(s)}\right]ds^{\prime}+B\right\\}$
(4)

is sufficient to satisfy the error bound

$|\braket{\psi(s)}{\psi_{U}(s)}|\geq 1-\delta,$ (5)

where $B$ is a boundary term that vanishes when $\partial_{s}H(0)$ and
$\partial_{s}H(1)$ both equal zero [50]. We see from Eq. (4) that, for any
fixed system, the required time $T$ is scaling inversely with the error
$\delta$.

The challenge with adiabatic state preparation for eigenstates of quantum
many-body systems is the fact that $\Delta(s)$ may be extremely small for
large systems. This is especially true when $H_{0}$ and $H_{1}$ have very
different eigenstates, and $H(s)$ must pass through one or more quantum phase
transitions. This motivates the search for initial Hamiltonians $H_{0}$ for
which the starting eigenstate can be prepared on a quantum computer, but the
eigenstate structure of $H_{0}$ is not completely trivial and has some
resemblance to that of $H_{1}$ [56]. Even in cases where $\ket{\psi_{U}(1)}$
is not a good approximation to the eigenstate $\ket{\psi(1)}$ of $H_{1}$, the
state $\ket{\psi_{U}(1)}$ can still be a useful starting vector for other
state preparation algorithms which converge more rapidly.

In order to perform the time evolution in Eq. (1), one usually uses some
version of the Trotter approximation. The conceptual starting point for the
Trotter approximation is the Baker-Campbell-Hausdorff formula, which states
that when $e^{A}e^{B}=e^{C}$, we have the formal series

$C=A+B+\frac{1}{2}[A,B]+\frac{1}{12}[A,[A,B]]-\frac{1}{12}[B,[A,B]]+\cdots.$
(6)

If our Hamiltonian has two non-commuting pieces,

$H=H_{A}+H_{B}.$ (7)

At first order in the Trotter-Suzuki expansion, we can use [57]

$\displaystyle e^{-iH\Delta t}$ $\displaystyle=e^{-iH_{A}\Delta
t}e^{-iH_{B}\Delta t}+O[(\Delta t)^{2}]$ $\displaystyle=e^{-iH_{B}\Delta
t}e^{-iH_{A}\Delta t}+O[(\Delta t)^{2}].$ (8)

At second order we have

$\displaystyle e^{-iH\Delta t}$ $\displaystyle=e^{-iH_{B}\Delta
t/2}e^{-iH_{A}\Delta t}e^{-iH_{B}\Delta t/2}+O[(\Delta t)^{3}]$
$\displaystyle=e^{-iH_{A}\Delta t/2}e^{-iH_{B}\Delta t}e^{-iH_{A}\Delta
t/2}+O[(\Delta t)^{3}].$ (9)

The generalization to higher-order expressions can be found in Ref. [58]. The
performance of the Trotter-Suzuki expansion can be improved in numerous ways,
such as using random orderings [59, 60], sums of Trotter products at different
orders [61, 62, 63], extrapolation methods [64], and renormalization [65].

##  3 Variational Methods

Variational quantum algorithms encompass a broad class of methods that are
among the most popular approaches to the preparation of eigenstates using
current and near-term quantum hardware. While the examples we consider here
are optimizing a single vector, there are also many different variational
methods that use subspaces [66, 67, 68, 69, 70]. The typical strategy is a
hybrid approach where the quantum device is used to prepare a parameterized
family of possible wave functions, and then a classical computation is
performed to minimize the associated cost function. Let $\boldsymbol{\theta}$
be an $L$-dimensional vector of parameters $\theta_{j}$. The most common
example is the search for the ground state of a quantum Hamiltonian $H$ by
minimizing a cost function $C({\boldsymbol{\theta}})$ given by the energy
expectation value
$C({\boldsymbol{\theta}})=\braket{{\boldsymbol{\theta}}}{H}{{\boldsymbol{\theta}}}$
[71, 72]. We consider a general ansatz for the wave function
$\ket{\boldsymbol{\theta}}$ that is a product of unitary operators acting upon
some simple initial state $\ket{\psi_{I}}$ [73, 74],

$\ket{{\boldsymbol{\theta}}}=V_{L}U_{L}(\theta_{L})\cdots
V_{1}U_{1}(\theta_{1})\ket{\psi_{I}}.$ (10)

where each $V_{j}$ is a fixed unitary operator. It is convenient to take each
$U_{j}(\theta_{j})$ as an exponential of a Hermitian operator $H_{j}$,

$U_{j}(\theta_{j})=\exp(-iH_{j}\theta_{j}/2),$ (11)

where we restrict $H_{j}$ to be its own inverse so that $H_{j}^{2}=I$. This
involutory condition is satisfied by any product of Pauli matrices on any
multi-qubit system. In such cases, we have the simple trigonometric relation,

$U_{j}(\theta_{j})=\cos(\theta_{j}/2)I-i\sin(\theta_{j}/2)H_{j}.$ (12)

For any operator $O$, we find that

$U_{j}^{\dagger}(\theta_{j})OU_{j}(\theta_{j})=O_{1}+O_{\sin}\sin(\theta_{j})+O_{\cos}\cos(\theta_{j}),$
(13)

for some operators $O_{1}$, $O_{\sin}$, and $O_{\cos}$ independent of
$\theta_{j}$. It follows that [74]

$\displaystyle\frac{\partial}{\partial\theta_{j}}U_{j}^{\dagger}(\theta_{j})OU_{j}(\theta_{j})$
$\displaystyle=O_{\sin}\cos(\theta_{j})-O_{\cos}\sin(\theta_{j})$
$\displaystyle=\frac{U_{j}^{\dagger}(\theta_{j}+\alpha_{j})OU_{j}(\theta_{j}+\alpha_{j})-U_{j}^{\dagger}(\theta_{j}-\alpha_{j})OU_{j}(\theta_{j}-\alpha_{j})}{2\sin(\alpha_{j})},$
(14)

for any $\alpha_{j}$ such that $\sin(\alpha_{j})\neq 0$. We note that this
parameter shift formula is exact and not simply a finite-difference
approximation. This allows values for $\alpha_{j}$ of $O(1)$, which is helpful
to measure gradient components in the presence of stochastic and systematic
errors. One can now compute the components of the gradient using

$\frac{\partial}{\partial\theta_{j}}C(\boldsymbol{\theta})=\frac{C(\boldsymbol{\theta}+\boldsymbol{\alpha}_{j})-C(\boldsymbol{\theta}-\boldsymbol{\alpha}_{j})}{2\sin(\alpha_{j})},$
(15)

where the vector $\boldsymbol{\alpha}_{j}$ has components
$[\boldsymbol{\alpha}_{j}]_{k}=\alpha_{j}\delta_{jk}$. These gradients can be
used to minimize the cost function $C(\boldsymbol{\theta})$.

Consider adiabatic evolution with initial Hamiltonian $H_{0}$, final
Hamiltonian $H_{1}$, and interpolating Hamiltonian $H(s)=sH_{1}+(1-s)H_{0}$.
We then have a string of exponentials for our adiabatic evolution operator,

$U(1)=e^{-iH(1)ds}\cdots e^{-iH(s)ds}\cdots e^{-iH(0)ds},$ (16)

which we apply to the ground state of $H_{0}$. Let $N$ be the number of time
steps and let $ds=1/(N+1)$. If we now use the Trotter approximation to write

$e^{-iH(s)ds}=e^{-i[sH_{1}+(1-s)H_{0}]ds}\approx
e^{-isH_{1}ds}e^{-i(1-s)H_{0}ds},$ (17)

then the adiabaic evolution operator has the form

$U(1)\approx e^{-i\gamma_{N}H_{1}ds}e^{-i\beta_{N}H_{0}ds}\cdots
e^{-i\gamma_{j}H_{1}ds}e^{-i\beta_{j}H_{0}ds}\cdots
e^{-i\gamma_{0}H_{1}ds}e^{-i\beta_{0}H_{0}ds},$ (18)

for $\beta_{j}=1-jds$ and $\gamma_{j}=jds$. This structure provides the
theoretical motivation for the quantum approximate optimization algorithm
(QAOA) [75]. Instead of using the values for $\gamma_{j}$ and $\beta_{j}$ as
prescribed by adiabatic evolution, they are treated as free variational
parameters optimized to minimize the energy expectation of the Hamiltonian
$H_{1}$.

For large quantum systems, the required number of variational parameters will
grow with system size. The number of variational parameters needed as a
function of the size of the system with fixed error tolerance remains an open
question. There are at least two major challenges that arise in quantum
variational algorithms for large systems. The first challenge is the problem
of barren plateaus. For parameterized random quantum circuits, the components
of the cost function gradient will become exponentially small in the number of
qubits of the quantum system [73].

The second challenge is the appearance of many local minima, making gradient
descent optimization difficult. Before discussing the problem of local minima,
we first review some terminology from computational complexity theory. A
decision problem is one where the two possible answers are yes or no. P refers
to the set of decision problems that can be solved using a deterministic
Turing machine in polynomial time. NP refers to the set of decision problems
whose solution, once given, can be confirmed by a deterministic Turing machine
in polynomial time. Equivalently, NP refers to decision problems that can be
solved using a non-deterministic Turing machine in polynomial time, where a
general non-deterministic Turing machine is endowed with the ability to branch
over all possible outcomes in parallel. A problem $p$ is NP-hard if all
problems in NP can be obtained in polynomial time from the solution of $p$. If
a decision problem in NP is NP-hard, then it is called NP-complete.

Consider a graph with $d$ vertices and an adjacency matrix $A_{i,j}$ marking
the edges of the graph that equal $0$ or $1$ for each pair of vertices
$\\{i,j\\}$. The MaxCut problem poses the task of finding the subset $S$ of
the vertices that maximizes the number of edges connecting $S$ and its
complement,

$\sum_{i\in S}\sum_{j\notin S}A_{i,j}.$ (19)

The MaxCut problem was shown to be NP-complete [76]. The continuous MaxCut
problem consists of finding the $d$-dimensional vector
$\boldsymbol{\phi}=[0,2\pi)^{d}$ that minimizes

$\mu({\boldsymbol{\phi}})=\frac{1}{4}\sum_{i=1}^{d}\sum_{j=1}^{d}A_{i,j}[\cos(\phi_{i})\cos(\phi_{j})-1].$
(20)

In Ref. [77], the continuous MaxCut problem is shown to be equivalent to the
MaxCut problem and therefore also NP-hard. Furthermore, the continuous MaxCut
problem can also be recast as a variational quantum optimization problem for
the Ising model Hamiltonian,

$\frac{1}{4}\sum_{i=1}^{d}\sum_{j=1}^{d}A_{i,j}(\sigma^{z}_{i}\sigma^{z}_{j}-1),$
(21)

with variational wave function

$e^{-i\sigma^{y}_{d}\phi_{d}/2}\cdots
e^{-i\sigma^{y}_{1}\phi_{1}/2}\ket{0}^{\otimes d}.$ (22)

Although there is no proof that NP contains problems outside of P, there is
much speculation that this is true. NP-hard problems would then belong to the
set of difficult problems outside of P, and this would include the problem of
minimizing the variational cost function for an Ising Hamiltonian.

Although the general performance of variational methods for large quantum
systems is challenging, there are many cases in which major simplifications
arise because of some simplification, such as the emergence of a mean-field
picture. There are many examples of such approaches for fermionic quantum
many-body systems [78, 79, 80, 81, 82, 83, 84, 85]. One popular example is the
unitary coupled cluster (UCC) method. In the UCC method, one starts with an
initial state $\ket{\psi_{I}}$, which is a mean-field reference state. For the
unitary transformation, $U$, we take the form

$U=e^{T(\boldsymbol{\theta})-T^{\dagger}(\boldsymbol{\theta})},$ (23)

where

$T(\boldsymbol{\theta})=\sum_{m}T_{m}(\boldsymbol{\theta}),$ (24)

and $T_{m}(\boldsymbol{\theta})$ is an $m$-body operator that produces
excitations. The singles excitation has the form,

$T_{1}(\boldsymbol{\theta})=\sum_{i}\sum_{a}\theta^{i}_{a}a^{\dagger}_{a}a_{i},$
(25)

where $a^{\dagger}_{a}$ and $a_{i}$ and fermionic creation and annihilation
operators for orbitals $a$ and $i$ respectively. The doubles terms has the
structure,

$T_{2}(\boldsymbol{\theta})=\frac{1}{4}\sum_{i<j}\sum_{a<b}\theta^{i,j}_{a,b}a^{\dagger}_{a}a^{\dagger}_{b}a_{j}a_{i}.$
(26)

For the general case, we have

$T_{m}(\boldsymbol{\theta})=\frac{1}{(m!)^{2}}\sum_{i<j<\cdots}\sum_{a<b<\cdots}\theta^{i,j,\cdots}_{a,b,\cdots}a^{\dagger}_{a}a^{\dagger}_{b}\cdots
a_{j}a_{i}.$ (27)

There are several ways to encode ferimonic antisymmetrization properties on a
quantum computer. Although often not the most efficient, the simplest approach
is the Jordan-Wigner transformation [86]. We define

$\displaystyle\sigma^{+}_{j}=$
$\displaystyle(\sigma^{x}_{j}+i\sigma^{y}_{j})/2,$
$\displaystyle\sigma^{-}_{j}=$
$\displaystyle(\sigma^{x}_{j}-i\sigma^{y}_{j})/2,$ (28)

and use the convention that $\ket{0}$ corresponds to occupation number $0$,
and $\ket{1}$ corresponds to occupation number $1$. We then have a faithful
representation of the algebra of creation and annihilation operators with the
mapping

$\displaystyle a^{\dagger}_{j}$
$\displaystyle=\sigma^{-}_{j}\otimes\sigma^{z}_{j-1}\otimes\cdots\otimes\sigma^{z}_{1},$
$\displaystyle a_{j}$
$\displaystyle=\sigma^{+}_{j}\otimes\sigma^{z}_{j-1}\otimes\cdots\otimes\sigma^{z}_{1}.$
(29)

This gives the required anticommutation relations,

$\\{a_{j},a^{\dagger}_{k}\\}=\delta_{j,k},\;\;\\{a_{j},a_{k}\\}=\\{a^{\dagger}_{j},a^{\dagger}_{k}\\}=0.$
(30)

Many other antisymmetrization techniques [87, 88, 89, 90, 91, 92] have been
designed that are computationally more efficient in cases where the products
of creation and annihilation operators in the Hamiltonian appear in
combinations with some locality restriction with respect to the orbital index.

A convenient choice for the mean-field reference state $\ket{\psi_{I}}$ is a
Hartree-Fock state, corresponding to a Slater determinant of single-particle
orbitals achieving the lowest energy expectation value. The Thouless theorem
[93] shows how to prepare any desired Slater determinant state starting from
any other Slater determinant state. Let $\alpha_{p}({\bf r})$ label the
original orbitals and let $\beta_{p}({\bf r})$ label the new orbitals. We take
$a^{\dagger}_{p},a_{p}$ to be the creation and annihilation operators for
$\alpha_{p}({\bf r})$, and $b^{\dagger}_{p},b_{p}$ to be the creation and
annihilation operators for $\beta_{p}({\bf r})$. Let $N$ be the number of
particles in our system of interest. The aim is to derive a simple relation
between $b^{\dagger}_{N}\cdots b^{\dagger}_{1}\ket{\rm vac}$ and
$a^{\dagger}_{N}\cdots a^{\dagger}_{1}\ket{\rm vac}$. Without loss of
generality, we use a linear transformation to redefine the orbitals
$\beta_{1}({\bf r}),\cdots,\beta_{N}({\bf r})$ so that for each
$p=1,\cdots,N$, we have

$b^{\dagger}_{p}=a^{\dagger}_{p}+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,p}$
(31)

for some coefficient matrix $u_{q,p}$. The linear transformation on
$\beta_{1}({\bf r}),\cdots,\beta_{N}({\bf r})$ has no effect on
$b^{\dagger}_{N}\cdots b^{\dagger}_{1}\ket{\rm vac}$ except for introducing an
overall normalization factor. Our convention will ensure that
$b^{\dagger}_{N}\cdots b^{\dagger}_{1}\ket{\rm vac}$ and
$a^{\dagger}_{N}\cdots a^{\dagger}_{1}\ket{\rm vac}$ have the same
normalization.

The Thouless theorem is based on the observation that for each $p=1,\cdots,N$,

$\displaystyle\left(a^{\dagger}_{p}+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,p}\right)F({\rm
no}\;a^{\dagger}_{p})\ket{\rm vac}$
$\displaystyle=\left(1+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,p}a_{p}\right)a^{\dagger}_{p}F({\rm
no}\;a^{\dagger}_{p})\ket{\rm vac},$ (32)

where $F({\rm no}\;a^{\dagger}_{p})$ is an arbitrary function of the creation
and annihilation operators where $a^{\dagger}_{p}$ does not appear. We then
have

$\displaystyle b^{\dagger}_{N}\cdots b^{\dagger}_{1}\ket{\rm vac}$
$\displaystyle=\left(a^{\dagger}_{N}+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,N}\right)\cdots\left(a^{\dagger}_{1}+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,1}\right)\ket{\rm
vac}$
$\displaystyle=\left(1+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,N}a_{N}\right)a^{\dagger}_{N}\cdots\left(1+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,1}a_{1}\right)a^{\dagger}_{1}\ket{\rm
vac}.$ (33)

This leads to simple relation,

$\displaystyle b^{\dagger}_{N}\cdots b^{\dagger}_{1}\ket{\rm vac}$
$\displaystyle=\left(1+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,N}a_{N}\right)\cdots\left(1+\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,1}a_{1}\right)a^{\dagger}_{N}\cdots
a^{\dagger}_{1}\ket{\rm vac}$
$\displaystyle=\exp\left(\sum_{p=1}^{N}\sum_{q=N+1}^{\infty}a^{\dagger}_{q}u_{q,p}a_{p}\right)a^{\dagger}_{N}\cdots
a^{\dagger}_{1}\ket{\rm vac}.$ (34)

Once the Hartree-Fock orbitals are determined using classical computing, one
can prepare a simple $N$-particle Slater determinant state with orbitals given
by the computational basis of the quantum computer and then apply the
transformation in Eq. (34) [94].

##  4 Phase Detection Algorithms

Quantum phase estimation [95] is a well-known example of a phase detection
algorithm that can be used to find energy eigenvalues and prepare energy
eigenstates of the quantum many-body problem [96, 97, 98, 24]. Suppose for the
moment that $\ket{\psi}$ is an eigenstate of the unitary operator $U$ with
eigenvalue $e^{2\pi i\theta}$. Of particular interest is the case where the
unitary operator $U$ is the time evolution operator for some Hamiltonian $H$
over some fixed time step $\Delta t$. The goal is to efficiently determine the
phase angle $\theta$. Since $U\ket{\psi}=e^{2\pi i\theta}\ket{\psi}$, we have
$U^{2^{j}}\ket{\psi}=e^{2\pi i\theta 2^{j}}\ket{\psi}$ for any nonnegative
integer $j$. Together with the state $\ket{\psi}$, we take $n$ ancilla qubits
with each initialized as $\ket{0}$. The resulting state is $\ket{0}^{\otimes
n}\otimes\ket{\psi}$. The Hadamard gate is a single qubit gate that maps
$\ket{0}$ to $1/\sqrt{2}(\ket{0}+\ket{1})$ and maps $\ket{1}$ to
$1/\sqrt{2}(\ket{0}-\ket{1})$. The action of the Hadamard gate for a general
linear combination of $\ket{0}$ and $\ket{1}$ is determined by linearity. We
apply Hadamard gates to each of the ancilla qubits so that we get

$\frac{1}{2^{n/2}}\left(\ket{0}+\ket{1}\right)^{\otimes n}\otimes\ket{\psi}.$
(35)

For each of the ancilla qubits $j=0,\cdots,n-1$, we use the ancilla qubit to
control the unitary gate $U^{2^{j}}$. This means that $U^{2^{j}}$ is applied
when the ancilla qubit $j$ is in state $\ket{1}$, but no operation is
performed if the ancilla qubit is in state $\ket{0}$. The result we get is
[99]

$\frac{1}{2^{n/2}}\left(\ket{0}+e^{2\pi i\theta
2^{n-1}}\ket{1}\right)\otimes\cdots\otimes\left(\ket{0}+e^{2\pi i\theta
2^{0}}\ket{1}\right)\otimes\ket{\psi}=\ket{f(\theta)}\otimes\ket{\psi},$ (36)

where

$\displaystyle\ket{f(\theta)}=$
$\displaystyle\frac{1}{2^{n/2}}\sum_{m=0}^{2^{n}-1}\left(e^{2\pi i\theta
2^{n-1}m_{n-1}}\ket{m_{n-1}}\right)\otimes\cdots\left(e^{2\pi i\theta
2^{0}m_{0}}\otimes\ket{m_{0}}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2^{n/2}}\sum_{m=0}^{2^{n}-1}e^{2\pi i\theta
m}\ket{m_{n-1}}\otimes\cdots\otimes\ket{m_{0}},$ (37)

and $m_{n-1}\cdots m_{0}$ are the binary digits of the integer $m$.

Let $k$ be an integer between $0$ and $2^{n}-1$ with binary representation
$k_{n-1}\cdots k_{0}$. We note that when $\theta$ equals $k$ divided by
$2^{n}$, then $\ket{f(k/2^{m})}$ is the quantum Fourier transform of the state
$\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}$,

$\ket{f(k/2^{m})}=\frac{1}{2^{n/2}}\sum_{m=0}^{2^{n}-1}e^{2\pi
ikm/2^{n}}\ket{m_{n-1}}\otimes\cdots\otimes\ket{m_{0}}.$ (38)

We can therefore extract information about the value of $\theta$ by applying
the inverse quantum Fourier transform to $\ket{f(\theta)}$,

$QFT^{-1}\ket{f(\theta)}=\frac{1}{2^{n}}\sum_{k=0}^{2^{n}-1}\sum_{m=0}^{2^{n}-1}e^{-2\pi
i(k/2^{n}-\theta)m}\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}.$ (39)

We see that if $\theta$ equals $k/2^{n}$ for some integer $k$ in the
summation, then $QFT^{-1}\ket{f(\theta)}$ equals
$\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}$. In the general case, we get a
superposition of such states $\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}$
that is highly peaked for integers $k$, where $k/2^{n}$ is close to $\theta$.
We simply measure each ancilla qubit and determine $k/2^{n}$ to obtain an
estimate of $\theta$. This is repeated over several trials to build a
probability distribution and refine the estimate of $\theta$.

Suppose now that $\ket{\psi}$ is not an eigenstate of $U$ but rather a general
superposition of eigenstates $\ket{\psi_{a}}$ with eigenvalues $e^{2\pi
i\theta_{a}}$,

$\ket{\psi}=\sum_{a}c_{a}\ket{\psi_{a}}.$ (40)

We can now apply phase estimation to the general state $\ket{\psi}$ in exactly
the same manner as before. Let us assume that the separation between each
$\theta_{a}$ is large compared to $1/2^{n}$. This ensures that the peaked
distributions we get for each eigenvector have negligible overlap. The outcome
after measuring the $n$ ancilla qubits will be

$\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}\otimes\ket{\psi_{a}},$ (41)

for some eigenstate $\ket{\psi_{a}}$. The probability of $\ket{\psi_{a}}$
being selected will equal $|c_{a}|^{2}$. The error of quantum phase estimation
in determining eigenvalues will scale inversely with $2^{n}$. This arises from
the discretization of energy values $k/2^{n}$, where $k$ is an integer from
$0$ to $2^{n}-1$. If we relate $U$ to the time evolution of a Hamiltonian $H$
for time step $\Delta t$, the error in the energy scales inversely with the
total time evolution required. This scaling of the uncertainty matches the
lower bound one expects from the Heisenberg uncertainty principle. The error
of phase estimation for eigenstate preparation arises from the admixture of
terms from different eigenstates,

$\frac{1}{2^{n}}\sum_{a}\sum_{m=0}^{2^{n}-1}c_{a}e^{-2\pi
i(k/2^{n}-\theta_{a})m}\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}\otimes\ket{\psi_{a}}.$
(42)

When the spacing between $\theta_{a}$ is much larger than $1/2^{n}$, then the
contamination of other eigenstates will be $O(2^{-n})$. For the case when $U$
is the time evolution of a Hamiltonian $H$ for time duration $\Delta t$, then
the error of eigenstate preparation scales inversely with the total amount of
time evolution needed.

We have mentioned the quantum Fourier transform, but have not yet discussed
how it is implemented. It suffices to describe its action on the state
$\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}$. We again use the notation that
$k_{n-1}\cdots k_{0}$ are the binary digits of the integer $k$. The desired
action of the quantum Fourier transform upon
$\ket{k_{n-1}}\otimes\cdots\otimes\ket{k_{0}}$ is

$\displaystyle\frac{1}{2^{n/2}}\left(\ket{0}+e^{2\pi
ik2^{n-1}/2^{n}}\ket{1}\right)$
$\displaystyle\otimes\cdots\otimes\left(\ket{0}+e^{2\pi
ik2^{0}/2^{n}}\ket{1}\right)$
$\displaystyle=\frac{1}{2^{n/2}}\sum_{m=0}^{2^{n}-1}e^{2\pi
ikm/2^{n}}\ket{m_{n-1}}\otimes\cdots\otimes\ket{m_{0}}.$ (43)

The first few steps of the quantum Fourier transform algorithm will actually
produce the desired result with the tensor product in the reverse order,

$\frac{1}{2^{n/2}}\left(\ket{0}+e^{2\pi
ik2^{0}/2^{n}}\ket{1}\right)\otimes\cdots\otimes\left(\ket{0}+e^{2\pi
ik2^{n-1}/2^{n}}\ket{1}\right).$ (44)

But this can be fixed by pairwise swap gates between qubits $0$ and $n-1$, $1$
and $n-2$, etc.

The quantum Fourier transform begins with the state

$\ket{k_{n-1}}\otimes\ket{k_{n-2}}\otimes\cdots\otimes\ket{k_{0}}.$ (45)

We first act upon qubit $n-1$ with a Hadamard gate and this gives

$\frac{1}{2^{1/2}}\left(\ket{0}+e^{2\pi
ik_{n-1}2^{n-1}/2^{n}}\ket{1}\right)\otimes\ket{k_{n-2}}\otimes\cdots\otimes\ket{k_{0}}.$
(46)

The coefficient in front of $\ket{1}$ equals $1$ if $k_{n-1}=0$ and equals
$-1$ if $k_{n-1}=1$. We use qubit $n-2$ to apply a controlled phase rotation
to qubit $n-1$ by a phase $e^{2\pi ik_{n-2}2^{n-2}/2^{n}}$. The result is

$\frac{1}{2^{1/2}}\left(\ket{0}+e^{2\pi
i(k_{n-1}2^{n-1}+k_{n-2}2^{n-2})/2^{n}}\ket{1}\right)\otimes\ket{k_{n-2}}\otimes\cdots\otimes\ket{k_{0}}.$
(47)

We continue in this manner with qubit $j$ applying a controlled phase rotation
on qubit $n-1$ by a phase $e^{2\pi ik_{j}2^{j}/2^{n}}$. After doing this for
all of the remaining qubits, we get

$\frac{1}{2^{1/2}}\left(\ket{0}+e^{2\pi
ik/2^{n}}\ket{1}\right)\otimes\ket{k_{n-2}}\otimes\cdots\otimes\ket{k_{0}}.$
(48)

We perform the analogous process for qubits $n-2,\cdots,1$. For the qubit $0$,
we simply apply the Hadamard gate. In the end, we get the desired result,

$\frac{1}{2^{n/2}}\left(\ket{0}+e^{2\pi
ik2^{0}/2^{n}}\ket{1}\right)\otimes\cdots\otimes\left(\ket{0}+e^{2\pi
ik2^{n-1}/2^{n}}\ket{1}\right).$ (49)

As described above, we now apply swap gates between pairs of qubits $0$ and
$n-1$, $1$ and $n-2$, etc. and then we obtain the desired quantum Fourier
transform.

Iterative phase estimation performs the determination of the binary digits of
$\theta$ one at a time [95, 100, 101, 102, 103, 104]. Let $\ket{\psi}$ again
be an eigenstate of $U$ with eigenvalue $e^{2\pi i\theta}$. We first consider
the case where $\theta$ is equal to $k/2^{n}$ where $k$ is an integer between
$0$ and $2^{n}-1$. We start with $\ket{0}\otimes\ket{\psi}$ and apply a
Hadamard gate to obtain

$\frac{1}{2^{1/2}}\left(\ket{0}+\ket{1}\right)\otimes\ket{\psi}.$ (50)

We now use the ancilla qubit to perform the controlled unitary operator
$U^{2^{n-1}}$. The result is then

$\ket{f_{0}(\theta)}\otimes\ket{\psi},$ (51)

where

$\ket{f_{0}(\theta)}=\frac{1}{2^{1/2}}\left(\ket{0}+e^{2\pi i\theta
2^{n-1}}\ket{1}\right).$ (52)

Applying a Hadamard gate to $\ket{f_{0}(\theta)}$ gives

$\displaystyle\left(\frac{1+e^{2\pi i\theta
2^{n-1}}}{2}\ket{0}+\frac{1-e^{2\pi i\theta
2^{n-1}}}{2}\ket{1}\right)=\delta_{k_{0},0}\ket{0}+\delta_{k_{0},1}\ket{1}.$
(53)

Therefore, we can determine the digit $k_{0}$. Let us assume that we have
determined the digits from $k_{0},\cdots,k_{j-1}$. We can determine $k_{j}$ by
taking

$\frac{1}{2^{1/2}}\left(\ket{0}+\ket{1}\right)\otimes\ket{\psi}$ (54)

and using the ancilla qubit to perform the controlled unitary operator
$U^{2^{n-j-1}}$ followed by the phase gate

$\ket{0}\bra{0}+e^{-2\pi
i(k_{j-1}2^{-2}+\;\cdots\;+k_{0}2^{-j-1})}\ket{1}\bra{1},\\\ $ (55)

on the ancilla qubit. This phase gate removes the complex phase associated
with the binary digits $k_{0},\cdots,k_{j-1}$ that have already been
determined. The net result is

$\ket{f_{j}(\theta)}\otimes\ket{\psi},$ (56)

where

$\ket{f_{j}(\theta)}=\frac{1}{2^{1/2}}\left(\ket{0}+e^{2\pi i\theta
2^{n-j-1}-2\pi i(k_{j-1}2^{-2}+\;\cdots\;+k_{0}2^{-j-1})}\ket{1}\right).\\\ $
(57)

Applying a Hadamard gate to $\ket{f_{j}(\theta)}$ gives

$\displaystyle\delta_{k_{j},0}\ket{0}+\delta_{k_{j},1}\ket{1}.$ (58)

For the general case where $\theta$ is not equal to $k/2^{n}$ for some integer
$k$ between $0$ and $2^{n}-1$, there will be some distribution of values
associated with the measurements of the binary digits $k_{n-1},\cdots,k_{0}$.
As with regular phase estimation, the error in energy resolution scales
inversely with $2^{n}$ and is therefore inversely proportional to the number
of operations of $U$ needed. If $U$ is the time evolution of a Hamiltonian $H$
over time step $\Delta t$, then the error in the energy scales inversely with
the total time evolution required. Iterative phase estimation is not designed
to perform eigenstate preparation. If we start from a general linear
combination of energy eigenstates, then the uncertainty in the sequential
measurements of the binary digits $k_{n-1},\cdots,k_{0}$ arising from the
different eigenvalues $e^{2\pi i\theta_{a}}$ will prevent the algorithm from
functioning as intended.

The rodeo algorithm is another phase detection algorithm [105, 106, 107] that
shares some structural similarities with iterative phase estimation. In
contrast to iterative phase estimation, however, the rodeo algorithm is
efficient in preparing energy eigenstates starting from a general initial
state. Let $H$ be the Hamiltonian for which we want to prepare energy
eigenstates. To explain the algorithm, we first consider the case where the
initial state is an eigenstate of $H$. We call it $\ket{\psi_{j}}$ with
eigenvalue $E_{j}$. We use one ancilla qubit and start with the state

$\ket{0}\otimes\ket{\psi_{j}},$ (59)

and apply the Hadamard gate on the ancilla qubit,

$\frac{1}{2^{1/2}}\left(\ket{0}+\ket{1}\right)\otimes\ket{\psi_{j}}.$ (60)

We then use the ancilla to perform the controlled unitary for $e^{-iHt_{1}}$
and apply the phase gate

$\ket{0}\bra{0}+e^{iEt_{1}}\ket{1}\bra{1},\\\ $ (61)

on the ancilla. This produces

$\frac{1}{2^{1/2}}\left(\ket{0}+e^{-i(E_{j}-E)t_{1}}\ket{1}\right)\otimes\ket{\psi_{j}}.$
(62)

We now apply a Hadamard gate to the ancilla qubit, which then gives

$\frac{1}{2}\left[\left(1+e^{-i(E_{j}-E)t_{1}}\right)\ket{0}+\left(1-e^{-i(E_{j}-E)t_{1}}\right)\ket{1}\right]\otimes\ket{\psi_{j}}.$
(63)

If we measure the ancilla qubit, the probability of measuring $\ket{0}$ is
$\cos^{2}[(E_{j}-E)t_{1}/2]$ and the probability of measuring $\ket{1}$ is
$\sin^{2}[(E_{j}-E)t_{1}/2]$. We call the measurement of $\ket{0}$ a success
and the measurement of $\ket{1}$ a failure. We repeat this process for $n$
cycles with times $t_{1},\cdots,t_{n}$. The probability of success for all $n$
cycles is

$\prod_{k=1}^{n}\cos^{2}[(E_{j}-E)t_{k}/2].$ (64)

If we take random times $t_{1},\cdots,t_{n}$ to be chosen from a Gaussian
normal distribution with zero mean and $\sigma$ standard deviation, then the
success probability averaged over many trials will equal

$P_{n}(E)=\frac{\left[1+e^{-(E_{j}-E)^{2}\sigma^{2}/2}\right]^{n}}{2^{n}}.$
(65)

We see that the peak value is equal to $1$ when $E_{j}=E$ and the width of the
peak is $O(\sigma^{-1}n^{-1/2}).$

Let us now consider a general linear combination of energy eigenstates

$\ket{\psi}=\sum_{j}c_{j}\ket{\psi_{j}}.$ (66)

For this case, the probability of success for $n$ cycles is

$P_{n}(E)=\sum_{j}\frac{\left[1+e^{-(E_{j}-E)^{2}\sigma^{2}/2}\right]^{n}\left|c_{j}\right|^{2}}{2^{n}}.$
(67)

When scanning over the input parameter $E$, peaks in $P_{n}(E)$ will appear at
places where there are eigenvalues $E_{j}$ and the overlap with the initial
state is not too small. For fixed $n$, the error of the energy determination
scales inversely with $\sigma$. Similarly to phase estimation and iterative
phase estimation, the rodeo algorithm saturates the Heisenberg bound, where
the error in the energy scales inversely with the total duration of time
evolution.

In contrast with both phase estimation and iterative phase estimation, the
rodeo algorithm is exponentially fast for eigenstate preparation. There are
several other energy projection and filtering methods with similar
characteristics [108, 109, 110]. Once the peak of the eigenstate energy in
$P_{n}(E)$ is located approximately, we set $E$ as the peak value. With $E$
fixed and $\sigma$ fixed, the error estimates for the eigenvector scale as
$1/2^{n}$ for small $n$ and accelerate to $1/4^{n}$ for asymptotically large
values of $n$ [105]. The $1/2^{n}$ comes from the fact that the arithmetic
mean of $\cos^{2}(\theta)$ equals $1/2$, while the $1/4^{n}$ comes from the
fact that the geometric mean of $\cos^{2}(\theta)$ equals $1/4$. In Ref.
[111], it was shown that the use of progressive smaller values for the time
evolution parameters $t_{j}$ accelerates the convergence of the rodeo
algorithm towards $1/4^{n}$. The main limitation of the rodeo algorithm for
large quantum many-body systems is the requirement that the initial state have
non-negligible overlap with the eigenstate of interest. This is a difficult
problem that is common to nearly all eigenstate preparation algorithms that
use measurement projection. Nevertheless, one can use techniques such as
adiabatic evolution, variational methods, or some other approach as a
preconditioner to significantly increase the overlap with the eigenstate of
interest [105].

##  5 Summary and Outlook

In this article, we have presented several methods that show the essential
features of adiabatic evolution, variational methods, and phase detection
algorithms. All of the algorithms have their strengths and limitations, and
one common theme is that the techniques can be combined with each other to
produce something that is potentially greater than the sum of its parts. For
example, adiabatic evolution provides a theoretical foundation for the QAOA
variational method. In turn, the variational method can be used to find a good
starting Hamiltonian for adiabatic evolution. Both adiabatic evolution and
variation methods can be used as an initial-state preconditioner for phase
detection algorithms.

There has been great interest by both scientists and the general public on the
question of quantum advantage, if and when quantum computers are able to
perform tasks exceeding the capabilities of classical computers. It is
generally believed that calculations of real-time dynamics and spectral
functions of quantum many-body systems are areas ripe for possible quantum
advantage. However, the dynamics of some quantum many-body system starting
from a trivial initial state is not something that connects directly with
real-world phenomena. To make connections with real-world experiments and
observations, one also needs the ability to prepare energy eigenstates. It is
not clear whether quantum advantage will be achievable for the task of
eigenstate preparation. However, this may not be necessary. It may be enough
for quantum eigenstate preparation to be competitive with classical computing
methods to achieve quantum advantage for calculating the real-time dynamics or
spectral functions for real-world applications. The algorithms described in
this article provide some of the tools needed, but much more work is needed to
address the remaining challenges.

Acknowledgments The author acknowledges support from the U.S. Department of
Energy (DE-SC0021152, DE-SC0013365, DE-SC0023658) and the SciDAC-4 and
SciDAC-5 NUCLEI Collaborations. This research used resources of the Oak Ridge
Leadership Computing Facility, which is a DOE Office of Science User Facility
supported under Contract DE-AC05-00OR22725.

## References

  *   * Loh Jr et al. [1990] Loh Jr, E., Gubernatis, J., Scalettar, R., White, S., Scalapino, D., Sugar, R.: Sign problem in the numerical simulation of many-electron systems. Physical Review B 41(13), 9301 (1990)
  * Henelius and Sandvik [2000] Henelius, P., Sandvik, A.W.: Sign problem in Monte Carlo simulations of frustrated quantum spin systems. Physical Review B 62(2), 1102 (2000)
  * Alhassid et al. [2014] Alhassid, Y., Gilbreth, C.N., Bertsch, G.F.: Nuclear deformation at finite temperature. Phys. Rev. Lett. 113(26), 262503 (2014) https://doi.org/10.1103/PhysRevLett.113.262503 arXiv:1408.0081 [nucl-th]
  * Lähde et al. [2015] Lähde, T.A., Luu, T., Lee, D., Meißner, U.-G., Epelbaum, E., Krebs, H., Rupak, G.: Nuclear Lattice Simulations using Symmetry-Sign Extrapolation. Eur. Phys. J. A 51(7), 92 (2015) https://doi.org/10.1140/epja/i2015-15092-1 arXiv:1502.06787 [nucl-th]
  * Frame et al. [2018] Frame, D., He, R., Ipsen, I., Lee, D., Lee, D., Rrapaj, E.: Eigenvector continuation with subspace learning. Phys. Rev. Lett. 121(3), 032501 (2018) https://doi.org/10.1103/PhysRevLett.121.032501 arXiv:1711.07090 [nucl-th]
  * Lonardoni et al. [2018] Lonardoni, D., Gandolfi, S., Lynn, J.E., Petrie, C., Carlson, J., Schmidt, K.E., Schwenk, A.: Auxiliary field diffusion Monte Carlo calculations of light and medium-mass nuclei with local chiral interactions. Phys. Rev. C 97(4), 044318 (2018) https://doi.org/10.1103/PhysRevC.97.044318 arXiv:1802.08932 [nucl-th]
  * Chow et al. [2021] Chow, J., Dial, O., Gambetta, J.: Ibm quantum breaks the 100-qubit processor barrier. IBM Research Blog 2 (2021)
  * Pelofske et al. [2022] Pelofske, E., Bärtschi, A., Eidenbenz, S.: Quantum volume in practice: What users can expect from nisq devices. IEEE Transactions on Quantum Engineering 3, 1–19 (2022)
  * Neyenhuis [2023] Neyenhuis, B.: Quantinuum’s h-series quantum computers. Bulletin of the American Physical Society (2023)
  * Ayral et al. [2023] Ayral, T., Besserve, P., Lacroix, D., Ruiz Guzman, E.A.: Quantum computing with and for many-body physics (2023) arXiv:2303.04850 [quant-ph]
  * Boykin et al. [2002] Boykin, P.O., Mor, T., Roychowdhury, V., Vatan, F., Vrijen, R.: Algorithmic cooling and scalable nmr quantum computers. Proceedings of the National Academy of Sciences 99(6), 3388–3393 (2002)
  * Rempp et al. [2007] Rempp, F., Michel, M., Mahler, G.: Cyclic cooling algorithm. Physical Review A 76(3), 032325 (2007)
  * Lee et al. [2020] Lee, D., Bonitati, J., Given, G., Hicks, C., Li, N., Lu, B.-N., Rai, A., Sarkar, A., Watkins, J.: Projected Cooling Algorithm for Quantum Computation. Phys. Lett. B 807, 135536 (2020) https://doi.org/10.1016/j.physletb.2020.135536 arXiv:1910.07708 [quant-ph]
  * Gustafson [2020] Gustafson, E.: Projective Cooling for the transverse Ising model. Phys. Rev. D 101(7), 071504 (2020) https://doi.org/10.1103/PhysRevD.101.071504 arXiv:2002.06222 [hep-lat]
  * Ball and Cohen [2022] Ball, C., Cohen, T.D.: Boltzmann Distributions on a Quantum Computer via Active Cooling (2022) arXiv:2212.06730 [quant-ph]
  * Boykin et al. [2002] Boykin, P.O., Mor, T., Roychowdhury, V., Vatan, F., R., V.: Algorithmic cooling and scalable NMR quantum computers. Proc. Natl. Acad. Sci. 99, 3388 (2002) https://doi.org/10.1073/pnas.241641898
  * Xu et al. [2014] Xu, J.-S., Yung, M.-H., Xu, X.-Y., Boixo, S., Zhou, Z.-W., Li, C.-F., Aspuru-Guzik, A., Guo, G.-C.: Demon-like algorithmic quantum cooling and its realization with quantum optics. Nature Photonics 8, 113 (2014) https://doi.org/10.1038/nphoton.2013.354
  * Kraus et al. [2008] Kraus, B., Büchler, H.P., Diehl, S., Kantian, A., Micheli, A., Zoller, P.: Preparation of entangled states by quantum Markov processes. Phys. Rev. A 78, 042307 (2008) https://doi.org/10.1103/PhysRevA.78.042307
  * Verstraete et al. [2009] Verstraete, F., Wolf, M.M., Ignacio Cirac, J.: Quantum computation and quantum-state engineering driven by dissipation. Nature Physics 5, 633 (2009) https://doi.org/10.1038/nphys1342
  * Sinayskiy and Petruccione [2012] Sinayskiy, I., Petruccione, F.: Efficiency of open quantum walk implementation of dissipative quantum computing algorithms. Quantum Information Processing 11, 1301–1309 (2012)
  * Kaplan et al. [2017] Kaplan, D.B., Klco, N., Roggero, A.: Ground States via Spectral Combing on a Quantum Computer (2017) arXiv:1709.08250 [quant-ph]
  * Yen et al. [2019] Yen, T.-C., Lang, R.A., Izmaylov, A.F.: Exact and approximate symmetry projectors for the electronic structure problem on a quantum computer. The Journal of chemical physics 151(16) (2019)
  * Siwach and Lacroix [2021] Siwach, P., Lacroix, D.: Filtering states with total spin on a quantum computer. Phys. Rev. A 104(6), 062435 (2021) https://doi.org/10.1103/PhysRevA.104.062435 arXiv:2106.10867 [quant-ph]
  * Ruiz Guzman and Lacroix [2022] Ruiz Guzman, E.A., Lacroix, D.: Accessing ground-state and excited-state energies in a many-body system after symmetry restoration using quantum computers. Phys. Rev. C 105(2), 024324 (2022) https://doi.org/10.1103/PhysRevC.105.024324 arXiv:2111.13080 [quant-ph]
  * Seki and Yunoki [2022] Seki, K., Yunoki, S.: Spatial, spin, and charge symmetry projections for a fermi-hubbard model on a quantum computer. Physical Review A 105(3), 032419 (2022)
  * Lacroix et al. [2023] Lacroix, D., Ruiz Guzman, E.A., Siwach, P.: Symmetry breaking/symmetry preserving circuits and symmetry restoration on quantum computers: A quantum many-body perspective. Eur. Phys. J. A 59(1), 3 (2023) https://doi.org/10.1140/epja/s10050-022-00911-7 arXiv:2208.11567 [quant-ph]
  * Childs and Wiebe [2012] Childs, A.M., Wiebe, N.: Hamiltonian Simulation Using Linear Combinations of Unitary Operations. arXiv e-prints, 1202–5822 (2012) https://doi.org/10.48550/arXiv.1202.5822 arXiv:1202.5822 [quant-ph]
  * Roggero et al. [2020] Roggero, A., Gu, C., Baroni, A., Papenbrock, T.: Preparation of excited states for nuclear dynamics on a quantum computer. Phys. Rev. C 102(6), 064624 (2020) https://doi.org/10.1103/PhysRevC.102.064624 arXiv:2009.13485 [quant-ph]
  * Lv et al. [2023] Lv, P., Wei, S.-J., Xie, H.-N., Long, G.-L.: QCSH: A full quantum computer nuclear shell-model package. Sci. China Phys. Mech. Astron. 66(4), 240311 (2023) https://doi.org/10.1007/s11433-022-2044-5 arXiv:2205.12087 [quant-ph]
  * Siwach and Arumugam [2022] Siwach, P., Arumugam, P.: Quantum computation of nuclear observables involving linear combinations of unitary operators. Phys. Rev. C 105(6), 064318 (2022) https://doi.org/10.1103/PhysRevC.105.064318 arXiv:2206.08510 [quant-ph]
  * Meister et al. [2022] Meister, R., Benjamin, S.C., Campbell, E.T.: Tailoring term truncations for electronic structure calculations using a linear combination of unitaries. Quantum 6, 637 (2022)
  * Motta et al. [2019] Motta, M., Sun, C., Tan, A.T.K., O’Rourke, M.J., Ye, E., Minnich, A.J., Brandão, F.G.S.L., Chan, G.K.-L.: Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution. Nature Physics 16(2), 205–210 (2019) https://doi.org/10.1038/s41567-019-0704-4
  * McArdle et al. [2019] McArdle, S., Jones, T., Endo, S., Li, Y., Benjamin, S.C., Yuan, X.: Variational ansatz-based quantum simulation of imaginary time evolution. npj Quantum Information 5(1), 75 (2019)
  * Yeter-Aydeniz et al. [2020] Yeter-Aydeniz, K., Pooser, R.C., Siopsis, G.: Practical quantum computation of chemical and nuclear energy levels using quantum imaginary time evolution and lanczos algorithms. npj Quantum Information 6(1), 63 (2020)
  * Turro et al. [2022] Turro, F., Roggero, A., Amitrano, V., Luchi, P., Wendt, K.A., DuBois, J.L., Quaglioni, S., Pederiva, F.: Imaginary-time propagation on a quantum chip. Phys. Rev. A 105(2), 022440 (2022) https://doi.org/10.1103/PhysRevA.105.022440 arXiv:2102.12260 [quant-ph]
  * Jouzdani et al. [2022] Jouzdani, P., Johnson, C.W., Mucciolo, E.R., Stetcu, I.: Alternative approach to quantum imaginary time evolution. Phys. Rev. A 106(6), 062435 (2022) https://doi.org/10.1103/PhysRevA.106.062435 arXiv:2208.10535 [quant-ph]
  * Born and Fock [1950] Born, M., Fock, V.: Beweis des adiabatensatzes. Zeitschrift für Physik 51, 435–439 (1950) https://doi.org/10.1007/BF01343193
  * Kato [1950] Kato, T.: On the adiabatic theorem of quantum mechanics. Journal of the Physical Society of Japan 5(6), 435–439 (1950) https://doi.org/10.1143/JPSJ.5.435 https://doi.org/10.1143/JPSJ.5.435
  * Farhi et al. [2000] Farhi, E., Goldstone, J., Gutmann, S., Sipser, M.: Quantum Computation by Adiabatic Evolution. arXiv e-prints, 0001106 (2000) arXiv:quant-ph/0001106 [quant-ph]
  * Farhi et al. [2001] Farhi, E., Goldstone, J., Gutmann, S., Lapan, J., Lundgren, A., Preda, D.: A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem. Science 292(5516), 472–475 (2001) https://doi.org/10.1126/science.1057726 https://www.science.org/doi/pdf/10.1126/science.1057726
  * Childs et al. [2001] Childs, A.M., Farhi, E., Preskill, J.: Robustness of adiabatic quantum computation. Physical Review A 65(1), 012322 (2001)
  * Van Dam et al. [2001] Van Dam, W., Mosca, M., Vazirani, U.: How powerful is adiabatic quantum computation? In: Proceedings 42nd IEEE Symposium on Foundations of Computer Science, pp. 279–287 (2001). IEEE
  * Farhi et al. [2002] Farhi, E., Goldstone, J., Gutmann, S.: Quantum adiabatic evolution algorithms with different paths. arXiv preprint quant-ph/0208135 (2002)
  * Roland and Cerf [2002] Roland, J., Cerf, N.J.: Quantum search by local adiabatic evolution. Physical Review A 65(4), 042308 (2002)
  * Santoro and Tosatti [2006] Santoro, G.E., Tosatti, E.: Optimization using quantum mechanics: quantum annealing through adiabatic evolution. Journal of Physics A: Mathematical and General 39(36), 393 (2006)
  * Wiebe and Babcock [2012] Wiebe, N., Babcock, N.S.: Improved error-scaling for adiabatic quantum evolutions. New Journal of Physics 14(1), 013024 (2012) https://doi.org/10.1088/1367-2630/14/1/013024 arXiv:1105.6268 [quant-ph]
  * Richerme et al. [2013] Richerme, P., Senko, C., Smith, J., Lee, A., Korenblit, S., Monroe, C.: Experimental performance of a quantum simulator: Optimizing adiabatic evolution and identifying many-body ground states. Physical Review A 88(1), 012334 (2013)
  * Albash and Lidar [2018] Albash, T., Lidar, D.A.: Adiabatic quantum computation. Rev. Mod. Phys. 90, 015002 (2018) https://doi.org/10.1103/RevModPhys.90.015002
  * Coello Pérez et al. [2022] Coello Pérez, E.A., Bonitati, J., Lee, D., Quaglioni, S., Wendt, K.A.: Quantum state preparation by adiabatic evolution with custom gates. Phys. Rev. A 105(3), 032403 (2022) https://doi.org/10.1103/PhysRevA.105.032403 arXiv:2111.12207 [quant-ph]
  * van Vreumingen and Schoutens [2023] Vreumingen, D., Schoutens, K.: Adiabatic ground state preparation of fermionic many-body systems from a two-body perspective (2023) arXiv:2305.01284 [quant-ph]
  * Elhatisari et al. [2022] Elhatisari, S., et al.: Wave function matching for the quantum many-body problem (2022) arXiv:2210.17488 [nucl-th]
  * Avron et al. [1987] Avron, J., Seiler, R., Yaffe, L.: Adiabatic theorems and applications to the quantum hall effect. Communications in Mathematical Physics 110, 33–49 (1987)
  * Avron et al. [1993] Avron, J., Seiler, R., Yaffe, L.: Adiabatic theorems and applications to the quantum hall effect. Communications in Mathematical Physics 156, 649–649 (1993)
  * Klein and Seiler [1990] Klein, M., Seiler, R.: Power-law corrections to the kubo formula vanish in quantum hall systems. Communications in mathematical physics 128, 141–160 (1990)
  * Jansen et al. [2007] Jansen, S., Ruskai, M.-B., Seiler, R.: Bounds for the adiabatic approximation with applications to quantum computation. Journal of Mathematical Physics 48(10) (2007) https://doi.org/10.1063/1.2798382 https://pubs.aip.org/aip/jmp/article-pdf/doi/10.1063/1.2798382/16055030/102111_1_online.pdf. 102111
  * Sarkar et al. [2023] Sarkar, A., Lee, D., Meißner, U.-G.: Floating block method for quantum Monte Carlo simulations (2023) arXiv:2306.11439 [nucl-th]
  * Trotter [1959] Trotter, H.F.: On the product of semi-groups of operators. Proceedings of the American Mathematical Society 10(4), 545–551 (1959)
  * Suzuki [1993] Suzuki, M.: General decomposition theory of ordered exponentials. Proceedings of the Japan Academy, Series B 69(7), 161–166 (1993) https://doi.org/10.2183/pjab.69.161
  * Childs et al. [2019] Childs, A.M., Ostrander, A., Su, Y.: Faster quantum simulation by randomization. Quantum 3, 182 (2019)
  * Campbell [2019] Campbell, E.: Random compiler for fast hamiltonian simulation. Physical review letters 123(7), 070503 (2019)
  * Low et al. [2019] Low, G.H., Kliuchnikov, V., Wiebe, N.: Well-conditioned multiproduct hamiltonian simulation. arXiv preprint arXiv:1907.11679 (2019)
  * Watkins et al. [2022] Watkins, J., Wiebe, N., Roggero, A., Lee, D.: Time-dependent hamiltonian simulation using discrete clock constructions. arXiv preprint arXiv:2203.11353 (2022)
  * Faehrmann et al. [2022] Faehrmann, P.K., Steudtner, M., Kueng, R., Kieferova, M., Eisert, J.: Randomizing multi-product formulas for hamiltonian simulation. Quantum 6, 806 (2022)
  * Rendon et al. [2022] Rendon, G., Watkins, J., Wiebe, N.: Improved error scaling for trotter simulations through extrapolation. arXiv preprint arXiv:2212.14144 (2022)
  * Carena et al. [2021] Carena, M., Lamm, H., Li, Y.-Y., Liu, W.: Lattice renormalization of quantum simulations. Physical Review D 104(9), 094519 (2021)
  * Parrish and McMahon [2019] Parrish, R.M., McMahon, P.L.: Quantum Filter Diagonalization: Quantum Eigendecomposition without Full Quantum Phase Estimation. arXiv e-prints, 1909–08925 (2019) https://doi.org/10.48550/arXiv.1909.08925 arXiv:1909.08925 [quant-ph]
  * Francis et al. [2022] Francis, A., Agrawal, A.A., Howard, J.H., Kökcü, E., Kemper, A.F.: Subspace Diagonalization on Quantum Computers using Eigenvector Continuation. arXiv e-prints, 2209–10571 (2022) https://doi.org/10.48550/arXiv.2209.10571 arXiv:2209.10571 [quant-ph]
  * Yang et al. [2022] Yang, Y., Zhang, Z., Xu, X., Lu, B., Li, Y.: Quantum algorithms for optimal effective theory of many-body systems. arXiv e-prints, 2211–14854 (2022) https://doi.org/10.48550/arXiv.2211.14854 arXiv:2211.14854 [quant-ph]
  * Cortes and Gray [2022] Cortes, C.L., Gray, S.K.: Quantum Krylov subspace algorithms for ground- and excited-state energy estimation. Phys. Rev. A 105(2), 022417 (2022) https://doi.org/10.1103/PhysRevA.105.022417 arXiv:2109.06868 [quant-ph]
  * Mejuto-Zaera and Kemper [2023] Mejuto-Zaera, C., Kemper, A.F.: Quantum Eigenvector Continuation for Chemistry Applications. arXiv e-prints, 2305–00060 (2023) https://doi.org/10.48550/arXiv.2305.00060 arXiv:2305.00060 [physics.chem-ph]
  * Peruzzo et al. [2014] Peruzzo, A., McClean, J., Shadbolt, P., Yung, M.-H., Zhou, X.-Q., Love, P.J., Aspuru-Guzik, A., O’brien, J.L.: A variational eigenvalue solver on a photonic quantum processor. Nature communications 5(1), 4213 (2014)
  * Cerezo et al. [2021] Cerezo, M., Arrasmith, A., Babbush, R., Benjamin, S.C., Endo, S., Fujii, K., McClean, J.R., Mitarai, K., Yuan, X., Cincio, L., et al.: Variational quantum algorithms. Nature Reviews Physics 3(9), 625–644 (2021)
  * McClean et al. [2018] McClean, J.R., Boixo, S., Smelyanskiy, V.N., Babbush, R., Neven, H.: Barren plateaus in quantum neural network training landscapes. Nature communications 9(1), 4812 (2018)
  * Mari et al. [2021] Mari, A., Bromley, T.R., Killoran, N.: Estimating the gradient and higher-order derivatives on quantum hardware. Physical Review A 103(1), 012405 (2021)
  * Farhi et al. [2014] Farhi, E., Goldstone, J., Gutmann, S.: A quantum approximate optimization algorithm. arXiv preprint arXiv:1411.4028 (2014)
  * Karp [2010] Karp, R.M.: Reducibility Among Combinatorial Problems. Springer, ??? (2010)
  * Bittel and Kliesch [2021] Bittel, L., Kliesch, M.: Training variational quantum algorithms is np-hard. Phys. Rev. Lett. 127, 120502 (2021) https://doi.org/10.1103/PhysRevLett.127.120502
  * Kandala et al. [2017] Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Brink, M., Chow, J.M., Gambetta, J.M.: Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets. nature 549(7671), 242–246 (2017)
  * Dumitrescu et al. [2018] Dumitrescu, E.F., McCaskey, A.J., Hagen, G., Jansen, G.R., Morris, T.D., Papenbrock, T., Pooser, R.C., Dean, D.J., Lougovski, P.: Cloud Quantum Computing of an Atomic Nucleus. Phys. Rev. Lett. 120, 210501 (2018) https://doi.org/10.1103/PhysRevLett.120.210501 arXiv:1801.03897 [quant-ph]
  * Qian et al. [2021] Qian, Z., Watkins, J., Given, G., Bonitati, J., Choi, K., Lee, D.: Demonstration of the Rodeo Algorithm on a Quantum Computer (2021) arXiv:2110.07747 [quant-ph]
  * Kiss et al. [2022] Kiss, O., Grossi, M., Lougovski, P., Sanchez, F., Vallecorsa, S., Papenbrock, T.: Quantum computing of the 6Li nucleus via ordered unitary coupled cluster. Phys. Rev. C 106(3), 034325 (2022) https://doi.org/10.1103/PhysRevC.106.034325 arXiv:2205.00864 [nucl-th]
  * Tilly et al. [2022] Tilly, J., Chen, H., Cao, S., Picozzi, D., Setia, K., Li, Y., Grant, E., Wossnig, L., Rungger, I., Booth, G.H., et al.: The variational quantum eigensolver: a review of methods and best practices. Physics Reports 986, 1–128 (2022)
  * Hlatshwayo et al. [2022] Hlatshwayo, M.Q., Zhang, Y., Wibowo, H., LaRose, R., Lacroix, D., Litvinova, E.: Simulating excited states of the Lipkin model on a quantum computer. Phys. Rev. C 106(2), 024319 (2022) https://doi.org/10.1103/PhysRevC.106.024319 arXiv:2203.01478 [nucl-th]
  * Pérez-Obiol et al. [2023] Pérez-Obiol, A., Romero, A.M., Menéndez, J., Rios, A., García-Sáez, A., Juliá-Díaz, B.: Nuclear shell-model simulation in digital quantum computers (2023) arXiv:2302.03641 [quant-ph]
  * Sarma et al. [2023] Sarma, C., Di Matteo, O., Abhishek, A., Srivastava, P.C.: Prediction of the neutron drip line in oxygen isotopes using quantum computation (2023) arXiv:2306.06432 [nucl-th]
  * Jordan and Wigner [1993] Jordan, P., Wigner, E.P.: Über Das Paulische Äquivalenzverbot. Springer, ??? (1993)
  * Bravyi and Kitaev [2002] Bravyi, S.B., Kitaev, A.Y.: Fermionic quantum computation. Annals of Physics 298(1), 210–226 (2002) https://doi.org/10.1006/aphy.2002.6254
  * Verstraete and Cirac [2005] Verstraete, F., Cirac, J.I.: Mapping local hamiltonians of fermions to local hamiltonians of spins. Journal of Statistical Mechanics: Theory and Experiment 2005(09), 09012 (2005)
  * Whitfield et al. [2016] Whitfield, J.D., Havlíček, V., Troyer, M.: Local spin operators for fermion simulations. Physical Review A 94(3), 030301 (2016)
  * Steudtner and Wehner [2018] Steudtner, M., Wehner, S.: Fermion-to-qubit mappings with varying resource requirements for quantum simulation. New Journal of Physics 20(6), 063010 (2018)
  * Derby et al. [2021] Derby, C., Klassen, J., Bausch, J., Cubitt, T.: Compact fermion to qubit mappings. Physical Review B 104(3), 035118 (2021)
  * Nys and Carleo [2023] Nys, J., Carleo, G.: Quantum circuits for solving local fermion-to-qubit mappings. Quantum 7, 930 (2023)
  * Thouless [1960] Thouless, D.J.: Stability conditions and nuclear rotations in the hartree-fock theory. Nuclear Physics 21, 225–232 (1960)
  * Quantum et al. [2020] Quantum, G.A., Collaborators*†, Arute, F., Arya, K., Babbush, R., Bacon, D., Bardin, J.C., Barends, R., Boixo, S., Broughton, M., Buckley, B.B., Buell, D.A., Burkett, B., Bushnell, N., Chen, Y., Chen, Z., Chiaro, B., Collins, R., Courtney, W., Demura, S., Dunsworth, A., Farhi, E., Fowler, A., Foxen, B., Gidney, C., Giustina, M., Graff, R., Habegger, S., Harrigan, M.P., Ho, A., Hong, S., Huang, T., Huggins, W.J., Ioffe, L., Isakov, S.V., Jeffrey, E., Jiang, Z., Jones, C., Kafri, D., Kechedzhi, K., Kelly, J., Kim, S., Klimov, P.V., Korotkov, A., Kostritsa, F., Landhuis, D., Laptev, P., Lindmark, M., Lucero, E., Martin, O., Martinis, J.M., McClean, J.R., McEwen, M., Megrant, A., Mi, X., Mohseni, M., Mruczkiewicz, W., Mutus, J., Naaman, O., Neeley, M., Neill, C., Neven, H., Niu, M.Y., O’Brien, T.E., Ostby, E., Petukhov, A., Putterman, H., Quintana, C., Roushan, P., Rubin, N.C., Sank, D., Satzinger, K.J., Smelyanskiy, V., Strain, D., Sung, K.J., Szalay, M., Takeshita, T.Y., Vainsencher, A., White, T., Wiebe, N., Yao, Z.J., Yeh, P., Zalcman, A.: Hartree-fock on a superconducting qubit quantum computer. Science 369(6507), 1084–1089 (2020) https://doi.org/10.1126/science.abb9811 https://www.science.org/doi/pdf/10.1126/science.abb9811
  * Kitaev [1995] Kitaev, A.Yu.: Quantum measurements and the Abelian stabilizer problem. Electronic Colloquium on Computational Complexity 3 (1995) arXiv:quant-ph/9511026 [quant-ph]
  * Abrams and Lloyd [1997] Abrams, D.S., Lloyd, S.: Simulation of many body Fermi systems on a universal quantum computer. Phys. Rev. Lett. 79, 2586–2589 (1997) https://doi.org/10.1103/PhysRevLett.79.2586 arXiv:quant-ph/9703054
  * Dorner et al. [2009] Dorner, U., Demkowicz-Dobrzanski, R., Smith, B.J., Lundeen, J.S., Wasilewski, W., Banaszek, K., Walmsley, I.A.: Optimal quantum phase estimation. Physical review letters 102(4), 040403 (2009)
  * Svore et al. [2013] Svore, K.M., Hastings, M.B., Freedman, M.: Faster phase estimation. arXiv preprint arXiv:1304.0741 (2013)
  * Nielsen and Chuang [2010] Nielsen, M.A., Chuang, I.L.: Quantum Computation and Quantum Information. Cambridge university press, ??? (2010)
  * Dobšíček et al. [2007] Dobšíček, M., Johansson, G., Shumeiko, V., Wendin, G.: Arbitrary accuracy iterative quantum phase estimation algorithm using a single ancillary qubit: A two-qubit benchmark. Physical Review A 76(3), 030306 (2007)
  * O’Loan [2009] O’Loan, C.J.: Iterative phase estimation. Journal of Physics A: Mathematical and Theoretical 43(1), 015301 (2009)
  * Lanyon et al. [2010] Lanyon, B.P., Whitfield, J.D., Gillett, G.G., Goggin, M.E., Almeida, M.P., Kassal, I., Biamonte, J.D., Mohseni, M., Powell, B.J., Barbieri, M., et al.: Towards quantum chemistry on a quantum computer. Nature chemistry 2(2), 106–111 (2010)
  * Wiebe and Granade [2016] Wiebe, N., Granade, C.: Efficient bayesian phase estimation. Physical review letters 117(1), 010503 (2016)
  * Ding and Lin [2023] Ding, Z., Lin, L.: Even shorter quantum circuit for phase estimation on early fault-tolerant quantum computers with applications to ground-state energy estimation. PRX Quantum 4(2), 020331 (2023)
  * Choi et al. [2021] Choi, K., Lee, D., Bonitati, J., Qian, Z., Watkins, J.: Rodeo Algorithm for Quantum Computing. Phys. Rev. Lett. 127(4), 040505 (2021) https://doi.org/10.1103/PhysRevLett.127.040505 arXiv:2009.04092 [quant-ph]
  * Bee-Lindgren et al. [2022] Bee-Lindgren, M., et al.: Rodeo algorithm with controlled reversal gates. work in progress (2022)
  * Qian et al. [2022] Qian, W., Basili, R., Pal, S., Luecke, G., Vary, J.P.: Solving hadron structures using the basis light-front quantization approach on quantum computers. Phys. Rev. Res. 4(4), 043193 (2022) https://doi.org/10.1103/PhysRevResearch.4.043193 arXiv:2112.01927 [quant-ph]
  * Ge et al. [2019] Ge, Y., Tura, J., Cirac, J.I.: Faster ground state preparation and high-precision ground energy estimation with fewer qubits. Journal of Mathematical Physics 60(2) (2019)
  * Lu et al. [2021] Lu, S., Banuls, M.C., Cirac, J.I.: Algorithms for quantum simulation at finite energies. PRX Quantum 2(2), 020321 (2021)
  * Stetcu et al. [2022] Stetcu, I., Baroni, A., Carlson, J.: Projection algorithm for state preparation on quantum computers (2022) arXiv:2211.10545 [quant-ph]
  * Cohen and Oh [2023] Cohen, T.D., Oh, H.: Optimizing rodeo projection (2023) arXiv:2305.19952 [quant-ph]


# Improved Power Decoding of
Algebraic Geometry Codes

Sven Puchinger, Johan Rosenkilde, Grigory Solomatov  Department of Applied
Mathematics and Computer Science,
Technical University of Denmark (DTU), Denmark

###### Abstract

Power decoding is a partial decoding paradigm for arbitrary algebraic geometry
codes for decoding beyond half the minimum distance, which usually returns the
unique closest codeword, but in rare cases fails to return anything. The
original version decodes roughly up to the Sudan radius, while an improved
version decodes up to the Johnson radius, but has so far been described only
for Reed–Solomon and one-point Hermitian codes. In this paper we show how the
improved version can be applied to any algebraic geometry code.

###### Index Terms:

Algebraic Geometry Codes, Power Decoding

##  I Introduction

###  I-A Related work

Power decoding was originally proposed in [1] and [2] as a method for decoding
Reed–Solomon codes beyond half the minimum distance, achieving roughly the
same decoding radius as Sudan’s list decoder [3]. It is a partial unique
decoder, meaning that it either returns a unique decoding solution, or
declares failure. It has been demonstrated through numerical simulations that
power decoding has a very low probability of failure for random errors; there
are proofs of this for small parameter choices [2, 4], but the general case
remains to be an open problem.

In [5] it was shown how to incorporate a multiplicity parameter into power
decoding, similar to that of Guruswami–Sudan (GS) algorithm [6], thereby
matching the latter’s decoding radius, called the _Johnson radius_. Both the
original and the improved power decoding with multiplicity has been applied to
one-point Hermitian codes, which belong to the large family of _algebraic
geometry_ (AG) _codes_ that were introduced by Goppa in 1981 [7, 8]. They are
a natural generalization of Reed–Solomon codes and some of their subfamilies
have been shown to have minimal distance that beats the Gilbert–Varshamov
bound [9].

At the heart of power decoding lies a set of “key equations”, which can be
solved by linear algebra, or by faster, more structured approaches, e.g. by
generalizing the Berlekamp–Massey or Euclidean algorithm, see e.g. [10, 11]
and their references. This step is algebraically connected with to the
“interpolation step” of the GS algorithm, which may be computed in similar
complexity [12, 11].

After interpolation, the GS algorithm proceeds with a “root-finding step”,
which is not required by power decoding. For Reed–Solomon and one-point
Hermitian codes, root-finding is asymptotically faster than the interpolation
step111Except if the code length is very short compared to the field size: a
sub-routine of the root-finding step is to find $\mathbb{F}_{q}$-roots in
$\mathbb{F}_{q}[x]$ polynomials.  [13, 12, 14], but in hardware-
implementations root-finding can take up significant circuit area [15]. For
more general AG codes, the root-finding step is not widely studied, and it is
unclear whether it is faster than interpolation; e.g. [16] describes root-
finding for one-point codes which is not superior to performing interpolation
using linear algebra. Further, both power decoding and the GS algorithm
generalize to interleaved codes, and here the root-finding step is currently
the bottleneck in the GS algorithm [17], and so power decoding is
asymptotically faster [18].

Recently, [19] generalized the original power decoding to apply to general AG
codes. Their work was principally motivated by their new variant of (the
original) power decoding, called power error-locating pairs, which may be
applied to attack AG-code-based cryptosystems in a setting where it seems the
GS decoder can not be used.

###  I-B Contributions

In this paper we:

  * •

formulate improved power decoding [5] in the language of function fields and
adapt it to general AG codes,

  * •

show how the resulting _key equations_ can be solved using linear algebra, and

  * •

derive the decoding radius for the proposed method under mild assumptions on
the code and verify it using numerical simulations.

The most important takeaway is that our decoding radius coincides with the one
from [20] for Reed–Solomon codes and the one from [21] for one-point Hermitian
codes, achieving the Johnson radius for suitable decoding parameters. Though
the computational complexity of our decoder can likely be improved by
replacing linear algebra with more advanced techniques, this remains outside
of the scope of this paper.

Our function field modelling revolves around “using an extra place”
$P_{\infty}$ as in Pellikaan [22], which allows controlling the “size” of
functions by their pole order at $P_{\infty}$. Our key equations are in the
style of Gao [23], as introduced for power decoding in [4].

We present evidence supporting that our decoder achieves a decoding radius
similar to the GS algorithm in two ways: we theoretically derive a radius at
which the decoder will surely fail, and then we see in simulations that for
fewer random errors, the decoder seems to almost always succeed. We have no
theoretic bound on this failure probability, which has so far proved very
challenging even in the case of Reed–Solomon codes, see e.g. [20]. On the
other hand, advantages of the power decoding compared to GS algorithm include:

  * •

The power decoder is structurally simpler than the GS algorithm since it does
not have a root-finding step, making it easier to implement, possibly faster,
and with less footprint in hardware.

  * •

It is reasonable that a generalization to interleaved AG codes is possible as
in [18, 21].

  * •

Potentially apply power decoding in a more efficient message attack on AG-
based cryptosystems than in [19], possibly by first reformulating it in terms
of power error-locating pairs.

##  II Preliminaries

In this section we briefly introduce AG codes as _evaluation codes_ and
formulate the decoding problem that we wish to solve. We assume that the
reader is familiar with the language of function fields, which is presented in
great detail in [24].

Let $F$ be a function field having genus $g$ over a base field $\mathbb{F}$.
Let $P_{1},\dots,P_{n},P_{\infty}$ be rational places of $F$, and fix two
divisors $D=P_{1}+\cdots+P_{n}$ and $G$ with $\deg G=\gamma$, such that
$\operatorname{supp}G\cap\operatorname{supp}D=\emptyset$. The code that we
will be considering is

$\mathcal{C}_{D,G}=\\{\big{(}f(P_{1}),\dots,f(P_{n})\big{)}\in\mathbb{F}^{n}\;|\;f\in\mathcal{L}(G)\\}\
,$

where $\mathcal{L}(G)=\\{f\in F\;|\;(f)\geq-G\\}\cup\\{0\\}$ is the _Riemann-
Roch space_ of $G$. The dimension of this code is $k=\mathfrak{l}(G)$, where
$\mathfrak{l}(A)$ denotes the dimension of $\mathcal{L}(A)$ for any divisor
$A$; and the code’s minimum distance is bounded from below by the _designed
minimum distance_ $d^{*}=n-\gamma$.

The problem that we wish to address is as follows:

######  Problem II.1.

Let $f\in\mathcal{L}(G)$ be a _message_ and let
$\boldsymbol{c}\in\mathcal{C}_{D,G}$ be the corresponding _codeword_. Given
the _received word_

$\boldsymbol{r}=(r_{1},\dots,r_{n})=\boldsymbol{c}+\boldsymbol{e}\ ,$

where $\boldsymbol{e}=(e_{1},\dots,e_{n})\in\mathbb{F}^{n}$ is some unknown
error, recover $f$.

In the next section, in the traditional spirit of power decoding, we show how
II.1 can be reformulated as a highly non-linear system of equations.

##  III The Key Equations

In this section we formulate the so called _key equations_ that lie at the
heart of power decoding using the language of function fields.

We begin by defining the _error locator_ and the _received word interpolator_
, following up with some results about their sizes.

######  Definition III.1.

If $\boldsymbol{r}=(r_{1},\dots,r_{n})$ is the received word, then an
$\boldsymbol{r}$-_interpolator_ of degree $\rho$ is an element

$\displaystyle R$ $\displaystyle\in\mathcal{L}(G+\rho
P_{\infty})\quad\text{such that}$ $\displaystyle R$
$\displaystyle\not\in\mathcal{L}\big{(}G+(\rho-1)P_{\infty}\big{)}\ ,$

and $R(P_{i})=r_{i}$ for $i=1,\dots,n$.

######  Lemma III.2.

There exists an ${\boldsymbol{r}}$-interpolator $R$ with degree $\rho$
satisfying $\rho\leq n-\gamma+2g-1$.

###### Proof.

Consider the $\mathbb{F}$-linear map

$\displaystyle\phi:\begin{cases}\mathcal{L}(G+\rho
P_{\infty})\rightarrow\mathbb{F}^{n}\\\
h\mapsto\big{(}h(P_{1}),\dots,h(P_{n})\big{)}\end{cases}$

We are guaranteed that $R$ exists if the dimension of the image of $\phi$,
which is given by the rank-nullity theorem as

$\dim\operatorname{im}\phi=\mathfrak{l}(G+\rho P_{\infty})-\mathfrak{l}(G+\rho
P_{\infty}-D)\ ,$

is $n$. Indeed, Riemann’s theorem theorem says that

$\mathfrak{l}(G+\rho P_{\infty})\geq\gamma+\rho+1-g\ ,$

and if $\rho\geq n-\gamma+2g-1$, then [24, Thm 1.5.17] promises that

$\mathfrak{l}(G+\rho P_{\infty}-D)=\gamma+\rho-n+1-g\ ,$

which shows that the image of $\phi$ has dimension $n$. ∎

######  Definition III.3.

Let $E=\sum_{i\in\mathcal{E}}P_{i}$, where $\mathcal{E}=\\{i\;|\;e_{i}\neq
0\\}$ is the set of error positions. An _error locator_ with multiplicity $s$
and degree $\lambda_{s}$ is an element

$\displaystyle\Lambda_{s}$
$\displaystyle\in\mathcal{L}(\lambda_{s}P_{\infty}-sE)\quad\text{such that}$
$\displaystyle\Lambda_{s}$
$\displaystyle\not\in\mathcal{L}\big{(}(\lambda_{s}-1)P_{\infty}-sE\big{)}\ .$

######  Lemma III.4.

For any $\lambda_{s}\geq s|\mathcal{E}|+g$ there exists an error locator with
degree at most $\lambda_{s}$, however, if $\lambda_{s}<s|\mathcal{E}|$, then
no such error locator exists.

###### Proof.

If $\lambda_{s}\geq s|\mathcal{E}|+g$, then the Riemann-Roch theorem promises
that

$\mathfrak{l}(\lambda_{s}P_{\infty}-sE)\geq\lambda_{s}-s|\mathcal{E}|+1-g\geq
0\ ,$

which guarantees existence of an error locator with degree at most
$\lambda_{s}$.

If on the other hand $\lambda_{s}<s|\mathcal{E}|$, then $\Lambda_{s}$ must
have more zeroes than poles, which is impossible. ∎

It will be convenient to define the following divisors:

$\displaystyle V_{t}$ $\displaystyle=\lambda_{s}P_{\infty}+tG\ ,$
$\displaystyle Q_{t}$ $\displaystyle=\lambda_{s}P_{\infty}+t(G+\rho
P_{\infty})\quad\text{and}$ $\displaystyle W_{j}$
$\displaystyle=\lambda_{s}P_{\infty}+j(G+\rho P_{\infty}-D),$

for $j=0,\dots,s-1$ and $t=1,\dots,\ell$.

The next lemma relates these divisors to $f$, $\Lambda_{s}$ and $R$.

######  Lemma III.5.

If $f$ is the sent message, $R$ is an $\boldsymbol{r}$-interpolator and
$\Lambda_{s}$ is an error locator, then for $t=1,\dots,\ell$ it holds that

$\Lambda_{s}f^{t}\in\mathcal{L}(V_{t})\ ,$

and for $j=0,\dots,s-1$ it holds that

$\Lambda_{s}(f-R)^{j}R^{t-j}\in\begin{cases}\mathcal{L}\big{(}W_{j}+(t-j)(G+\rho
P_{\infty})\big{)}&\text{ if }j<s\\\ \mathcal{L}(Q_{t}-sD)&\text{ if }j\geq
s\end{cases}\ .$

###### Proof.

The first claim is obvious. For the second claim observe that

$\displaystyle f-R$ $\displaystyle\in\mathcal{L}\big{(}G+\rho
P_{\infty}-D+E\big{)}\ .$

If $j<s$, then

$\displaystyle\Lambda_{s}(f-R)^{j}$
$\displaystyle\in\mathcal{L}\big{(}\lambda_{s}P_{\infty}+j(G+\rho
P_{\infty})-sE-j(D-E)\big{)}$
$\displaystyle\in\mathcal{L}\big{(}\lambda_{s}P_{\infty}+j(G+\rho
P_{\infty}-D)-(s-j)E\big{)}$ $\displaystyle\subset\mathcal{L}(W_{j})\ .$

If on the other hand $j\geq s$, then

$\displaystyle\Lambda_{s}(f-R)^{j}R^{t-j}$
$\displaystyle\in\mathcal{L}\big{(}\lambda_{s}P_{\infty}+t(G+\rho
P_{\infty})-sE-j(D-E)\big{)}$
$\displaystyle\subset\mathcal{L}\big{(}\lambda_{s}P_{\infty}+t(G+\rho
P_{\infty})-sE-s(D-E)\big{)}$ $\displaystyle=\mathcal{L}(Q_{t}-sD)\ .$

∎

We are now ready to present the key equations, which express an algebraic
relationship between the sent message $f$, an ${\boldsymbol{r}}$-interpolator
$R$ and a corresponding error locator $\Lambda_{s}$.

######  Theorem III.6 (The Key Equations).

For $t=1,\dots,\ell$

$\displaystyle\begin{split}\Lambda_{s}f^{t}-\sum_{j=0}^{\min\\{t,s-1\\}}\binom{t}{j}\Lambda_{s}(f-R)^{j}R^{t-j}\\\
\in\begin{cases}\\{0\\}&\text{ if }1\leq t\leq s-1\\\
\mathcal{L}(Q_{t}-sD)&\text{ if }s\leq t\leq\ell\end{cases}\ .\end{split}$ (1)

###### Proof.

We have that

$\displaystyle\Lambda_{s}f^{t}$
$\displaystyle=\Lambda_{s}\big{(}(f-R)+R\big{)}^{t}$
$\displaystyle=\sum_{j=0}^{t}\binom{t}{j}\Lambda_{s}(f-R)^{j}R^{t-j}\ ,$

which proves the claim when $1\leq t\leq s-1$, while if $s\leq t\leq\ell$,
then it follows from Lemma III.5 that

$\displaystyle\Lambda_{s}f^{t}-\sum_{j=0}^{\min\\{t,s-1\\}}\binom{t}{j}\Lambda_{s}(f-R)^{j}R^{t-j}$
$\displaystyle=\sum_{j=s}^{t}\binom{t}{j}\Lambda_{s}(f-R)^{j}R^{t-j}\in\mathcal{L}(Q_{t}-sD)$

∎

Since (1) is satisfied by the sent message, it makes sense to look for a
solution of (1) in order to address II.1. We explain how this can be done in
the next section.

##  IV Solving the Key Equations

In this section we show how to solve the following problem, which is motivated
by Theorem III.6 and II.1:

######  Problem IV.1.

Consider an instance of II.1. Fix $\lambda_{s}\geq 0$ and let $R$ be an
${\boldsymbol{r}}$-interpolator. Find $\phi_{t}\in\mathcal{L}(V_{t})$ for
$t=1,\dots,\ell$ and $\psi_{j}\in\mathcal{L}(W_{j})$ for $j=0,\dots,s-1$, not
all zero, such that

$\displaystyle\phi_{t}-\sum_{j=0}^{\min\\{t,s-1\\}}\binom{t}{j}\psi_{j}R^{t-j}$
$\displaystyle\quad\quad\in\begin{cases}\\{0\\}&1\leq t\leq s-1\\\
\mathcal{L}(Q_{t}-sD)&s\leq t\leq\ell\end{cases}\ .$

According to Theorem III.6, if $f$ is the sent message and there exists an
error locator $\Lambda_{s}$ with degree at most $\lambda_{s}$, then among the
possible solutions of IV.1 we would find

$\displaystyle\phi_{t}$ $\displaystyle=\Lambda_{s}f^{t}$
$\displaystyle\psi_{j}$ $\displaystyle=\Lambda_{s}(f-R)^{j}.$ (2)

Since $\deg\Lambda_{s}$ is roughly proportional to the number of errors (see
Lemma III.4), this means that a solution of IV.1 of the type (2) for the
smallest-possible $\lambda_{s}$ corresponds to a codeword that is among the
closest codewords to our received word.222Note that it might happen that for
two errors of similar weight, the error locator degree of the one of smaller
weight is larger or equal to the one of larger weight. For large enough
parameter $s$, this case becomes less likely for random errors. It might be
that IV.1 for small $\lambda_{s}$ has solutions that do not correspond to an
error locator/message pair. This might happen since IV.1 poses weaker
constraints on its solutions compared to (1). Previous works on (improved)
power decoding indicate that such “spurious” solutions only occur with very
small probability. Also our simulation results in Section VI indicate this. We
will get back to types of failure in the next section.

These arguments motivate our decoding strategy: Find a solution of IV.1 for
the smallest $\lambda_{s}$ for which the problem has a solution. If we are
successful, we get a solution of IV.1 of the form in (2) (or a scalar multiple
thereof), where $\Lambda_{s}$ and $f$ correspond to a close codeword to
${\boldsymbol{r}}$. We can extract $f$, and thus ${\boldsymbol{c}}$ from this
solution by division
$\tfrac{\phi_{1}}{\psi_{0}}=\tfrac{\Lambda_{s}f^{1}}{\Lambda_{s}}=f$.

Since computational complexity is not a primary concern here, we describe in
this section how to solve the problem for a fixed $\lambda_{s}$. The “smallest
solution” can then be found by iterating $\lambda_{s}$ from $0$ to the maximal
possible value for which we expect to be able to decode.

The entire decoding strategy is outlined in Algorithm 1.

Input: Received word ${\boldsymbol{r}}\in\mathbb{F}^{n}$ and positive integers
$s\leq\ell$

Output: $f\in\mathcal{L}(G)$ such that
${\boldsymbol{c}}=(f(P_{1}),\dots,f(P_{n}))$ is a codeword with corresponding
minimal error locator degree; or “decoding failure”.

1 for  _$\lambda_{s}=0,\dots,sn+g$_ do

2 if  _IV.1 with parameter $\lambda_{s}$ has a unique solution up to scalar
multiples _ then

3 $(\phi_{t})_{t=1,\dots,\ell},(\psi_{j})_{j=0,\dots,s-1}\leftarrow$ a
solution of IV.1

4 $\hat{f}\leftarrow\tfrac{\phi_{1}}{\psi_{0}}$

5 if  _$\hat{f}\in\mathcal{L}(G)$ and $\psi_{0}$ is a valid error locator for
the error $\hat{{\boldsymbol{e}}}={\boldsymbol{r}}-\hat{{\boldsymbol{c}}}$,
where $\hat{c}$ is the codeword corresponding to $\hat{f}$_ then

6 return _$\hat{f}$_

7 else

8 return _“decoding failure”_

9

10

return _“decoding failure”_

Algorithm 1 Decoder

In the remainder of this section we give an explicit construction of a matrix
$\boldsymbol{U}$ over $\mathbb{F}$ whose right kernel contains vectors from
which solutions to IV.1 can be obtained, i.e., how to implement Lines 1 and 1
of Algorithm 1.

For any divisor $A$ let
$\boldsymbol{\beta}_{A}\in\mathcal{L}(A)^{\mathfrak{l}(A)}$ denote a vector
whose entries form _any_ $\mathbb{F}$-basis of $\mathcal{L}(A)$ such that for
any two divisors $A$ and $B$ it holds that

$\mathcal{L}(A)=\mathcal{L}(B)\iff\boldsymbol{\beta}_{A}=\boldsymbol{\beta}_{B}\
.$

For any $h\in\mathcal{L}(A)$ let $h_{A}\in\mathbb{F}^{\mathfrak{l}(A)}$ be the
unique vector such that $\boldsymbol{\beta}_{A}\cdot h_{A}=h$, and for any
$p\in\mathcal{L}(B-A)$ let
$p_{B,A}\in\mathbb{F}^{\mathfrak{l}(B)\times\mathfrak{l}(A)}$ be the matrix
whose $i$-th column is $(pa_{i})_{B}$, where
$(a_{1},\dots,a_{\mathfrak{l}(A)})=\boldsymbol{\beta}_{A}$. It is easy to see
that

$\boldsymbol{\beta}_{B}p_{B,A}h_{A}=p\boldsymbol{\beta}_{A}\cdot
h_{A}=ph\in\mathcal{L}(B)\ .$

In order to construct $\boldsymbol{U}$ we will need a satisfactory received
word interpolator $R$ which can be for example be obtained as follows: If
$\boldsymbol{M}\in\mathbb{F}^{n\times\mathfrak{l}(H)}$, where $H=\rho
P_{\infty}+G$, is the matrix whose $i$-th column is

$[b_{i}(P_{1})\ \cdots\ b_{i}(P_{n})]^{T}\ ,$

where $(b_{1},\dots,b_{\mathfrak{l}(H)})=\boldsymbol{\beta}_{H}$, then
$R=\boldsymbol{\beta}_{H}\cdot\boldsymbol{z}$, where $\boldsymbol{z}$ is any
solution to the linear system $\boldsymbol{M}\boldsymbol{z}=\boldsymbol{r}$.

The next lemma together with the the subsequent Corollary IV.3 give an
explicit construction of $\boldsymbol{U}$ and show its relation to IV.1.

######  Lemma IV.2.

Let

$\boldsymbol{u}=(\boldsymbol{v}_{1}|\dots|\boldsymbol{v}_{\ell}|\boldsymbol{w}_{0}|\dots|\boldsymbol{w}_{s-1})\in\mathbb{F}^{\nu}\
,$

where

$\displaystyle\boldsymbol{v}_{r}$ $\displaystyle=(\phi_{r})_{V_{r}}\ ,$
$\displaystyle\boldsymbol{w}_{j}$
$\displaystyle=(\psi_{j})_{W_{j}}\quad\text{and}$ $\displaystyle\nu$
$\displaystyle=\sum_{r}\mathfrak{l}(V_{r})+\sum_{j}\mathfrak{l}(W_{j})$

for $r=1,\dots,\ell$ and $j=0,\dots,s-1$, where the $\phi_{r},\psi_{j}$ are a
solution to IV.1. If
$\boldsymbol{U}_{t}\in\mathbb{F}^{\mathfrak{l}(Q_{t})\times\nu}$ is the matrix
defined by

$\boldsymbol{U}_{t}=\begin{bmatrix}\boldsymbol{V}_{t,1}|\cdots|\boldsymbol{V}_{t,\ell}|\boldsymbol{W}_{t,0}|\cdots|\boldsymbol{W}_{t,s-1}\end{bmatrix}\
,$

where

$\displaystyle\boldsymbol{V}_{t,r}$
$\displaystyle=\begin{cases}1_{Q_{t},V_{r}}&\text{ if }t=r\\\
0_{Q_{t},V_{r}}&\text{ if }t\neq r\\\ \end{cases}\quad\text{and}$
$\displaystyle\boldsymbol{W}_{t,j}$
$\displaystyle=\begin{cases}-\binom{t}{j}(R^{t-j})_{Q_{t},W_{j}}&\text{ if
}j\leq t\\\ 0_{Q_{t},W_{j}}&\text{ if }j>t\\\ \end{cases}\ ,$

for $t=1,\dots,\ell$, then

$\boldsymbol{\beta}_{Q_{t}}\boldsymbol{U}_{t}\boldsymbol{u}\in\begin{cases}\\{0\\}&\text{
if }1\leq t\leq s-1\\\ \mathcal{L}(Q_{t}-sD)&\text{ if }s\leq
t\leq\ell\end{cases}\ .$

###### Proof.

Simply observe that

$\displaystyle\boldsymbol{\beta}_{Q_{t}}\boldsymbol{U}_{t}\boldsymbol{u}$
$\displaystyle=\boldsymbol{\beta}_{Q_{t}}\sum_{r=1}^{\ell}\boldsymbol{V}_{t,r}\boldsymbol{v}_{r}+\boldsymbol{\beta}_{Q_{t}}\sum_{j=0}^{s-1}\boldsymbol{W}_{t,j}\boldsymbol{w}_{j}$
$\displaystyle=\phi_{t}-\sum_{j=0}^{\min\\{t,s-1\\}}\binom{t}{j}\psi_{j}R^{t-j}\
.$

∎

######  Corollary IV.3 (of Lemma IV.2).

If

$\boldsymbol{U}=\left[\begin{array}[]{cc}\boldsymbol{K}_{1}\boldsymbol{U}_{1}\\\
\hline\cr\vdots\\\
\hline\cr\boldsymbol{K}_{\ell}\boldsymbol{U}_{\ell}\end{array}\right]\ ,$

where

$\boldsymbol{K}_{t}=\begin{cases}1_{Q_{t},Q_{t}}&\text{ if }1\leq t\leq s-1\\\
\text{left kernel matrix of }1_{Q_{t},Q_{t}-sD}&\text{ if }s\leq
t\leq\ell\end{cases}\ $

and $\boldsymbol{u}$ is as in Lemma IV.2, then

$\displaystyle\boldsymbol{U}\boldsymbol{u}=\boldsymbol{0}\ .$ (3)

In the next section we describe the conditions under which decoding is
expected to fail.

##  V Decoding Radius

We turn to investigating the decoding performance of the proposed decoder.
Since the decoder returns at most one codeword while attempting to decode
beyond half the minimum distance, it must obviously fail for certain received
words. Theoretically bounding this probability seems a difficult problem even
for the simplest case of Reed–Solomon codes, and only partial results are
known [1, 5]. We follow the approach of several previous papers on power
decoding: we derive theoretically a number of errors at which our decoder is
always guaranteed to fail, and we _define_ this to be the “decoding radius”;
this name will be supported by simulations in Section VI which indicates that
the decoding algorithm succeeds with high probability whenever fewer, random
errors occur.

It can be seen from Algorithm 1 that our decoder declares a decoding failure
if and only if there is no value of $\lambda_{s}$ for which the matrix
$\boldsymbol{U}$ in Corollary IV.3 has a right kernel of
$\mathbb{F}_{q}$-dimension exactly 1.

On the other hand, when we do _not_ declare decoding failure, the returned
message candidate $\hat{f}\in\mathcal{L}(G)$ must correspond to the closest
codeword to the received word, as in the discussion in the preceding section.

This notion of decoding failure motivates a natural bound for when we – simply
by linear algebraic arguments – would know a priori that the decoder will fail
in the generic case when the smallest error locator corresponding to a
codeword has degree $\lambda_{s}=s|\mathcal{E}|+g$ (and not smaller), where
$|\mathcal{E}|$ is the distance of the codeword to the received word:

######  Definition V.1.

For a given code $\mathcal{C}$, the _decoding radius_ of power decoding,
denoted $\tau_{\mathsf{max}}$, is the greatest value of $\tau$ such that for
$\lambda_{s}=s\tau+g$, if $\epsilon,\nu$ are the row resp. column dimensions
of the matrix $\boldsymbol{U}$ of Lemma IV.2, then $\nu>\epsilon+1$.

Since the proposed decoder will typically return no solutions (and therefore
fail) when $|\mathcal{E}|>\tau$, the decoding radius tells us how many errors
we should at most expect to correct.333It might succeed in rare cases in which
the error locator has surprisingly low degree. Our simulation results suggest
that this happens with very small probability. We reiterate that, at this
point we have given no indications that the decoder should usually succeed up
to $\tau_{\mathsf{max}}$ errors, but we address this by simulations in Section
VI.

For a given code $\mathcal{C}$, one may compute $\tau_{\mathsf{max}}$ exactly,
but since the equations involve dimensions of sequences of Riemann–Roch
spaces, it seems impossible to give a general precise closed form. However, we
may lower-bound these numbers using common function field tools:

######  Lemma V.2.

For $t=1,\dots,\ell$ and $j=0,\dots,s-1$ it holds that

$\displaystyle\mathfrak{l}(V_{t})$ $\displaystyle\geq\lambda_{s}+t-g+1\ ,$
$\displaystyle\mathfrak{l}(W_{j})$
$\displaystyle\geq\lambda_{s}+j(\gamma+\rho-n)-g+1\ ,$
$\displaystyle\mathfrak{l}(Q_{t}-sD)$
$\displaystyle\geq\lambda_{s}+t(\gamma+\rho)-sn-g+1\ ,$

and if $\gamma\geq 2g-2$, then

$\mathfrak{l}(Q_{t})=\lambda_{s}+t(\gamma+\rho)-g+1\ .$

###### Proof.

The lower bounds on $\mathfrak{l}(V_{t})$, $\mathfrak{l}(W_{j})$ and
$\mathfrak{l}(Q_{t}-sD)$ are given by the Riemann-Roch theorem, while the
exact value of $\mathfrak{l}(Q_{t})$ is given by [24, Thm 1.5.17] since
$\gamma\geq 2g-2$. ∎

The next lemma gives bounds on the dimensions of $\boldsymbol{U}$.

######  Lemma V.3.

If $\gamma\geq 2g-2$, then $\boldsymbol{U}\in\mathbb{F}^{\epsilon\times\nu}$,
where

$\displaystyle\epsilon$
$\displaystyle\leq(s-1)(\lambda_{s}-g+1)+\tfrac{s(s-1)}{2}(\gamma+\rho)+(\ell-s+1)sn\
,$ $\displaystyle\nu$
$\displaystyle\geq(\ell+s)(\lambda_{s}-g+1)+\tfrac{\ell(\ell+1)}{2}\gamma+\tfrac{s(s-1)}{2}(\gamma+\rho-n)\
.$

###### Proof.

From Lemma V.2 it follows that

$\displaystyle\nu$
$\displaystyle=\sum_{t=1}^{\ell}\mathfrak{l}(V_{t})+\sum_{j=0}^{s-1}\mathfrak{l}(W_{j})$
$\displaystyle\geq\sum_{t=1}^{\ell}(\lambda_{s}+t\gamma-g+1)+\sum_{j=0}^{s-1}(\lambda_{s}+j(\gamma+\rho)-jn-g+1)$
$\displaystyle=\ell(\lambda_{s}-g+1)+\tfrac{\ell(\ell+1)}{2}\gamma$
$\displaystyle+s(\lambda_{s}-g+1)+\tfrac{(s-1)s}{2}(\gamma+\rho-n)\ .$

Similarly, observing that $\boldsymbol{K}_{t}$ must have exactly

$\mathfrak{l}(Q_{t})-\mathfrak{l}(Q_{t}-sD)$

rows for $t=s-1,\dots,\ell$, it holds that

$\displaystyle\epsilon$
$\displaystyle=\sum_{t=1}^{s-1}\mathfrak{l}(Q_{t})+\sum_{t=s}^{\ell}\big{(}\mathfrak{l}(Q_{t})-\mathfrak{l}(Q_{t}-sD)\big{)}$
$\displaystyle\leq\sum_{t=1}^{s-1}\big{(}\lambda_{s}+t(\gamma+\rho)-g+1\big{)}+\sum_{t=s}^{\ell}sn$
$\displaystyle=(s-1)(\lambda_{s}-g+1)+\frac{(s-1)s}{2}(\gamma+\rho)+(l-s+1)sn\
.$

∎

With the aid of Lemma V.3 we can deduce a lower bound on values of
$\lambda_{s}$ for which decoding must fail.

######  Lemma V.4.

It holds that (3) has at least two linearly independent solutions (and
decoding fails) if

$\lambda_{s}>s\tfrac{2\ell-s+1}{2(\ell+1)}n-\tfrac{\ell}{s}\gamma+\tfrac{\ell}{\ell+1}+g\
.$ (4)

###### Proof.

Keeping the notation from Lemma V.3 we have that (3) has at least two linearly
independent solutions if

$\nu>\epsilon+1\ ,$

which is equivalent to (4). ∎

######  Corollary V.5.

The decoding radius as defined in Definition V.1 is given by

$\displaystyle\tau_{\mathsf{max}}(\ell,s):=\left\lfloor\tfrac{2\ell-s+1}{2(\ell+1)}n-\tfrac{\ell}{2s}\gamma+\tfrac{\ell}{s(\ell+1)}\right\rfloor\
.$

######  Remark V.6.

The decoding radius $\tau_{\mathsf{max}}(\ell,s)$ coincides with the one in
[20] in the special case of Reed–Solomon, and the one in [21] for one-point
Hermitian codes.

######  Remark V.7.

Although this definition of failure is convenient for theoretical analysis it
means our decoder must try every sensible (sufficiently small) value of $\tau$
before declaring failure, which is computationally inefficient. In this paper,
our primary concern is not computational complexity, but we will remark that
in practice, one only needs to use the largest sensible value of $\tau$: we
then pick _any_ non-zero vector in the right kernel of $\boldsymbol{U}$ and
check if the corresponding message $f$ is in $\mathcal{L}(G)$. This is how the
simulations in Section VI have been conducted.

###  V-A Parameter Choice & Asymptotic Behavior

The following theorem shows that the decoding radius achieves the Johnson
radius asymptotically for large parameters $\ell$ and $s$. It also implies a
good practical choice of the parameters $\ell$ and $s$.

######  Theorem V.8.

Define the sequence $(\ell_{i},s_{i})_{i\in\mathbb{Z}_{>0}}$ as $\ell_{i}=i$
and $s_{i}=\left\lfloor\sqrt{\tfrac{\gamma}{n}}i\right\rfloor+1$. Then, for
$i\to\infty$, we have

$\displaystyle\tau_{\mathsf{max}}(\ell_{i},s_{i})=n\left(1-\sqrt{\tfrac{\gamma}{n}}-O(\tfrac{1}{i})\right).$

###### Proof.

This is a special case of [21, Theorem 5]. ∎

##  VI Numerical Results

In this section we present Monte-Carlo simulation results of the proposed
decoder for a few different AG codes in order to experimentally verify our
hypothesis that Algorithm 1 corrects up to $\tau_{\mathsf{max}}(\ell,s)$
errors with high probability and fails with high probability above.

To be precise, for a given code and radius $\tau$, we test the following:

  * •

Draw a random message $f$ and encode it into ${\boldsymbol{c}}$

  * •

Draw an error ${\boldsymbol{e}}$ uniformly at random from the set of vectors
in $\mathbb{F}^{n}$ of Hamming weight $\tau$

  * •

Decode ${\boldsymbol{r}}={\boldsymbol{c}}+{\boldsymbol{e}}$ using Algorithm 1

  * •

If the algorithm does not return a decoding failure and the returned message
polynomial is exactly $f$, then count it as success.

Note that the simulation might be unsuccessful also in cases in which
Algorithm 1 returns a valid close codeword, but not the transmitted one. This
is called a miscorrection, and typically occurs only rarely (see, e.g., [25]
for RS codes).

The above described simulation was tested on a few one- and two-point codes
over two well known function fields:

  * •

The _Hermitian_ function field is defined by the equation

$H_{q}:y^{q}+y=x^{q+1}$

over $\mathbb{F}_{q^{2}}$, i.e. the finite field of cardinality $q^{2}$. It
has genus $\frac{1}{2}q(q-1)$ and $q^{3}+1$ rational places.

  * •

The _Suzuki_ function field is defined by the equation

$S_{q}:y^{q}+y=x^{q_{0}}(x^{q}+x)$

over $\mathbb{F}_{q^{4}}$, where $q=2q_{0}^{2}>2$. It has genus $q_{0}(q-1)$
and $q^{2}+1$ rational places.

  * •

The function from [26, Lemma 3.2] field is defined by

$T_{q}:y^{q}+y=\frac{x^{q}}{x^{q-1}+1}$

over $\mathbb{F}_{q^{2}}$. It has genus $(q-1)^{2}$ and $q^{3}-q^{2}+2q$
rational places.

Table I contains simulation results for various function field family, code,
and decoder parameters. All results confirm our hypothesis.

Table I: Simulation results

Curve | $|\mathbb{F}|$ | $\gamma$ | $n$ | $k$ | $d^{*}$ | $\ell$ | $s$ | $\tau$ | OFR | $N\geq$
---|---|---|---|---|---|---|---|---|---|---
$H_{4}$ | $4^{2}$ | $15$ | $64$ | $10$ | $49$ | $4$ | $2$ | $29^{+}$ | $0.00$ | $10^{2}$
|  |  |  |  |  |  |  | $30$ | $1.00$ | $10^{2}$
$H_{4}$ | $4^{2}$ | $10+5$ | $63$ | $10$ | $48$ | $4$ | $2$ | $28^{+}$ | $0.00$ | $10^{2}$
|  |  |  |  |  |  |  | $29$ | $1.00$ | $10^{2}$
$H_{5}$ | $5^{2}$ | $55$ | $125$ | $46$ | $70$ | $5$ | $2$ | $36^{+}$ | $0.00$ | $10^{2}$
|  |  |  |  |  |  |  | $37$ | $1.00$ | $10^{2}$
$H_{5}$ | $5^{2}$ | $30+25$ | $124$ | $46$ | $69$ | $3$ | $2$ | $35^{+}$ | $0.00$ | $10^{2}$
|  |  |  |  |  |  |  | $36$ | $0.94$ | $10^{2}$
$S_{1}$ | $2^{4}$ | $12$ | $24$ | $12$ | $12$ | $2$ | $2$ | $5^{+}$ | $0.00$ | $10^{4}$
|  |  |  |  |  |  |  | $6$ | $>0.99$ | $10^{4}$
$S_{1}$ | $2^{4}$ | $6+6$ | $23$ | $12$ | $11$ | $2$ | $2$ | $5^{+}$ | $0.00$ | $10^{4}$
|  |  |  |  |  |  |  | $6$ | $1.00$ | $10^{4}$
$S_{1}$ | $2^{4}$ | $4$ | $24$ | $4$ | $20$ | $6$ | $2$ | $12^{+}$ | $<6.98\cdot 10^{-4}$ | $10^{3}$
|  |  |  |  |  |  |  | $13$ | $1.00$ | $10^{3}$
$S_{1}$ | $2^{4}$ | $2+2$ | $23$ | $4$ | $19$ | $6$ | $2$ | $11^{+}$ | $0.00$ | $10^{3}$
|  |  |  |  |  |  |  | $12$ | $>0.99$ | $10^{3}$
$T_{4}$ | $4^{2}$ | $15$ | $55$ | $7$ | $40$ | $4$ | $2$ | $23^{+}$ | $<1.74\cdot 10^{-3}$ | $10^{3}$
|  |  |  |  |  |  |  | $24$ | $1.00$ | $10^{3}$

  * •

Code parameters $\gamma,n,k,d^{*}$. $\gamma=a+b$ means that $G=aP+bP^{\prime}$
for some rational places $P$ and $P^{\prime}$. Decoder parameters $\ell,s$.
Number of errors $\tau$, + means that $\tau=\tau_{\max}$. Observed failure
rate OFR. Each simulation was repeated at least $N$ times.

## References

  * [1] G. Schmidt, V. Sidorenko, and M. Bossert, “Decoding Reed-Solomon Codes Beyond Half the Minimum Distance Using Shift-Register Synthesis,” in _IEEE International Symposium on Information Theory_ , 2006, pp. 459–463.
  * [2] ——, “Syndrome Decoding of Reed-Solomon Codes Beyond Half the Minimum Distance Based on Shift-Register Synthesis,” _IEEE Transactions on Information Theory_ , vol. 56, no. 10, pp. 5245–5252, 2010.
  * [3] M. Sudan, “Decoding of Reed–Solomon Codes beyond the Error-Correction Bound,” _Journal of Complexity_ , vol. 13, no. 1, pp. 180–193, 1997.
  * [4] J. S. R. Nielsen, “Power Decoding of Reed–Solomon Codes Revisited,” in _International Castle Meeting on Coding Theory and Applications_ , Sep. 2014. [Online]. Available: http://jsrn.dk/publications.html
  * [5] J. Rosenkilde, “Power Decoding of Reed–Solomon Up to the Johnson Radius,” _Advances in Mathematics of Communications_ , vol. 12, no. 1, pp. 81–106, Feb. 2018.
  * [6] V. Guruswami and M. Sudan, “Improved Decoding of Reed–Solomon and Algebraic-Geometric Codes,” in _IEEE Annual Symposium on Foundations of Computer Science_ , 1998, pp. 28–37.
  * [7] V. D. Goppa, “Codes on algebraic curves,” _Dokl. Akad. Nauk SSSR_ , vol. 259, no. 6, pp. 1289–1290, 1981.
  * [8] ——, “Algebraico-Geometric Codes,” _Izvestiya: Mathematics_ , vol. 21, no. 1, pp. 75–91, 1983.
  * [9] M. A. Tsfasman, S. G. Vladut, and T. Zink, “Modular curves, Shimura curves, and Goppa codes, better than Varshamov-Gilbert bound,” _Mathematische Nachrichten_ , vol. 109, no. 1, pp. 21–28, 1982.
  * [10] V. Sidorenko and G. Schmidt, “A Linear Algebraic Approach to Multisequence Shift-Register Synthesis,” _Problems of Information Transmission_ , vol. 47, no. 2, pp. 149–165, 2011.
  * [11] J. Rosenkilde and A. Storjohann, “Algorithms for simultaneous Hermite–Padé approximations,” _Journal of Symbolic Computation_ , vol. In press, Oct. 2019. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0747717119301300
  * [12] M. Chowdhury, C.-P. Jeannerod, V. Neiger, E. Schost, and G. Villard, “Faster Algorithms for Multivariate Interpolation With Multiplicities and Simultaneous Polynomial Approximations,” _IEEE Transactions on Information Theory_ , vol. 61, no. 5, pp. 2370–2387, May 2015.
  * [13] V. Neiger, J. Rosenkilde, and E. Schost, “Fast Computation of the Roots of Polynomials Over the Ring of Power Series,” in _International Symposium on Symbolic and Algebraic Computation_ , Jul. 2017. [Online]. Available: https://hal.inria.fr/hal-01457954/document
  * [14] J. Nielsen and P. Beelen, “Sub-Quadratic Decoding of One-Point Hermitian Codes,” _IEEE Transactions on Information Theory_ , vol. 61, no. 6, pp. 3225–3240, Jun. 2015.
  * [15] A. Ahmed, R. Koetter, and N. R. Shanbhag, “VLSI Architectures for Soft-Decision Decoding of Reed-Solomon Codes,” _IEEE Transactions on Information Theory_ , vol. 57, no. 2, pp. 648–667, Feb. 2011.
  * [16] X.-W. Wu and P. H. Siegel, “Efficient Root-Finding Algorithm With Application to List Decoding of Algebraic-Geometric Codes,” _IEEE Transactions on Information Theory_ , vol. 47, no. 6, pp. 2579–2587, 2001.
  * [17] H. Cohn and N. Heninger, “Approximate Common Divisors via Lattices,” _The Open Book Series_ , vol. 1, no. 1, pp. 271–293, 2013.
  * [18] S. Puchinger and J. Rosenkilde né Nielsen, “Decoding of Interleaved Reed-Solomon Codes Using Improved Power Decoding,” _IEEE International Symposium on Information Theory_ , 2017.
  * [19] A. Couvreur and I. Panaccione, “Power error locating pairs,” _Designs, Codes and Cryptography_ , vol. 88, no. 8, pp. 1561–1593, 2020.
  * [20] J. Rosenkilde, “Power decoding Reed-Solomon codes up to the Johnson radius,” _Advances in Mathematics of Communications_ , vol. 12, no. 1, p. 81, 2018.
  * [21] S. Puchinger, J. Rosenkilde, and I. Bouw, “Improved Power Decoding of Interleaved One-Point Hermitian Codes,” _Designs, Codes and Cryptography_ , vol. 87, no. 2-3, pp. 589–607, 2019.
  * [22] S. C. Porter, B.-Z. Shen, and R. Pellikaan, “Decoding geometric goppa codes using an extra place,” _IEEE transactions on information theory_ , vol. 38, no. 6, pp. 1663–1676, 1992.
  * [23] S. Gao, “A New Algorithm for Decoding Reed-Solomon Codes,” in _Communications, Information and Network Security_ , ser. The Springer International Series in Engineering and Computer Science. Springer, Jan. 2003, no. 712, pp. 55–68.
  * [24] H. Stichtenoth, _Algebraic Function Fields and Codes_ , 2nd ed. Springer, 2009.
  * [25] G. Schmidt, V. R. Sidorenko, and M. Bossert, “Collaborative Decoding of Interleaved Reed–Solomon Codes and Concatenated Code Designs,” _IEEE Transactions on Information Theory_ , vol. 55, no. 7, pp. 2991–3012, 2009.
  * [26] A. Garcia and H. Stichtenoth, “On the asymptotic behaviour of some towers of function fields over finite fields,” _Journal of number theory_ , vol. 61, no. 2, pp. 248–273, 1996.