Hybrid classical-quantum linear solver using Noisy Intermediate-Scale Quantum machines.

We propose a realistic hybrid classical-quantum linear solver to solve systems of linear equations of a specific type, and demonstrate its feasibility with Qiskit on IBM Q systems. This algorithm makes use of quantum random walk that runs in [Formula: see text](N log(N)) time on a quantum circuit made of [Formula: see text](log(N)) qubits. The input and output are classical data, and so can be easily accessed. It is robust against noise, and ready for implementation in applications such as machine learning.

Introduction

Algorithms that run on quantum computers hold promise to perform important computational tasks more efficiently than what can ever be achieved on classical computers, most notably Grover’s search algorithm and Shor’s integer factorization[1]. One computational task indispensable for many problems in science, engineering, mathematics, finance, and machine learning, is solving systems of linear equations . Classical direct and iterative algorithms take and time[2,3]. Interestingly, the Harrow-Hassidim-Lloyd (HHL) quantum algorithm[4-13], which is based on the quantum circuit model[14], takes only to solve a sparse system of linear equations, while for dense systems it requires [11]. Linear solvers and experimental realizations that use quantum annealing and adiabatic quantum computing machines[15-17] are also reported[18-20]. Most recently, methods[21,22] inspired by adiabatic quantum computing are proposed to be implemented on circuit-based quantum computers. Whether substantial quantum speedup exists in these algorithms remains unknown.

In practice, the applicability of quantum algorithms to classical systems are limited by the short coherence time of noisy quantum hardware in the so-called Noisy Intermediate-Scale Quantum (NISQ) era[23] and the difficulty in executing the input and output of classical data. Other roadblocks toward practical implementation include limited number of qubits, limited connectivity between qubits, and large error correction overhead. At present, experiments demonstrating the HHL linear solver on circuit quantum computers are limited to matrices[24-29], while linear solvers inspired by adiabatic quantum computing are limited to matrices[21,22]. For quantum annealers, the state-of-the-art linear solvers can solve up to matrices[20].

In addition to the problems of limited available entangled qubits and short coherence time, the HHL-type algorithms for the so-called Quantum Linear Systems Problem (QLSP) are designed to work only when input and output are quantum states[30]. This condition imposes severe restriction to practical applications in the NISQ era[23,30,31]. It has been shown that the HHL algorithm can not extract information about the norm of the solution vector [4]. A state preparation algorithm for inputting a classical vector would take time[30,32-34], with large overhead for current hardware. In addition, quantum state tomography is required to read out the classical solution vector , which is a demanding task[35,36], except for special cases like one-dimensional entangled qubits[37]. Inputting the matrix A is also a challenge that may kill the quantum speedup[1,24-29].

In this work, we propose a hybrid classical-quantum linear solver that uses circuit-based quantum computer to perform quantum random walks. In contrast to the HHL-type linear solvers, the solution vector and the constant vector in this hybrid algorithm stay as classical data in the classical registers. Only the matrix A is encoded in quantum registers. The idea is similar to that of variational quantum eigensolvers[38-41], where quantum speedup is exploited only for sampling exponentially large state Hilbert spaces, while the rest of computational task is done by classical computer. This makes it easy to perform data input and output: the vector can be arbitrary, and the components and the norm of the vector can be easily accessed.

We consider matrices that are useful for Markov decision problems such as in reinforcement learning[42]. We show that these matrices can be efficiently encoded by introducing the Hamming cube structure: a square matrix of size N requires quantum bits only. The quantum random walk algorithm we here propose takes time to obtain one component of the vector. We also show that in the quantum random walk algorithm the matrices produced as a result of qubit-qubit correlation are inherently complex, which can be an advantage for performing difficult tasks. For the same amount of time, the matrices the classical random walk algorithm can solve are limited to factorisable ones only.

We have tested the quantum random walk algorithm using software development kit Qiskit on IBM Q systems[43,44]. Numerical results show that this linear solver works on ideal quantum computer, and most importantly, also on noisy quantum computer having a short coherence time, provided the quantum circuit that encodes the A matrix is not too long. The limitation due to machine errors is discussed.

Results

We consider a system of linear equations of real numbers , where A is a matrix to be solved, vectors and are, respectively, the solution vector and a vector of constants. Without loss of generality, we rewrite A aswhere 1 is the identity matrix, and is a real number. We take P as a (stochastic) Markov-chain transition matrix, such that and , where refers to the P matrix element in the J-th column of the I-th row. This type of linear systems appears in value estimation for reinforcement learning[42,45,46], and radiosity equation in computer graphics[47]. In reinforcement learning algorithms, given a fixed policy of the learning agency, the vector is the value function that determines the long-term cumulative reward, and efficient estimation of this function is key to successful learning[42]. Note that the matrix A given in Eq. (1) used as model Hamiltonian matrix belongs to the so-called stoquastic Hamiltonians[48,49].

To solve , we expand the solution vector as Neumann series, that is, . Let us define the I0 component of truncated up to terms as

This expression for can be evaluated by random walks on a graph of N nodes, with the probability of going from node I and node J of the graph given by the matrix element , which we set as symmetric (undirected), namely . An example of a four-node graph is shown in Fig. 1(a). By performing a series of random walks starting from node I0, walking c steps according to the transition probability matrix P, and ending at some node I, Eq. (2) can be readily calculated to get the value, which is close to the solution for some large c steps. Truncating the series introduces an error . So, for a given , the number of steps necessary to meet a given tolerance is equal to .

Figure 1

(a) Quantum (or classical) random walk on an undirected graph. The transition probability of going from node to node or vice versa is equal to , these elements forming a matrix. (b) The four nodes on this Hamming cube are labeled by integers ; they are encoded as four different states , , , , respectively.

The above procedure can be extended to general matrices A by setting where for real matrix elements (see Methods 0.4). The calculation converges[50,51] provided that the spectral radius where the matrix B* is defined by . The matrices we here consider is a special case where is a constant, and this simplification guarantees convergence of the calculation.

For classical Monte Carlo methods to compute Eq. (2), it takes time to calculate the cumulative distribution function that is used to determine the next walking step. So, these linear systems can be solved by classical Monte Carlo methods within time[52-56]. Similar Monte Carlo methods have been extended to more general matrices for applications in Green’s function Monte Carlo method for many-body physics[57-59].

Encoding state spaces on Hamming cubes

As for material resources, in general it takes at least classical bits to store a row of a stochastic transition matrix P (or A). However, for the classical and quantum random walks we here consider, it is possible to reduce significantly the number of classical or quantum bits necessary to encode the corresponding transition probability matrix P to by introducing the Hamming cube (HC) structure[60]. To do it, we first associate each graph node with a bit string. As shown in Fig. 1(b), the four nodes of the graph are fully represented by two bits. Node states , , , and represent binary string states , , , and , respectively. For a N-node graph, only (to base 2) bits are needed to encode the integers , each representing the n-bit binary string state, namely , where is 0 or 1.

Classical random walk

Before we introduce our quantum random walk algorithm, let us first consider classical random walks.

To perform random walks on a N-node graph, we use a simple coin-flipping process with time steps. The -th bit flips with probability or does not flip with probability , the total probability being equal to 1. The transition probability matrix elements are given bywhere the n-bit binary string state is determined by , where denotes the bitwise exclusive or (XOR) operation, and the subscript c denotes classical states. The total number of , given by , is the Hamming weight of , and so corresponds to the Hamming distance between and states. This metric measures the number of steps that a walker needs to go from to on the Hamming cube.

For the four-node graph shown in Fig. 1, the transition probability matrix P for classical random walks readswhere denotes the Kronecker product. The lower triangular part of the matrix is omitted due to symmetry. This simple case demonstrates a general feature for classical transition probability matrix : the probability of flipping both bits is simply a product of the probabilities of flipping the 0-th bit and the 1-th bit in arbitrary order. For instance, ; similarly for the other ’s. The fact that can be factorized into a Kronecker product of the matrices of each individual bit indicates that each bit flips independently, as for a Markovian process.

Quantum random walk

We can simulate quantum walks[61-67] on a N-node graph to obtain the solution vector from Eq. (2). To do it, we use discrete-time coined quantum walk circuit[68,69]. The circuit for the four-node graph in Fig. 1 is shown in Fig. 2. The first two qubits and are state registers that will be initialized to encode the four-node graph, while the third qubit is the coin register.

Figure 2

Discrete-time coined quantum walk circuit for the transition matrix given in Eq. (10). Qubits and are state register qubits to represent the four-node graph in Fig. 1, first set as 0 before initialization, while the qubit is the coin register qubit. The measured registers and are fed back to initialize the next iteration. The classical-step is repeated c times to obtain the Neumann expansion up to order c.

To derive the quantum transition probability matrix on a graph of N nodes, we consider the state space of the -qubit circuit as spanned by with : the -th qubit registers the coin state , and the other n qubits encode the N-node graph. We take the convention that the rightmost bit is i0. Given a n-bit string , the initialized quantum state reads

Next we let the state evolve in random walk: in each walking step, we toss the coin by rotating the coin qubit, and then flip a graph qubit by applying the CNOT gate. This process is repeated on all the n qubits in the state, starting with the 0-th qubit. The corresponding evolution operator readswhere the prime (′) on the Π denotes that the operator applies first to the right, followed by the operator, and so on; the 1 operator is an identity map on the n-qubit state , is a Pauli X gate (the Pauli matrix ) that acts on the -th qubit, and is a single-qubit rotation operatorthat acts on the coin qubit state. Note that the first parentheses in Eq. (6) represents a CNOT gate. It is important to note that here we use one quantum coin only to decide on the Pauli X gate operation over all the n qubits, so the order of qubit operations plays a role in the determination of the transition probability matrix P.

The first step is to project on , which leads towith . By tracing out the coin degree of freedom, we obtain the reduced density matrix for the graph and hence the probability matrix . The resulting quantum transition probability matrix elements then readwhere is determined by . For one quantum evolution, the complex phase factors and play no role. We will see later that these phases come into play in the case of multiple evolutions .

To understand the transition probability matrix produced by the quantum walk circuit (Fig. 2), let us again consider the four-node graph in Fig. 1, where

Unlike the above classical random walk, this matrix cannot be factorized into a Kronecker product of the matrices of each individual qubit. The probability of one qubit flipping depends on the other, indicating that the two qubits are correlated, or in quantum information theory entangled.

In comparison to Eq. (3) obtained from the classical random walk, we see that additional XOR operations are required for classical computer to obtain the same quantum transition probability matrix, as can be seen from Eq. (9). In the case of , the classical and quantum transition probability matrices given by Eqs (4) and (10) are related by a permutation . The quantum version of the Hamming distance between and is given by , which clearly shows the temporal correlation between the -th and -th qubits. We attribute this correlation to the fact that only one quantum coin is used to decide on the Pauli X gate over all the n qubits, thus creating some connection between qubits, and to the non-Markovian nature of quantum walk dynamics[70,71], in which the quantum circuit memorizes the qubit state when it is walking in the direction that has the qubit state in the Hamming cube.

It can be of interest to note that the circuit given in Eq. (6) is just one possible design leading to a particular correlation between qubits. In general, there are numerous ways to rearrange the walking steps to obtain different kinds of correlation, and it is possible to design the circuit for specific purposes. A simple way is to perform the walking steps in Eq. (6) in a reverse order, operating the operator to the right first, followed by the operator, and so on. This leads to a different metric with . It turns out that this corresponds to the Hamming distance in the Gray code representation.

The Gray code uses single-distance coding for integer sequence , where adjacent integers differ by single bit flipping. In the case of the four-node graph in Fig. 1, the integers in the Gray code representation correspond to the , , , states, respectively. It is obvious that this Gray code representation can be obtained from the natural binary code representation by a permutation . There also exists a permutation that transforms to in the Gray code basis. The proof of this correspondence for arbitrary N is given in Methods 0.1. Both the transform and inverse transform between the natural binary code and Gray code representations take operations using classical computer[72]. This again shows that the quantum random walk algorithm gains improvement over the classical one.

As the change of the Hamming distance for each walking step in the Gray code representation is , a quantum walker in a geodesic of a Hamming cube automatically walks with the least action, that is, with the minimum change of the Hamming distance. This geodesic is a Hamiltonian path on hypercubes[73].

It is possible to increase the level of correlation in the probability matrix by performing multiple quantum evolutions, , where q is the number of quantum walk evolutions. The probability matrix produced by two quantum walk evolutions, , is given by (see Methods 0.2 for derivation)where, for and ,and

The fact that the summation over I in Eq. (11) runs over state configurations before the square is taken points to the complicated mixing of negative signs and complex phases ’s and ’s. The sign problem makes it difficult for pure classical Monte Carlo methods to simulate this transition.

In general, the dependence of the two-evolution quantum probability matrix on ’s, ’s and ’s, is not trivial. Its explicit expression for the graph is given in Methods 0.3. The phases ’s and ’s enter into play for graph sizes . On the other hand, the two-evolution probability matrix for classical random walk is given bywhich is still factorisable.

Numerical results

Figure 3 shows the performance of our hybrid quantum random walk algorithm on linear systems of dimension and . Their relative errors decrease with increasing sampling number. The relative error is defined as for the I-th component of the solution vector , where is the exact result obtained with the NumPy package. To demonstrate, we use randomly generated vectors and matrices A with a uniform distribution, and . We choose and c such that the error introduced by the Neumann expansion is within . See Table 1 for the relevant parameters of the two matrices. The program is written and compiled with Qiskit version 0.7.2. The simulation results (upper figure) are obtained using QASM simulator[43], while the quantum machine results (lower figure) are obtained using IBM Q 20 Tokyo device or Poughkeepsie device[74,75].

Figure 3

Relative errors as a function of the sampling number n for and matrices. The relevant parameters and estimated errors for these two matrices can be found in Table 1. Black solid lines represent the error reduction expected for Monte Carlo calculations. (Upper figure) Red dashed line and green dash-dotted line are the results computed by the QASM simulator. (Lower figure) Blue dash-dotted line and red dotted line are data for the same matrices computed by the IBM Q 20 Tokyo machine or Poughkeepsie machine. Cyan and magenta horizontal dashed lines depict the estimated errors.

Table 1

Relevant parameters for the matrices A of various sizes used for numerical experiments.

N	c	q	γ	Condition number	Estimated error ε0
64	6	2	0.3	1.457
128	6	2	0.3	1.599
256	6	1	0.3	1.857	0.1255
1024	10	1	0.5	2.973	0.2010

Estimated error is defined in the text.

The curves obtained by the QASM simulator are results averaged over ten runs. Their relative errors decrease as , where n is the number of random walk samplings. This reduction is typical of Monte Carlo simulations, because the hybrid quantum walk algorithm has essentially the same structure as classical Monte Carlo methods. So, we do not gain any speedup in sampling number. Yet, this result substantiates the fact that our proposed algorithm works on ideal quantum computers.

For real IBM Q quantum devices, the accuracy stops improving after a certain number of samplings (see the plateau (blue dash-dotted curve) and oscillation (red dotted curve) in Fig. 3). This hardware limitation can be estimated using an error formula , where is the condition number for the matrix A and E is the readout error of real machines. The condition number gauges the ratio of the relative error in the solution vector to the relative error in the A matrix[3]: some perturbation in the matrix, , can cause an error in the solution vector, , such that . By taking E as an estimate for , we obtain the above error for the solution vector as . The condition numbers given in Table 1 are computed by using Eq. (9) to construct the A matrices. For the average readout error of IBM Q 20 Tokyo device, we use [74]. The estimated errors are given in Table 1. We see that the relative errors fall below the respective errors, indicating that the precision limit is due to the readout error of the current NISQ hardware. Note that the machines are calibrated several times during data collection, so the hardware error varies and the E value is only an estimate.

Figure 4 shows the results for linear systems of dimension and , obtained by the QASM simulator that performs two quantum walk evolutions with uniformly distributed . The relevant parameters for these two matrices are given in Table 1. The results again evidence that the algorithm works well, even in the presence of complex phases ’s and ’s. Note that we here take (, , ) as random variables to demonstrate the efficiency of our algorithm, but in real applications, these variables must be provided by other algorithms to generate a proper P matrix Fig. 4.

Figure 4

Relative errors as a function of the sampling number n for and matrices, obtained by performing two quantum walk evolutions, . Black solid lines represent the error reduction expected for Monte Carlo calculations. Red dashed line and blue dotted line are the results computed by the QASM simulator.

The communication latency between classical and quantum computer is the most time-consuming part, containing communications. Fortunately, this number does not scale as N. For users with direct access to the quantum processors, communication bottleneck should be less severe.

Discussion

A comparison of computational resources is given in Table 2. For hybrid quantum walk algorithm, we need qubits, CNOT gates, and U3 gates, where q is the number of evolutions. The initialization takes log(N) X gates; but since they can be executed simultaneously, the initialization occupies one time slot only. Totally time slots are required for one quantum walk evolution to obtain one component of the solution vector . This can be an advantage when one is interested in partial information about . The same amount of time slots can be similarly derived for the classical random walk algorithm. Yet, we stress that these two algorithms deal with different transition probability matrices: factorisable matrices for classical random walk, and more complex correlated matrices for quantum random walk. The qubit-qubit correlation built into the correlated matrix can potentially be harnessed to perform complex tasks.

Table 2

Comparison of various algorithms for solving N × N linear systems , with respect to time and space complexities, and Input/Output issues.

Algorithm	Time	Space for A	Input/Output
Classical Direct[2,3]	O(N^3)	O(N^2)	efficient for any A,x→,b→
Classical Iterative[2,3]	O(N^2)	O(N^2)	efficient for any A,x→,b→
Quantum HHL[4]	O(log(N))	O(log(N)) qubits	norm ∥x→∥ not available difficult for A,x→,b→
Classical MC[45,53,55] (for one component xI)	O(N)	O(N)	efficient for any x→,b→ limited A (stochastic P)
Classical RW on HC (for one component xI)	O(log(N))	O(log(N))	efficient for any x→,b→ limited A (factorisable P)
Hybrid QW on HC (for one component xI)	O(log(N))	O(log(N)) qubits	efficient for any x→,b→ limited A (correlated P)

Note that for classical Monte Carlo (MC) method, classical random walk (RW) and hybrid quantum random walk (QW), the time complexities in the table are per sampling time. It takes samplings to achieve the desired accuracy (see the text).

Other advantages of the algorithms we propose are:

By restricting the matrices A to those that can be encoded in Hamming cubes, we can sample both classical and quantum random walk spaces that scale exponentially with the number of bits/qubits, and hence gain space complexity.

Classical Monte Carlo methods have time complexity of for general P matrices. For the matrices here considered, our algorithms have .

It is easier to access input and output than the HHL-type algorithm.

Random processes in a quantum computer are fundamental, and so are not plagued by various problems associated with pseudo-random number generators[76], like periods and unwanted correlations.

Our quantum algorithm can run on noisy quantum computers whose coherence time is short.

We propose a hybrid quantum algorithm suitable for NISQ quantum computers to solve systems of linear equations. The solution vector and constant vector we consider here are classical data, so the input and readout can be executed easily. Numerical simulations using IBM Q systems support the feasibility of this algorithm. We demonstrate that, by performing two quantum walk evolutions, the resulting probability matrix become more correlated in the parameter space. As long as the quantum circuit in this framework produces highly correlated probability matrix with a relatively short circuit depth, we can always gain quantum advantages over classical circuits.

Methods

Gray code basis

The natural binary code is transformed to the Gray code basis[72] according to

with . The probability matrix in the Gray code basis is given bywith .

Lemma 1

Let S be the set of all possible n-bit strings with , and be a permutation of the set S. If there exists a function such that for ,

, and if is bitwise XOR homomorphic, then we have .

Proof 1 Since is bitwise XOR homomorphic, Eq. (A.3) leads to.

Lemma 2

Let be represented by . Let be a function that transforms from natural bit string to Gray code according to , with . Then g is a bitwise XOR homomorphism.

Proof 2 Let be represented by bit strings and , respectively. Usingwith , we get.

Using Lemma 1 and Lemma 2, the following theorem is clear.

Theorem 1

There exists a permutation that maps the probability matrix produced by classical random walk to the probability matrix given in Eq. (A.2) produced by the quantum random walk circuit in a reverse order, that is, in Gray code basis.

Derivation of Eq. (11)

We use the evolution operator given in Eq. (6),to compute the two-evolution operator

Next we project the operator on the state,where is given in Eq. (12) and

We then project on the final state which leads to the probability matrix elements as

Two-evolution quantum walk on N = 4 graph

The probability matrix elements for two quantum evolutions on the four-node graph read

Surprisingly, in this case the matrix elements do not depend on the and phases. However, the matrix elements do depend on complex phases when , as can be numerically checked. Note that depend on only: the destructive interference between configurations totally eliminates the dependence, which is difficult to do by simple classical random walks.

Solving for general matrices

Here we discuss the applicability of our quantum random walk algorithm to general matrices[50,51,77,78]. Given an arbitrary matrix A, we can obtain and . Then the linear system can be solved by performing random walks according to the transition probabilities and by multiplying the factor at each walking step, provided that the linear solver converges to a solution. In classical random walk algorithms, it has been shown[50] that the convergence of the linear solver depends on the spectral radius of the matrix B* where , that is, the necessary and sufficient condition for convergence is . We expect a similar condition for quantum random walk algorithms. However, one should consider the hybrid solver presented in this work as a special-purpose solver, in which the quantum circuit is designed for a specific matrix problem. The quantum circuits demonstrated in this work show that there are probability transition matrices that are easy to sample using quantum circuits but difficult using classical circuits. How to tailor a circuit design along with the relevant parameters suitable for the kind of application we are looking for is beyond the scope of this work.