Quantum Linear System Algorithm for General Matrices in System Identification

Solving linear systems of equations is one of the most common and basic problems in classical identification systems. Given a coefficient matrix A and a vector b, the ultimate task is to find the solution x such that Ax=b. Based on the technique of the singular value estimation, the paper proposes a modified quantum scheme to obtain the quantum state |x〉 corresponding to the solution of the linear system of equations in O(κ2rpolylog(mn)/ϵ) time for a general m×n dimensional A, which is superior to existing quantum algorithms, where κ is the condition number, r is the rank of matrix A and ϵ is the precision parameter. Meanwhile, we also design a quantum circuit for the homogeneous linear equations and achieve an exponential improvement. The coefficient matrix A in our scheme is a sparsity-independent and non-square matrix, which can be applied in more general situations. Our research provides a universal quantum linear system solver and can enrich the research scope of quantum computation.


Introduction
System identification [1][2][3] is a common method to determine the mathematical model describing the behavior of classical systems. Thus, the future evolution of the system can be predicted through the identified system model, which is widely applied to common weather forecast, flood forecast, market trend, etc. The traditional system identification method, namely the classical identification method, mainly includes least squares method [4], impulse response method and maximum likelihood method [5,6]. Existing studies [2,3] found that solving linear systems of equations is the basis of system identification problems. In fact, not only system identification problems, the application of linear equations involves various fields of science and engineering, including machine learning [7], partial differential equations [8], classic control system, and so on. Therefore, solving linear systems of equations for general matrices is of great significance.
Due to the importance of linear systems of equations in various fields, the solution of linear equations has become an enduring issue, and many algorithms derived therefrom. The classical solvers mainly include: matrix elimination method [9] and Kaczmarz method [10]. The most famous one of the former is the Gaussian elimination method, which is often used to solve small linear systems of equations and is suitable for a general coefficient matrix. The Kaczmarz method is generally more practical in the field of largescale linear equations. The running time for these classical solvers scales as O(n 3 ), where n is the size of the matrix, which will cost a lot of computing resources in solving large-scale linear systems. However, quantum computation [11][12][13] is capable of greatly reducing the time complexity for matrix operation and numerical calculation, which can be regarded as a promising attempt as a computing tool to improve the identification efficiency.
Quantum computation is an emerging computing technology that regulates quantum information units to perform high-efficiency calculations based on the laws of quantum mechanics, including coherent superposition and entanglement [14]. In 1994, Shor proposed the algorithm for prime factorization [15] with exponential acceleration over classical algorithms, which shows the potential of quantum computation for the first time. Since then, quantum computation has reached an era of rapid development. In recent years, scholars have also made significant progress in quantum algorithm research, including Grover algorithm [16], quantum simulation [17][18][19], duality algorithm [20][21][22], linear systems of equations solver [23][24][25], matrix multiplication algorithm [26,27], and so on. For the high-dimensional linear systems of equations, there have been breakthroughs in the field of quantum computation. In 2009, Harrow, Hassidim and Lloyd [23] proposed the quantum linear system algorithm (HHL) to obtain the quantum state |x = |A −1 b corresponding to the solution of Ax = b in time O(poly log(n)), where the sparse matrix A ∈ R n×n and x, b ∈ R n , which can improve the computational efficiency with an exponential speedup over classical algorithms. The HHL algorithm is of great significance in the field of quantum information processing and has a wide range of applications in big data, machine learning, numerical computing and other scenarios. In 2018, Wossnig et al. [28] proposed a sparsity-independent quantum linear system algorithm (QLSA) based on a quantum singular value estimation algorithm (QSVE). After that, Shao and Xiang [29] modified the QSVE algorithm to adapt to the non-Hermitian case. Current algorithms for linear systems have been widely applied in the emerging research area of quantum information processing. However, existing quantum algorithms have different restrictions on matrix A, such as the most typical one of HHL algorithm, which requires A to be a sparse Hermitian matrix so that the unitary transformation e iAt [30,31] can be realized in a constant time. At present, the quantum algorithm suitable for arbitrary linear system of equations has not been fully studied.
Without loss of generality, existing quantum algorithms assumed that the coefficient matrix A is Hermitian as it is well known that the general case can be reduced to the Hermitian case by embedding a general rectangular matrix M into a block antidiagonal Hermitian matrix with the elements of M † and M in the lower and upper half, respectively [28]. Different from previous algorithms, we proposed a modified quantum scheme to solve the cases of general matrices directly, which can reduce the time complexity of solving the linear system of equations. Moreover, it may not be easy to expand A into a Hermitian matrix when A is given as quantum information. However, our scheme does not need such expansion and works well on the original non-Hermitian matrix, and hence it can be implemented more efficiently. Based on this idea, this paper considers three cases of the solution of linear systems and proposes a quantum linear system algorithm for general matrices, where A is not required to be sparse or square, which can effectively improve the computational efficiency and expand the application range of quantum computation. For the homogeneous linear equations, we design the corresponding quantum circuit to ensure the completeness of the solution, which supplies exponential speed-up over classical algorithms. Meanwhile, we modify the quantum phase estimation (QPE) circuit to determine the sign of the phase by setting a sign qubit, which can be widely applied to various quantum algorithms.
The rest of our paper is organized as follows. Section 2 analyzes a general model of classical identification system based on semi-tensor product and shows the detailed process of our quantum algorithms. In Section 3, we make a time complexity comparison between existing algorithms and our algorithms. Then, we perform a numerical simulation to clarify the process of quantum algorithm in Section 4. Finally, we conclude in Section 5.

The Classical System Identification Problem
Consider a general discrete model of system identification as follows: where x(i) is an n dimensional system state of the i-th sampling, u(i) is the input with dimension m, A is an n × n system matrix and B is an n × m matrix. The goal of system identification is to estimate the matrices A and B from a set of inputs {u(i)} and states {x(i)}. System identification problems can be expressed in terms of the semi-tensor product method [32]. As a kind of special matrix multiplication, the semi-tensor product generalizes the ordinary matrix multiplication to the general case. T ⊗ S denotes the Kronecker product of matrices T m×n and S p×q , which is expressed as Just as a computational tool for solving the model, T S denotes the semi-tensor product of matrices T and S: where l =lcm(n, p) is the least common multiple of n and p. The semi-tensor product is the generalization of matrix multiplication. When n = p, there are l = n = p and T S = TS.
, where s i is the ith column vector of the matrix S. Therefore, we may estimate A and B from a set of u(i) and x(i).
where x(i) T is a 1 × n matrix and V C (A) is an n 2 × 1 matrix. According to Equation (3), the least common multiple lcm(n, n 2 ) = n 2 , and x(i) T Suppose there are N + 1 observed samples, and . . .
Then we have The Equation (6) is a linear system of equations, and the task is to find the solution Y. In the Equation (6), H ⊗ I n is an Nn × (m + n)n matrix, where n is the dimension of system states, m is the dimension of the input, and N is the number of samples. For the high-dimensional identification system, the time complexity of classical algorithms is enormous and existing quantum algorithms can not directly solve the non-square linear systems of equations. In order to reduce the cost of computing resources, it is necessary to propose a quantum algorithm for general linear equations.

The Quantum Linear System Algorithm for General Matrices
Inspired by the singular value estimation algorithm [23,28,29], we propose a quantum algorithm for general linear systems of equations as follows.
Given a general linear equation Ax = b, the singular value decomposition is where A ∈ R m×n , x ∈ R n and b ∈ R m , σ i is the singular value of A, µ i ∈ R m and ν i ∈ R n are the left and right singular vectors, and . Let the rank of A be r(r ≤ m, n) and the rank of [A b] be q; the relation between the solution vector x of Ax = b and the r, q is: The linear system of equations Ax = b can be solved by a mathematical optimization technique of minimizing the sum of squares of errors between the solution and the actual data, which is the so-called least squares method In the Equation (7), {µ i } ∈ R m and {ν i } ∈ R n are a set of basis in m and n dimensional spaces. Therefore, x and b can be expressed as When Note that when r < n, α i (i = r + 1, . . . , n) is not assigned, and the equation Ax = b has infinitely many solutions. In engineering, we usually want to find out the lowest energy solution state x with x|x minimality, that is The goal is to convert whose detailed quantum process of our scheme is described as follows.
The following mappings to access to the data structure can be performed in O(polylog(mn)) time.
The data structure is based on an array of binary trees, each binary tree contains enough leaves that store the squared amplitudes of the corresponding matrix entry, which can be found in [28,33] with a detailed description of such a binary tree memory structure. In order to facilitate mathematical operation, we define two degenerate operators P and Q that operate only on valid input information |ξ , where the dimension of the input state |ξ |0 is reduced to the dimension of valid information |ξ , so P and Q are called degenerate operators. The maps P and Q append an arbitrary input state |ξ to a register that encodes: Similarly, it follows that P and Q have orthonormal columns and thus P † P = I m and The eigenvalues of S are e ±2πiϕ i , and the corresponding eigenvectors are ω ± i |w ± i = −P|µ i + e ∓2πiϕ i Q|ν i , where ϕ i is the phase of eigenvalues and ω ± i is the norm of eigenvectors. Then, it can be obtained Through phase rotation, the process of [28,29]. It is worth noting that the above step is avoidable for the case of the coefficient matrix A being Hermitian. At this point, the singular value decomposition is A = ∑ i σ i µ i µ T i , and the task is to convert |b = ∑ i β i |µ i to the solution |x = ∑ i β i /σ i |µ i such that Ax = b. However, for the case of non-Hermitian, it is necessary to realize the transformation of quantum states |µ i to |ν i .
For a general linear system of equations with m = n, the above derivation will have some changes. For i ∈ [1, r], Equations (14) and (15) are valid. While i > r, A|ν i = 0, A † |µ i = 0, and we can obtain and At this point, e 2πiϕ i = −1 is the eigenvalue of S, that is, ϕ i = ±1/2, and the corresponding eigenvectors are Q|ν i and P|µ i . In order to achieve |b we first need to eliminate the formula ∑ m i=r+1 β i |µ i . Based on these definitions, we show the basic procedure of our algorithm: 1.
Preparing the initial quantum state |b = Σ m i=1 b i |i , which can be represented as: 2.
Apply P in the initial state |b 3.

4.
Apply a phase shift operator controlled by the phase ϕ i , then we obtain

5.
Perform a controlled rotation on the ancillary qubit based on the register storing phase value ϕ i and will obtain where σ i = cos(π ϕ i ) A F and t = min i |σ i |, i ∈ [1, r].

6.
Apply the inverse transformation of step 3 to obtain

7.
Measure the ancillary register. When the measurement result is |0 , the quantum state will collapse to 8.
Apply the inverse of Q and we will obtain the desired state which is the particular solution of the equation Ax = b, that is, the lowest energy solution state.
The quantum gate circuit of our quantum algorithm is shown in Figure 2.  The operator R is a quantum controlled rotation gate. When the phase ϕ i = ± 1 2 , R is a NOT gate, otherwise Note that the actual phase value is ϕ i ∈ (−1, 1), which serves as the control qubits of the phase shift operation in the step 4, while the previous quantum phase estimation algorithm outputs phase value ϕ i ∈ (0, 1). Therefore, we design a modified quantum phase estimation circuit to determine the sign of the phase in Figure 1. In the modified QPE circuit, we can estimate the phase value in the range ϕ ∈ (−1, 1).

The Quantum Algorithm for Homogeneous Linear Equations
For the condition of r < n, the equation Ax = b has infinitely many solutions. Therefore, in order to obtain general solutions of the equation Ax = b, we need to solve the homogeneous linear equation Ax = 0. Since The description of solving homogeneous linear equations is essentially just finding the projection of a state onto the ground state for an operator [34]. Through the quantum circuit shown in Figure 3, we obtain the combination of the eigenvectors corresponding to σ i = 0 and make the solution of homogeneous linear equations complete. The input state is |c = ∑ n i=1 c i |ν i and the output state is |x = ∑ n i=r+1 c i |ν i , where |x is the combination of right singular vectors corresponding to the singular value 0 of A contained in |c . When the phase ϕ i = ±1/2, the output of the controlled-NOT gate R N is |1 , otherwise it outputs |0 . We can obtain the solution of the Ax = 0 when an arbitrary input |c contains right singular vectors of A. In order to ensure that the output |x is complete, we input n linearly independent |c = |0 , |1 . . . |n . In addition, the arbitrary r linearly independent x i can form the solution vector basis of the homogeneous linear equation.
Based on QSVE, our quantum algorithms are sparsity-independent and may be applied to non-square dense matrices.

Algorithms Complexity Analysis
Then, we analyse the time complexity of our quantum algorithms. The time complexity of our scheme includes the following two parts: quantum data generation and the quantum algorithm process. On the one hand, relying on a binary tree memory structure detailed as described in [28,33], where the matrix entries associated with matrix A m×n are stored as suitable data structure, the oracle from classical data to quantum data can be implemented efficiently in time O(log 2 mn) and the data structure size is O(w log mn) where w is the number of non zero entries in A. On the other hand, based on the quantum singular value estimation algorithm, our algorithm achieves a runtime O(κ polylog(mn)/δ), where κ is the condition number of the coefficient matrix A, and δ denotes the precision parameter.
Define that the additive error achieved in output statex is , which means if x is the exact result andx is the result obtained from quantum algorithms, then x −x ≤ . In order to achieve accuracy , the precision parameter of our algorithm needs to reach δ = /(κ A F ). Assuming the spectral norm A * is bounded by a constant, since where r is the rank of matrix A. Therefore, our quantum algorithm has the time complexity of O(κ 2 √ rpolylog(mn)/ ). Remarkably, when A is sparse, exponential acceleration is achievable.
In the quantum algorithm of homogeneous linear equations, the time complexity of QSVE is O( polylog(mn)/ ). In view of the success probability of the ancillary register collapses to |1 , we need to repeat the coherent computation n/(n − r) times on average. Therefore, the runtime of the quantum homogeneous linear equation algorithm is given by O(n polylog(mn)/((n − r) )).
After obtaining the quantum state |x corresponding to the solution of Ax = b, we need to simulate the subsequent system states through the identified system, where . According to Equations (4)-(6), we can obtain where The inner product between pairs of states can be implemented in time O(ploylog(m + n)) by the swap text algorithm [35]. Therefore, we can predict the system state at the next moment based on the known system state x(i) and input u(i).
For the case of general matrix A m×n , previous quantum algorithms generally convert A to a Hermitian matrix:

Numerical Simulation
To clarify the process of our algorithm and prove the feasibility of algorithms, we perform simulation on an illustrative example.
For simplicity, we consider a first-order discrete model of classical system as follows: where x(i) is the system state of the i-th sampling and u(i) is the input state. The goal of system identification is to estimate coefficients a and d from a set of u(i) and x(i). Assuming that the initial system state x(1) = 3, the input states u = {4, 3, 0} and the evolved system states x(2) = −4, x(3) = x(4) = 0, the mathematical model can be transformed into a general linear systems of equations The maps U P and U Q append an arbitrary input state to a register that encodes the Equation (12), which can be realized by quantum gate circuits in Figure 4, and matrix forms of these maps The following shows the detailed procedure of the numerical simulation: 1.

Perform phase estimation on
then we obtain the following state where the eigenvalues of S are λ i = i, i, −i, −i, 1, −1, and ω i |w i is the corresponding eigenvector. 4.
Change the phase, then we obtain
Apply the inverse transformation of step 3 to obtain 7.
Apply the inverse of Q and we will obtain the desired state 8.
Measure the ancillary register. When the result is |0 , the quantum state will collapse to that is proportional to the solution of the equation Ax = b, so we obtain a = 3 5 C and d = 4 5 C. Substituting a and d into Equation (29), we obtain C = − 4 5 . So far, the first-order discrete identification model is: We simulated a 6-qubit quantum circuit diagram on the Origin Cloud, as shown in Figure 5.  [5] is the ancillary qubit and q [4] is the register storing input and output information. From left to right, the controlled rotation gate is 3 that generates the initial state P|b , The controlled gates According to the simulation result in Figure 6, when the ancillary qubit q [5] = 0 and the register storing phase information is restored to q[0] = q[1] = q[2] = 0, the probabilities of the output qubit q [4] are P{|0 } = 0.086 and P{|1 } = 0.16. Therefore, the amplitudes of q [4] are A{|0 } = P{|0 }/(P{|0 } + P{|1 }) = 0.59 and A{|1 } = P{|1 }/(P{|0 } + P{|1 }) = 0.81, and the solution quantum state is q [4] = 0.59|0 + 0.81|1 , which is consistent with the expected quantum state based on our algorithm.  As a comparison, we simulated a second-order discrete system identification model When the ancillary qubit q[6] = 0 and the register q [5] is restored to 0, the probabilities of the output qubits q [4] and q [3] are shown in the Figure 7. Thus, the solution quantum state is |x = 0.493|00 + 0.5|01 + 0.516|10 + 0.485|11 , which is the lowest energy solution among infinitely many solutions. So far, according to the existing samples, the secondorder discrete identification model with the lowest energy is: x(i + 1) = 0.493 0 0.5 0 x(i) + 0.516 0.485 u(i).

Conclusions
This paper develops a quantum algorithm of general linear equations for solving classical system identification problems. Our scheme can be finished in time O(κ 2 √ rpolylog(mn)/ ) for an m × n dimensional linear systems of equations Ax = b, where κ is the condition number of the linear equation, r is the rank of the matrix A and is the precision parameter, which is superior to existing algorithms. For the linear equation with non-square coefficient matrix, we discuss three cases of solutions, including the unique solution, approximate solution and infinitely many solutions. Our algorithm can obtain the unique solution, the approximate solution with the minimum error and the lowest energy solution among infinitely many solutions, which adapts to all cases of linear systems of equations. For the case of infinitely many solutions, we design a quantum circuit to obtain general solutions in time O(n polylog(mn)/((n − r) )), which can achieve an exponential improvement over classical algorithms. In addition, we design a modified QPE circuit to obtain a wider range of phase values, which can expand the application range of quantum phase estimation.
Based on QSVE, our algorithms is sparsity-independent compared with HHL algorithm. Meanwhile, we have extended the existing quantum linear system algorithms to general equations, which can effectively enrich the application area of linear systems of equations. For large-scale linear systems, such as machine learning, numerical calculation of partial differential equations, etc., our algorithms will have a wider range of applications and is of research significance. In the future work, we will focus on how our algorithms are implemented on quantum computers and how to apply on them to real practical problems.