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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{O}}$$\end{document}O(N log(N)) time on a quantum circuit made of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{O}}$$\end{document}O(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.

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 o(N log(N)) time on a quantum circuit made of o(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.
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 → = → x b A . Classical direct and iterative algorithms take  N ( ) 3 and  N ( ) 2 time 2,3 . Interestingly, the Harrow-Hassidim-Lloyd (HHL) quantum algorithm [4][5][6][7][8][9][10][11][12][13] , which is based on the quantum circuit model 14 , takes only  N (log( )) to solve a sparse × N N system of linear equations, while for dense systems it requires N N ( log( ))  11 . Linear solvers and experimental realizations that use quantum annealing and adiabatic quantum computing machines [15][16][17] are also reported [18][19][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 × 2 2 matrices [24][25][26][27][28][29] , while linear solvers inspired by adiabatic quantum computing are limited to × 8 8 matrices 21,22 . For quantum annealers, the state-of-the-art linear solvers can solve up to 12 12 × 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 x → 4 . A state preparation algorithm for inputting a classical vector → b would take N ( )  time 30,[32][33][34] , with large overhead for current hardware. In addition, quantum state tomography is required to read out the classical solution vector x → , which is a demanding task 35,36 , except for special cases like one-dimensional entangled qubits 37 .
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 x → and the constant vector → b 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][39][40][41] , where quantum speedup is exploited only for sampling exponentially large state Hilbert spaces, while the rest of computational 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 N (log( ))  quantum bits only. The quantum random walk algorithm we here propose takes  N (log( )) time to obtain one component of the x → 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.

We consider a system of linear equations of real numbers
, where A is a N N × matrix to be solved, × N 1 vectors x → and b → are, respectively, the solution vector and a vector of constants. Without loss of generality, we rewrite A as where 1 is the identity matrix, and 0 1 γ < < is a real number. We take P as a (stochastic) Markov-chain transition matrix, such that P 0 , where P I J , 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 → x 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 , we expand the solution vector as Neumann series, that is,  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 P I J , , which we set as symmetric (undirected), namely = P P I J J I , , . An example of a four-node graph is shown in Fig. 1(a). By performing a series of random walks starting from node I 0 , walking c steps according to the transition probability matrix P, and ending at some node I c , Eq. (2) can be readily calculated to get the x I c ( ) 0 value, which is close to the solution x I 0 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 c log(1/ )/ log(1/ ) ε γ ∼ . The above procedure can be extended to general matrices A by setting  www.nature.com/scientificreports www.nature.com/scientificreports/ For classical Monte Carlo methods to compute Eq. (2), it takes N ( )  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 N ( ) 2  time [52][53][54][55][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 N ( )  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 N (log( ))  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 N 4 = graph are fully represented by two bits. Node states 0 | 〉, 1 | 〉, | 〉 2 , and | 〉 3 represent binary string states | 〉 00 , | 〉 01 , 10 | 〉, and | 〉 11 , respectively. For a N-node graph, only N n log ( ) where ⊗ 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 P classical : 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, Quantum random walk. We can simulate quantum walks 61-67 on a N-node graph to obtain the solution vector → x 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 j 0 and j 1 are state registers that will be initialized to encode the four-node graph, while the third qubit j 2 is the coin register.
To derive the quantum transition probability matrix on a graph of N nodes, we consider the state space of the : the n ( 1) + -th qubit registers the coin state i n | 〉 ◊ , and the other n qubits encode the N-node graph. We take the convention that the rightmost bit is i 0 . Given a n-bit string … www.nature.com/scientificreports www.nature.com/scientificreports/ where the prime (′) on the Π denotes that the = k 0 operator applies first to the right, followed by the = k 1 operator, and so on; the 1 q operator is an identity map on the n-qubit state | 〉 J q , X k is a Pauli X gate (the Pauli matrix σ x ) that acts on the k-th qubit, and U u ( ) 3 is a single-qubit rotation operator that 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 J 0, ψ | 〉, which leads to By tracing out the coin degree of freedom, we obtain the reduced density matrix for the graph and hence the probability matrix P J J The resulting quantum transition probability matrix elements then read For one  quantum evolution, the complex phase factors e iφ  and λ  e i play no role. We will see later that these phases come into play in the case of multiple evolutions q  . 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  www.nature.com/scientificreports www.nature.com/scientificreports/ 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  N (log( )) 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 = N 4, the classical and quantum transition probability matrices given by Eqs (4) and (10)  -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  | 〉 − i 1 when it is walking in the direction that has the qubit state i | 〉  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 = − k n 1 operator to the right first, followed by the k n 2 = − operator, and so on. This leads to a different metric It turns out that this d quantum 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 (0, 1, 2, 3) in the Gray code representation correspond to the 00 | 〉, | 〉 01 , 11 | 〉, 10 | 〉 states, respectively. It is obvious that this Gray code representation can be obtained from the natural binary code representation by a permutation ( ) 0 1 2 3 0 1 3 2 . There also exists a permutation that transforms P classical to P quantum 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 N (log( ))  operations using classical computer 72 . This again shows that the quantum random walk algorithm gains N (log( ))  improvement over the classical one. As the change of the Hamming distance for each walking step in the Gray code representation is d 1 δ = , 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, q  , where q is the number of quantum walk evolutions. The probability matrix produced by two quantum walk evolutions, 2  , is given by (see Methods 0.2 for derivation)

( )
The fact that the summation over I in Eq. (11) runs over (2 ) n  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 N 4 = graph is given in Methods 0.3. The phases φ  's and λ  's enter into play for graph sizes ≥ N 8. On the other hand, the two-evolution probability matrix for classical random walk is given by which is still factorisable.
Numerical results. Figure 3 shows the performance of our hybrid quantum random walk algorithm on lin-  74,75 . The curves obtained by the QASM simulator are results averaged over ten runs. Their relative errors decrease as n 1/ s , where n s is the number of random walk samplings. This n 1/ s 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 E r 0 ε κ ∼ × , where κ is the condition number for the matrix A and E r is the  www.nature.com/scientificreports www.nature.com/scientificreports/ readout error of real machines. The condition number κ gauges the ratio of the relative error in the solution vector → x to the relative error in the A matrix 3 : some perturbation in the matrix, A A δ + , can cause an error in the solution vector, . By taking E r as an estimate for A δ , 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 = . × − E 6 76 10 r 274 . The estimated errors ε 0 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 r value is only an estimate. Figure 4 shows the results for linear systems of dimension = N 64 and N 128 = , 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 (  θ , l φ , λ l ) 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.
The communication latency between classical and quantum computer is the most time-consuming part, containing cn ( ) s  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 N 1 log( ) + qubits, q N log( ) CNOT gates, and q N log( ) U 3 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 q N 1 2 log( ) + time slots are required for one quantum walk evolution to obtain one component of the solution vector x → . This can be an advantage when one is interested in partial information about x → . 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. Other advantages of the algorithms we propose are: (i) 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. (ii) Classical Monte Carlo methods have time complexity of N ( )  for general P matrices. For the matrices here considered, our algorithms have N (log( ))  . (iii) It is easier to access input and output than the HHL-type algorithm. (iv) 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. (v) 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 x → and constant vector b → we consider here are classical data, so the input and www.nature.com/scientificreports www.nature.com/scientificreports/ 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.

Gray code basis. The natural binary code
n n 1 2 1 0 is transformed to the Gray code basis 72 according to The probability matrix in the Gray code basis is given by  ∀ ∈ … − i n 0, , 1 . 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. Eq. (11). We use the evolution operator given in Eq. (6),  Surprisingly, in this case the matrix elements do not depend on the ( , ) 0 1 φ φ and λ λ ( , ) 0 1 phases. However, the matrix elements do depend on complex phases when N 8 ≥ , as can be numerically checked. Note that P P P P ( , , , ) 01 02 23 13 depend on θ 1 only: the destructive interference between configurations totally eliminates the θ 0 dependence, which is difficult to do by simple classical random walks.

Derivation of
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 = − B 1 A and B P v I J I J I J , , , = . Then the linear system x b A → = → can be solved by performing random walks according to the P I J , transition probabilities and by multiplying the factor v I J , 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 , that is, the necessary and sufficient condition for convergence is B ( ) 1 ⁎ ρ < . 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.