Finding spin-glass ground states using quantum walks

TL;DR

This paper explores the use of continuous-time quantum walks to find spin glass ground states, revealing practical parameter settings and potential for real-world quantum optimization applications.

Contribution

It demonstrates that quantum walks can effectively find spin glass ground states without precise parameter tuning, unlike in search problems.

Findings

Quantum walks outperform classical methods in finding spin glass ground states.

Parameter tuning is less critical for spin glass problems compared to search problems.

Quantum walks achieve better than square-root scaling in this context.

Abstract

Quantum computation using continuous-time evolution under a natural hardware Hamiltonian is a promising near- and mid-term direction toward powerful quantum computing hardware. We investigate the performance of continuous-time quantum walks as a tool for finding spin glass ground states, a problem that serves as a useful model for realistic optimization problems. By performing detailed numerics, we uncover significant ways in which solving spin glass problems differs from applying quantum walks to the search problem. Importantly, unlike for the search problem, parameters such as the hopping rate of the quantum walk do not need to be set precisely for the spin glass ground state problem. Heuristic values of the hopping rate determined from the energy scales in the problem Hamiltonian are sufficient for obtaining a better than square-root scaling. This makes it practical to use quantum walks for solving such problems, and opens the door for a range of applications on suitable quantum hardware.

Equations

$${\hat{H}}_I=-\sum _{(j\ne k)=0}^{n-1}{J}_{jk}{\hat{Z}}_j{\hat{Z}}_k-\sum _{j=0}^{n-1}{h}_j{\hat{Z}}_j,$$

$${\hat{Z}}_j=\left(\underset{r=0}{\overset{j-1}{⨂}}{\hat{1}}_2\right)\otimes \hat{Z}\otimes \left(\underset{r=j+1}{\overset{n-1}{⨂}}{\hat{1}}_2\right),$$

$${\hat{H}}_T=-\Gamma \sum _{j=0}^{n-1}{\hat{X}}_j,$$

$${\hat{H}}_{\mathrm{TI}}(t)=A(t){\hat{H}}_T+B(t){\hat{H}}_I,$$

$$⟨j|{\hat{H}}_G|k⟩=-\gamma L_{jk},$$

$$|\psi (t)⟩=\exp \{-\text{rmi}{\hat{H}}_{G}t\}|\psi (0)⟩,$$

$$|\psi (0)⟩=N^{-1/2}\sum _{j=0}^{N-1}|j⟩,$$

$$\hat{H}(\gamma )\equiv {\hat{H}}_G+{\hat{H}}_P,$$

$$\bar{P}(t,\Delta t)\equiv \frac{1}{\Delta t}\underset{t}{\overset{t+\Delta t}{\int }}\mathrm{d}{t}_{f}P(t_f).$$

$$P_{\infty }\equiv \bar{P}(0,\infty )\equiv \underset{\Delta t\to \infty }{\lim }\bar{P}(0,\Delta t),$$

$${\hat{H}}_K$$

$${\hat{H}}_h=\gamma \left[n1-\sum _{j=0}^{n-1}{\hat{X}}_j\right].$$

$${\hat{H}}_{\mathrm{S}}=-|m⟩⟨m|.$$

$$2{\gamma }_{\mathrm{opt}}^{(h)}=\frac{1}{N}\sum _{r=1}^n\left(\genfrac{}{}{0ex}{}{n}{r}\right)\frac{1}{r},$$

$$H_{\mathrm{SK}}=-\frac{1}{2}\sum _{(j\ne k)=0}^{n-1}{J}_{jk}{S}_j{S}_k$$

$${\hat{H}}_{\mathrm{SK}}=-\frac{1}{2}\sum _{(j\ne k)=0}^{n-1}{J}_{jk}{\hat{Z}}_j{\hat{Z}}_k-\sum _{j=0}^{n-1}{h}_j{\hat{Z}}_j,$$

$${\hat{H}}_{\mathrm{REM}}$$

$$\hat{H}(\gamma )=\sum _{a=0}^{N-1}{E}_a(\gamma )|E_a(\gamma )⟩⟨E_a(\gamma )|,$$

$$P(t)$$

$$P_{\infty }=\sum _{a=0}^{N-1}{\left|⟨E_0^{(P)}|E_a(\gamma )⟩\right|}^2{\left|⟨E_a(\gamma )|\psi (0)⟩\right|}^2.$$

$${\gamma }_{\mathrm{heur}}^{(h)}\equiv \frac{1}{2n}⟨E_{N-1}^{(P)}-E_0^{(P)}⟩.$$

$$\left\langle E_{N-1}^{(P)}-E_{0}^{(P)}\right\rangle\simeq-(2^{\frac{3}{2}}\sigma_{P}^{(\mathrm{energy})})\,\mathrm{erf}^{-1}\Big{(}\frac{1}{N}-1\Big{)},$$

$$\langle P_{\infty}\rangle=\begin{array}[]{ll}\tilde{O}(N^{-0.402\pm 0.001})&\mbox{with }\gamma^{(h)}_{\mathrm{opt}}\\ \tilde{O}(N^{-0.417\pm 0.002})&\mbox{with }\gamma^{(h)}_{\mathrm{heur}}\end{array},$$

$${\tau }_{\mathrm{mix}}^{(ϵ)}$$

$$⟨{\tau }_{\mathrm{mix}}^{(0.05)}{\omega }_{\mathrm{SK}}⟩$$

$$\displaystyle\Big{\langle}\bar{P}(t_{\mathrm{short}},\Delta t_{\mathrm{short}})\Big{\rangle}$$

$$⟨\hat{H}(\gamma ){⟩}_{\psi (t)}$$

$$⟨{\hat{H}}_G{⟩}_{\psi (t)}-⟨{\hat{H}}_G{⟩}_{\psi (0)}$$

$$\frac{1}{t_f}\underset{t=0}{\overset{t_f}{\int }}\mathrm{d}t⟨{\hat{H}}_P{⟩}_{\psi (t)}$$

Full text

Keywords: quantum computing, quantum walks, spin glasses

1. 1 Introduction

2. 2 Computing with quantum walks

3. 2.1 Continuous-time quantum walks

4. 2.2 Computing using a quantum walk

5. 2.3 Graph choice for quantum walk computing

6. 2.4 Solving the search problem using quantum walks

7. 3 Spin glass problem definitions

8. 3.1 Sherrington-Kirkpatrick spin glass

9. 3.2 Random energy model

10. 4 Numerical methods

11. 5 Quantum walks with spin glasses

12. 5.1 Setting the hopping rate

13. 5.2 Success probability

14. 5.3 Mixing times

15. 6 Computational mechanisms

16. 6.1 Role of correlations in SK

17. 6.2 Energy conservation dynamics

18. 7 Summary and outlook

1 Introduction

Optimization problems need to be solved in a broad range of areas, such as scheduling, route planning, supply chains, finance. This is often computationally intensive, so the prospect of quantum enhanced solution methods is an important research direction for practical quantum computing. One way to tackle optimization in a quantum setting is to use a device which realises an Ising Hamiltonian with a transverse field. Computing using the Ising Hamiltonian works as follows: The optimization problem is encoded into the Ising Hamiltonian $\hat{H}_{I}$

[TABLE]

on $n$ qubits, such that the solution corresponds to the ground state of $\hat{H}_{I}$. In our notation, the operator $\hat{Z}_{j}$ on the full Hilbert space applies the single qubit Pauli-$Z$ operator $\hat{Z}$ to the $j$th qubit,

[TABLE]

where $\hat{\mathbb{1}}_{2}$ is the identity operator on a single qubit. The (real) values of the coupling strengths $J_{jk}$ and fields $h_{j}$ define the optimization problem, and efficient methods are known for expressing optimization problems in terms of these coupling and field strengths [e.g., Choi, 2010]. The transverse field term $\hat{H}_{T}$

[TABLE]

drives transitions between states, where $\Gamma$ is a real-valued transverse field strength, and $\hat{X}_{j}$ is the operator on the full Hilbert space that applies the single qubit Pauli-$X$ operator to the $j$th qubit, defined by analogy with $\hat{Z}_{j}$ in (2). The qubits are initialised in the ground state of $\hat{H}_{T}$, this is easy to do by applying a strong transverse field to align all the qubits in the state $\left|+\right\rangle=2^{-1/2}(\left|0\right\rangle+\left|1\right\rangle)$. Then, the computation is carried out by applying the full transverse Ising Hamiltonian

[TABLE]

where $t$ is time and $A(t)$, $B(t)$ are real-valued control functions. To obtain a candidate solution to the optimization problem, the qubit register is measured after a time $t_{f}$. For some problems, sampling from the distribution of low energy states provides the required solution – this can be done by repeating the computation, which will in general not produce the lowest energy state with certainty.

The Ising Hamiltonian is a natural choice for encoding problems for two reasons. First, it is proven to be universal for classical problems [De las Cuevas and Cubitt, 2016]. There are efficient methods for mapping NP-hard optimization problems to the Ising model [Lucas, 2014, Choi, 2010], providing a practical route to quantum algorithms. Since many optimization problems are NP-hard, an exponential speed up is not expected, but even modest polynomial improvements are useful for practical applications. There is increasing interest in how to obtain polynomial advantages through quantum algorithms [Moylett et al., 2017, Montanaro, 2018, Ambainis et al., 2019]. Interesting results have been presented for a wide range of applications, such as mathematics [Bian et al., 2013, Li et al., 2017], computer science [Chancellor et al., 2016], computational biology [Perdomo-Ortiz et al., 2012], finance [Marzec, 2016], and aerospace [Coxson et al., 2014]. Second, the Ising Hamiltonian can be implemented in a range of different physical systems. The quantum Ising Hamiltonian is the basic interaction Hamiltonian in the D-Wave Systems Inc. programmable superconducting devices [D-Wave, 1999–, Boixo et al., 2013, Johnson et al., 2011]. Implementations in other promising architectures include Rydberg systems [Bernien et al., 2017] and trapped ions [Kim et al., 2011]. The Ising Hamiltonian is also the basic tool for specialised optimization hardware, such as coherent Ising machines [Inagaki et al., 2016, McMahon et al., 2016]. Optimization using the Ising Hamiltonian can be implemented in digital quantum architectures by using the quantum approximate optimization algorithm (QAOA) [Farhi et al., 2014a, b, Marsh and Wang, 2019] or quantum alternating operator ansatz [Hadfield et al., 2019]. Studies by Zhou et al. [2018] show how to exploit non-adiabatic effects in QAOA on early quantum hardware.

There are several known methods for driving the quantum system from its initial state into the ground state of a Hamiltonian defining the problem to be solved. These methods correspond to different choices for the control functions $A(t)$ and $B(t)$ in (4). Adiabatic quantum computing [Kadowaki and Nishimori, 1998, Farhi et al., 2000, 2001] keeps the quantum system in the ground state while the initial Hamiltonian is slowly changed into the problem Hamiltonian. Quantum annealing [Finnila et al., 1994] takes advantage of open quantum systems effects to cool the system towards the ground state. Continuous-time quantum walks evolve the system under a time-independent Hamiltonian for a suitable time before measurement of the final state. Computation by continuous-time quantum walk and adiabatic quantum computing are end points of a family of continuous-time protocols that use the same Hamiltonian terms but are applied with different time dependent modulation [Morley et al., 2019]. In this work, we focus on computation by quantum walk using time-independent transverse Ising Hamiltonians.

Quantum walks can solve the search problem [Childs and Goldstone, 2004], achieving the same quadratic $O(N^{1/2})$ quantum speed up as is obtained by Grover’s algorithm [Grover, 1996]. We describe the search problem further in Subsection 2.4. For particular graphs, quantum walks can solve problems exponentially faster [e.g., Childs et al., 2003], and quantum walks are now widely used as subroutines in more complex quantum algorithms. However, in the continuous-time setting, the application of quantum walks to optimization problems has not been studied in detail. There is increasing interest in quenches [Amin et al., 2018] or pauses [Marshall et al., 2019, Passarelli et al., 2019] in quantum annealing, which effectively run an open-system version of a quantum walk during part of the computation. Thermal relaxation effects dominate in the regime currently accessible by flux qubit quantum annealers, which is the focus of these works. An algorithm which is essentially a quantum walk on a spin glass, although presented using different terminology, has been analysed by Hastings [2019]. Along with the same energy conservation arguments we describe in section 6.2, Hastings’ findings suggest that quantum walks on spin glasses will be interesting to explore. Given that quantum walks provide a better performance for searching than adiabatic quantum computing, especially when limited coherence time and other practical factors, such as precision of control settings, are considered [Morley et al., 2019], it is important to understand how they perform for a wider range of problems.

In this work, we tackle the question of if, and how, a quantum walk can be useful for practical quantum optimization. We present a detailed numerical investigation of continuous-time quantum walks applied to solving combinatorial optimization problems, using the Sherrington-Kirkpatrick spin glass ground state problem as a prototypical example. Finding the ground state of a frustrated Sherrington-Kirkpatrick spin glass [Sherrington and Kirkpatrick, 1975] is known to be not only NP-hard, but also uniformly-hard, as suggested by its finite-temperature spin glass transition. Without a finite temperature spin glass transition, a problem cannot be uniformly hard, since the lack of a transition implies that typical cases will be easy for the Monte Carlo family of algorithms, as discussed in [Katzgraber et al., 2014]. As has been shown for a random problem type used in early benchmarks of quantum annealing hardware [Katzgraber et al., 2014], uniform hardness is crucial: without this property, randomly generated instances of NP-hard problems are not necessarily hard to solve [Beier and Vöcking, 2004, Krivelevich and Vilenchik, 2006, Lucas, 2014].

We use a random energy model [Derrida, 1980] for comparisons, to draw out the effects of the correlations between energy difference and Hamming distance in the spin glass. A problem with perfect correlations is easy to solve, like finding the ground state of a spin system with only local fields, no couplings. A completely random problem, such as finding the ground state of a random energy model instance, has no correlation to exploit and so is very hard to solve, essentially requiring random guessing. However, a completely random model is fully characterised by average values of its properties, and finding exact ground states of specific instances is typically not interesting. Intermediate problems with some correlations are both hard and interesting, with complex behaviour and phase diagrams, like spin models with frustration and spin glass phases. Real optimization problems typically have correlations; they are often hard to solve but also produce interesting solutions. The inherent complexity of a problem comes from the structures of the problem and its correlations, not the structure of the solution itself. One illustration of this is the construction of hard benchmarking problems with ‘planted’ solutions defined at the time of construction, which therefore have no special structure related to the problem’s hardness, see for example [Hen, 2019, Hamze et al., 2019].

The paper is structured as follows: In section 2, we review the setting for computation by continuous-time quantum walk encoded into qubits, including application to the search problem. In section 3, we introduce the Sherrington-Kirkpatrick spin glass model, and the random energy model we use for comparison. In section 4, we describe the numerical methods used in this investigation. In section 5, we present the main results showing how quantum walks can find spin glass ground states more effectively than a quantum search algorithm. In section 6, we identify the computational mechanisms and important aspects of the problem structure that contribute to the effectiveness of quantum walk computation. Finally, in section 7, we summarize and conclude.

2 Computing with quantum walks

Both discrete (coined) quantum walks [Aharonov et al., 2001, Shenvi et al., 2003] and continuous-time quantum walks [Farhi and Gutmann, 1998, Childs et al., 2003] are used for computation. This work only uses the continuous-time quantum walk, and also only as an encoded quantum walk, in which qubits are used to store the binary labels of the positions of the quantum walker (see figure 1 for a simple example).

2.1 Continuous-time quantum walks

A continuous-time quantum walk is defined on an undirected graph $G(V,E)$, with $V=\{j\}_{j=0}^{N-1}$ the set of $N$ vertex labels and $E$ the set of label-pairs $(j,k)$ associated with edges. The vertices correspond to the positions of the walker, and the edges indicate the allowed transitions between vertices. This is conveniently encoded in the adjacency matrix $A$ of the graph, which has entries $A_{jk}=1$ for $(j,k)\in E$ and $A_{jk}=0$ otherwise. The Laplacian of $G$ is $L=A-D$, where $D$ is a diagonal matrix formed from the degree of each vertex, $D_{jj}=\mbox{deg}(j)$, where $\mbox{deg}(j)$ is the number of edges connected to vertex $j$. Both the adjacency matrix $A$ and Laplacian $L$ are symmetric matrices which can thus be used to define a quantum Hamiltonian for the dynamics of the continuous-time quantum walk on the graph. In this work, we only need regular graphs, for which $\mbox{deg}(j)$ is constant with respect to $j$. For regular graphs, the only difference between using the adjacency matrix $A$ or Laplacian $L$ is an irrelevant global phase [Childs and Goldstone, 2004]. We use the Laplacian form of the Hamiltonian for consistency with prior work. We thus define the quantum walk Hamiltonian $\hat{H}_{G}$ for a quantum walk on graph $G$ by

[TABLE]

where $\gamma$ is the hopping rate between connected vertices per unit time. The states $\left|j\right\rangle,\left|k\right\rangle$ for $j,k\in V$ are associated with the vertices of $G$ and form a basis for a Hilbert space of dimension $N$. In the Ising model context, the dimension of the Hilbert space is $N=2^{n}$ where $n$ is the number of qubits, and $\{|j\rangle\}_{j=0}^{N-1}$ is the computational basis. For a quantum walk starting in state $\left|\psi(0)\right\rangle$, the state of the walker evolves according to the Schrödinger equation, with formal solution

[TABLE]

using units in which $\hbar=1$.

2.2 Computing using a quantum walk

The task is to solve an optimization problem whose $N=2^{n}$ candidate solutions $j$ are represented in the computational basis $\{\left|j\right\rangle\}_{j=0}^{N-1}$, where $j$ is a bit string corresponding to the state of $n$ qubits. The problem is encoded in an Ising Hamiltonian $\hat{H}_{P}$, of the form described by $\hat{H}_{I}$ in (1) and whose eigenbasis is the computational basis. We write the basis state with eigenvalue $E^{(P)}_{a}$ as $\left|E^{(P)}_{a}\right\rangle$, with $a\in\{0\dots N-1\}$, and adopt the convention that $E^{(P)}_{a}\leq E^{(P)}_{a+1}$. In other words, \Big{\{}\left|E^{(P)}_{a}\right\rangle\Big{\}}_{a=0}^{N-1} is a reordering of $\{\left|j\right\rangle\}_{j=0}^{N-1}$ based on the corresponding eigenenergies of $\hat{H}_{P}$. The encoding is chosen such that the solution corresponds to the ground state $\left|E^{(P)}_{0}\right\rangle$ of the problem Hamiltonian $\hat{H}_{P}$.

To use a quantum walk to solve the problem, we must first choose a suitable state in which to initialize the system. With no prior knowledge of the solution, the equal superposition of all basis states

[TABLE]

is a sensible choice that avoids bias. More generally, the initial state can be prepared as weighted or biased superposition, to incorporate prior knowledge about the solution [Perdomo-Ortiz et al., 2011, Duan et al., 2013, Chancellor, 2017, Graß and Lewenstein, 2017, Baldwin and Laumann, 2018, Kechedzhi et al., 2018, Graß, 2019]. Next, we choose a suitable walk graph $G$. The main requirement is that the ground state of the quantum walk Hamiltonian $\hat{H}_{G}$ coincides with the initial state, either biased or unbiased (see section 6.2). A simple way to achieve a biased starting state would be to ‘tilt’ the driver fields so they are no longer completely transverse. We only treat the unbiased case in this work, so our initial state will be $\left|\psi(0)\right\rangle$ throughout. The full Hamiltonian $\hat{H}(\gamma)$ is defined by adding the quantum walk Hamiltonian $\hat{H}_{G}$ to the problem Hamiltonian $\hat{H}_{P}$

[TABLE]

where the key parameter is the hopping rate $\gamma$ in $\hat{H}_{G}$, see (5). The computation is performed by evolving the initial state (7) under the full Hamiltonian $\hat{H}(\gamma)$ for a time $t_{f}$, then measuring the qubit register in the computational basis. The intuition, based on the faster spreading of quantum walks over classical found in prior work [Farhi and Gutmann, 1998], is that the quantum walk dynamics provide rapid exploration of the basis states, while the energy structure of the problem Hamiltonian $\hat{H}_{P}$ causes localisation around low-energy states.

The success probability $P(t_{f})=\left|\left\langle E^{(P)}_{0}\middle|\psi(t_{f})\right\rangle\right|^{2}$ of finding the solution state when measuring will not in general be unity. It will typically be necessary to repeat the protocol multiple times to obtain a high probability of success over all the repeats. In general, it will be best to use different measurement times $t_{f}$ for each repeat. Different measurement times will produce different success probabilities $P(t_{f})$, and varying the measurement time avoids repeatedly measuring at a time for which the probability $P(t_{f})$ happens to be atypically small. More precisely, we choose the measurement time $t_{f}$ uniformly at random in an interval $[t,t+\Delta t]$, and define an average single run success probability

[TABLE]

Operationally, choosing the measurement time $t_{f}$ randomly in the interval $[t,t+\Delta t]$ samples success probabilities from the distribution with $\bar{P}(t,\Delta t)$ as its mean. Sampling measurement times in this way means that the protocol typically needs to be repeated $M_{\mathrm{rep}}\sim 1/\bar{P}(t,\Delta t)$ times to achieve an overall $O(1)$ success probability. Note that it is not generally possible to check whether the state measured is indeed the ground state of $\hat{H}_{P}$. However, it is easy to calculate the energy of the state measured in each repeat. If only the lowest energy state is accepted, it is only necessary for the ground state of $\hat{H}_{P}$ to be measured once out of all the repeats. The more repeats, the more confidence is gained that the lowest energy state found is the ground state. And studying the distribution of the sampled energies can provide more information about the problem.

The procedure described in this subsection does not in general provide an optimal quantum algorithm, because the repeats do not use information gained from the outcomes of previous runs. We will discuss this further in section 7; for most of this paper we are concerned with understanding the average single run success probability, as an essential prerequisite to building optimal algorithms.

In the limit of small interval width $\Delta t$, the average success probability defined in (9) reduces to the single time probability $P(t_{f})=\lim_{\Delta t\rightarrow 0}\bar{P}(t_{f},\Delta t)$. The long time limit of this average,

[TABLE]

is particularly useful, because it can be calculated via a numerical diagonalization of the Hamiltonian (see section 4) and it predicts the short time average well (see subsection 5.3). In this paper, we will often use the long time average $P_{\infty}$ as an indication of the success probability achievable in a single run, and thus the number of repeats required to achieve $O(1)$ success probability overall. We will separately address the timescale required to reach this probability in each run.

2.3 Graph choice for quantum walk computing

There are many graph-based Hamiltonians with the initial state $\left|\psi(0)\right\rangle$ defined in (7) as the ground state. A common choice is the complete graph $K$, in which every vertex is connected to every other. This graph has the quantum walk Hamiltonian $\hat{H}_{K}$ that couples every computational basis state $\left|j\right\rangle$ state to every other,

[TABLE]

The complete graph is useful because it makes some algorithms analytically tractable [see, e.g., Childs and Goldstone, 2004]. However, for implementation on qubit-based hardware, the complete graph is not in general practical, requiring higher order interaction terms than the transverse Ising term (3). In this qubit setting, an implementation of the complete graph requires a sum over every one-body term (e.g $\hat{X}_{j}$), every two-body term (e.g $\hat{X}_{j}\hat{X}_{k}$), every three-body term (e.g $\hat{X}_{j}\hat{X}_{k}\hat{X}_{l}$) … up to the $n$-body term $\prod_{j=0}^{n-1}\hat{X}_{j}$, a total of $N$ terms. One- and two-body terms are relatively easy to implement, since they correspond to Hamiltonians found naturally. Terms in three or more Pauli-$X$ operators are much more difficult and generally require extra qubits to engineer in real physical systems.

A more natural choice of graph for qubits is the hypercube. The $n$-bit labels are associated with the vertices of the graph such that the edges correspond to flipping one bit, as illustrated in figure 1. The hypercube quantum walk Hamiltonian $\hat{H}_{h}$ on $n$ qubits is composed of single-body terms

[TABLE]

With $\hat{H}_{h}$ as the graph Hamiltonian, the full quantum walk computational Hamiltonian $\hat{H}(\gamma)$ defined in (8) is a transverse Ising Hamiltonian in the form of $\hat{H}_{\mathrm{TI}}$ in (4), with the control functions $A(t)$ and $B(t)$ kept constant throughout the computation.

In this work, we predominantly use the hypercube graph, with some comparisons made with the same problems on the complete graph.

2.4 Solving the search problem using quantum walks

The simplest example of an algorithm in this continuous-time quantum walk setting is the search problem. The problem is to find the marked state, a single bit-string $m\in\{0,1\}^{n}$ out of $N=2^{n}$ possible bit strings. Finding a marked state was shown to have a quantum algorithm with a speed up over classical algorithms by Grover [1996]. To map this problem to the continuous-time Hamiltonian setting, the marked basis state $\left|m\right\rangle$ is given one less unit of energy than all the rest of the basis states, by defining the problem Hamiltonian $\hat{H}_{\mathrm{S}}$ as

[TABLE]

By construction, the problem Hamiltonian $\hat{H}_{\mathrm{S}}$ has the marked state $|m\rangle$ as its ground state.

The continuous-time quantum walk search problem has been analytically solved [Childs and Goldstone, 2004] for several different walk graphs. For the complete graph and the hypercube graph, a quantum speed up is obtained for carefully chosen optimal values of the hopping rate $\gamma$. For the complete graph Hamiltonian $\hat{H}_{K}$, the optimal value is $\gamma^{(K)}_{\mathrm{opt}}=1/N$, while for the hypercube Hamiltonian, $\hat{H}_{h}$, the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ is given by

[TABLE]

where ${{n}\choose{r}}=\frac{n!}{r!(n-r!)}$ is the binomial coefficient. For a quantum speed up, the hopping rate must be set to $\gamma^{(h)}_{\mathrm{opt}}$ as defined by (14) with high precision. It has been shown [Morley et al., 2019] that the fractional tolerance to misspecification of the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ falls as $O(N^{-1/2})$.

The measurement time must also be chosen appropriately. In the limit of large problem size $N$, the marked state can be found with unit success probability, \lim_{N\rightarrow\infty}\Big{[}P(t^{\mathrm{(opt)}}_{f})\Big{]}=1, by measuring in the computational basis at an optimal measurement time $t^{\mathrm{(opt)}}_{f}$. For both the hypercube and complete graphs, the optimal time $t^{\mathrm{(opt)}}_{f}$ scales with the square-root of the problem size $N$ as $t^{\mathrm{(opt)}}_{f}\simeq\frac{\pi}{2}N^{1/2}$. This corresponds to a quadratic speed up compared to the best classical algorithm. Due to the absence of structure in the search problem specifically, such a quadratic speed up has been proven to the best possible quantum speed up [Bennett et al., 1997].

The variation of $P(t_{f})$ with $t_{f}$ is shown in figure 2 for search on hypercube graphs of size $N=2^{30}$ (i.e., $n=30$ qubits) and $N=2^{11}$ (i.e., $n=11$ qubits), using the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$.

The sinusoidal oscillations of the probability $P(t_{f})$ occur because the quantum walk is performing Rabi oscillations between the initial state and the marked state. The two lowest energy levels of the full Hamiltonian $\hat{H}(\gamma)$ with varying $\gamma$ undergo an avoided level crossing at $\gamma^{(h)}_{\mathrm{opt}}$ and the associated eigenstates $\left|E_{0}(\gamma^{(h)}_{\mathrm{opt}}\right\rangle)$ and $\left|E_{1}(\gamma^{(h)}_{\mathrm{opt}}\right\rangle)$ are approximately the orthogonal equal superpositions of the starting state and marked state, $\left|E_{0,1}(\gamma^{(h)}_{\mathrm{opt}}\right\rangle)\simeq(\left|\psi(0)\right\rangle\pm\left|m\right\rangle)/2^{\frac{1}{2}}$. The gap $E_{1}(\gamma^{(h)}_{\mathrm{opt}})-E_{0}(\gamma^{(h)}_{\mathrm{opt}})$ scales with the problem size $N$ as $O(N^{-1/2})$ [Childs and Goldstone, 2004].

These simple, two-level dynamics describe the quantum walk solution to the search problem well for large problem size $N$: the oscillations in the $N=2^{30}$ case have no visible irregularities. For smaller sizes, finite-size effects due to population of higher energy levels are apparent: the oscillations in the $N=2^{11}$ case have lower probability peaks and show some irregular behaviour, such as the small dip on the first peak. These finite-size effects are further illustrated in figure 2, which shows the instantaneous success probability $P(t_{f})$ at the asymptotically optimal and numerically determined best times, as well as the infinite-time average success probability $P_{\infty}$ defined in (10). All three probabilities show a pronounced dip around $n=8$ qubits, with smooth behaviour only settling in for $n>12$ qubits. Figure 2 also shows that the infinite-time probability $P_{\infty}$ asymptotes to a half. Hence, a quantum walk search with a random measurement time should on average only need to be repeated twice to locate the marked state; knowing the exact time to measure for the optimal success probability is not necessary for the success of the algorithm. Fixed point quantum search algorithms [Yoder et al., 2014, Dalzell et al., 2017] are another approach that avoids the need to know how long to run the algorithm for.

The search problem in the continuous-time quantum computing setting has two important drawbacks. Firstly, implementing the problem Hamiltonian $\hat{H}_{S}$ directly on $n$ qubits requires $O(2^{n})$ terms of products of up to $n$ Pauli-$Z$ operators, similar to the problem with implementing the complete-graph Hamiltonian $\hat{H}_{K}$, defined in (11), on qubits. Implementing higher order Pauli-$Z$ terms can be done using extra qubits as “gadgets”, e.g., [Jordan and Farhi, 2008]. An alternative type of gadget, specifically for permutation-symmetric problems like search, is given in [Dodds et al., 2019], building on classical problem mapping techniques in [Chancellor et al., 2016, 2017]. Secondly, it is impossible to map the problem Hamiltonian to qubits without specifying the solution outright. Hence, the search problem serves as a useful toy problem, especially in contexts where having analytic, computational, and physical implementations available for comparisons facilitates benchmarking and other testbed procedures.

3 Spin glass problem definitions

In this work we focus on spin glass problems that have features in common with real life hard optimizations problems and, unlike the search problem, do not admit analytic solutions. The search problem solved by quantum walk provides useful comparisons with these spin glass problems.

3.1 Sherrington-Kirkpatrick spin glass

The Sherrington-Kirkpatrick (SK) spin glass Hamiltonian $H_{\mathrm{SK}}$ [Sherrington and Kirkpatrick, 1975] is defined on $n$ spins as

[TABLE]

where $S_{j}$ are the classical spins ($S_{j}\in\{-1,1\}$) and the couplings $J_{jk}$ are drawn independently from the normal distribution $\mathcal{N}(\mu,\sigma_{\mathrm{SK}}^{2})$ with mean $\mu$ and variance $\sigma_{\mathrm{SK}}^{2}$. Finding the ground state of this Hamiltonian is NP-hard [Choi, 2010], and uniformly hard, due to its finite-temperature phase transition [Sherrington and Kirkpatrick, 1975].

It is computationally convenient to break the spin inversion symmetry by adding single-body field terms of the form $\sum_{j=0}^{n-1}h_{j}S_{k}$, where $h_{j}$ are the field strength values. Like the couplings $J_{jk}$, the fields $h_{j}$ are also drawn independently from $\mathcal{N}(\mu,\sigma_{\mathrm{SK}}^{2})$. When the fields strengths $h_{j}$ are drawn from the same distribution as the coupling strengths $J_{jk}$, the hardness of finding the ground state follows directly from the hardness of the $h_{j}=0$ case. The SK spin glass with such fields is mathematically equivalent to a zero field spin glass with one more spin which is “fixed” in one orientation. This is not true in general for different distributions of field strength $h_{j}$. There are known examples in which fields can destroy spin glass behaviour [see, e.g., Young and Katzgraber, 2004, Feng et al., 2014]. In particular, if the field strengths are much larger than the coupling strengths ($|h_{j}|\gg|J_{jk}|$ for all $j,k$), then the energy is minimized trivially when all the spins each minimize the energy with respect to their individual fields. While the distribution of field strengths could be used to tune the problem hardness, we do not use it in this way here, and only consider cases where the field and coupling strengths are drawn from the same distribution.

An astute reader will notice that if one effectively un-fixes the spin which corresponds to the fields (thus making all states two fold degenerate and converting the system to a double cover of the orignal system), these couplings will effectively be on average stronger by a factor of $\sqrt{2}$. As this increase in coupling strength does not scale with the number of spins, it is going to become less and less significant as the size of the system is scaled up the hardness will be preserved.

The mapping into the quantum Ising model is almost trivial: the classical spin variables $S_{j}$ are simply mapped to Pauli-$Z$ operators. Thus, the problem Hamiltonian $\hat{H}_{\mathrm{SK}}$ becomes

[TABLE]

The SK problem Hamiltonian differs from the search problem by having structure, produced by the $\hat{Z}_{j}\hat{Z}_{k}$ terms. As a result, the covariances between the energies of two basis states depends on the Hamming-distance between them [Baldwin and Laumann, 2018]. Knowing the energy of one state gives some information about the energy of states that differ by a small number of bit-flips. This results in a distribution of the eigenenergies that is almost normal (as can be seen by plotting the distributions and numerically calculating moments), but which deviates from normal in the tails of the distribution.

3.2 Random energy model

To isolate the effect of the correlations in the SK problem, we compare it with the random energy model (REM) [Derrida, 1980], in which the eigenenergies themselves are independently drawn from a normal distribution. The problem Hamiltonian $\hat{H}_{\mathrm{REM}}$ for REM is

[TABLE]

with $\{\left|j\right\rangle\}_{j=0}^{N-1}$ the computational ($Z$) basis and the energies $F_{j}$ drawn independently from the normal distribution $\mathcal{N}(0,\sigma_{\mathrm{REM}}^{2})$.

REM has a similar energy level distribution to that of SK, apart from the tails. By definition it lacks the correlations: knowing the energy of one state gives no information about the energies of other states. Comparison between these two models highlights the effect of the pairwise structure in the SK model.

4 Numerical methods

The main tool used for the investigations in this work is numerical simulation. We are studying computationally hard problems for which there are no tractable analytical solutions except in special cases.

For each number of qubits $5\leq n\leq 20$ we generated 10,000 random instances of the SK spin glass Hamiltonian, defined in (16), with the couplings $J_{jk}$ and fields $h_{j}$ drawn with a standard deviation $\sigma_{\mathrm{SK}}=\omega_{\mathrm{SK}}$, where $\omega_{\mathrm{SK}}$ is an arbitrary energy unit. The value $\omega_{\mathrm{SK}}=5$ was used for computational convenience. We also generated 10,000 random instances of the REM Hamiltonian, defined in (17), for each number of qubits $5\leq n\leq 15$, with normally-distributed energies $F_{j}$ drawn with a standard deviation $\sigma_{\mathrm{REM}}=\omega_{\mathrm{REM}}$. The value $\omega_{\mathrm{REM}}=1$ was used for computational convenience. Note that choosing any arbitrary constant for $\omega$ will only affect overall time and energy scales by a constant factor, and the energy unit $\omega_{\mathrm{SK}}$ has been scaled out of the plots where relevant.

The key quantity to determine numerically is the probability that the ground state is found by running a quantum walk computation on each spin glass instance. It is particularly convenient to compute the infinite-time probability $P_{\infty}$ given by (20), for sizes where full diagonalization is possible. Writing the spectral expansion of the full computational quantum walk Hamiltonian as

[TABLE]

with indices ordered such that $E_{a}(\gamma)\leq E_{a+1}(\gamma)$ and $\left|E_{a}(\gamma)\right\rangle$ the eigenstate with eigenvalue $E_{a}(\gamma)$, we can write the instantaneous probability in terms of the spectral expansions as

[TABLE]

Assuming no degeneracy (that is, all gaps $E_{a}-E_{b}$ are nonzero), which is justified for the randomized nature of the SK and REM problems, the oscillatory terms cancel in the infinite limit (because $\intop_{0}^{\infty}\mathrm{d}t\exp{(-it\theta)}=0$ for nonzero $\theta$) to leave the infinite-time average probability $P_{\infty}$ given by

[TABLE]

All of the numerical simulation in this work has been performed using the Python3 language [Van Rossum and Drake, 2003], aided extensively by the IPython [Perez and Granger, 2007] interpreter and the Jupyter Notebook [Kluyver et al., 2016] system. The numerical heavy-lifting has been done using NumPy [Oliphant, 2006], SciPy [Jones et al., 2001–], and pandas [McKinney, 2010], and the plotting has been done using matplotlib [Hunter, 2007]. The dynamical simulations have been performed by computing the action of the propagator $\exp{(-it\hat{H}(\gamma))}$ on the initial state $\left|\psi(0)\right\rangle$, using the sparse matrix functions within SciPy when possible. For the more computationally demanding analyses, we were limited to $n\leq 11$ by the computational resources available. Where relevant, figures in this paper have error bars included. However, in most cases the error bars are much smaller than the size of the marker symbols used and so are not visible. This is due to the size of the data sets (10k instances per value of $n$), which provides a good level of accuracy for the average quantities.

Simulations were run on the Imperial and Durham University high performance computing facilities. The data for all the instances used is available on a permanent data archive [Chancellor et al., 2019].

5 Quantum walks with spin glasses

In order to implement a quantum walk algorithm for finding the ground states of the spin glasses defined in section 3, we follow the procedure described in section 2.2: Choose a quantum walk graph $G$ and associated Hamiltonian $\hat{H}_{G}$, and add the spin glass Hamiltonian to get the full computational quantum walk Hamiltonian $\hat{H}(\gamma)=\hat{H}_{G}+\hat{H}_{P}$, where $\hat{H}_{P}$ refers to $\hat{H}_{\mathrm{SK}}$ or $\hat{H}_{\mathrm{REM}}$ as appropriate. Since the hypercube is the natural choice of graph for qubit implementations, we use this graph, with quantum walk Hamiltonian $\hat{H}_{h}$ defined in (12), unless otherwise indicated. For the initial state $|\psi(0)\rangle$, we use the equal superposition (7), which is the ground state of the hypercube Hamiltonian $\hat{H}_{h}$.

5.1 Setting the hopping rate

In contrast to the search problem, for SK and REM it is impossible to efficiently calculate the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ that maximizes the success probability. It is not even clear which measure of success probability should be maximized because, unlike the search problem, there will be no efficient way to find the optimal measurement time $t^{\mathrm{(opt)}}_{f}$ for any choice of hopping rate $\gamma$. To bootstrap the investigation, we choose to define the optimal hopping-rate $\gamma^{(h)}_{\mathrm{opt}}$ with respect to one of the average probabilities defined in (9); in particular, we choose the hopping rate that maximizes the infinite-time average probability $P_{\infty}$ defined in (10). We make this choice because the infinite-time average probability $P_{\infty}$ is numerically convenient to calculate, and because it has been seen to be a relevant measure of probability in the search example, see figure 2. We will see in Subsection 5.3 that the probability $P_{\infty}$ typically agrees well with probabilities averaged over shorter and more practical time windows.

Some plots of the infinite-time probability $P_{\infty}$ against hopping rate $\gamma$ for typical 11-qubit examples of the SK and REM are shown in figure 3. Note that the maximal success probability varies by an order of magnitude between the two problem-types, with REM highest and SK lowest. While the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ is instance-dependent, these plots show that the dependence of infinite-time probability $P_{\infty}$ on hopping rate $\gamma$ is typically characterised by broad, bumpy peaks for SK, and by narrow, well-defined peaks for REM. This implies that a precise value of the hopping rate $\gamma$ is needed for REM, while there is some tolerance to non-optimal values of the hopping rate $\gamma$ for SK for the sizes that we have studied.

To investigate the success probability more systematically, we performed a brute-force numerical search to find the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ that maximizes the success probability $P_{\infty}$ for each spin glass instance from the data sets of 10k random instances for $5\leq n\leq 11$. This gives a baseline maximum average single run success probability for the quantum walk algorithm.

The optimal hopping rates $\gamma^{(h)}_{\mathrm{opt}}$ correspond to the best a quantum walk algorithm on the hypercube can possibly do in a single run. For practical algorithms, we need a heuristic method for choosing the hopping rate that can be calculated from the known parameters. For the quantum walk search algorithm, the optimal hopping rate balances the energy between the two components of the Hamiltonian, $\hat{H}_{P}$ and $\hat{H}_{G}$. Guided by this, we define the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ for SK and REM such that it balances these overall energy-scales on average. We match the energy-spread $E_{N-1}^{(h)}-E_{0}^{(h)}$ of the hypercube quantum walk Hamiltonian $\hat{H}_{h}$ with the average energy-spread $\langle E_{N-1}^{(P)}-E_{0}^{(P)}\rangle$ of the problem Hamiltonian $\hat{H}_{P}$. For the hypercube Hamiltonian $\hat{H}_{h}$ defined in (12), we have the energy spread $E_{N-1}^{(h)}-E_{0}^{(h)}=2n\gamma$; hence, we define the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ by

[TABLE]

To demonstrate that this heuristic is sensible, we compare in figure 4 the distributions of optimal hopping rates $\gamma^{(h)}_{\mathrm{opt}}$ for SK (blue) and REM (red), as well as the heuristic hopping rates (black, dashed and dotted lines for SK and REM respectively) calculated according to (21), for the 11-qubit data set. For both SK and REM, the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ falls in the centre of the $\gamma^{(h)}_{\mathrm{opt}}$ distributions. Note that the SK distribution is much broader than for REM: not only are the individual peaks for $\gamma^{(h)}_{\mathrm{opt}}$ for SK much broader than for REM (figure 3), but the distribution of the maxima of those peaks is also much broader (figure 4). This may seem to be a problem for specifying a heuristic value for $\gamma$ for SK from average energies, but as we will show, it is actually REM that fails for the heuristic $\gamma$, while SK works well.

For a normal distribution of energy levels, the average problem energy spread can be estimated as

[TABLE]

where $\sigma_{P}^{(\mathrm{energy})}$ is the standard deviation of the energy eigenvalues of the problem Hamiltonian. For REM, the standard deviation $\sigma_{\mathrm{REM}}^{(\mathrm{energy})}$ is equal to the energy unit $\omega_{\mathrm{REM}}$ by definition (see section 3.2). For SK, the standard deviation $\sigma_{\mathrm{SK}}^{(\mathrm{energy})}$ can be shown to be equal to $\frac{\omega_{\mathrm{SK}}}{2}[n(n+3)]^{\frac{1}{2}}$. Equation (22) is accurate for REM (which has normally-distributed energy levels by definition) but, as already noted, the distribution of the eigenenergies in SK deviates from normal, especially in the tails. Numerically, we find that there is a multiplicative constant factor of approximately 0.887 that corrects the formula in (22) for SK for the effects of the non-normal tails. For the numerical analysis, we use the numerically calculated average energy-spread at each number of qubits $n$.

Figure 5 compares the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ and average optimal hopping rate $\langle\gamma^{(h)}_{\mathrm{opt}}\rangle$ at different numbers of qubits $5\leq n\leq 11$. The full width at half maximum (FWHM) has also been calculated for each instance, to estimate the tolerance $\Delta\gamma^{(h)}_{\mathrm{opt}}$ to deviations from the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ (illustrated in figure 3). The width of the shaded regions in figure 5 corresponds to the average tolerance range $\langle\Delta\gamma^{(h)}_{\mathrm{opt}}\rangle$ at each number $n$ of qubits. While the heuristic hopping rate differs slightly from the the average optimal hopping rate for SK, the average tolerance range $\langle\Delta\gamma^{(h)}_{\mathrm{opt}}\rangle$ is much broader, and does not shrink with increasing number of qubits $n$. For REM, however, while we see close agreement on average, the tolerance range shrinks quickly with the number of qubits $n$ as the peaks (as in figure 3, right) become narrower. This means that the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ is more likely to lie further than $2\Delta\gamma^{(h)}_{\mathrm{opt}}$ outside of the actual probability peak for each instance, even though it agrees well with the average optimal hopping rate $\langle\gamma^{(h)}_{\mathrm{opt}}\rangle$. Consequently, a quantum walk with the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ does not perform well for most REM instances.

It is instructive to quantify this sensitivity to deviations from the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$. Figure 5 shows log-linear and log-log plots of the average fractional tolerance range $\langle\Delta\gamma^{(h)}_{\mathrm{opt}}/\gamma^{(h)}_{\mathrm{opt}}\rangle$ against number $n$ of qubits for SK (blue circles), REM (red squares) and search (green triangles) on the hypercube. For SK, the fractional tolerance range $\langle\Delta\gamma^{(h)}_{\mathrm{opt}}/\gamma^{(h)}_{\mathrm{opt}}\rangle$ decreases as approximately $1/n$, while for REM and search the decrease is approximately $N^{-0.5}$. This decrease is expected theoretically for search [Childs and Goldstone, 2004]. The fitted lines do not show exactly a square-root dependence (exponent of $-0.5$) due to the finite size effects for small numbers of qubits $n\leq 12$.

Thus, we see that REM behaves like the search problem in a quantum walk setting. For a precisely optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$, the success probability is high, but this instance-dependent hopping rate is hard to predict, unlike for the analytically tractable quantum walk search algorithm. Without this precise hopping rate, quantum walks perform no better than guessing for the search problem and for REM. In contrast, quantum walks applied to SK give a better-than-guessing success probability $P_{\infty}>1/N$ for the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ calculated according to (21).

With the conditions under which we can achieve a better-than-guessing success probability characterised for the three problem types, SK, REM, and search, we turn to the scaling of this success probability with problem size $N$.

5.2 Success probability

Figure 6 shows how the single-time success probability $P(t_{f})$ varies with the measurement time $t_{f}$ for two typical 11-qubit examples of SK and REM. In the REM case, the behaviour is similar to that shown in figure 2 for search: an oscillatory nature indicating the dominance of a two-level avoided-crossing feature, but with evidence of the population of other energy-levels that lead to finite-size effects in search. For REM, these finite-size effects are more pronounced, and are instance-dependent. The random nature of the REM problems means there is not such a clear cut off size, as there is for the search problem, above which finite size effects are negligible. In any case, based on search, we expect finite size effects to be significant at $n=11$. For SK, the behaviour is quite different from search or REM. There is no indication of dominant oscillatory behaviour; instead, these plots show unpredictable, highly instance-dependent fluctuating dynamics for all the sizes we are using. This indicates that for SK, the behaviour is determined by the excitation of many energy levels.

As with finding a suitable hopping-rate $\gamma$, both REM and SK differ from the search problem in that there is no practical way to find the optimal measurement time $t^{\mathrm{(opt)}}_{f}$; a different approach must be taken instead. As already noted for the search problem, this can be handled by using the time averaged probabilities defined in (9). We first consider the infinite-time probability $P_{\infty}$, as defined in (10), since it is easy to calculate (see section 4).

Figure 7 shows the average infinite-time success probability $\langle P_{\infty}\rangle$ against the number $n$ of qubits for the two problems using both the optimal $\gamma^{(h)}_{\mathrm{opt}}$ and heuristic $\gamma^{(h)}_{\mathrm{heur}}$ hopping rates. For SK, this gives exponential decay with the number of qubits $n$ in both cases: the average probability $\langle P_{\infty}\rangle$ changes with $n$ according to

[TABLE]

where $\tilde{O}$ may neglect factors logarithmic in its argument. That is, using the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ instead of the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ has only a minor impact on the average success probability $\langle P_{\infty}\rangle$.

For REM, the behaviour is quite different. With the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ we see a success probability $P_{\infty}$ of constant order but with a pronounced dip. This behaviour is similar to that seen for the search problem, where the dip seen in figure 2 is a finite-size effect. This similarity is expected, given the similarity between the dynamical behaviour shown in figure 2 for search and in figure 6 for REM. With the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ for REM, we see a significantly reduced success probability $P_{\infty}$ compared to the optimal case. That is, the heuristic is performing poorly, despite the good agreement shown in figure 5.

The clear difference in behaviour between SK and REM can be explained by the different tolerances $\Delta\gamma^{(h)}_{\mathrm{opt}}$ to deviations from the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ shown in figure 5 and figure 5. For SK, the tolerance range is broad enough for the heuristic to lie within it, while for REM the heuristic hopping rate $\gamma_{\mathrm{heur}}^{(h)}$ almost always misses this range entirely even though it is close to the average optimal hopping rate $\left\langle\gamma_{\mathrm{opt}}^{(h)}\right\rangle$.

5.3 Mixing times

We have thus numerically determined an average success probability scaling with problem size of ${\sim\tilde{O}(N^{-0.42})}$ for a quantum walk finding SK spin glass ground states, using the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$. This is based on the infinite time-success probability $P_{\infty}$, i.e., uniform sampling from the distribution of all possible run times. We now investigate the time dependence in more detail: can we sample from a finite run time and still obtain the same speed up? Since $P(0)=1/N$ corresponds to random guessing, there must be a minimum time before which it is not effective to measure.

We define a mixing-time $\tau_{\mathrm{mix}}^{(\epsilon)}$ to be the latest time, $t$, for which the time averaged probabilities $\bar{P}(0,t)$ and $\bar{P}(0,2t)$ at the two times $t$ and $2t$ differ by a fraction greater than the fluctuation parameter $\epsilon$,

[TABLE]

This definition of $\tau_{\mathrm{mix}}^{(\epsilon)}$ is based on similar definitions found in prior work [Aharonov et al., 2001], with modifications for computational convenience. We numerically estimated the mixing-time $\tau_{\mathrm{mix}}^{(0.05)}$ for each SK instance up to $n=11$ qubits, using the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ for each instance. We simulated the quantum walk computation dynamics for a successively-doubling duration until a time at which the condition is met was reached. The fluctuation parameter $\epsilon=0.05$ corresponds to a deviation of $5\%$. To verify that the mixing-time $\tau_{\mathrm{mix}}^{(0.05)}$ correctly captures the relevant dynamical timescale, we also numerically estimated it for the search problem at each system size from $n=5$ to $n=30$ qubits. The search problem using continuous-time quantum walks can be mapped to the symmetric subspace, allowing larger sizes to be analysed. The mixing-time $\tau_{\mathrm{mix}}^{(0.05)}$ for search exhibits the expected exponential timescale: the solid green line of best fit in figure 8 has the expected scaling with problem size $N$ of $\tau_{\mathrm{mix}}^{(0.05)}=\tilde{O}(N^{1/2})$.

For search, the scaling is dominated by the run time, the success probability is $O(1)$. However, this behaviour only emerges clearly above $n\sim 20$. Below this, the behaviour is influenced by the finite-size effects that arise due to population of higher energy levels. This means it is not useful to analyse the behaviour of the REM time scaling, finite size effects mask the scaling behaviour for computationally tractable sizes. However, unlike search and REM, the SK behaviour is influenced by higher energy levels at all sizes, through the frustration provided by the random couplings between the spins. Hence, we do not expect to see such finite-size effects in SK; the behaviour is already dominated by the frustration at small sizes. Figure 8 shows a log-log plot of the mixing-time (scaled by $\omega_{\mathrm{SK}}$) averaged over the ensemble $\langle\tau_{\mathrm{mix}}^{(0.05)}\omega_{\mathrm{SK}}\rangle$. The solid blue line of best fit has a logarithmic scaling with problem size $N$ of

[TABLE]

Thus it contributes a logarithmic factor to the overall scaling. We emphasise that while this single-run timescale is polynomial in the number of spins $n$, the overall timescale is still exponential in $n$ due to the exponential number of repeats required to achieve $O(1)$ success probability.

To confirm the subsidiary nature of the time scaling for each SK run, we show in figure 9 a log-plot comparing, for the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$, the success probability $P_{\infty}$ in the infinite-time case (as in figure 7) and in the case of an early, logarithmically-scaling (with respect to $N$) measurement window $12.5n^{\frac{1}{2}}/\omega_{\mathrm{SK}}\leq t\leq 17.5n^{\frac{1}{2}}/\omega_{\mathrm{SK}}\equiv(t_{\mathrm{short}},\Delta t_{\mathrm{short}})$. This $n^{0.5}$ scaling of the window is even shorter than the fitted scaling of $n^{0.75}$, although at these sizes the difference is not significant. This finite-time probability $\bar{P}(t_{\mathrm{short}},\Delta t_{\mathrm{short}})$ is similar to the infinite-time probability $P_{\infty}$: the solid purple line of best fit in figure 9 has a scaling with problem size $N$ of

[TABLE]

This should be compared with (23), where the value of the exponent for the average infinite time success probability $P_{\infty}$ with the heuristic hopping rate $\gamma^{(h)}_{\mathrm{heur}}$ is given by $-0.417\pm 0.002$.

As the dominant factor in the total runtime comes from the required number of repeats, and because the single-run timescale contributes only a logarithmic factor, these results constitute good numerical evidence for an average total runtime which scales with problem size $N$ as ${\sim\tilde{O}(N^{0.41})}$ for using quantum walks to find spin glass ground states, over the range of $N$ in our data sets. This scaling is a better than the best possible (quadratic) speed up achievable for quantum walk search algorithms. Moreover, it comes without the requirement for exponential precision in setting the hopping rate that renders practical use of quantum walk searching difficult for large problems. We now present some insights into where the improvement over search comes from.

6 Computational mechanisms

6.1 Role of correlations in SK

To investigate whether the energy correlations with Hamming distance in SK play a significant role in the computational process of finding the ground state with a quantum walk, we performed three additional sets of numerical tests.

Firstly, we used the same SK instances but performed the quantum walk using a complete graph Hamiltonian $\hat{H}_{K}$, defined in (11), instead of the hypercube graph Hamiltonian $\hat{H}_{h}$. This removes the correspondence of Hamming-distance between classical states with the distance between those states on the graph – for the complete graph, every state is one unit (edge) away from every other state. In terms of the Hamiltonian, the transverse Ising term is replaced by sums of products of up to $n$ Pauli-$X$ operators that flip up to $n$ qubits at the same time, in all possible combinations. For each SK instance up to $n=11$, we estimated the optimal hopping rate $\gamma^{(K)}_{\mathrm{opt}}$ for the complete graph, and then used it to calculate the infinite-time probability $P_{\infty}$.

Secondly, we constructed ‘scrambled SK’ instances, denoted sSK, by randomizing which state corresponds to which energy in the SK instances. In doing so, we arrive at Hamiltonians with identical energy spectra to the SK instances, but without the correlations between energy difference and Hamming distance on the hypercube graph. This approach has similarities with previous work [Farhi et al., 2008, 2011, Hen, 2014]. For each sSK instance, we estimated the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$, which is different from that used for the ordinary SK versions. This hopping rate was then used to calculate $P_{\infty}$.

Thirdly, we sorted the eigenenergies of each REM instance in increasing size and assigned them to the computational basis states in the order of a binary-reflected Gray code on their bitstrings, to arrive at a problem denoted REMGC. In doing so, we added some amount of Hamming-distance structure by ensuring that the closest energies are assigned to states that differ by only a single bit-flip. For each REMGC instance, we estimated an optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$, which is different from that used for the ordinary REM problem. This was used to calculate the infinite-time probability $P_{\infty}$. While REMGC is not a hard problem as defined, it provides a useful example to compare with how the quantum walk finds the ground state of a SK spin glass.

These three variants provide separate tests of the influence of the graph structure (choice of quantum walk Hamiltonian) and problem structure (pairwise correlations in SK). Figure 10 shows how the infinite-time probability $P_{\infty}$ varies with the number of qubits $n$ for these three variants, alongside SK and REM on a hypercube graph from figure 7.

The variation of $P_{\infty}$ with the number of qubits $n$ for the five variants is clearly split into two groups, behaviour like REM and search on the one hand, and behaviour like SK on the other. Removing the correlations from SK by scrambling the energies (sSK) results in behaviour like REM and search. Moreover, removing the correspondence between distance and Hamming weight by using the complete graph instead of the hypercube also changes the SK problem behaviour to be like REM and search. In the opposite direction, inserting pairwise correlations into REM via a Gray code (REMGC) results in problems that are much more like SK than like the REM problems on a hypercube graph.

From this, we infer that the problem structure – in this case the pairwise correlations in SK – needs to be matched by a compatible driver Hamiltonian – in this case the hypercube/transverse Ising – to obtain better than quadratic scaling. This type of local structure in the solution space is exploited in many classical algorithms. For example, classical Monte Carlo optimizations that use a single bit flip update rule are naturally using this hypercube structure. Using a complete graph instead would correspond to flipping a random number of bits, which is equivalent to guessing at each step.

6.2 Energy conservation dynamics

Continuous-time quantum walk time evolution is unitary, and there is no time dependence in the Hamiltonian that can lead to energy gain or loss by the system. Hence, it is important to consider how it can find a lower energy state than it starts in (with respect to $\hat{H}_{P}$) with any better-than-guessing probability. For the search problem, this happens through an analog of Rabi flopping (see figure 2), cycling between the initial and solution states. However, the dominant avoided level crossing structure is not present in the spin glasses to provide this mechanism.

We now show that there is a very generic mechanism (also described independently by Hastings [2019]) that relies on starting in the ground state of the quantum walk part of Hamiltonian $\hat{H}_{G}$. Let $\langle\hat{O}\rangle_{\psi(t)}$ for operator $\hat{O}$ be defined by $\langle\psi(t)|\hat{O}|\psi(t)\rangle=\langle\hat{O}\rangle_{\psi(t)}$. Then, by linearity, and the definition of $\hat{H}(\gamma)$ in (8), the energy expectation at time $t$ is

[TABLE]

Due to the unitarity of the evolution under a time-independent Hamiltonian, this expectation energy will not change over time, giving

[TABLE]

which yields

[TABLE]

As $|\psi(0)\rangle$ is chosen to be the ground state of $\hat{H}_{G}$, the LHS must be non-negative. Furthermore, as $|\psi(0)\rangle$ is not an eigenstate of $\hat{H}(\gamma)$, some dynamics are guaranteed to occur and so the LHS must become positive at early times. Therefore, the RHS must also be non-negative always and positive at early times. Thus, taking any final time $t_{f}$, we get the inequality

[TABLE]

Equation (30) shows that performing time evolution under the computational quantum walk Hamiltonian from the initial state $\left|\psi(0)\right\rangle$ is guaranteed to lower the energy of the system with respect to $\hat{H}_{P}$ (the expectation value $\langle\hat{H}_{P}\rangle_{\psi(t)}$). This implies that the overlap with low energy eigenstates of $\hat{H}_{P}$ will increase, at least for short times. A measurement in the computational basis will thus be on average more likely than a random guess to produce a low energy state.

Starting in a low energy state is thus important for the success of the quantum walk algorithm (we have checked this numerically). It also implies that encoding prior information into the initial state will help, provided this is given in the form of a lower energy state than the uniform superposition state. This could be the final state from a previous run, for example, which will be explored further in Nita et al. [2020]. It is also necessary to bias the quantum walk Hamiltonian so that its ground state matches this biased initial state. Since this starting state is a known computational basis state, it is possible to do this biasing for suitably designed hardware.

For many optimization problem applications, it is helpful to find a low energy state, even if it is not actually the true ground state. From this point of view, that quantum walks necessarily lower the expectation energy with respect to the problem Hamiltonian is very appealing as a computational mechanism. This argument by itself does not provide a guaranteed scaling or quantum speed up, but it does explain how the quantum walk dynamics work in this setting, where there is no way to lose (or gain) energy. It is possible to generalise these arguments beyond time-independent Hamiltonians [Callison et al., 2020b], to include monotonic functions $A(t)$ and $B(t)$ in (4).

To illustrate this energy redistribution mechanism, the plots in figure 11 show how the expectation value $\langle\hat{H}_{G}\rangle_{\psi(t)}$ of the quantum walk Hamiltonian (green solid-line) and the expectation value $\langle\hat{H}_{P}\rangle_{\psi(t)}$ of the problem Hamiltonian (red solid-line) vary during a quantum walk. We have included the instantaneous success probability $P(t)$ (faint grey) to show that the timescale used is long enough for significant dynamics to take place. A typical 10-qubit SK example is shown in figure 11 and a typical 10-qubit REM example is shown in figure 11, both on the hypercube using their respective optimal hopping rates $\gamma^{(h)}_{\mathrm{opt}}$. Also shown is the ground state eigenvalue $\langle\hat{H}_{P}\rangle_{E^{(P)}_{0}}$ of the problem Hamiltonian (red, dash-dotted line) and the ground state eigenvalue $\langle\hat{H}_{G}\rangle_{\psi(0)}$ of the quantum walk Hamiltonian (green, dashed line). In both SK and REM, the initial evolution takes the state away from the $\hat{H}_{G}$ ground state, raising the $\hat{H}_{G}$ expectation value, and thereby lowering the $\hat{H}_{P}$ expectation value to a point around which it fluctuates for the duration simulated. This clearly shows the energy redistribution mechanism at work, and the short time scale over which it appears.

7 Summary and outlook

In this work, we have shown numerically that continuous-time quantum walks are a viable computational method for finding ground states of hard spin glass problems. We have produced strong numerical evidence for a better-than-search polynomial quantum speed up over random guessing, with a scaling of the average single run success probability ${\sim\tilde{O}(N^{-0.41})}$, using data sets of size $5\leq n\leq 20$ spins ($32\leq N\leq 1,048,576$). Moreover, and importantly, this is obtained without the need to set parameters exponentially precisely, as is required for quantum walk search algorithms. The hopping rate $\gamma$, that determines the relative strengths of the quantum walk and problem Hamiltonians, can be estimated from the overall energy scales, which are determined by the hardware and encoding of the problem.

To explain why quantum walks are able to do better than quantum searching in this case, we compared variants on the spin glass problems that remove or add pairwise correlations, and compared the hypercube graph quantum walk Hamiltonian with the complete graph quantum walk Hamiltonian. This showed that the combination of pairwise correlations in the encoding of the problem and a matching single spin flip quantum walk Hamiltonian is required to exploit the correlations. The single spin-flips driven by the transverse field terms $\hat{X}_{j}$ in the hypercube quantum walk Hamiltonian are the correct operators for the pairwise interaction terms $\hat{Z}_{j}\hat{Z}_{k}$ in the spin glass Hamiltonian. A single spin flip on either qubit $j$ or $k$ changes the energy for that term from high to low, or vice versa. Since we can choose how to encode the problems into the Hamiltonians, and there are known methods to convert higher order terms to pairwise terms [Bremner et al., 2002, Dattani, 2019], we can arrange to use this mechanism both for its computational advantages and practicality for hardware implementation as the transverse Ising Hamiltonian.

To explain how quantum walks are able to find low energy states when the closed quantum dynamics have no mechanism for losing energy, we showed how starting in the ground state of the quantum walk part of the Hamiltonian guarantees dynamics that decrease the expectation value of the energy with respect to the problem Hamiltonian. This also ensures that prior information can be provided by starting in lower energy states, from which improved solutions can be found. Exploiting this process will allow an optimal quantum algorithm to be built from multiple quantum walk runs that use the information gained from prior runs. Performing multiple quantum walk runs in early, noisy quantum hardware is a more viable approach than maintaining coherence for sufficiently accurate adiabatic algorithms. Quantum walks may also be simpler to implement since they do not require time dependent controls. This work thus provides a significant advance in understanding how to exploit quantum walks in practical hardware for optimization problems.

It is likely that further insights into the computational effectiveness of quantum walks in this transverse Ising Hamiltonian setting are to be found in current knowledge of spin glass phases in the presence of transverse fields. The spin glass transition itself is not fully understood, in neither the quantum nor classical case [see, e.g., Parisi, 1980, Fisher and Huse, 1987, 1988, Thirumalai et al., 1989, Larson et al., 2013, Young, 2017, Magalhaes et al., 2017]. However, the phases of interest for computation are not the spin glass phases themselves, but the phases where transitions between states are still occurring at a rapid enough rate to find solution states. Extremely long equilibration timescales are a defining property of all glass phases, including spin glasses [Bouchaud et al., 1998, Cugliandolo, 2002]. Since the equilibration (mixing) times $\tau_{\mathrm{mix}}^{(\epsilon)}$ we find in section 5.3 for the SK spin glass only scale polynomially with the number of spins, it is most likely that at the optimal hopping rates $\gamma^{(h)}_{\mathrm{opt}}$, our quantum walks are not in a finite size precursor to a spin glass phase, but rather in a precursor to a paramagnetic phase, for which equilibration times can be fast. Given that the system should localize more in lower energy states for smaller transverse fields, it is reasonable that our optimal hopping rates $\gamma^{(h)}_{\mathrm{opt}}$ occur near the edge of the precursor to the spin glass phase. Furthermore, the mild scaling of the width $\Delta\gamma^{(h)}_{\mathrm{opt}}$ of the peak around the optimal hopping rate $\gamma^{(h)}_{\mathrm{opt}}$ suggests that the regime where quantum walks performs well may correspond to the second paramagnetic phase observed in Magalhaes et al. [2017]. Polynomial gaps have been found around the spin glass–paramagnetic phase transition in a related model in [Knysh, 2016].

A numerical study such as this inevitably leaves open questions regarding the asymptotic scaling of the problems. In particular, we observed a range of hardness in the SK data sets and future work will investigate what fraction of the instances are actually hard for classical algorithms. Forthcoming work applying similar techniques to Max2SAT [Callison et al., 2020a] will characterise the hardness of small random instances in more detail, and establish quantum walks as an effective tool for hard optimization problems more generally. While general methods are known to speed up the best classical algorithms [Hartwig et al., 1984] for this type of problem [Montanaro, 2018, 2019], further work is required to determine whether an optimal continuous-time quantum walk algorithm can be devised that fully leverages the advantage from the correlations. Nonetheless, our work represents a significant advance in developing continuous-time quantum walk computation for hard optimization problems, and provides key insights into the computational mechanisms that can be exploited over short timescales, well-suited to the limited coherence times of noisy, intermediate scale quantum hardware.

VK and NC funded by UK EPSRC fellowship EP/L022303/1 and NC funded by EPSRC grant EP/S00114X/1. AC funded by EPSRC grant EP/L016524/1 via the Imperial College London CDT in Controlled Quantum Dynamics. We thank Prof. Ifan G. Hughes and Dr Ashley Montanaro for helpful discussions.