Efficient Quantum Walk Circuits for Metropolis-Hastings Algorithm
Abstract
We present a detailed circuit implementation of Szegedy’s quantization of the Metropolis-Hastings walk. This quantum walk is usually defined with respect to an oracle. We find that a direct implementation of this oracle requires costly arithmetic operations and thus reformulate the quantum walk in a way that circumvents the implementation of that specific oracle and which closely follows the classical Metropolis-Hastings walk. We also present heuristic quantum algorithms that use the quantum walk in the context of discrete optimization problems and numerically study their performances. Our numerical results indicate polynomial quantum speedups in heuristic settings.
pacs:
I Introduction
Markov chain Monte Carlo (MCMC) methods are a cornerstone of modern computation, with applications ranging from computational science to machine learning. The key idea is to sample a distribution by constructing a random walk which reaches this distribution at equilibrium . One important characteristic of a Markov chain is its mixing time, the time it requires to reach equilibrium. This mixing time is governed by the inverse spectral gap of , where the spectral gap is defined as the difference between the two largest eigenvalues. The runtime of a MCMC algorithm is thus determined by the product of the mixing time and the time required to implement a single step of the walk.
Szegedy Szegedy (2004) presented a general method to quantize reversible walks, resulting in a unitary transformation . The eigenvalues of a unitary matrix all lie on the unit complex circle, and we choose The steady state of the quantum walk has and is essentially a coherent version of the classical equilibrium distribution . The main feature of the quantum walk is that its spectral gap is quadratically larger than the one of the corresponding classical walk. Combined with the quantum adiabatic algorithm Farhi et al. (2000); Aharonov and Ta-Shma (2003); Boixo et al. (2010), this yields a quantum algorithm to reach the steady state that scales quadratically faster with than the classical MCMC algorithm Somma et al. (2008).
While at first glance this is an important advantage with far reaching applications, additional considerations must be taken into account to determine if quantum walks offer a significant speedup for any specific application. One of the reasons is that the time required to implement a single step of the quantum walk can be significantly larger than the time required to implement a single step of the classical walk. Thus, quantum walks are more likely to offer advantages in situations with extremely long equilibration times. Moreover, classical walks are often used heuristically out of equilibrium. When training a neural network for instance, where a MCMC method called stochastic gradient is used to minimize a cost function, it is in practice often not necessary to reach the true minimum, and thus the MCMC runs in time less than its mixing time. Similarly, simulated annealing is typically used heuristically with cooling schedules far faster than prescribed by provable bounds – and combined with repeated restarts. Such heuristic applications further motivate the constructions of efficient implementations of , and the development of heuristic methods for quantum computers.
This article addresses these two points. First, we present a detail realization and cost analysis of the quantum walk operator for the special case of a Metropolis-Hastings walk Metropolis et al. (1953); Hastings (1970). This is a widespread reversible walk, whose implementation only requires knowledge of the relative populations of the equilibrium distribution. While Szegedy’s formulation of the quantum walk builds on a classical walk oracle, our implementation circumvents the direct implementation of the classical walk, which would require costly arithmetic operations. Instead, we directly construct a related but different quantum unitary walk operator with an effort to minimize circuit depth. Second, we suggest heuristic uses of this oracle inspired by the adiabatic algorithm, and study their performances numerically.
Ii Preliminaries
ii.1 Quantum Walk
We define a classical walk on a -dimensional state space by a transition matrix where the transition probability is given by matrix element . Thus, the walk maps the distribution to the distribution , or in matrix form . An aperiodic walk is irreducible if every state in is accessible from every other state in , which implies the existence of a unique equilibrium distribution . Finally, a walk is reversible if it obeys the detailed balance condition
(1) |
We now explain how to quantize an reversible classical walk .
Szegedy’s quantum walk Szegedy (2004) is formulated in an oracle setting. For a classical walk , it assumes a unitary transformation acting on a Hilbert space with the following action
(2) |
where . Define as the projector onto the subspace spanned by states . Combining to the reflection and the swap operator , we can construct the quantum walk defined by
(3) | ||||
(4) |
Szegedy’s walk is defined as , so it is essentially the square of the operator we have defined, but this will have no consequence on what follows aside from a minor simplification.
To analyze the quantum walk , let us define the state and consider the operator
(5) | ||||
(6) | ||||
(7) |
At this point, we must assume that the walk is reversible and make use of the detailed balance condition of Eq. (1) to obtain
(8) |
or, if we restrict the operator to its support , we get in matrix notation and are thus similar so they have the same eigenvalues. Define its eigenvectors . The matrices
(9) |
where are the eigenvalues of . Because the operator is obtained by projecting the operator onto the subspace , its eigenvectors with non-zero engenvalues in the full Hilbert space clearly have the form .
If we consider the action of without those projections, we get
(10) |
where is orthogonal to the subspace , so in particular it is orthogonal to all the vectors . Finally, because is a unitary, we also obtain that the are orthogonal to each other and that . This implies that the vectors are all mutually orthogonal and that is block diagonal in that basis.
Given the above observations, it is straightforward to verify that
(11) | ||||
(12) |
so the eigenvalues of on the subspace spanned by are where with corresponding eigenvectors .
ii.2 Adiabatic state preparation
We can use quantum phase estimation Kitaev (1995) to measure the eigenvalue of . In particular, we want this measurement to be sufficiently accurate to resolve the eigenvalue , or equivalently , from the rest of the spectrum. Assuming that the initial state is supported on the subspace , the spectral gap of is , so we only need about applications of to realize that measurement. This is quadratically faster than the classical mixing-time , which is at the heart of the quadratic quantum speed-up.
A measurement outcome corresponding to would produce the coherent stationary distribution . Indeed, first note that for any such that , Eq. (10) implies that . We can verify that this condition holds for :
(13) | ||||
(14) | ||||
(15) |
where we have used detailed balance Eq. (1) in the second step and in the last step.
From an initial state , the probability of that measurement outcome is , so the initial state must be chosen with a large overlap with the fixed point to ensure that this measurement outcome has a non-negligible chance of success. If no such state can be efficiently prepared, one can use adiabatic state preparation Farhi et al. (2000); Aharonov and Ta-Shma (2003) to increase the success probability. In its discrete formulation Somma et al. (2008) inspired by the quantum Zeno effect, we can choose a sequence of random walks with coherent stationary distributions . The walks are chosen such that is easy to prepare and consecutive walks are nearly identical, so that . Thus, the sequence of measurements of the eigenstate of the corresponding quantum walk operators all yield the outcomes with probability , which results in the desired state. The overall complexity of this algorithm is
(16) |
where is the spectral gap of the -th quantized walk and is the time required to implement a single quantum walk operator.
ii.3 Metropolis-Hastings Algorithm
The Metropolis-Hastings algorithm Metropolis et al. (1953); Hastings (1970) is a special class of Markov chains which obey detailed balance Eq. (1) by construction. The basic idea is to break the transition probability into two steps. In a first step, a transition from to is proposed with probability . This transition is accepted with probability , and if the transition is rejected the state remains . The overall transition probability is thus
(17) |
The detailed balance condition Eq. (1) becomes
(18) |
which in the Metropolis-Hastings algorithm is solved with the choice
(19) |
We note that our quantum algorithm can also be applied to the Glauber, or heat-bath choice
(20) |
The Metropolis-Hastings algorithm is widely used to generate a Boltzmann distribution as used in statistical physics and machine learning. Given a real energy function on the configuration space , the Boltzmann distribution at inverse temperature is defined as where the partition function ensures normalization. In this setting, it is common practice to choose a symmetrical proposed transition , so the acceptance probability depends only on the energy difference
(21) |
Iii Circuit for Walk operator
Quantum algorithms built from quantization of classical walks Ambainis (2004); Szegedy (2004); Somma et al. (2008); Magniez et al. (2011) usually assume an oracle formulation of the walk operator, where the ability to implement the transformation of Eq. (2) is taken for granted. As we discuss below in Appendix A, this transformation requires costly arithmetic operations. One of the key innovation of this article is to provide a detailed and simplified implementation of walk operator along with a detailed cost analysis for Metropolis-Hastings walks. As it will become apparent, our implementation circumvents the use of altogether.
For concreteness, we will assume a -local Ising model, where , and the energy function takes some simple form
(22) |
where are subsets of at most Ising spins and the are real coupling constants, and each spin interacts with at most other spins.
As it is always the case for Ising models, we will assume that the proposed transition of the Metropolis-Hastings walk are obtained by choosing a random set of spins and inverting their signs. In other words, where the product is taken bit by bit and where is some simple probability distribution on (it does not contain a trivial move), so is clearly symmetrical. The distribution is sparse, in the sense that it has only non-zero entries.
For concreteness, we will suppose that is uniform over some set of moves, with :
(23) |
The most common example consists of single-spin moves, where a single spin is chosen uniformly at random to be flipped. More generally, we will suppose that moves are sparse in the sense that each move flips a constant-bounded number of spins and that each spin belongs to a constant-bounded number of different moves . For , we use as a shorthand for where the are the elements of . With a further abuse of notation, we will view both as Ising spin configurations and as subsets of , where the correspondance is given by the locations of spins in .
A direct implementation of the oracle generally requires costly quantum circuits using arithmetic operations. The complexity arises from the need to uncompute a move register and a Boltzmann coin when implementing . This turns out to be non-trivial and costly if a move is rejected. We thus do not implement , but instead present a circuit which is isometric to the entire walk operator , thus avoiding the problem. In other words, we construct a circuit for where maps
(24) |
To minimize circuit depth, the second register above is encoded in a unary representation, so it contains qubits and is encoded as with a at the -th position. In addition to these two registers, the circuit acts on an additional coin qubit. Thus, we will denote the System, Move, and Coin registers with corresponding subscripts , and they contain , , and qubits respectively.
Our implementation of the walk operator combines 4 components:
(25) |
where
(26) | ||||
(27) | ||||
(28) | ||||
(29) |
While these definitions differ slightly from the ones of Sec. II.1, it can be verified straightforwardly, similar to our discussion in Sec. II.1, that these realize the desired walk operator. In what follows, we provide a complete description of each of these components, and their complexity is summarized in Table 1.
Gate | 3L depth | 3L count | Total depth | Qubits |
---|---|---|---|---|
iii.1 Move preparation
Recall that the Move register is encoded in unary. For a general distribution , the method of Rudolph and Grover (2002) can be adapted to realize the transformation Eq. (26). Here, we focus on the case of a uniform distribution. To begin, suppose that is a power of 2. Starting in the state , the state (in unary) is obtained by applying a sequence of gates in a binary-tree fashion. Recall indeed that . The gate is in the third level of the Clifford hierarchy, so it can be implemented exactly using a constant number of gates. This represents a substantial savings compared to the method of Rudolph and Grover (2002) for general distributions which requires arbitrary rotations obtained from costly gate synthesis.
To avoid such costly rotations, even when is not a power of 2, we choose to pad the distribution with additional states and prepare a distribution where . The states encode the moves of the classical walk , while the additional states correspond to trivial moves . This padding has the effect of slowing down the classical walk by a factor , and hence the quantum walk by a factor less than , which is less than the additional cost of preparing a uniform distribution over a range which is not a power of .
iii.2 Spin flip
The operator of Eq. (28) flips a set of system spins conditioned on the coin qubit and on the -th qubit of the move register being in state 1. This can be implemented with at most Toffoli gates (controlled-controlled-NOT), where the constant upper-bounds the number of spins that are flipped by a single move of . The coin register acts as one control for each gate, the -th bit of the move register acts as the other control, and the targets are the system register qubits that are in , for . No gate is applied to the padding qubits .
This implementation has the disadvantage of being purely sequential. An alternative implementation uses additional scratchpad qubits but is entirely parallel. The details of the implementation depends on the sparsity of the moves , and in general there is a tradeoff between the scratchpad size and the circuit depth. When the moves consist of single-spin flips for instance, this uses CNOTs in a binary-tree fashion (depth ) to make copies of the coin qubit. The Toffoli gates can then be applied in parallel for each move, and lastly the CNOTs are undone.
iii.3 Reflection
The transformation of Eq. (29) is a reflection about the state . Using standard phase kickback methods, it can be implemented with a single additional qubit in state and a control-NOT gate. The latter can be realized from serial Toffoli gates Barenco et al. (1995) and linear depth.
Since our goal is to minimize circuit depth, we use a different circuit layout that uses at most ancillary qubits and Toffoli gates to realize the -fold controlled-not. The circuit once again proceeds in a binary tree fashion, dividing the set of qubits into pairs and applying a Toffoli gate between every pair with a fresh ancilla in state 0 as the target. The ancillary qubit associated to a given pair is in state if and only if both qubits of the pair are in state 0. The procedure is repeated for the ancillary qubits, until a single bit indicates if all qubits are in state 0. The ancillary bits are then uncomputed.
iii.4 Boltzmann coin
The Boltzmann coin given in Eq. (27) is the most expensive component of the algorithm, simply because it is the only component which requires rotations by arbitrary angles. Specifically, conditioned on move qubit being 1 and the system register being in state , the coin register undergoes a rotation by an angle
(30) |
for Metropolis-Hastings or
(31) |
for Glauber dynamics, where . Given the sparsity constraints of the function and of the moves , the quantity can actually be evaluated from a subset of qubits of the system register, namely . For single-spin flips on a -local Hamiltonian, by definition. For multi-spin flips , we get .
Thus, the Boltzmann coin consists of a sequence of conditional gates , where itself is a single-qubit rotation by an angle determined by the qubits in the set . Since each are of constant-bounded sizes, each can be realized from a constant number of gates, so the entire Boltzmann coin requires gates, where is the desired accuracy.
Because all gates act on the Coin register, they must be applied sequentially. An alternative to this consists of copying the coin register in the conjugate basis , i.e. since a sequence of rotations is equivalent to a tensor product of these rotations under this mapping. Moreover, any set of gates with non-overlapping can be executed in parallel, so the total depth can be bounded by a constant at the expense of additional qubits.
The complexity of the Boltzmann coin does scale exponentially with the sparsity parameters of the model however, namely as . A circuit that achieves consists of a sequence of single-qubit rotations by an angle given by Eq. (30) or Eq. (31), conditioned on the bits in taking some fixed value. Each of these multi-controlled rotations require Toffoli gates along with gates, for an overall circuit depth of to realize .
Perhaps a more efficient way to realize the Boltzmann coin uses quantum signal processing methods Low et al. (2016); Low and Chuang (2016, 2017); Haah (2018). This is a method to constructs a unitary transformation from a controlled version of . In the current setting, and we choose where is given at Eq. (30). Applying a Hadamard to the Coin qubit, followed by a controlled with the Coin acting as control, and followed by a Hadamard on the Coin qubit again results in the transformation that we called above and that build up the Boltzmann coin transformation .
Above, the constant is chosen in such a way that the argument of the exponential is restricted to some finite interval which does not span the entire unit circle, say in the range . The exponential can be further decomposed as a product
(32) |
Each of these components consists of a rotation by an angle , whose sign is conditioned on the parity of the bits in . The parity bit can be computed using CNOTs, and the rotation is implemented using gate synthesis, with a -gate count per transformation of , which is dictated by the accuracy . The complexity of quantum signal processing depends on the targeted accuracy. More precisely, it scales with the number of Fourier coefficients required to approximate the function or to some constant accuracy on the domain .
Iv Heuristic use
The Metropolis-Hastings algorithm is widely used heuristically to solve minimization problems using simulated annealing or related algorithms Kirkpatrick et al. (1983). The objective function is the energy . Starting from a random configuration or an informed guess, the random walk is applied until some low-energy configuration is reached. The parameter can be varied in time, with an initial low value enabling large energy fluctuations to prevent the algorithm from getting trapped in local minimums, and large final value to reach a good (perhaps local) minimum.
In this section, we propose heuristic ways to use the quantum walk in the context of a minimization problem. We first recall the concept of total time to solution Ronnow et al. (2014) which we use to benchmark and compare different heuristics. We then present two quantum heuristics which we compare using numerical simulations on small instances.
Since the purpose of our study is to compare a classical walk to its different quantum incarnations – as opposed to optimizing a classical walk – we will use a schedule with a linearly increasing value of in time up to a fixed final value of in our comparison and expect our conclusions to hold if an optimized schedule was used instead in both the classical and the quantum walks.
iv.1 Total time to solution
When a random walk is used to minimize some function , the minimum is only reached with some finite probability . Starting from some distribution and applying the walk sequentially times, the success probability is . To boost this probability to some constant value , it is necessary to repeat the procedure times the number of repetitions , times. The total time to solution is then defined as the duration of the walk
(33) |
There is a compromise to be reached between the duration of the walk and the success probability – longer walks can reach a higher success probability and therefore be repeated fewer times, but increasing the duration of the walk beyond a certain point has a negligible impact on its success probability . We thus define the minimum total time to solution as min(TTS) .
iv.2 Zeno with rewind
In Sec. II.2, we explained how to prepare the eigenstate of with eigenvalue 1 using a sequence of walks . In the setting of Metropolis-Hastings where is the walk with parameter , a natural choice of are given by . An optimized schedule is also possible, but for a systematic comparison with the classical walk, we choose this fixed schedule, whose only parameter is the number of steps .
Let us revisit the argument of Sec. II.2 to establish some notation. Define the binary projective measurement . This binary measurement can be realized from uses of , where denotes the spectral gap of . Starting from the state , the Zeno algorithm consists in performing the sequence of binary measurements in increasing value of . The outcome on state yields state and occurs with probability . The sequence of measurements succeeds if they all yield this outcome, which occurs with probability and requires applications of quantum walk operator. For the algorithm to be successful, the final measurement in the computational basis must also yield the optimal outcome , which occurs with probability . Thus, the total time to solution for an the -step algorithm is
(34) |
In the method outlined above, a measurement outcome requires a complete restart of the algorithm to . There exists an alternative to a complete restart which we call rewind. It was first described in the context of Zeno state preparation in Ref. Lemieux and Poulin (2019), but originates from Refs. Marriott and Watrous (2005); Temme et al. (2011). It consists of iterating between the measurements and until the measurement is obtained. It can easily be shown that a transition between outcomes or is while the probability of a transition between outcomes or is . Given the cost of each of these measurements, we obtain a simple recursion relation for the expected cost of a successful transition with rewind, and thus for the total time to solution for a -step Zeno protocol with rewind. The minimal total time to solution is obtained by minimizing over . In Ref. Lemieux and Poulin (2019), it was found that rewinding yields substantial savings compared to the regular Zeno strategy for the preparation of quantum many-body ground states.
iv.3 Unitary implementation
We propose another heuristic use of the quantum walk which does not use measurement. Starting from state , it consists in applying the quantum walk operators sequentially, resulting in the state
(35) |
and ending with a computational basis measurement. The algorithm is successful if a computational basis measurement yields the outcome on state (rewind could be used otherwise), so the total time to solution for the -step algorithm is
(36) |
While we do not have a solid justification for this heuristic use, in Ref. Boixo et al. (2009), a protocol was proposed which used a similar sequence of unitaries, but where each unitary was applied a random number of times. The motivation for these randomized transformations was to phase randomize in the eigenbasis of the instantaneous unitary operator. When the spectral gap of a unitary operator is and that unitary is applied a random number of times in the interval , then the relative phase between the eigenstate with eigenvalue 1 and the other eigenstates is randomized over the unit circle, thus mimicking the effect of a measurement (but with an unknown outcome). From this analogy, we could expect that the unitary implementation yields a minimal total time to solution roughly equal to the Zeno-based algorithm with no rewind. But as we will see in the next section, its behavior is much better than anticipated – this method is more efficient than the Zeno algorithm with rewind, which itself is more efficient than Zeno without rewind.
iv.4 Numerical results
We have numerically benchmarked three heuristic algorithms: the classical walk with a variable-length linear interpolation between and and starting from a uniform distribution; the discrete, or Zeno-based adiabatic algorithm with rewind; and the unitary algorithm of the last subsection. The first system considered is a one dimensional Ising model. Figure 1 shows the quantum versus classical minimal total time to solution. The results clearly indicate a polynomial advantages of the quantum algorithms over the classical algorithm. Surprisingly, both quantum approaches show a similar improvement over the classical approach that exceed the expected quadratic speedup, with a power law fit of 0.42 using the unitary algorithm and 0.39 using the Zeno algorithm.
The second system considered is a sparse random Ising model: it has gaussian random couplings of variance 1, and the interactions sets (c.f. Eq. (32)) consist of a random subsets of of all the pairs of sites. Figure 2 shows quantum versus classical minimum total time to solution for a random ensemble of 100 systems of each sizes to 14. We observe that the unitary algorithm is consistently faster than the classical algorithm, with an average polynomial speedup of degree 0.75, less than the expected quadratic gain. Moreover, the different problem instances are all quite clustered around this average behavior, suggesting that the quantum speedup is fairly general and consistent. In contrast the Zeno algorithm shows large fluctuations about its average, particularly on very small problem instances. The average polynomial speedup is of degree 0.92, far worse than the unitary algorithm. Overall, the results indicate a polynomial advantage of the quantum methods over the classical method, but these advantages are much less pronounced than for the 1D Ising model.
In both the one-dimensional and the random graph Ising model, the unitary quantum algorithm achieves very similar and sometimes superior scaling to the Zeno with rewind algorithm. This is surprising given the observed improvement obtained from rewind in Ref. Lemieux and Poulin (2019) and our expectation that the unitary algorithm behaves essentially like Zeno without rewind.
V Discussion
In Appendix A, we discuss the difficulty of implementing the quantum walk oracle . Our conclusion, and perhaps one of the key messages of this Article, is that even though the quantum walk is traditionally defined with the help of a walk oracle, its circuit implementation does not necessarily require this walk oracle, and this can lead to substantial savings. Appendix B presents an improved parallelized heuristic classical walk for discrete sparse optimization problems which could potentially lead to significant improvements on a quantum computer. Unfortunately this walk is not reversible, which motivates further generalization of Szegedy’s quantization to include irreversible classical walks. In the rest of this section, we discuss the prospect of using the quantum walk to outperform a classical supercomputer.
We have proposed heuristic quantum algorithms based on the Szegedy walk for solving discrete optimization problems. Theoretical bounds show that the quantum algorithm can benefit from a quadratic speed-up () over its classical counterpart. Our numerical simulations on small problem instances indicate a super-quadratic speed-up () for the Ising chain and sub-quadratic speed-up () for random sparse Ising graphs. It remains an interesting question to understand more broadly what type of problems can benefit from what range of speed-up, and with these crude estimates in hand we can already look into the feasibility of a quantum speed-up on realistic devices.
To do so we will compare performances to the special-purpose supercomputer “Janus” Janus Collaboration (2009, 2012) which consists of a massive parallel field programable gate array (FPGA). This system is capable of performing Monte Carlo spin updates per second on a three-dimensional Ising spin glass of size . A calculation that lasts a bit less than a month will thus realize Monte Carlo steps. On the one hand, assuming that the theoretically predicted quadratic speed-up holds, the quantum computer must realize at least steps per month in order to keep up with the classical computer. This requires that a single step of the quantum walk be realized in a few milliseconds. On the other hand, the super-quadratic speed-up we have observed would allow almost a tenth of a second to realize a single quantum step, while the sub-quadratic speed-up would require that a single step be realized within 0.1 microseconds.
Taking the circuit depth reported in Table 1 as reference with for a three-dimensional lattice leads to a circuit depth of . To avoid harmful error accumulation, the gate synthesis accuracy should be chosen as the inverse volume of the quantum circuit, roughly , so on the order of logical gates are required per fine-tuned rotation Ross and Selinger (2016); Bocharov et al. (2015), for a total logical circuit depth of 200,000. With these estimates, the three scenarios described above require logical gate speeds ranging from an unrealistically short 0.5 picoseconds (sub-quadratic speed-up), to an extremely challenging 1 nanosecond (quadratic speed-up), and allow 0.5 microseconds for super-quadratic speed-up.
We could instead compile the rotations offline and teleport them in the computation Jones et al. (2012), which requires at least more qubits, but increase the time available for a logical gate by the same factor. Under this scenario, the time required for each logical gate would range from 0.1 nanoseconds (sub-quadratic speed-up), to 20 microseconds (quadratic speed-up), and to 1 miliseconds (super-quadratic speed-up). These estimates are summarized in Table 2. The latter is a realistic logical gate time for many qubit architectures, while there is no current path to achieve nanosecond logical gate times. It is thus important to understand the class of problems for which heuristic super-quadratic speed-ups can be achieved and to continue to optimize circuit implementations.
Quantum speedup | Synthesis online | Synthesis offline |
---|---|---|
Sub-quadratic | 0.5ps | 0.1ns |
Quadratic | 1ns | 20s |
Super-quadratic | 0.5s | 1ms |
Vi Acknowledgements
We thank Jeongwan Haah, Thomas Häner, Matt Hastings, and Guang Hao Low for stimulating discussions.
Appendix A Walk oracle
Our implementation of the walk operator does not make use of the walk oracle of Eq. (2). Since the transition matrix elements can be computed efficiently, we know that can be implemented in polynomial time. But this requires costly arithmetics which would yield a substantially larger complexity than the approach presented above.
To see how this complexity arises, consider the following implementation of , which uses much of the same elements as introduced above. The computer comprises two copies of the system register, which we now label Left and Right. As before, it also comprises a Move register and a Coin register. Begin with the Left register in state and all other registers in state . Use the transformation to prepare the state . Using CNOTs, copy the state of the Left register onto the Right register, resulting in . Apply the move proposed by the Move register to the Right register. If the Move register is encoded in unary representation as above, this requires CNOTs, and results in the state
(37) |
Using a version of the Boltzmann coin transformation on the Left, Right and Coin register yields
(38) | |||
(39) |
At this point, we swap the Left and Right registers conditioned on the Coin qubit being in state 1, resulting in the state
(40) |
Finally, reset the move register to using CNOTS with controls from the Left and Right registers. At this point, the move register is disentangled and discarded, resulting in the state
(41) |
This is quite similar to the state that would result from the quantum walk operator of Eq. (2), save for one detail. When the acceptance register is in state 0, the state of the right register needs to be mapped to the state . Such a rotation clearly depends on all the coefficients , and all implementations we could envision used arithmetic operations that compute .
Appendix B Irreversible parallel walk
Note that the Boltzmann operator has a total number of gates that scales with the system size , even though it is used to implement a single step of the quantum walk and that on average, a single spin is modified per step of the walk. This contrasts with the classical walk where in a single step of , a spin transition is chosen with probability , the acceptance probability is computed, and the move is either accepted or rejected. Each transition typically involves only a few spins (one in the setting we are currently considering), so implementing such a transition in the classical walk does not require an extensive number of gates. The complexity in that case is actually dominated by the generation of a pseudo-random number selecting the location of the spin to be flipped. As a consequence, the quantum algorithm suffers an -fold complexity increase compared to the quantum algorithm.
This motivates the construction of a modified classical walk for the lattice spin model which also affects every spin of the lattice, putting the classical and quantum walks on equal footing in terms of gate count. For simplicity, suppose that the set of moves consist in single-spin flips. We define a parallel classical walk with transition matrix
(42) |
where and is the transition which consist of flipping the th spin only, so only spin differ in and . The variable is a tunable parameter of the walk. In other words, a single step of this walk can be decomposed into a sequence over spins , and consists of flipping with probability and accepting the flip with probability . Importantly, even if the moves are applied sequentially, the acceptance probability is always evaluated relative to the state at the beginning of the step, even though other spins could have become flipped during the sequence.
If instead the acceptance probability was evaluated conditioned on the previously accepted moves – i.e. where is the total transition accumulated up to step – then this acceptance probability would be the same as used in the Metropolis-Hastings algorithm. Note that for a local spin model with, e.g., nearest-neighbor interactions, the two acceptance probabilities only differ if a neighbor of site has been flipped prior to attempting to flip spin . Because a transition on each spin is proposed with probability , the probability of having two neighboring spins flipped is . Thus, we essentially expect a single step of this modified walk to behave like steps of the original Metropolis-Hastings walk, with a systematic error that scales like . Moreover, this systematic error is expected to decrease over time since once the walk settles in a low-energy configuration, very few spin transitions will turn out to be accepted, thus further decreasing the probability of a neighboring pair of spin flips.
To verify the above expectation, we have performed numerical simulations on a Ising model
where were randomly chosen from . Results are shown on Fig. 3. What we observe is that, for an equal amount of computational resources, the parallelized walk outperforms the original walk. This is true both in terms in reaching a quick pseudo minimum configuration at short times and in terms of reaching the true minimum at longer times. Thus, while this parallelization was introduced to ease the quantization procedure, it appears to be of interest on its own.
In this case, the quantum walk oracle can easily be applied. We first proceed as in the previous subsection and use CNOTs to copy the Right register onto the Left register, yielding state . Then, sequentially over all spins , apply a rotation to spin of the Left register conditioned on the state of the spin and its neighbors on the Right register. This rotation transforms . Note that the function only depends on the bits of that are adjacent to site , so this rotation acts on a constant number of spins so requires a constant number of gates. Thus, the cost of the classical and the quantum parallel walks have the same scaling in . Combined to its observed advantages over the original classical walk, the parallel walk thus appears as the ideal version for a quantum implementation.
Unfortunately, the parallel walk is not reversible – it does not obey the detail-balance condition Eq. (1). Thus, it is not directly suitable to quantization à la Szegedy. While quantization of non-reversible walks were considered in Magniez et al. (2011), they require an implementation of time-reversed Markov chain defined from and its fixed point as
(43) |
Unfortunately, we do not know how to efficiently implement a quantum circuit for the time-reversed walk , so at present we are unable to quantize this parallel walk.
References
- Szegedy (2004) M. Szegedy, in Proceedings of the 45th Annual IEEE Symposium on Foundations of Computer Science (2004).
- Farhi et al. (2000) E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, “Quantum computation by adiabatic evolution,” (2000), quant-ph/0001106 .
- Aharonov and Ta-Shma (2003) D. Aharonov and A. Ta-Shma, Proc. 35th Annual ACM Symp. on Theo. Comp., 20 (2003).
- Boixo et al. (2010) S. Boixo, E. Knill, and R. D. Somma, “Fast quantum algorithms for traversing paths of eigenstates,” (2010), 1005.3034 .
- Somma et al. (2008) R. Somma, S. Boixo, H. Barnum, and E. Knill, Phys. Rev. Lett., 101, 130504 (2008), arXiv:0804.1571 .
- Metropolis et al. (1953) N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, J. Chem. Phys., 21, 1087 (1953).
- Hastings (1970) W. Hastings, Biometrika, 57, 97 (1970).
- Kitaev (1995) A. Kitaev, “Quantum measurements and the abelian stabilizer problem,” (1995), quant-ph/9511026 .
- Ambainis (2004) A. Ambainis, in Proc. 45th Annual IEEE Symp. on Foundations of Computer Science (FOCS) (2004) quant-ph/0311001 .
- Magniez et al. (2011) F. Magniez, A. Nayak, J. Roland, and M. Santha, SIAM journal on computing, 40, 142 (2011).
- Rudolph and Grover (2002) T. Rudolph and L. Grover, “A –rebit gate universal for quantum computing,” (2002), quant-ph/0210187 .
- Barenco et al. (1995) D. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, N. Margolus, P. Shor, T. Sleator, J. Smolin, and H. Weinfurter, Phys. Rev. A, 52, 3457 (1995), quant-ph/9503016 .
- Low et al. (2016) G. Low, T. Yoder, and I. L. Chuang, Phys. Rev. X, 6, 041067 (2016).
- Low and Chuang (2016) G. H. Low and I. L. Chuang, “Hamiltonian simulation by qubitization,” (2016), arXiv:1610.06546 .
- Low and Chuang (2017) G. H. Low and I. L. Chuang, Phys. Rev. Lett., 118, 010501 (2017).
- Haah (2018) J. Haah, “Product decomposition of periodic functions in quantum signal processing,” (2018), arXiv:1806.10236 .
- Kirkpatrick et al. (1983) S. Kirkpatrick, C. J. Gelatt, and M. Vecchi, Science, 220, 4598 (1983).
- Ronnow et al. (2014) T. Ronnow, S. Wang, J. Job, S. Boixo, S. Isakov, D. Wecker, J. Martinis, D. Lidar, and M. Troyer, Science, 345, 420 (2014), arXiv:1401.2910 .
- Lemieux and Poulin (2019) J. Lemieux and D. Poulin, “Resource estimate for quantum many-body ground state preparation on a quantum computer,” In preparation (2019).
- Marriott and Watrous (2005) C. Marriott and J. Watrous, Computational Complexity, 14, 122 (2005), cs/0506068 .
- Temme et al. (2011) K. Temme, T. Osborne, K. Vollbrecht, D. Poulin, and F. Verstraete, Nature, 471, 87 (2011), arXiv:0911.3635 .
- Boixo et al. (2009) S. Boixo, E. Knill, and R. Somma, Quant. Info. and Comp., 9, 833 (2009).
- Janus Collaboration (2009) Janus Collaboration, Comp. Sci. Eng., 11, 48 (2009).
- Janus Collaboration (2012) Janus Collaboration, Eur. J. Phys., 210, 33 (2012).
- Ross and Selinger (2016) N. Ross and P. Selinger, Quant. Info. and Comp., 16, 901 (2016), arXiv:1403.2975 .
- Bocharov et al. (2015) A. Bocharov, M. Rotteler, and K. Svore, Phys. Rev. A, 91 (2015).
- Jones et al. (2012) N. C. Jones, J. D. Whitfield, P. L. McMahon, M.-H. Yung, R. Van Meter, A. Aspuru-Guzik, and Y. Yamamoto, New Journal of Physics, 14, 115023 (2012).