Quantum Algorithms for Calculating Determinant and Inverse of Matrix and Solving Linear Algebraic Systems

Quantum algorithms exploiting principles of quantum physics can be run on quantum computers and promise to offer considerable advantages over classical counterparts. Algorithms for quantum calculations are an increasingly prominent field in quantum information research. Notable examples include Shor’s factoring algorithm [1,2], Deutsch parallelism [3], Grover’s search algorithm [4], Quantum Fourier transform [5–7], and Quantum phase estimation [7,8]. Further advancements include the HHL algorithm for solving systems of linear equations [9–16], in which the crucial point is the conditional rotation of the state of the one-qubit ancilla. This operation, where inversion of the eigenvalue requires involving classical computing, is placed between the blocks of the direct and inverse phase estimation algorithms. HHL algorithm exponentiates a Hermitian matrix of the considered linear system, which can be performed via the Trotterization formula [17,18]. There is extensive literature devoted to improvement of this algorithm, especially its scaling and complexity [19–28]. However, the principal problem of eigenvalue inversion is not resolved via quantum methods yet. For this purpose, for instance, quantum state tomography can be used to obtain the eigenvalues of the exponent of Hermitian operator (the matrix of the linear system) [17]. Nevertheless, the HHL algorithm combining the phase estimation and classical inversion of eigenvalues is very popular in implementation of various algorithms of quantum computation, and it is widely applicable in quantum machine learning [11,29,30]. In particular, it is applied in the least-square linear-regression algorithms [31–33] working with large data sets, in quantum algorithms for singular-value decomposition [21,34], in algorithms for solving linear [35] and nonlinear [36] differential equations.

Another method of matrix inversion is presented in [37,38]. Both algorithms are based on the singular value decomposition and apply the function evaluation problem [21] to approximate the inverse of singular values by oddpolynomial, which requires introducing a special scale parameter whose value depends on the length of the interval including nonzero eigenvalues.

We also refer to algorithms for matrix manipulations based on the Trotterization method and the Baker-Champbell-Hausdorff formula for exponentiating matrices [39]. In [21], algebraic manipulations with matrices are implemented using unitary transformations that include the matrices as specific blocks (block-encoding model). A similar blockencoding approach was used in the algorithm of [40] to embed the inverse of the matrix of an algebraic system into the unitary transformation. In [41], an algorithm for matrix multiplication using binary encoding of matrix elements in the states of quantum systems followed by binary multiplication using Toffoli gates was proposed. Matrix multiplication over rings, including the multiplication of Boolean matrices, was studied in [42].

In our recent work [43], we presented alternative algorithms for matrix manipulations, such as sum and product of two matrices, calculation of determinant, matrix inversion, and solving linear systems, based on special unitary transformations of a quantum system. The main distinction of these algorithms from those mentioned above is the encoding of matrix elements into the pure states of specific (sub)systems as the probability amplitudes. However, the realization of those unitary operators in terms of the basic quantum operations was not explored in that study.

Later, in [44], we proposed the implementation of unitary transformations corresponding to vector inner product, matrix sum, and product using multi-qubit conditional operators (reducible to Toffoli gates) and Hadamard transformations, including the ancilla measurement for garbage removal. The depth of these algorithms increases logarithmically with the matrix dimension, or even remains independent of it (in the case of matrix sum), demonstrating the advantage of these algorithms compared to their classical counterparts.

In this paper, we propose quantum algorithms for determinant and inversion of matrix with their application to solving systems of linear equations. The principal part of these algorithms is the algorithm for computing the determinant of a square matrix. It is built on the ideas of Refs. [43,44]: -

Encode matrix elements into the states of certain quantum subsystems through probability amplitudes in such a way that all necessary multiplications of matrix elements are performed automatically and appear in the probability amplitudes of the initial state of the whole system;

We follow the above steps to propose an algorithm for calculating the determinant of a square matrix whose N rows are encoded into the pure states of N separate quantum subsystems. Notice that creating such initial states is a special problem in itself, but we do not discuss this problem here. The state of the whole system results in all possible products of N elements from different rows. One only needs to select the proper products to calculate the determinant. We propose such an algorithm and construct the quantum circuit realizing it. While developing the algorithm, our focus is on the space and depths of its specific blocks. Essentially, our algorithm is a quantum realization of the cofactor expansion of the classical algorithm, whose depth exponentially increases with matrix dimension. However, the resulting depth of our algorithm is O(N log² N) (or O(N² log N) in the worst case) with space of O(N log N) qubits. The final step of this algorithm is the one-qubit ancilla measurement which removes all the garbage. The probability of obtaining the required result of the ancilla measurement (which is 1) diminishes exponentially as the matrix dimension N increases, necessitating an exponentially larger number of algorithm evaluations. This is a central problem that we plan to investigate in future work.

It is remarkable that the minor modification of the algorithm for calculating the determinant yields an alternative quantum algorithm for matrix inversion using its definition in terms of algebraic complements. We again use the row-wise encoding of the matrix elements into the state of a quantum system so that all necessary products of matrix elements appear automatically in the superposition state. One thus only needs to properly sort these products. The depth of this algorithm is O(N² log N) (or O(N² log² N) in the worst case) and the space is O(N log N) qubits. Based on such inversion algorithm we develop an alternative algorithm for solving systems of linear algebraic equations. Thus, we do not use Trotterization methods and classical subroutine for inversion of eigenvalues. Therefore, the proposed algorithm is purely quantum in nature, and it gives an alternative avenue to the inversion algorithm considered in [9,17,28]. The estimations for the depth and space of this algorithm remain the same.

We also emphasize that, in comparison with inversion algorithms in [37,38], which approximately calculate the inverse eigenvalues, our algorithm constructs exact (up to the available accuracy of realization of quantum operations) inverse matrix without appealing to its eigenvalues which provides certain advantages of the matrix-encoding approach.

Regarding the question of algorithmic complexity, we remark that the specialized classical techniques of matrix manipulations facilitate a reduction in the depth of the algorithms for determinant calculation to O(log² N) [45,46] with the size O(N^2.496). In our study, we employ the Leibniz formula as a definition for the matrix determinant and define the inverse matrix in terms of the algebraic complements. Our quantum algorithms based on these definitions demonstrate clear advantages over their classical counterparts. It is plausible that the adaptation of specific classical algorithms for matrix operations into their quantum equivalents could yield further enhancements in the efficiency and complexity of quantum algorithms designed for calculating determinants and inverse matrices. However, this particular aspect is not addressed in the present paper. The primary contribution of our work lies in promoting the representation of matrix elements as probability amplitudes of pure states, and we elucidate the benefits of this representation in relation to standard algorithms for matrix manipulations.

The paper is organized as follows. The algorithm for calculating the matrix determinant is proposed in Section 2. This algorithm is taken as a basis for the algorithm calculating the inverse matrix presented in Section 3. Based on those results, the algorithm for solving systems of linear algebraic equations is proposed in Section 4. The main results are discussed in Section 6. The special subroutine for the controlled measurement is presented in Section 5. Appendix A contains the modified algorithms for matrix multiplication and examples of realization of the algorithms developed in this paper.

To calculate the determinant of N × N matrix à (we assume N = 2ⁿ) with elements aj_k(j, k = 0, 1, …, N – 1), 1A˜=(a00a01⋯a0(N−1)a10a11⋯a1(N−1)⋯⋯⋯⋯a(N−2)0a(N−2)1⋯a(N−2)(N−1)a(N−1)0a(N−1)1⋯a(N−1)(N−1)), we encode each row consisting of N elements into the pure state |Ψ_j〉 of the n-qubit subsystem S_j, 2|Ψj〉=∑k=0N−1ajk|k〉Sj, j=0,1,…,N−1, where k in the state |k〉_Sj is the binary encoding of the integer k, therefore all states |k〉_Sj with fixed j are mutually orthogonal. We shall note that encoding the elements of the matrix into the probability amplitudes of a quantum state is not a simple problem, but we do not consider aspects of creating such states in this paper. Representing rows of the matrix à in the form of pure states (2) creates N obvious constraints for the matrix elements a_ij because of the normalization of pure states: 3∑k=0N−1|ajk|2=1, j=0,…,N−1.

These constraints can be circumvented by incorporating additional terms with different probability amplitudes into the states |Ψ_i〉, but this issue will be addressed in future study. The whole system is the union of the subsystems, S=∪jSj , and the pure state of the whole system reads 4|Φ0〉=|Ψ0〉⋯|ΨN−1〉=∑k0,…,kN−1=0N−1a0k0…a(N−1)kN−1× |k0〉S0⋯|kN−1〉SN−1, which includes N^(S) = nN qubits. In computing the determinant, we need only those terms in which all k_j are different. We select these terms from the whole pure state of our system: 5|Φ0〉=∑k0,…,kN−1=0k0≠…≠kN−1N−1a0k0…a(N−1)kN−1×|k0〉S0⋯|kN−1〉SN−1+|g1〉S, where |g₁〉_S collects all extra terms which form the garbage to be removed later. We emphasize that the selection of the needed terms in (5) is quite a formal procedure, at this stage we don’t do anything to label these terms. We will do this later by introducing the ancillae B and removing all the garbage via the ancilla measurement, see Eqs. (19), (20) and Figure 2. We have to arrange the appropriate signs for each term of the first part on the right-hand side of (5). This can be achieved by setting the subsystem S_j (in each term) into the state |j〉 using the SWAP operators together with sign-switching. This pair of operators performs permutation and changes the sign of the term.

Thus, our purpose is to exchange the states of subsystems S_j(j = 0, …, N – 1) so that the state of the subsystem S_k becomes |k〉S_k in each term of the first part of superposition state (5). For that we introduce the projectors PSj(k) , 6PSj(k)=|k〉Sj Sj〈k|,k=0,…N−2, j=k+1,…N−1, and the ñ_k-qubit ancillae A_k, (ñ_k = ⌈log(N – k)⌉, k = 0, … N – 2). The ancilla A_k serves to save the index j of the subsystem S_j whose state is |k〉 with j > k. Therefore, the number of qubits ñ_k in A_k is less than n = log N unless k = 0, 1: ñ_k≤ n.

For instance, in Figure 2 for the case N = 4, there are three ancillae A₀, A₁, A₂. Ancilla A₂has to encode either only one state |2〉_S3 or nothing, therefore one qubit is enough: we can encode |2〉_S3 by the state |1〉_A2 of ancilla A₂. Next, A₁ has to encode two states |1〉_S2 and |1〉_S2 or nothing. Two qubits are enough for that: |1〉_S2 → |1〉_A1, |1〉_S3 → |2〉_A1. Similarly, A₀ requires 2 qubits as well: |0〉_S1 → |1〉_A0, |0〉_S2 → |2〉_A0, |0〉_S3 → |3〉_A0. Thus, we have the map 7|k〉Sj→|j˜〉Ak, which can be realized as follows.

Using the projector PSj(k) and ancilla A_k we can prepare the operator V_SjAk, which relates the index j of S_j to the appropriate state |j˜〉Ak of the ancilla A_k: 8VSjAk=PSj(k)⊗Zk(j)σkj(x)+(ISj−PSj(k))⊗IAk, where σkj(x) includes the set of σ^(x)-operators which change the ground state of the ancilla A_k to the state |j˜〉Ak , subscripts k and j are related to k and j in the map (7). For instance, to map the state |0〉_S2 into the state |2〉_A0 of the two-qubit ancilla A₀ we need to apply σ^(x) to the second qubits of A₀, i.e., σ02(x)=I⊗σ(x) , where I is the identity operator applied to the first qubit of A₀; Z_k(j) changes the sign of the ancilla state |j˜〉Ak (essentially, it counts swaps and can be represented by either the σ^(z)-operator applied to one of the excited qubits, orby the σ^(x)σ^(z)σ^(x)-operator applied to one of the ground-state qubits (only to one!) of A_k).

It is hard to derive the explicit formula for Z_k(j) in general case, but it can be simply done in any particular case. For instance, in the above example of the transformation |0〉_S2 → |2〉_A0 we can set Z₀ (2) = I ⊗σ^(z). Thus, the operator V_SjAk acts as follows. If the state of S_j is |k〉_Sj, then the state of A_k becomes |j˜〉Ak , otherwise the state of A_k remains |0〉_Ak.

To arrange the exchange of states between S_j and S_k, we introduce the projectors PAj(j) , 9PAk(j)=|j˜〉Ak Ak〈j˜|, k=0,…N−2,j=k+1,…N−1, j˜=0,…N−k−2, j˜=j−k−1, and construct the operators U_SkSj: 10USkSj =(SWAPSk,Sj⊗PAk(j)+ISkSj⊗(IAk−PAk(j))), k≠j.

Here, the operator SWAP_Sk,Sj swaps the states of the subsystems S_k and S_j by swapping states of appropriate qubits of these subsystems. Thus, the operator U_SkSj swaps the states of S_k and S_j if only the state of A_k is |j˜〉Ak . Otherwise the states of the above subsystems remain unchanged.

The circuits for the operators V_SjAk and USkSj(k) are given in Figure 1a (the general matrix case) and in Figure 1b (the particular 4 × 4-matrix case).

Figure 1.

V and U operators for (a) the general case and (b) the 4 × 4-matrix case. In the latter case, ñ₀ = ñ₁ = 2, ñ₂ = 1, ñ_k < n. Hence, A₁ is a two-qubit ancilla.

The depth of each operator V_SjAk is O(n) (due to the projectors PSj(k) with n control qubits and we assume that all σ^(x)-operators of σkj(x) are applied in parallel to different qubits of A_k), while the depth of U_SkSj is the depth of the control-SWAP-operator, which includes n swaps each controlled by n qubits, i.e., the total depth is nO(n) = O(n²). Theoretically, it can be implemented with less operations but this would not change the asymptotic complexity of the circuit.

Applying the operators 11Wk=∏j=k+1N−1USkSjVSjAk, k=0,…,N−2, to the state |Φ₀〉|0〉_A, we put the subsystem S_k into the state |k〉_Sk in each term of the first part of (5). Although W_k has the product of N – k – 1 pairs of U- and V-operators, only one pair of such operators is applied to each particular term in |Φ₀〉 which is due to the structure of the projector operators PSj(k) . We assume that this means the UV-products can be applied in parallel and thus the depth of W_k is O(n) + O(n²) = O(n²). Otherwise we have to multiply the obtained depth by O(N). The operator ∏k=0N−2Wk provides the needed sign of the probability amplitudes which will compose the expression for the determinant. To collect all of them in a single probability amplitude we apply the Hadamard transformation to each qubit of ancillae A_k, k = 0, …, N – 2, and denote this transformation H_A, i.e. 12HA=HA0…HAN−2′ HAk=⊗k=0n˜k−1H.

There are Ñ Hadamard operators in H_A,

Applying the operator WSA(1) , 14WSA(1)=HA∏k=0N−2Wk, to the state |Φ₀〉|0〉_A, A=∪k=0N−2Ak , and selecting the terms which carry the information about the determinant (these terms correspond to the ground state of A) we obtain 15|Φ1(D)〉=WSA(1)|Φ0〉|0〉A=det(A˜)2N˜/2 |0〉S0|1〉S1…|N−1〉SN−1|0〉A+|g2〉SA.

Here we introduce the determinant det(Ã) represented by the Leibniz formula, 16det(A˜)=∑k0≠…≠kN−1ϵk0…kN−1ak00…akN−1(N−1), where ϵ_k0…kN–1 is the permutation tensor: ϵk0…kN−1={ (−1)perm(k0…kN−1),k0≠k1≠…≠kN−10otherwise.

The depth of W⁽¹⁾ is O(Nn²). Next, we can set states of all S_j in the first term of |Φ1(D)〉 into the |0〉_Sj states. To this end we apply the σ^x-operators to the excited qubits of S, i.e. prepare the operator X_S. 17XS=∏j=1N−1XSj(j), with X_Sj(j) being the set of σ^(x) operators acting on the excited qubits of S_j in the state |j〉_Sj. The depth of X_S is O(1). Applying X_S to |Φ1(D)〉 we obtain: 18|Φ2(D)〉=XS|Φ1(D)〉=det(A˜)2N˜/2|0〉S|0〉A+|g3〉SA.

We note that the operator H_A in WSA(1) and operator X_S can be applied simultaneously because they affect different qubits. We separate them to underline the distinction in their functions. Next, we remove the garbage. For that we introduce the projector P_SA = |0〉_S|0〉_{A A} 〉0|_S〉0|, include the one-qubit ancilla B in the ground state, and construct the operator WSAB(2) , 19WSAB(2)=PSAσB(x)+(ISA−PSA)IB, whose depth is O(Nn + Nn) = O(Nn). Applying it to |Φ2(D)〉|0〉B we obtain 20|Φ3(D)〉=WSAB(2)|Φ2(D)〉|0〉B=det(A˜)2N˜/2|0〉S|0〉A|1〉B+|g3〉SA|0〉B.

Now we measure the ancilla B with the output |1〉_B. The probability of this result of measurement is 2^−Ñ|det(Ã)|². After the measurement we get the following state: 21|Φ4(D)〉=|Ψout(D)〉|0〉A, |Ψout(D)〉=arg(det(A˜))|0〉S,

The particular circuit for calculating the 4 × 4 determinant using the proposed algorithm is given in Figure 2. The circuit for the determinant of N × N matrix includes N log N qubits of the subsystem S, Ñ = O(N log N) qubits of subsystem A, and one qubit of ancilla B, therefore the space of this circuit is O(N log N) qubits.

Figure 2.

Quantum circuit illustrating computation of determinant of a 4 × 4 matrix, see notations in Figure 1. We omit superscript (D) in |Φ_k〉(k = 1, …,4) and |Ψ_out〉 for brevity. Operators H_A and X_S can be applied in parallel since they affect different qubits. The lower circuit is the notation for multi-qubit controlled σ^(x) operation.

We remark here that the minor modification of the algorithm for calculating the determinant of a square N × N matrix in Section 2 allows to calculate the inverse of a square (N – 1) × (N – 1) matrix. Moreover, this modification leads to a minor increase in the depth of the algorithm which becomes O(N² log N) instead of O(N log² N) in the algorithm for calculating the determinant. These facts motivate working out the quantum inversion algorithm based on the definition of the inverse matrix in terms of algebraic complements, although quantum analogs of other algorithms of matrix inversion might be more promising.

3.2.

Basic Algorithm

The circuit for the algorithm is presented in Figure 3. According to (2), the pure state |Φ₀〉 of the whole system is given in (4) with normalization (3) and constraints (22) for the elements in the first row and first column of the matrix Ã, see (23). From this state, we select only those terms that compose the determinant of à and therefore contain the terms that appear in the resulting inverse matrix A−1 , i.e., we consider |Φ₀〉 in form (5).

Figure 3.

(Color online) (a) Structure of quantum circuit for matrix inversion (ancillae are shown in green). (b) The structure of the block P.

However, not all terms in the first part of |Φ₀〉 are needed in our algorithm. We have to keep only those of them that include the 1st power of the factor q, (see constraints (22)), i.e., exactly one state |0〉_Sj, 1 ≤ j ≤ N – 1, can appear in the product |k₁〉_S1 … |k_N–1〉_SN–1. We also have to associate each of the probability amplitudes in (5) with the appropriate element Aji−1 of the inverse matrix. To do so, we introduce two n-qubit subsystems R and C whose basis states serve to enumerate, respectively, rows and columns of the inverse matrix A⁻¹.

To select the terms that compose Aji−1 from the superposition state we use the projector 28Pij=|0〉Si Si〈0|⊗|j〉S0 S0〈j|, i,j=1,…,N−1.

The idea is that this projector allows labeling the 1st and (j + 1)th (j = 0, …, N – 1) columns and the 1st and (i + 1)th (0 ≤ i ≤ N – 1) rows in the matrix Ã(23). As a result, we select those terms from the superposition state whose probability amplitudes compose M_ij.

Now we can construct the controlled operator 29WSRC(ij)=Pij⊗σR(j)C(i)(x)+(IS−Pij)⊗IRC.

This operator conditionally applies σR(j)C(i)(x) which is the product of σ^(x) (acting on the subsystems R and C) such that σR(j)C(i)(x)|0〉R|0〉C=|j〉R|i〉C . The depth of WSRC(ij) is O(n). To establish it we use the fact that the controlled operator WSRC(ij) with n-qubit controlling register can be represented in terms of Toffoli gates by a circuit of depth O(n) [47] and require that all σ^(x) operators in σR(j)C(i)(x) can be applied simultaneously. Applying the operator W^SRC=∏i=1N−1∏j=1N−1WSRC(ij) (whose depth is O(N²n)) to |Φ₀〉|0〉_R|0〉_C we obtain 30|Φ1〉=W^SRC|Φ0〉|0〉R|0〉C =1N∑i=1N−1∑j=1N−1∑k1…ki−1,ki+1…kN−1=00≠k1≠…≠ki−1≠ki+1≠…≠kN−1≠jN−1a1k1…a(i−1)ki−1ai0a(i+1)ki+1…a(N−1)kN−1×|j〉S0|k1〉S1…|ki−1〉Si−1|0〉Si|ki+1〉Si+1…|kN−1〉SN−1|j〉R|i〉C+|g2〉SRC, where a_i0 = q. In the RHS of (30), the first part collects those terms of the state |Φ₁〉 that are labeled by the states |j〉_R and |i〉_C of the subsystems R and C, respectively, with j > 0 and i > 0. That is why the states of S₀ and S_i are, respectively, |j〉_S0 and |0〉_Si in the selected part of |Φ₁〉. Now we have to arrange the proper sign of each term in the first part of (30). According to Section 3.1, this can be done by calculating the determinant of à using the algorithm (1) proposed in Section 2. Applying the transformation WSA(1) given in (14) to |Φ₁〉 we obtain 31|Φ2〉=WSA(1)|Φ1〉|0〉A=−qdet(A)2(N˜+n)/2∑i=1N−1∑j=1N−1Aji−1|0〉S0|1〉S1⋯⊗|N−1〉SN−1|j〉R|i〉C|0〉A+|g3〉SRCA.

Next, applying the operator XS=∏j=1N−1XSj(j) given in (17) to |Φ₂〉 we obtain (we note that the subsystem S=∪j=0N−1Sj has the same number of qubits in both the algorithm of determinant calculation and the algorithm of matrix inversion, therefore we can apply the same operator X_S in both cases) 32|Φ3〉=XS|Φ2〉=−qdet(A)2(N˜+n)/2∑i=1N−1∑j=1N−1Aji−1|0〉S|j〉R|i〉C|0〉A+|g4〉SRCA,

Finally, applying WSAB(2) (19) to |Φ₃〉|0〉_B we obtain 33|Φ4〉=WSAB(2)|Φ3〉|0〉B=∑i=1N−1∑j=1N−1a˜ji|0〉S|j〉R|i〉C|0〉A|1〉B+|g5〉SRCA|0〉B, where 34a˜ji=−qdet(A)Aji−12(N˜+n)/2.

The depth of the operator DSAB=WSAB(2)XSWSA(1) (under certain assumptions on the physical realization of the operator W_k, see Eq. (11)) is O(Nn²) < O(N²n) (the depth of W^SRC ). After measuring the ancilla B to get the output |1〉 with the probability q2|det(A)|22N˜+nG2~2−O(N˜)=2−O(Nn)) , G=∑i,j=1N−1|Aij−1|2 , we get the desired result 35|Φ5〉=|Ψout〉|0〉S|0〉A, where 36|Ψout〉=G−1∑i=1N−1∑j=1N−1Aji−1|j〉R|i〉C=G−1|A−1〉.

The state |Ψ_out〉 encodes A⁻¹ in the normalized form, i.e., 〈Ψ_out|Ψ_out〉 = 1, while the state |A−1〉=∑i=1N−1∑j=1N−1Aji−1|j〉R|i〉C is not normalized. The depth of the discussed algorithm for matrix inversion is determined by the operator W^SRC and equals O(N² log N). The circuit includes N log N qubits of the subsystem S, Ñ = O(N log N) qubits of the ancilla A, one qubit of the ancilla B, and 2 log N qubits of R and C, thus the space of this circuit is approximated by the space of the circuit for calculating the determinant which is O(N log N).

We note that this depth is obtained under the assumption that the depth of the operator WSA(1) in (14) is N log² N. If this depth is N² log² N, then the depth of the matrix inversion algorithm is also N² log² N.

Having constructed the inverse matrix A⁻¹ encoded into the state of the subsystem R ∪ C, we obtain a tool for solving a system of linear algebraic equations represented in the matrix form, Ax = b, whose solution reads x = A⁻¹b.

Therefore, we have to apply the multiplication algorithm proposed in Ref. [44] (see also Appendix A.1, for its modified version) to the matrices A⁻¹ and b. The structure of the circuit is shown in Figure 4a. A more detailed circuit is presented in Figure 4b.

Figure 4.

(Color online) (a) Structure of quantum circuit for solving a system of linear algebraic equations. Both blocks share the same ancillae (not shown in figure). (b) Circuit for quantum algorithm solving a system of linear algebraic equations Ax = b. Circuit of matrix inversion (up to the operator WSAB(2) ) is denoted by a block A⁻¹. All input qubits of the matrix-inversion algorithm (system S) become ancillae qubits. The output of the matrix-inversion algorithm (systems R and C) together with system b encoding the vector b becomes the input of the matrix-multiplication block. The output of the whole algorithm is encoded in the state |x〉_R of the system R. All ancillae qubits are shown in green color. We omit the superscript (L) in |Φ_k〉(k = 5,6,7) and |Ψ_out〉 for brevity.

Notice that we can use the same ancillae in both blocks of this circuit (ancillae are not shown in Figure 4a). Moreover, the subsystem S encoding the matrix A can be used as ancilla in further calculations.

In this case we do not measure ancillae B (unlike Section 3.2) and thus start with the state |Φ₄〉 given in (33). The vector b is encoded as 37|b〉=∑j=0N−1bj|j〉b, b0=0, ∑j=1N−1|bj|2=1.

The zero value for b₀ is prescribed by the encoding of A⁻¹ in (36), where the 0-states of R and C are not used. Out of system S ∪ A (which forms an ancilla in Figure 4b), we take the one-qubit ancilla B˜ needed in the multiplication algorithm and neglect other qubits. We also apply the operator X_B = σ^(x) to invert the state of B. Thus, the state of the whole system reads 38|Φ5(L)〉=XB|Φ4〉|b〉|0〉B˜=∑i,j,l=1N−1a˜jlbi|j〉R|l〉C|i〉b|0〉B|0〉B˜+|g1〉RCb|1〉B|0〉B˜.

Now we refer to the Appendix, Section A.1, where the matrix multiplication algorithm is described. Applying the result of this Appendix we identify K = N, M = 1 and replace the subsystems R₁, C₁, R₂ with R, C, b, while the subsystem C₂ is not used because the second matrix in our multiplication is just a column. Thus, WC1R2(1) in (A5) becomes WCb(1) , WC1(2) in (A7) becomes WC(2) , while the operator WC1R2B˜(3) in (A8) must be modified including the projector acting on B, i.e., the new operator reads: 39PCBb=|0〉C|0〉B|0〉b C〈0|B〈0|b〈0|,WCBbB˜(3)=PCBb⊗σB˜(x)+(ICBb−PCBb)⊗IB˜.

Applying the operator WCBbB˜(3)WC(2)WCb(1) to |Φ5(L)〉 (38) we obtain 40|Φ6(L)〉=WCBbB(3)WC(2)WCb(1)|Φ5(L)〉=12n/2(∑j1=1N−1∑j=1N−1a˜j1jbj|j1〉R|0〉C|0〉b)|0〉B|1〉B^+|g2〉RCb|1〉B|0〉B˜.

Performing measurement over B˜ with the output |1〉B˜ we remove garbage and obtain 41|Φ7(L)〉=|Ψout(L)〉|0〉C|0〉b|0〉B,|Ψout(L)〉=G−1∑j1=1N−1∑j=1N−1Aj1j−1bj|j1〉R=G−1|x〉R, where |x〉R=∑j=1N−1Aj1j−1bj|j1〉R is an unnormalized vector, the normalization G=(∑i|xi|2)1/2 , xi=∑j=1(N−1)Aij−1bj , is known from the probability of the above measurement which equals (qGdet(A))²/2^Ñ+2n due to the expression (34) for a˜jl . Thus, the probability is ∼ 1/2^Ñ+2n = 2^−O(Nn). We note that, although we use two algorithms (matrix inversion and matrix multiplication), only one measurement is performed to remove all the garbage. The result is stored in the register R. The depth of the whole algorithm is defined by the depth of the subroutine inverting A and thus equals O(N²n). In comparison with the circuit for inverse matrix, the circuit in Figure 4b includes log N qubits encoding the column b and thus the space remains essentially O(N log N).

In all algorithms presented above the weak point is the small probability of success in the measurement of the ancilla state. In fact, this probability is ∼ 2^{−O(N log N)} in all cases, in other words, it exponentially decreases with an increase in the dimensionality N of the matrix à (1). To smooth this effect we suggest to replace the usual ancilla measurement by the controlled measurement [48], which can be presented as a special subroutine.

We call B the one-qubit ancilla whose state we have to measure and write the state of the whole system as 42|χ0〉=G−1(|E〉|1〉B+|F〉|0〉), 〈E|F〉=0, where only the state |E〉 contains the useful information, G=〈E|E〉+〈F|F〉 is the normalization, and the states |E〉 and |F〉 are not normalized. First, before measurement, we double the state of ancilla adding one more one-qubit ancilla B˜ in the ground state and apply the C-NOT CNOTBB˜=|1〉BB〈1|⊗σB˜(x)+|0〉BB〈0|⊗IB˜, where σB(x) and IB˜ are, respectively, σ^(x) and the identity operator applied to the ancilla B˜ . We have 43|χ1〉=CNOTBB˜|χ0〉|0〉B˜= G(|E〉|1〉B|1〉B˜+|F〉|0〉B|0〉B˜).

Now we introduce the controlled measurement 44MBB˜=|1〉BB〈1|⊗MB˜+|0〉BB〈0|⊗IB˜, where MB˜ is the measurement operator applied to B˜ . Applying the operator MBB˜ to |X₁〉 we have |χ2〉=MBB˜|χ1〉=|Ψout〉|1〉B, |Ψout〉=(〈E|E〉)−1/2|E〉.

In this way we avoid necessity of performing multiple runs of the algorithm to get the needed result of the ancilla measurement. However, we also lose the possibility to define the normalization (〈E|E〉)^1/2, which might be important information. For instance, in the algorithm for determinant calculation this normalization equals the absolute value of the matrix determinant, which is the half of useful information (another half is the phase). But in the other two algorithms this normalization doesn’t include the principal information. So, getting information about that normalization is the problem to be resolved. Another problem is the practical realization of the above controlled measurement that is also beyond the scope of this paper. The circuit for the subroutine of controlled measurement is shown in Figure 5.

Figure 5.

(Color online) Replacement of ordinary measurement of ancilla state (left circuit) with the subroutine of controlled measurement (right circuit).

We propose two quantum algorithms for matrix manipulation, which are calculating the determinant and matrix inversion. Both are based on special row-wise encoding of the matrix elements in the pure states of quantum subsystems.

The algorithm for computing the determinant of a square matrix is discussed in detail. It uses multiqubit control-SWAP and control-Z operations. The depth of this algorithm is O(N log² N). This algorithm utilized the subroutine realizing the completely antisymmetric tensor ε (operator WSA(1) ) that may have relevance in other algorithms. The outcome of this algorithm is split into two components: the absolute value of the determinant of the matrix A, denoted as |det(A)|, and the argument of the determinant of the matrix A, denoted as arg(det(A)). The first part is obtained from the measurement of an ancilla (although the probability of the needed outcome in such measurement decreases exponentially with N), while the argument of the determinant remains embedded in the probability amplitude of the output state |Ψout(D)〉 . The size of the subsystem S encoding the matrix A is N log N, the size of ancillae A_k, k = 0, …, N – 2, is Ñ = O(N log N). Both these sizes serve to estimate the space of the whole algorithm which is O(N log N) qubits.

It is remarkable that the (N – 1) × (N – 1)-matrix inversion algorithm is essentially the algorithm for calculating the determinant of the N × N matrix obtained by replacing the first row by the row of 1/N and the first column by the column of q’s (up to the diagonal element). This special encoding provides the superposition of all products of matrix elements which are necessary for calculating minors included into the definition of the inverse matrix. The algorithm distributes these products with proper signs among the appropriate inverse-matrix elements which are encoded into the quantum state of certain subsystem thus forming the output state. The depth of the circuit is O(N² log N).

As an application of inversion we develop the algorithm for solving systems of linear algebraic equations. This algorithm consists of two blocks: matrix inversion and matrix multiplication. At that, the matrix inversion is the most resource-consuming block which finally determines the depth of the algorithm (which is O(N² log N)), and probability of the required ancilla state after measurement, ∼ 2^{−O(N log N)}. We have to note that the depths of the algorithms for calculating the determinant and matrix inverse were obtained under the assumption that the U-, V- operators in (11) can be applied in parallel. Otherwise we have to change the depths O(N log² N) and O(N² log N) to O(N² log ²N) in both cases.

We emphasize that the matrix-multiplication algorithm is modified in comparison with that proposed in [44]. The detailed description of the modified version is presented in Appendix, Section A.1.

One notable aspect is the use of ancillae to gather all the garbage generated during the calculation, which can then be discarded with a measurement of a one-qubit ancilla. The weak point of algorithms is the probability of obtaining the required state of ancilla after measuring. This probability is ∼ 2^{−O(N log N)} in all three algorithms.

Combining the algorithm of matrix inversion with the algorithm of matrix multiplication we demonstrate the usage of two types of matrix encoding into the state of a quantum system. The first one is the row-wise encoding used in the input of both algorithm for calculating the determinant of the N × N matrix and matrix inversion algorithm (for simplicity, we assume that N = 2ⁿ). It requires N subsystems of log N qubits each, and the state of ith subsystem encodes the elements of the ith row of the matrix. The result of matrix inversion is stored in the two subsystems of log N qubits each. These subsystems enumerate rows and columns, respectively. This form of encoding of the inverse matrix allows to perform its multiplication with column b properly encoded in the state of log N-qubit subsystem to obtain the solution x of the linear system. This solution is a column whose elements are encoded into the state of a log N-qubit subsystem R.