Hybrid classical-quantum linear solver using Noisy Intermediate-Scale Quantum machines
Chih-Chieh Chen Shiue-Yuan Shiau Ming-Feng Wu Yuh-Renn Wu























[]itriA70068@itri.org.tw
Electronic and Optoelectronic System Research Laboratories, Industrial Technology Research Institute, Hsinchu 31057, Taiwan



Physics Division, National Center for Theoretical Sciences, Hsinchu 30013, Taiwan


Electronic and Optoelectronic System Research Laboratories, Industrial Technology Research Institute, Hsinchu 31057, Taiwan







Graduate Institute of Photonics and Optoelectronics and Department of Electrical Engineering, National Taiwan University, Taipei 10617, Taiwan
Electronic and Optoelectronic System Research Laboratories, Industrial Technology Research Institute, Hsinchu 31057, Taiwan















 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 𝒪(Nlog(N)) time on a quantum circuit made of 𝒪(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.  








Hybrid classical-quantum linear solver using Noisy Intermediate-Scale Quantum machines
    Yuh-Renn Wu
    December 30, 2023
======================================================================================














§ 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 <cit.>.  One   computational task indispensable for many problems in science, engineering, mathematics, finance, and machine learning, is solving systems of linear equations Ax⃗=b⃗. Classical direct and iterative algorithms take 𝒪(N^3) and 𝒪(N^2) time <cit.>. Interestingly, the Harrow-Hassidim-Lloyd (HHL) quantum algorithm <cit.>, which is based on the quantum circuit model <cit.>, takes only  𝒪(log(N))  to solve a sparse N× N system of linear equations, while for dense systems it requires 𝒪(√(N)log(N)) <cit.>. Linear solvers and experimental realizations that use  quantum annealing and adiabatic quantum computing machines <cit.> are also reported <cit.>. Most recently, methods <cit.> 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 <cit.> 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 <cit.>,  while  linear solvers   inspired by adiabatic quantum computing are limited to 8× 8 matrices <cit.>. For quantum annealers, the state-of-the-art linear solvers can solve up to 12× 12 matrices <cit.>.  





In addition to the problems of limited available entangled qubits and short coherence time, the HHL-type algorithms are designed to work  only when  input and output are  quantum states <cit.>. This condition imposes severe restriction to practical applications in the NISQ era <cit.>. It has been shown that the HHL algorithm can not extract information about the norm of the solution vector x⃗ <cit.>. A state preparation algorithm for inputting a classical vector b⃗ would take 𝒪(N) time <cit.>, 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 <cit.>, except  for special cases like  one-dimensional entangled qubits <cit.>. Inputting the matrix A is also a challenge that may kill the quantum speedup <cit.>.


 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 <cit.>, 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 b⃗ vector can be arbitrary, and the components and the norm of  the x⃗ vector can be easily accessed. 

We consider matrices  that are useful for Markov decision problems such as in reinforcement learning <cit.>. We show that  these matrices can be efficiently encoded by introducing the Hamming cube structure: a square matrix of size N requires 𝒪(log(N)) quantum bits only. The quantum random walk algorithm we here propose takes  𝒪(log(N)) 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 <cit.>. 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.



§ METHODS




We consider a system of linear equations of real numbers Ax⃗=b⃗, 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 

    A=1-γP ,

 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_I,J≥ 0 and ∑_J P_I,J=1, 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 <cit.>, and radiosity equation in computer graphics <cit.>. 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 <cit.>.  Note that the matrix A given in Eq. (<ref>) used as model Hamiltonian matrix  belongs to the so-called  stoquastic Hamiltonians <cit.>.   





To solve Ax⃗=b⃗, we expand the solution vector as Neumann series, that is, x⃗=A^-1b⃗=(1-γP)^-1b⃗=∑_s=0^∞γ^s P^s b⃗. Let us define the I_0 component of x⃗ truncated up to γ^c terms as 

    x^(c)_I_0=∑_s=0^c γ^s ∑_I_1,...,I_s=1^N P_I_0,I_1...P_I_s-1,I_s b_I_s .

 This expression for x^(c)_I_0  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_I,J=P_J,I. An example of a four-node graph is shown in Fig. <ref>(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. (<ref>) can be readily calculated to get the x^(c)_I_0 value, which is close to the solution x_I_0 for some large c steps.  Truncating the series introduces an error ϵ∼𝒪(γ^c). So, for a given γ, the number of steps necessary to meet a given tolerance ϵ is equal to  c∼log (1/ϵ)/log (1/γ).
 
 
The above expansion procedure can be extended to more general matrices A by setting A=1-B where B_I,J=P_I,Jv_I,J for  real matrix elements v_I,J provided that the eigenvalues of B are bounded by -1<(B)_I<1.  
 
 
For classical Monte Carlo methods to compute Eq. (<ref>), 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 <cit.>. Similar Monte Carlo  methods have been extended to more general matrices for applications in Green's function Monte Carlo method for many-body physics <cit.>.

 





 §.§ 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 𝒪(log(N)) by introducing the Hamming cube (HC) structure <cit.>.
To do it, we first associate each graph node with a bit string. As shown in Fig. <ref>(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 log_2(N)= n (to base 2)  bits are needed to encode the integers J∈{0,1,..., N-1}, each representing the n-bit binary string state, namely  |J⟩=|j_n-1,...j_1,j_0⟩, where j_ℓ 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 𝒪(log(N)) time steps. The ℓ-th bit flips with probability sin^2(θ_ℓ/2) or does not flip with probability cos^2(θ_ℓ/2), the total probability being equal to 1. The transition probability matrix elements are given by

    P^classical_J',J=∏_ℓ=0^n-1|cos^2 (θ_ℓ/2)|^1-i_ℓ|sin^2 (θ_ℓ/2)|^i_ℓ ,

where the  n-bit binary string state |I⟩_c=|i_n-1,...,i_1,i_0⟩_c  is determined by |J'⟩_c= |I⟩_c ⊕ |J⟩_c, where ⊕ denotes the bitwise exclusive or (XOR) operation, and the subscript c denotes classical states. The total number of |sin^2(θ_ℓ/2)|, given by d^classical=∑_ℓ=0^n-1i_ℓ, is the Hamming weight of |I⟩_c, and so corresponds to the Hamming distance between |J'⟩_c and |J⟩_c states. This metric measures the number of steps that a walker needs to go from |J⟩_c  to |J'⟩_c  on the Hamming cube.


 For the four-node graph shown in Fig. <ref>, the transition probability matrix P for classical random walks reads


    P^classical   =   [ cos^2(θ_0/2 )cos^2(θ_1/2) sin^2(θ_0/2 )cos^2(θ_1/2) cos^2(θ_0/2 )sin^2(θ_1/2) sin^2(θ_0/2 )sin^2(θ_1/2);                           cos^2(θ_0/2 )cos^2(θ_1/2) sin^2(θ_0/2 )sin^2(θ_1/2) cos^2(θ_0/2 )sin^2(θ_1/2);                           cos^2(θ_0/2 )cos^2(θ_1/2) sin^2(θ_0/2 )cos^2(θ_1/2);                           cos^2(θ_0/2 )cos^2(θ_1/2);                           ]
       =   [ cos^2(θ_1/2 ) sin^2(θ_1/2 ); sin^2(θ_1/2 ) cos^2(θ_1/2 ) ]⊗[ cos^2(θ_0/2 ) sin^2(θ_0/2 ); sin^2(θ_0/2 ) cos^2(θ_0/2 ) ] ,


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, P^classical_0,3=P^classical_|00⟩,|11⟩=sin^2(θ_0/2 )sin^2(θ_1/2)=P^classical_|00⟩,|01⟩P^classical_|00⟩,|10⟩; similarly for the other P^classical_I,J's. The fact that  P^classical  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 <cit.> on a N-node graph to obtain the solution vector x⃗ from Eq. (<ref>).  To do it, we use discrete-time coined quantum walk circuit <cit.>. The circuit for the four-node graph in Fig. <ref> is  shown in Fig. <ref>.  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 (n+1)-qubit circuit as spanned by {|i_n⟩_♢⊗ |i_n-1,...,i_1,i_0 ⟩_q} with n=log_2(N): 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  (j_n-1,...,j_1,j_0), the initialized quantum state reads 

    |ψ_0,J⟩   =    |0 ⟩_♢⊗ |j_n-1,j_n-2,...,j_2,j_1,j_0 ⟩_q
       =    |0 ⟩_♢⊗|J⟩_q .



Next we let the |ψ_0,J⟩ 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 |j_n-1,j_n-2,...,j_2,j_1,j_0 ⟩_q state, starting with the 0-th qubit.
The corresponding evolution operator reads 

    𝒰=∏_k=0^n-1 ' (|0⟩_♢⟨ 0|_♢⊗ 1_q + |1⟩_♢⟨ 1|_♢⊗ X_k) ·(U_3( u_k)⊗ 1_q) ,

 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_3( u) is a single-qubit rotation operator

    U_3( u)=U_3(θ,ϕ,λ)= 
    [         cos(θ/2)    -e^iλsin(θ/2);     e^iϕsin(θ/2) e^i(λ+ϕ)cos(θ/2);                  ]

that acts on the coin qubit state.  Note that the first parentheses in Eq. (<ref>) 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 |ψ_0,J⟩, which leads to

    𝒰 |ψ_0,J⟩= ∑^1_i_n-1,...,i_0=0∏_ℓ=0^n-1 U_3( u_ℓ)_i_ℓ, i_ℓ-1
       ×   |i_n-1⟩_♢⊗ |i_n-1⊕ j_n-1,...,i_1⊕ j_1,i_0⊕ j_0⟩_q ,

with i_-1= 0. 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 =  ⟨ J'| _♢ [𝒰 |ψ_0J⟩⟨ψ_0J | 𝒰^† ]  |J'⟩. The resulting quantum transition probability matrix  elements then read

    P^quantum_J',J   =   ∏_ℓ=0^n-1|U_3( u_ℓ)_i_ℓ, i_ℓ-1|^2
       =   ∏_ℓ=0^n-1|cos^2 (θ_ℓ/2)|^1-(i_ℓ⊕ i_ℓ-1)|sin^2 (θ_ℓ/2)|^i_ℓ⊕ i_ℓ-1 ,

where |I⟩_q=|i_n-1,...,i_1,i_0 ⟩_q  is determined by |J'⟩_q= |I⟩_q ⊕ |J⟩_q.  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. <ref>), let us again consider the four-node graph in  Fig. <ref>, where 


    P^quantum=[ cos^2(θ_0/2 )cos^2(θ_1/2) sin^2(θ_0/2 )sin^2(θ_1/2) cos^2(θ_0/2 )sin^2(θ_1/2) sin^2(θ_0/2 )cos^2(θ_1/2);                           cos^2(θ_0/2 )cos^2(θ_1/2) sin^2(θ_0/2 )cos^2(θ_1/2) cos^2(θ_0/2 )sin^2(θ_1/2);                           cos^2(θ_0/2 )cos^2(θ_1/2) sin^2(θ_0/2 )sin^2(θ_1/2);                           cos^2(θ_0/2 )cos^2(θ_1/2);                           ].


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. (<ref>) obtained from the classical random walk, we see that  additional 𝒪(log(N)) XOR operations are required for classical computer to obtain the same quantum transition probability matrix, as can be seen from Eq. (<ref>). In the case of N=4, the classical and quantum transition probability  matrices given by Eqs. (<ref>) and (<ref>) are related by a permutation 
[ 0 1 2 3; 0 3 2 1 ]. The quantum  version of the Hamming distance  between |J⟩_q and |J'⟩_q is  given by d^quantum=∑_ℓ=0^n-1 i_ℓ⊕ i_ℓ-1, which clearly shows the temporal correlation between the  ℓ-th and (ℓ-1)-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 <cit.>, 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. (<ref>) 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. (<ref>) 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 d^quantum=∑_ℓ=0^n-1 i_ℓ+1⊕ i_ℓ with i_n = 0. 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 0 → 1 →⋯→ N-1, where adjacent integers differ  by single bit flipping.  In the case of the four-node graph  in Fig. <ref>, 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  Appendix <ref>. Both the transform and inverse transform between the natural binary code and Gray code representations take 𝒪(log(N)) operations using classical computer <cit.>. This again shows that the quantum random walk algorithm gains  𝒪(log(N)) 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 <cit.>. 
 
 

 



















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 Appendix <ref> for derivation)

    P_J',J^quantum   =   ∑_k=0^1 |∑_I   f(I, J'⊕ J ⊕ I) δ_i_n-1,k|^2  ,

where, for I = (i_n-1,...,i_0) and K = (k_n-1,...,k_0), 

    f(I, K)   =    [U_3( u_n-1)]_i_n-1,i_n-2⋯ [U_3( u_0)]_i_0,k_n-1
          × [U_3( u_n-1)]_k_n-1,k_n-2⋯[U_3( u_0)]_k_0,0  ,

and

    [U_3(θ,ϕ,λ)]_μ,ν
          =  e^i[μϕ + νλ] (-1)^ (1-μ) ν(cos(θ/2))^1-(μ⊕ν)(sin(θ/2))^μ⊕ν  .

The fact that the summation over  I in Eq. (<ref>) 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 Appendix <ref>. 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

    P_J',J^classical = 
          ∏_ℓ=0^n-1|cos^4(θ_ℓ/2)+sin^4(θ_ℓ/2)|^1-i_ℓ|2cos^2(θ_ℓ/2) sin^2(θ_ℓ/2) |^i_ℓ ,

which is still factorisable.











§ NUMERICAL RESULTS


Figure <ref> shows  the performance of our hybrid quantum random walk algorithm on  linear systems  of dimension N=256 and N=1024. Their relative errors decrease with increasing sampling number. The relative error is defined as ϵ=|x^exact_I-x_I|/|x^exact_I| for the I-th component of the solution vector x⃗, where x⃗^exact is the exact result obtained with the NumPy package. To demonstrate, we use randomly generated   vectors b⃗ and matrices A with a uniform distribution, b_I ∈ [-1,1] and  θ_ℓ∈ [0,π].  We choose γ and c such that the error introduced by the Neumann expansion is within 𝒪(10^-4). See Table <ref> 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 <cit.>, while the quantum machine results (lower figure) are obtained using IBM Q 20 Tokyo device or Poughkeepsie device <cit.>. 



The  curves obtained by the QASM simulator are results averaged over ten runs. Their relative errors decrease as 1/√(n_s), where  n_s is the number of random walk samplings. This 1/√(n_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. <ref>). This hardware limitation can be estimated using an error  formula ϵ_0 ∼κ× E_r, where κ is the condition number for the matrix A and E_r is the 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 <cit.>: some perturbation in the matrix, A+δA, can cause an error in the solution vector, x⃗+δx⃗, such that ||δx⃗||∼κ× ||δA||.  By taking E_r as an estimate for ||δA||, we obtain the above error  for the solution vector as ϵ_0 = ||δx⃗||∼κ× E_r.
The condition numbers given  in Table <ref> are computed by using Eq. (<ref>) to construct the A matrices. For the  average readout error of IBM Q 20 Tokyo device, we use E_r=6.76× 10^-2 <cit.>. The estimated errors ϵ_0 are given in Table <ref>. 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 <ref> 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 (θ_ℓ,ϕ_ℓ,λ_ℓ)∈[0,π]. The relevant parameters for these two matrices are given in Table <ref>. The results again evidence that the algorithm works well, even in the presence of complex phases ϕ_ℓ's and λ_ℓ's.




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. 














§ DISCUSSIONS




A comparison of  computational  resources is given in Table <ref>. For  hybrid quantum walk algorithm, we need 1+log(N) qubits, qlog(N) CNOT gates, and qlog(N) 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 1+2qlog(N) 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 𝒪(log (N)).



(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 <cit.>, like periods and unwanted correlations.

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








§ CONCLUSION


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 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. 


























































































































































We thank Chia-Cheng Chang, Yu-Cheng Su, and Rudy Raymond for discussions. Access to IBM Q systems is provided by IBM Q Hub at  National Taiwan University. This work is supported in part by Ministry of Science and Technology, Taiwan, under grant No. MOST 107-2627-E-002 -001 -MY3, MOST 106-2221-E-002 -164 -MY3, and MOST 105-2221-E-002 -098 -MY3.


















 
 

lemmaLemma
theoremTheorem







§ GRAY CODE BASIS

The natural binary code B=(B_n-1,B_n-2,...,B_1,B_0) is transformed to the Gray code basis <cit.> according to

    g(B)_i=B_i+1⊕ B_i ,

∀ i ∈{0,...,n-1 } with B_n=0. The probability matrix in the Gray code basis is given by
	
    P^quantum_J',J   =   ∏_ℓ=0^n-1|U_3(θ_ℓ)_i_ℓ, i_ℓ+1|^2
       =   ∏_ℓ=0^n-1|cos^2 (θ_ℓ/2)|^1-(i_ℓ⊕ i_ℓ+1)|sin^2 (θ_ℓ/2)|^i_ℓ⊕ i_ℓ+1

 with i_N=0. 


 
	Let S_N be the set of all possible n-bit strings { (S_n-1,S_n-2,...,S_1, S_0) | S_i ∈{0,1} ∀  i  ∈{0,1,...,n-1 }} with n=log_2 N, and π be a permutation of the set S_N. If there exists a function f:S_N ↦ℝ such that for A ∈ℝ^N× N,
	
    A_I⊕ J,J=f(I)

	∀  I,J ∈ S_N, and if π is bitwise XOR homomorphic, then we have
	A_π(I⊕ J),π(J)=f(π(I)).



Since π is bitwise XOR homomorphic, 
 Eq. (<ref>)  leads to

    A_π(I⊕ J),π(J)   =    A_π(I)⊕π(J),π(J)
       =   f(π(I))

∀  I,J ∈ S_N. 



 
	Let B ∈ S_N be represented by (B_n-1,...,B_0).
	Let g:S_N ↦ S_N be a function that transforms from natural bit string  to Gray code according to g(B)_i=B_i+1⊕ B_i, ∀  i ∈{0,1,...,n-1 } with B_n= 0. Then g is a bitwise XOR homomorphism.



	Let I,J ∈ S_N be represented by bit strings (I_n-1,...,I_0) and (J_n-1,...,J_0), respectively. Using
	
    g(I)_i   =   I_i+1⊕ I_i
    
    	g(J)_i   =   J_i+1⊕ J_i

	with I_n=J_n=0, we get
	
    [g(I)⊕ g(J)]_i     =   g(I)_i ⊕ g(J)_i
       =   (I_i+1⊕ I_i) ⊕ (J_i+1⊕ J_i)
       =   (I_i+1⊕ J_i+1) ⊕ (I_i⊕ J_i)
       =   g(I⊕ J)_i .
 
	∀  i ∈0,...,n-1.


Using Lemma <ref> and Lemma <ref>, the following theorem is clear.

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





 
 



§ DERIVATION OF EQ. (<REF>)


We use the evolution operator given in Eq. (<ref>),

    𝒰 
      = ∑_i_n-1,...,i_0,i_-1∏_ℓ=0^n-1  [U_3( u_ℓ)]_i_ℓ,i_ℓ-1 (X_ℓ)^i_ℓ ( | i_n-1⟩_♢⟨ i_-1 |_♢ )

to compute the two-evolution operator

    𝒰^2    =   ∑_i_n-1,...,i_0,i_-1
    k_n-1,...,k_0,k_-1∏_ℓ=0^n-1  [U_3( u_ℓ )]_i_ℓ,i_ℓ-1  [U_3( u_ℓ )]_k_ℓ,k_ℓ-1
          × (X_ℓ)^i_ℓ + k_ℓδ_i_-1,k_n-1 ( | i_n-1⟩_♢⟨ k_-1 |_♢ )   .




Next we project the 𝒰^2 operator  on the  |ψ_0J⟩ state,

    𝒰^2|ψ_0J⟩   =   ∑_i_n-1,...,i_0,i_-1
    k_n-1,...,k_0,k_-1∏_ℓ=0^n-1  [U_3( u_ℓ )]_i_ℓ,i_ℓ-1  [U_3( u_ℓ )]_k_ℓ,k_ℓ-1
          ×δ_i_-1,k_n-1δ_0,k_-1  | i_n-1⟩_♢ | I ⊕ K ⊕ J ⟩_q    
       =   ∑_I,K   f(I,K)   | i_n-1⟩_♢     | I ⊕ K ⊕ J ⟩_q  ,

where f(I,K) is given in Eq. (<ref>) and 

    | I ⊕ K ⊕ J ⟩_q= | i_n-1⊕ k_n-1⊕ j_n-1,..., i_0⊕ k_0⊕ j_0⟩_q .

We then project 𝒰^2|ψ_0J⟩ on the final state | k ⟩_♢ | J' ⟩_q

    ⟨ J'|_q ⟨ k|_♢𝒰^2|ψ_0J⟩    =   ∑_I,K   f(I,K)   δ_k, i_n-1δ_J', I ⊕ K ⊕ J 
        =    ∑_I   f(I,I⊕ J' ⊕ J)   δ_k, i_n-1  ,

which leads to the  probability matrix  elements as

    P_J',J      =  ⟨ J'| _♢ [𝒰^2 |ψ_0J⟩⟨ψ_0J | (𝒰^†)^2 ]  |J'⟩
          = ∑_k |⟨ J'|_q ⟨ k|_♢𝒰^2 |ψ_0J⟩|^2 
          = ∑_k=0^1 |∑_I   f(I, J'⊕ J ⊕ I) δ_i_n-1,k|^2 .


 
 



§ TWO-EVOLUTION QUANTUM WALK ON N=4 GRAPH




The probability matrix elements P^quantum_J',J for two quantum evolutions 𝒰^2 on the four-node graph read 

    P_00   =   P_11=P_22=P_33
       =   1/4sin^2θ_0+1/8(1+cos^2θ_1 +cos^2θ_0 
           +4 cosθ_1 cosθ_0 +cos^2 θ_1 cos^2 θ_0 -sin^2 θ_1 sin^2 θ_0) ,
    
    P_01   =   P_23=1/4sin^2θ_1,
    
    P_02   =   P_13=1/4sin^2θ_1,
    
    P_03   =   P_12=1/4(1-2cosθ_1cosθ_0+cos^2θ_1 )   .
 
Surprisingly, in this case the matrix elements do not depend on the (ϕ_1,ϕ_2) and (λ_1,λ_2) phases. However,  the matrix elements do depend on complex phases when N≥ 8, as can be numerically  checked. Note that (P_01,P_02,P_23,P_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.














































