Fundamental matrix operations and solving linear systems of equations are ubiquitous in scientific investigations. Using the 'sender-receiver' model, we propose quantum algorithms for matrix operations such as matrix-vector product, matrix-matrix product, the sum of two matrices, and the calculation of determinant and inverse matrix. We encode the matrix entries into the probability amplitudes of the pure initial states of senders. After applying proper unitary transformation to the complete quantum system, the desired result can be found in certain blocks of the receiver's density matrix. These quantum protocols can be used as subroutines in other quantum schemes. Furthermore, we present an alternative quantum algorithm for solving linear systems of equations.
Wentao Qi, Alexandr I Zenchuk, Asutosh Kumar, Junde Wu. Quantum algorithms for matrix operations and linear systems of equations[J]. Communications in Theoretical Physics, 2024, 76(3): 035103. DOI: 10.1088/1572-9494/ad2366
1. Introduction
Feynman and Deutsch have pointed out that it would be impossible to accurately and efficiently simulate quantum mechanical systems on a classical computer. They also conceived the idea that machines 'built of quantum mechanical elements which obey quantum mechanical laws' [1, 2] might be able to process information fundamentally more efficiently than typical classical computers. Such computational power would then have applications to a multitude of problems within and outside quantum mechanics, including information theory, cryptography, mathematics and statistics. Foreseeing such a situation, arose the need of building 'quantum computers' and devising 'quantum algorithms' that could be run on quantum computers. Over the past three decades there has been a substantial growth in research in both theoretical and experimental aspects of quantum computing. Quantum algorithms exploit principles of quantum physics, and are known to offer significant advantages over classical counterparts for a number of problems. The first evidence of quantum advantage was the Deutsch algorithm [2], a simple quantum algorithm to quickly verify if a Boolean function is constant or balanced. Henceforth, quantum algorithms were paid much attention and had been significantly developed. Namely Simon's algorithm [3, 4], Shor's factoring algorithm [5, 6] based on the quantum Fourier transform [7-9], and Grover's algorithm [10] are much celebrated.
It is possible to reduce many scientific problems for matrix formalism and for the problem of solving linear systems of equations. Fundamental matrix operations and solving linear systems of equations frequently arise on their own, as well as subroutines in more complex systems, and are ubiquitous in all realms of science and engineering including machine learning and optimization. However, with huge data sets and large dimensions of physical systems, such tasks can be practically intractable for classical computers in terms of data processing speed, storage and time consumption. A quantum algorithm from Harrow, Hassidim and Lloyd (HHL) for solving linear systems of equations [11] was devised recently, which provides an exponential speedup over the best known classical algorithms, and has led to further remarkable developments [12-14]. Moreover, the HHL algorithm has been experimentally realized for some simple instances [15-18]. An alternative protocol for solving systems of linear algebraic equations with particular realization on superconducting the quantum processor of the IBM Quantum Experience was proposed in [19]. After the advent of the HHL algorithm, more efficient quantum algorithms and many substantial results for matrix operations have been obtained [20-24]. In particular, Zhao et al [21] have presented a novel method for performing elementary linear algebraic operations on a quantum computer for complex matrices. To simulate these matrix operations, they constructed a Quantum Matrix Algebra Toolbox, a set of unitary matrices, and used the Trotter product formula together with a phase estimation circuit. These results, added with quantum computation, can be used for solving many problems in machine learning and data processing. The scalar product of arbitrary vectors has also been investigated [21, 22] using an ancilla and Hadamard operator. Stolze and Zenchuk, using a 'sender-receiver' model, proposed an alternative scheme for the scalar product via a two-terminal quantum transmission line [25]. They encoded the given vectors into the pure initial states of two different senders and obtained the result, as an element of the two-qubit receiver's density matrix, after time evolution and a proper unitary transformation. In the case that the model does not consider time evolution, design of the unitary transformation becomes very significant for the whole quantum computation.
In this paper, based on the 'sender-receiver' model, we present quantum algorithms for matrix operations such as the matrix-vector product, the matrix-matrix product, the sum of two matrices, and the calculation of determinant and inverse of a matrix. We stress here that the protocols proposed in this paper do not use Trotterization which is a basic tool of algebraic calculations in [21] and is an operation-consuming process requiring multiple repetition of a set of certain operators. Instead, for each matrix operation, we prepare the unitary transformation of a quantum system which is universal for a given matrix dimension and does not depend on particular matrices to be operated. We recall that known quantum algorithms of matrix inversion [26, 27] include classical arithmetics for inverting eigenvalues and therefore, they are not completely quantum algorithms. Here, we propose a truly quantum protocol for matrix inversion based on its definition in terms of algebraic complements. This is to avoid inverting eigenvalues of a matrix. Further to this, using these quantum algorithms, we propose an alternative quantum scheme for solving linear systems of equations. The proposed quantum algorithms are interesting quantum schemes for matrix operations, especially for solving linear systems of equations. In each case, the unitary operation W that acts on the initial quantum state ρ(0) of the whole system, i.e. ρ = Wρ(0)W† is physically realizable. For instance, the multiple-quantum NMR technology can be used to achieve these operations. Regarding the running time, except for the calculation of the determinant and inverse of a matrix, the unitary transformation W, according to the Solovay-Kitaev theorem [9, 28], can be approximated by the standard set of universal gates in polynomial time.
The paper is organized as follows. In section 2, we describe the sender-receiver model which serves as a general structural scheme for various matrix operations. In section 3, the detailed protocols for matrix operations including matrix-vector product, matrix-matrix product, matrix-matrix sum, calculation of matrix determinant and matrix inverse and solving systems of linear equations are proposed. A brief complexity analysis is presented in section 4, and the conclusion follows in section 5. Some additional calculations are presented in appendix A.
2. Sender-receiver model
In this model, the N-node quantum system S consists of two or more subsystems called senders. We denote the ith sender by Si, so that $S\,:=\,{\bigcup }_{i=1}^{n}{S}_{i}$. A node is actually a qubit and ${N}^{({S}_{i})}$ is the number of qubits with the sender Si, so that $N={\sum }_{i=1}^{n}{N}^{({S}_{i})}$. The receiver (denoted by R) doesn't involve any additional node into the system S; it consists of some nodes of the senders as shown schematically by the big dashed circle in figure 1(a). We denote the ith sender without the receiver's nodes by $S{{\prime} }_{i}$ and join all $S{{\prime} }_{i}$ into the subsystem $S^{\prime} \,:=\,{\bigcup }_{i=1}^{n}{S}_{i}^{\prime} $. We also employ the multi-index notation [29] as follows. Let X be a string of N-entries 0 and 1 indicating, respectively, ground or excited state of the corresponding spin. Let ∣X∣ be the sum of the entries of X which equals to the number of excited spins X (we call it the excitation number). In all our protocols, each sender carries a single excitation, i.e. the excitation number equals to the number of senders n. The multi-index MX is assotiated with the appropriate subsystem X of N(X) nodes. We conceive different 'sender-receiver' models for different matrix operations (see figure 1).
Figure 1. The Sender (S) - Receiver (R) model for various matrix operations: (a) General schematic. (b) Matrix-vector product. S1 has m rows of k-qubits and S2 has only one k-qubit row. R has $\left(m+1\right)$ qubits where the first m qubits are part of S1 and the last one belongs to S2. (c) Matrix-matrix product. S1 and S2 have m and n rows of k-qubits, respectively. R has $\left(m+n\right)$ qubits where the first m qubits belong to S1 and the remaining n qubits belong to S2. (d) Sum of two matrices. Both S1 and S2 have m rows. The first row of each sender has n + 1 qubits, and the remaining m − 1 rows have n qubits. R includes $\left(m+n\right)$ qubits where the first m qubits belong to the last column of S1 and the rest n qubits belong to the last row of S2. The columns in S1 are enumerated from left to right, while the columns in S2 are enumerated from right to left. (e) Determinant of an n × n square matrix. There are n senders, and each sender Si is a single n-qubit row. R includes the last qubit of every Si. (f) Inverse of a non-degenerate n × n square matrix. There are n senders, and each sender Si is an (n + 1)-qubit row. R includes the last two qubits of every Si.
For the matrix-vector product, the matrix-matrix product, and the sum of two matrices, all entries of a particular matrix (vector) are allocated to a single sender. However, for the determinant and inverse of a square matrix, we allocate elements of each row of the matrix to separate senders. Thus, there are n senders for an n × n matrix. Moreover, for calculating the inverse matrix, n auxiliary qubits are added to the model which renders assistance for reserving n2 algebraic complements.
In our treatment, matrices and vectors are considered in the complex number field. The initial states of senders are all single-excitation pure normalized states ∣ψi〉, and the elements of matrix A (a vector is also a matrix) in a particular matrix operation are encoded as the probability amplitudes in those pure states. Hence, ${\parallel A\parallel }_{2}^{2}\lt 1$ (see6(6If a matrix M does not satisfy the condition ${\parallel M\parallel }_{2}^{2}\lt 1$, it can be multiplied by a small enough number on a classical computer so that the Frobenius norm of the resulting matrix is less than one. Without the loss of generality, all matrices are constrained by this requirement.)), where ${\parallel A\parallel }_{2}:=\sqrt{{\rm{Tr}}({A}^{\dagger }A)}$ is the Frobenius norm. The initial state of the whole system is
Note that it is the structure of the initial state that provides us all the necessary multiplications of the matrix elements placing the products into the appropriate position in the density matrix of a quantum state. All we need to get the final result is to calculate the sum of the certain elements of ρ(0) (maybe multiplied by ±1 in calculating the determinant and inverse matrix) and move the result into the desired position in the receiver's density matrix. This can be reached after evolving ρ(0) by an appropriate unitary transform W, ρ := Wρ(0)W†, with subsequent tracing over $S^{\prime} $: ${\rho }^{(R)}:={\mathrm{Tr}}_{S^{\prime} }\rho $. The unitary transform W satisfies the commutation relation
where ${I}_{z}={\sum }_{i=1}^{N}{I}_{i,z}$ is the z-projection operator of the total spin-momentum with every Ii,z having eigenvalues ±1/2. This commutational restriction on W is an important condition for our algorithms. It ensures that W is block-diagonal with respect to the excitation number, and that instead of the whole state space we only need to consider some subspace constrained by a certain excitation number. If the maximal number of excitations is n then W = diag(W(0), W(1), ⋯ ,W(n)). Thus, W promises that the result can be found in a certain block of the receiver's density matrix ρ(R). At that, the highest-order (by absolute value) coherence matrices [30],7(7The density matrix of an N-spin quantum system can be written as a sum [30] $\rho ={\sum }_{n=-N}^{N}{\rho }^{(n)}$, where each submatrix ρ(n) consists of the elements of ρ responsible for spin-state transitions that change the number of excited spins by n. This number is positive, negative or zero respectively when the number of excited spins increases, decreases or remains the same after such a transition. Matrices ρ(n) are called multi-quantum coherence matrices and the number n is called the order of coherence matrices. The above representation is especially evident if we go to the basis of N-qubit quantum states sorted by the excitation number.) included into ρ(R) are of the orders ±n, which we denote by ρ(R;±n). Of these, the (−n)-order coherence matrix is of our significance. The (−n)-order coherence matrix of the receiver's state, in the multi-index notation, can be written as
We can simplify equation (3) as follows. Since n =∣MR∣ − ∣NR∣ and $\max | {M}_{R}| =\max | {N}_{R}| =n$, then ∣NR∣ = n and ∣MR∣ = 0. Moreover, due to the commutation relation equation (2), ${\sum }_{i=1}^{m+1}| N{{\prime} }_{i}| +| {N}_{R}| =| I| =n$ and ${\sum }_{i=1}^{m+1}| N{{\prime} }_{i}| +| {M}_{R}| =| J| =0$. Therefore $| {J}_{i}| =| N{{\prime} }_{i}| =| {M}_{R}| =0$. Finally, since ${W}_{{0}_{S};{0}^{{\prime} }{0}_{R}}\equiv {W}_{{0}_{S};{0}_{S}}$ is the only element of the 0-excitation block of W, we set ${W}_{{0}_{S};{0}_{S}}=1$. Hereafter 0A means the multi-index of zeros assotiated with the subsystem A, $0{{\prime} }_{{S}_{i}}$ is the multi-index of zeros assotiated with the subsystem ${S}_{i}^{\prime} $, and $0^{\prime} =0{{\prime} }_{{S}_{1}}\cdots 0{{\prime} }_{{S}_{n}}$. All in all, equation (3) reduces to
This is the general formula for the elements of the (−n)-order coherence matrix of the receiver. Our aim is to find the appropriate unitary transform for every matrix operation mentioned earlier, and present the corresponding receiver's state in a convenient form.
2.1. Pure-state analysis
The density-matrix formulation is useful for theoretical study of the matrix-operation protocols. For practical computation, however, we turn to pure-state analysis. Let ∣ψi〉 be the pure initial state of the ith sender, i = 1,…,n. Then, after applying W, we can write the pure state of the whole system as follows:
where $| {{\rm{\Phi }}}_{i}^{(n)}{\rangle }_{R}$, i = 1,…,DΦ, is the set of n-excitation basis elements of the receiver's state space (not all n-excitation basis elements of the receiver's state space are required in general), the amplitudes Ai, i = 1,…,DΦ, carry the useful information (result) of a matrix operation and ${\sum }_{i=1}^{{D}_{{\rm{\Phi }}}}| {A}_{i}{| }^{2}+| a{| }^{2}=1$ is the normalization. The second term in equation (5) is not of our interest and includes other states of the quantum system,
Now we can adopt one of the following two strategies.
1.Strategy of amplitudes. It is the probabilistic method of calculating amplitudes Ai, see figure 2(a). Performing measurements on the receiver qubits with outputs $| {{\rm{\Phi }}}_{i}^{(n)}{\rangle }_{R}$, we obtain the quantities ∣Ai∣2 as the probabilities of measuring $| {{\rm{\Phi }}}_{i}^{(n)}{\rangle }_{R}$.
2.Strategy of state. It is the probabilistic method of obtaining the output state of the receiver R with the amplitudes Ai encoded into the n-excitation receiver state, see figure 2(b). Performing measurements on qubits of the subsystem $S^{\prime} $ with output $| 0^{\prime} {\rangle }_{S^{\prime} }$ results in the following state of the receiver R:
Figure 2. The schemes for probabilistic realization of matrix operations. (a) Probabilistic method of calculating the amplitudes Ai. (b) Obtaining the output state of the receiver with Ai encoded into the (part of) n-excitation states, equation (7).
3. Matrix operations
In this section, different strategies are adopted for different matrix operations to design the unitary transformations.
3.1. Matrix-vector product
We encode the elements of matrix $A={\left({a}_{{ij}}\right)}_{m\,\times \,k}$ and vector ${\boldsymbol{v}}={\left({v}_{j}\right)}_{k\times 1}$ as the probability amplitudes into the pure normalized states
where $| {{\rm{\Psi }}}_{{S}_{1}}^{({ij})}\rangle \equiv | 1{\rangle }_{{ij}}\otimes | 0{\rangle }^{\otimes {mk}-1}$ is a single excited spin state of S1 in the ith row and jth column, and $| {{\rm{\Psi }}}_{{S}_{2}}^{(j)}\rangle \equiv | 1{\rangle }_{j}\otimes | 0{\rangle }^{\otimes k-1}$ is the state of S2 having one excited spin at the jth position.
Let ${\rho }_{{N}_{R};{0}_{R}}^{(R;-2)}$ be the elements of the (−2)-order coherence matrix. Expression (4) with n = 2 reduces to
Here 0X is the multi-index of zeros related to the subsystem X. Since the receiver R is an (m + 1)-qubit subsystem, its (−2)-order coherence matrix ρ(R;−2) has $\left(\displaystyle \genfrac{}{}{0em}{}{m+1}{2}\right)$ elements (two excitations out of m + 1 receiver spins). But we only consider the matrix elements with index ${N}_{{R}_{i}}=\{{\underbrace{0\cdots 0}}_{i-1}1{\underbrace{0\cdots 0}}_{m-i}1\}$ with 1 at the ith and (m + 1)th positions. Equation (12) for these elements yields
where for the sake of notational convenience the 0-parts in W are omitted here, and also in later discussions unless stated otherwise. That is ${W}_{{N}_{{R}_{i}};{I}_{{S}_{1}}{I}_{{S}_{2}}}^{(2)}\equiv {W}_{{0}_{{S}_{1}}^{{\prime} }{0}_{{S}_{2}}^{{\prime} }{N}_{{R}_{i}};{I}_{{S}_{1}}{I}_{{S}_{2}}}^{(2)}.$ Let
where ${\theta }_{1}=\tfrac{1}{\sqrt{k}}$. This involves m rows of W(2), and the other rows must fulfill the unitarity. Then with ${\rho }_{{I}_{{S}_{1}}^{({ij})};{0}_{{S}_{1}}}^{({S}_{1})}={a}_{{ij}}{a}_{00}^{* }$, ${\rho }_{{I}_{{S}_{2}}^{(j)};{0}_{{S}_{2}}}^{({S}_{2})}={v}_{j}{v}_{0}^{* }$ and $\alpha ={\theta }_{1}{a}_{00}^{* }{v}_{0}^{* }$, we have ${\rho }_{{N}_{{R}_{i}};{0}_{R}}^{(R;-2)}\,=\alpha {\sum }_{j=1}^{k}{a}_{{ij}}{v}_{j}$. Hence, the matrix-vector product is
Pure-state analysis of matrix-vector product.-Here, we provide the pure-state analysis as discussed in section 2. With DΦ = m, ${A}_{i}={\theta }_{1}{\sum }_{j=1}^{k}{a}_{{ij}}{v}_{j}$ and the initial state ∣ψ1〉 ⨂ ∣ψ2〉, equation (5) takes the form
To register the result via the probabilistic method, figure 2(a), we perform measurements on the receiver qubits with the output $| {{\rm{\Phi }}}_{i}^{(2)}{\rangle }_{R}$ and obtain the quantity ${\theta }_{1}^{2}| {(A{\boldsymbol{v}})}_{i}{| }^{2}={\theta }_{1}^{2}| {\sum }_{j\,=\,1}^{K}{a}_{{ij}}{v}_{j}{| }^{2}$, which is proportional to the square of the absolute value of the ith element of Av.
To produce the output state of the receiver R with the vector Av encoded into the 2-excitation state, figure 2(b), we perform the measurements on the qubits of the subsystem $S^{\prime} $ with the output $| 0^{\prime} {\rangle }_{S^{\prime} }$. Then the receiver's state equation (7) takes the form:
where γ is the normalization constant determined in equation (9).
3.2. Matrix-matrix product
We consider the matrix-matrix product AB of two matrices $A={({a}_{{ij}})}_{m\times k}$ and $B={({b}_{{ij}})}_{k\times s}$, see figure 1(c). The initial states of the two senders are
Here $| {{\rm{\Psi }}}_{{S}_{1}}^{({ij})}\rangle \equiv | 1{\rangle }_{{ij}}\otimes | 0{\rangle }^{\otimes {mk}-1}$ and $| {{\rm{\Psi }}}_{{S}_{2}}^{({ij})}\rangle \equiv | 1{\rangle }_{{ij}}\otimes | 0{\rangle }^{\otimes {sk}-1}$ are the single-excitation states of S1 and S2 respectively. Note that the transpose of B is encoded into the quantum state ∣ψ2〉 of S2. That is, the probability amplitude of $| {{\rm{\Psi }}}_{{S}_{2}}^{({ij})}\rangle $ is bji (rather than bij) to satisfy the multiplication rule. In the following, we carry out similar treatment as in the matrix-vector product. We consider the multi-index NR in which the entries corresponding to the ith and the (m + j)th spin, respectively, in the first and second columns of R are 1, and others are 0. That is,
and with ${\rho }_{{I}_{{S}_{1}}^{({il})};{0}_{{S}_{1}}}^{({S}_{1})}={a}_{{il}}{a}_{00}^{* }$, ${\rho }_{{I}_{{S}_{2}}^{({jl})};{0}_{{S}_{2}}}^{({S}_{2})}={b}_{{lj}}{b}_{00}^{* }$ and $\beta ={\theta }_{1}{a}_{00}^{* }{b}_{00}^{* }$, we have ${\rho }_{{N}_{{R}_{{ij}}};{0}_{R}}^{(R;-2)}=\beta {\sum }_{l=1}^{k}{a}_{{il}}{b}_{{lj}}$. Thus, the matrix product is
Pure-state analysis of the matrix-matrix product.-With DΦ = ms, ${A}_{{il}}={\theta }_{1}{\sum }_{j=1}^{k}{a}_{{ij}}{b}_{{jl}}$ and the initial state ∣ψ1〉 ⨂ ∣ψ2〉, equation (5) can be written as
To register the elements of the matrix AB via the probabilistic method, figure 2(a), we perform the measurements on the receiver qubits with the output $| {{\rm{\Phi }}}_{{il}}^{(2)}{\rangle }_{R}$ to obtain the quantity ${\theta }_{1}^{2}| {({AB})}_{{il}}{| }^{2}={\theta }_{1}^{2}| {\sum }_{j=1}^{k}{a}_{{ij}}{b}_{{jl}}{| }^{2}$, which is proportional to the square of the absolute value of the (il)th element of AB.
To obtain the output state of the receiver R with the matrix AB encoded into the 2-excitation state, figure 2(b), we perform the measurements on the qubits of the subsystem $S^{\prime} $ with the output $| 0{\rangle }_{S^{\prime} }$. Then equation (7) yields the following state of the receiver R:
where γ is the normalization constant determined in equation (9).
3.3. Sum of matrices
Consider two m × s matrices C and D. The first row of each sender has s + 1 elements, while the remaining rows have n elements. The initial states of two subsystems read
where c00 ≠ 0, d00 ≠ 0, λ ≠ 0 and the extra node in the first row of each sender is marked by 0. Here $| {{\rm{\Lambda }}}_{{S}_{k}}^{({ij})}\rangle =| 1{\rangle }_{{ij}}\otimes | 0{\rangle }^{\otimes {mn}}$, k = 1, 2. The receiver consists of the last column of S1 and the last row of S2 as seen in figure 1(d).
Here, elements of the (−2)-order coherence matrix of the receiver density matrix ρ(R) can be written as ${\rho }_{{N}_{R};{0}_{R}}^{(R;-2)}={\sum }_{{I}_{{S}_{1}}{I}_{{S}_{2}}}{W}_{0{{\prime} }_{{S}_{1}}0{{\prime} }_{{S}_{2}}{N}_{R};{I}_{{S}_{1}}{I}_{{S}_{2}}}^{(2)}{\rho }_{{I}_{{S}_{1}};{0}_{{S}_{1}}}^{({S}_{1})}{\rho }_{{I}_{{S}_{2}};{0}_{{S}_{2}}}^{({S}_{2})}$. We choose elements of the receiver density matrix ρ(R) marked by the index ${N}_{{R}_{{ij}}}$ in equation (18) corresponding to the ith excited spin of the column and jth excited spin of the row of the receiver. Hence,
There are ms such elements which coincide with the total number of elements in either matrix C or D. Therefore, the elements ${\rho }_{{N}_{{R}_{{ij}}},{0}_{R}}^{(R;-2)}$ can be used to store the sum C + D. Let
is proportional to the sum of two proper elements of C and D, where ${I}_{{S}_{1}}^{(10)}$ and ${I}_{{S}_{2}}^{(10)}$ are extra elements of the 1st row of S1 and S2, respectively, from the figure 1(d). Hence,
Pure-state analysis of sum of two matrices.-Equation (5), with DΦ = ms, Aij = λθ2(aij + bij) and the pure initial state ∣ψ1〉 ⨂ ∣ψ2〉, yields
Using the probabilistic method for calculating elements of the matrix C + D, figure 2(a), we perform measurements on the receiver qubits with the output $| {{\rm{\Phi }}}_{{ij}}^{(2)}{\rangle }_{R}$ and obtain the quantity ${\left({\theta }_{2}\lambda \right)}^{2}| {a}_{{ij}}+{b}_{{ij}}{| }^{2}$, which is proportional to the square of the absolute value of the (ij) th element of C + D.
To obtain the output state of the receiver R with the matrix C + D encoded into the 2-excitation state, figure 2(b), we perform measurements on the qubits of the subsystem $S^{\prime} $ with output $| 0^{\prime} {\rangle }_{S^{\prime} }$. Equation (7) yields the following state of the receiver R:
where γ is the normalization constant determined in equation (9).
3.4. Determinant
To compute the determinant of a square matrix $E={\left({e}_{{ij}}\right)}_{n\times n}$, we encode every ith row of E into the pure initial state of the senders Si, i = 1,…,n,
where $| {{\rm{\Phi }}}_{{S}_{i}}^{(j)}\rangle \equiv | 1{\rangle }_{j}\otimes | 0{\rangle }^{\otimes n-1}$ represents the state of Si with only one excited spin at the jth position.
Here, contrary to the previous investigations, there are n senders and one receiver. Consequently, a slight modification in the treatment is in order.
We only consider NR = {11 ⋯ 1} and MR = {00 ⋯ 0}, so that ${\rho }_{{N}_{R};{M}_{R}}^{(R)}\equiv {\rho }_{{1}_{R};{0}_{R}}^{(R)}$ is exactly the unique element of the receiver's (−n)-order coherence matrix ρ(R;−n) (in the rest of this paper we adopt notation ${X}_{i}\equiv {X}_{{S}_{i}}$, for any index X),
where ${\theta }_{3}=\tfrac{1}{\sqrt{n!}}$ and ${\epsilon }_{{i}_{1}{i}_{2}\cdots {i}_{n}}$ is the permutation symbol which equals 1 (if the inversion number of the permutation i1i2 ⋯ in is even), −1 (if the inversion number of the permutation i1i2 ⋯ in is odd) or 0 (if at least two indices coincide). Then with ${\rho }_{{I}_{l}^{({i}_{l})};{0}_{l}}^{({S}_{l})}={e}_{{{li}}_{l}}{e}_{l0}^{* }$ and $\gamma ={\theta }_{3}{\prod }_{i=1}^{n}{{e}_{i0}}^{* }$,
Using the probabilistic method for calculating the determinant of the matrix E, figure 2(a), we perform measurements on the receiver qubits with the output ∣Φ(n)〉R and obtain the quantity ${\theta }_{3}^{2}| {\rm{\det }}(E){| }^{2}$, which is proportional to the square of the determinant of E.
To produce the output state of the receiver R with the ${\rm{\det }}(E)$ encoded into the n-excitation state, figure 2(b), we perform measurements on the qubits of the subsystem $S^{\prime} $ with the output $| 0^{\prime} {\rangle }_{S^{\prime} }$. Equation (7) yields the following state of the receiver R:
where γ is the normalization constant determined in equation (9).
3.5. Inverse of matrix
Like the determinant, to compute the inverse of a non-degenerate square matrix $E={({e}_{{ij}})}_{n\times n}$, we encode every ith row of E into the pure state of sender Si. Moreover, unlike previous 'sender-receiver' models, we add the auxiliary module called Aux to the senders (the last qubits of the senders). These auxiliary qubits appear in R which includes the two last spins of each sender. See figure 1(f). The normalized initial state of each sender Si is given by
where ${\hat{e}}_{i0}\ne 0,\,\sigma \ne 0$, ${\hat{e}}_{{ij}}$ are the elements of the non-degenerate matrix E whose inverse is intended, and $| {{\rm{\Xi }}}_{{S}_{i}}^{(j)}\rangle \equiv | 1{\rangle }_{j}\otimes | 0{\rangle }^{\otimes n}$ is the state of sender Si corresponding to the jth excited spin (qubit). In particular, $| {{\rm{\Xi }}}_{{S}_{i}}^{(n+1)}\rangle $ is the state associated with the excited spin of Si included in the Aux module.
We have already obtained the determinant. To find the inverse ${E}^{-1}={E}^{* }/\det (E)$, where E* is the adjoint matrix, our goal is to find every algebraic complement Eij of eij (the matrix element of E).
We consider the (−n)-order coherence matrix with elements ${\rho }_{{N}_{R};{0}_{R}}^{(R;-n)}$ and select n2 elements as follows. The receiver R can be viewed as a two-column subsystem. Let the entries of the multi-index NR associated with the first column (the nth spins of each sender) be 1 except the ith entry which is 0. Similarly, let the entries of the multi-index NR associated with the second column (i.e. the entries corresponding to the module Aux) be 0 except the jth one which is 1. We denote such multi-index by ${\hat{N}}_{{R}_{{ij}}}$, i, j = 1,…,n. Thus, there are n2 different elements in selected part of ρ(R;−n). These are used for storing the elements of the inverse matrix E−1.
The elements ${\rho }_{{\hat{N}}_{{R}_{{ij}}};{0}_{R}}^{(R;-n)}$ are
where ${\theta }_{4}=\tfrac{1}{\sqrt{(n-1)!}}$ and ${\epsilon }_{j;{l}_{1}\cdots {l}_{i-1}{l}_{i+1}\cdots {l}_{n}}$ is +1 (− 1) if {l1 ⋯ li−1li+1 ⋯ ln} is an even (odd) permutation of {1, ⋯ ,j − 1, j + 1, ⋯ n}, and 0 otherwise. Thus, we obtain
We can also find the determinant of E from ρ(R) by continuing to design the unitary transformation W. Let ${\hat{N}}_{R}=\{1010\cdots 10\}$ be the element ${\rho }_{{\hat{N}}_{R};{0}_{R}}^{(R;-n)}$ which corresponds to the states ∣0〉 of all spins from Aux and states ∣1〉 of all spins from the first column of R. Let the elements of W(n) satisfy the constraint:
Pure-state analysis of matrix inverse.-Equation (5), with DΦ = n2 + 1 and the pure initial state ∣ψ1〉 ⨂ ⋯ ⨂ ∣ψn〉, reads
where the subscript last means the set of nth qubits of each sender and the subscript aux means the set of auxiliary qubits of each sender.
To register the inverse of the matrix E via the probabilistic method, figure 2(a), we perform measurements on the receiver qubits with the output $| {{\rm{\Phi }}}_{{ij}}^{(n)}{\rangle }_{R}$ to obtain either the quantities $| {\theta }_{4}\sigma \det (E){E}_{{ji}}^{-1}{| }^{2}$ proportional to the square of elements of E−1 (if (ij) ≠ (00)) or ${\theta }_{3}^{2}{\det }^{2}(E)$ proportional to the square of the determinant (if (ij) = (00)).
To form the output state of the receiver R with the elements of E−1 encoded into the n-excitation state, figure 2(b), we perform measurements on the qubits of the subsystem $S^{\prime} $ with the output $| 0^{\prime} {\rangle }_{S^{\prime} }$. Then equation (7) yields the following state of the receiver R:
where γ is the normalization constant determined in equation (9).
3.6. Linear system of equations
We apply the quantum algorithms investigated in this paper to solve the linear system of equations, Ex = b, where $E\,={({e}_{{ij}})}_{n\times n}$ is a non-degenerate matrix and ${\boldsymbol{b}}={\left({b}_{1},{b}_{2},\cdots ,{b}_{n}\right)}^{{\rm{T}}}$ is a unit vector. To obtain the solution x = E−1b, we first construct the unitary transformation V satisfying the condition (2) by encoding b into V, and then apply it to the receiver's state ρ(R) of the inverse matrix operation: ξ(R) = Vρ(R)V†. Let
where ${L}_{{R}_{j}}={\hat{N}}_{{R}_{{jj}}}=\{{\underbrace{10\cdots 10}}_{2(j-1)}01{\underbrace{10\cdots 10}}_{2(n-j)}\}$, and other entries need to fulfill the unitarity of V. Now with ${V}_{{0}_{R};{0}_{R}}^{\dagger }=1$ and using equation (41), we have
Pure-state analysis of linear system of equations.-Equation (5), with DΦ = n and pure initial state ∣ψ1〉 ⨂ ⋯ ⨂ ∣ψn〉 after replacing W with $({I}_{S^{\prime} }\otimes V)W$, reads:
where $| {{\rm{\Phi }}}_{{ii}}^{(n)}{\rangle }_{R}$ and $| {{\rm{\Phi }}}_{00}^{(n)}{\rangle }_{R}$ are determined using equation (42).
Using the probabilistic method for registration of the elements of the inverse of the matrix E, figure 2(a), we perform measurements on the receiver qubits with the output $| {{\rm{\Phi }}}_{{ii}}^{(n)}{\rangle }_{R}$ to obtain the quantities $| {\theta }_{4}\sigma {\rm{\det }}(E){x}_{i}{| }^{2}$, which are proportional to ∣xi∣2.
To obtain the output state of the receiver R with xi encoded into the n-excitation state, figure 2(b), we perform measurements on the qubits of the subsystem $S^{\prime} $ with the output $| 0^{\prime} {\rangle }_{S^{\prime} }$. Then, equation (7) yields the following state of the receiver R:
where γ is the normalization constant determined in equation (9).
For a summary of the proposed quantum algorithms and the illustrations of the quantum algorithms for the determinant and solving linear systems of equations, see Table B1 in Appendix B.
4. Complexity analysis and realization of W
In this section, we consider some features of our protocols. First, due to the tensor-product structure of the initial state of the sender S, all products between the elements of the matrices in the algorithms of the matrix-vector product and the matrix-matrix product, computation of the determinant and the inverse of matrix are obtained automatically and appear in the (−n)-order coherence matrix of the sender S (here n is the excitation number). Thus, to obtain the result, we have to perform the necessary number of addition and maybe multiplication by −1. This is realized by the block W(n) of the unitary transformation W. The subsequent partial trace over $S^{\prime} $ serves just to locate the result in the density matrix of the receiver R; this operation does not perform any algebraic manipulation.
The fact that multiplication of the matrix elements is performed automatically without special efforts reduces the number of algebraic operations imposed on the block W(n) of the unitary transformation W. However, there is a significant problem in the realization of such a block. Apparently, the representation of W(n) in terms of the two-level matrices [9] creates an obstacle in harnessing the above advantage. In fact, the number of 'active' elements of W(n) involved into the algebraic operation is predicted by two parameters: (i) the number of n-excitation basis elements in the initial senders state,
(this number is determined by the number of senders and equals ${K}_{1}={\prod }_{i=1}^{n}{N}^{({S}_{i})};$ do not mistake it for the total number of n-excitation basis elements which is $\left.\left(\displaystyle \genfrac{}{}{0em}{}{N}{n}\right)\right)$ and (ii) the number of entries of the receiver density matrix where the result of a matrix operation is finally placed. The latter fixes the number K2 of the involved rows of W(n), while the former fixes the number of involved entries of each above row. Then the number of required two-level unitary operations [9] is
For illustration, let us consider the calculation of the inverse of an n × n matrix. In this case, we need K2 = n2 rows of W(n) corresponding to n excited spins of the 2n-qubit receiver, according to equation (35) of section 3.5. Each of these rows has K1 = n! nonzero entries which serve to combine all necessary products of the elements of the matrix E whose inverse is to be calculated. The result is placed in the probability amplitudes of certain n-excitation receiver states. Thus, all in all, we only need n!n2 exchange-operations instead of n!n3 operations of a classical algorithm. This is an advantage of a quantum algorithm.
However, the n-excitation subspace of the state-space of the sender S includes ${K}_{1}={\left(n+1\right)}^{n}$ basis elements ∣i1〉 ⨂ ... ⨂ ∣in〉 (here ∣ik〉 means ikth qubit excited in (n + 1)-qubit sender Sk), and we need K2 = n2 rows of W(n). Equation (50) then yields $K={n}^{2}(2{\left(n+1\right)}^{n}-{n}^{2}\,+\,1)/2$. This is the number of two-level matrices required to encode such entries. The values of the other entries in those rows are not important for us. We have listed values of K1 and K2 for various matrix operations in table 1.
Table 1. Parameters K1 and K2 for the matrix-vector product (Av), the matrix-matrix product (AB), the sum of matrices (C + D), the determinant ($\det (E)$), and the inverse of matrix (E−1).
Av
AB
C + D
$\det (E)$
E−1
K1
mk
msk2
${\left({ms}+1\right)}^{2}$
nn
${\left(n+1\right)}^{n}$
K2
k
ms
mn
1
n2
Thus, the proposed unitary transformation W(n) is an example of a multiqubit transformation which cannot be effectively written in terms of elementary unitary transformations. Therefore, the development of special techniques for generating such transformations becomes relevant. For instance, Floquet engineering of Hamiltonians might be effective [31-33].
4.1. Unitary transformation via Hamiltonian evolution
In [19], the evolution of the 4-qubit chain under XX-Hamiltonian was used to generate the proper unitary transformation for solving a system of linear equations. We build on this idea and introduce the Hamiltonian H that conserves the excitation number in the system. The XX-Hamiltonian can serve this purpose:
where Iαi is an α-projection of the ith spin momentum, α = x, y, z, and Dij are coupling constants. In the case of dipole-dipole interaction in the strong magnetic field directed along z-axis, we have
where φij is the angle between vector ${\vec{r}}_{{ij}}$ and the magnetic field, γ is the gyromagnetic ratio and ℏ is the Planck's constant. Then, at some fixed time instant t0, we can assign
Notice that there are N(N − 1)/2 different coupling constants Dij parameterized by three N coordinates xj, yj and zj. Of course, this number of free parameters is insufficient to parameterize all needed entries of W in the protocol of each matrix algorithm of sections 3.1-3.6. Partially, this obstacle can be overcome by taking into consideration the set of additional nodes as qubits of the extended receiver [29]. This allows us to increase the number of coupling constants (and consequently the number of free parameters) in the quantum system. See the construction of W(2) in equation (53) for calculating the determinant of a 2 × 2 matrix in Appendix A and the linear system of 2 equations with 2 variables in Appendix B.
5. Conclusions
Based on the 'sender-receiver' model, we have constructed the quantum computational framework for calculating the matrix-vector product, the matrix-matrix product, the sum of two matrices, the determinant and the inverse of a matrix over a complex field. By encoding the information of the matrices into the pure initial states of the senders and performing the appropriate unitary transformation W, the final results are elements of the receiver's density matrix which, in turn, can be considered as input for other quantum algorithms. The basic feature of our algorithms is that they don't use Trotterization. In addition, the algorithm for calculating the inverse matrix does not use inversion of eigenvalues λi, i.e. transformation λi → 1/λi. Therefore, the construction of the inverse matrix is completely quantum in nature.
Moreover, the unitary transformation W in each protocol is a universal transformation, i.e. it can be used to perform the appropriate operation with any matrices. Finally, as an application of the proposed quantum algorithms, in particular, the inverse of a matrix, we have presented a quantum scheme for solving linear systems of equations different from the HHL algorithm. The proposed algorithms also provide new insights to develop other quantum algorithms. Although the representation of W in terms of elementary transformations is not effective, the development of methods for constructing appropriate evolution Hamiltonians might overcome this obstacle. The construction of unitary transformation via Hamiltonian evolution is discussed.
In future work, we intend to study quantum algorithms implementing the matrix elementary transformation, as it is a significant tool for studying fundamental matrix algebra problems, which can complete the calculation of the determinant and inverse matrix from a different standpoint within this paper. Further, declining the number of initial qubits and improving the efficiency of quantum algorithms are also valuable research problems that need to be addressed.
Acknowledgments
This project is supported by the National Natural Science Foundation of China (Grant No. 12031004 and Grant No. 12271474, 61877054), the Fundamental Research Foundation for the Central Universities (Project No. K20210337) and the Zhejiang University Global Partnership Fund, 188170+194452119/003. The work was partially funded by a state task of Russian Fundamental Investigations (State Registration No. FFSG-2024-0002).
Appendix A. Illustration for obtaining the determinant
For the square matrix $E=\left(\begin{array}{cc}3/4 & 1/4\\ 1/4 & 3/4\end{array}\right)$, the pure initial states of two senders can be written as
Let ${W}_{0101;1001}=\tfrac{1}{\sqrt{2}}$ and ${W}_{0101;0110}=-\tfrac{1}{\sqrt{2}}$, where $\sqrt{2}$ is designed to ensure that W is a unitary matrix. The commutation relation [W, Iz] = 0 shows that the W is block-diagonal with respect to the number of excitations. Since the 4-qubit system contains no more than four excitations, the unitary transformation W includes five blocks: W = diag(W(0), W(1), W(2), W(3), W(4)). We let W(0) = W(4) = 1 and W(1) = W(3) = I4. The block of major concern here is W(2) written in the basis {permut ∣0⨂21⨂2〉} = {∣0011〉, ∣0101〉, ⋯ ,∣1100〉}:
Hence, ${\rho }_{11;00}^{(R;-2)}=\tfrac{3\sqrt{2}}{32}$. Also, ${e}_{10}={e}_{20}=\tfrac{\sqrt{6}}{4}$ and ${\theta }_{3}=\tfrac{1}{\sqrt{2}}$ yield $\gamma ={\theta }_{3}{e}_{10}{e}_{20}=\tfrac{3\sqrt{2}}{16}$. Therefore, the determinant is $\det (E)=\tfrac{1}{\gamma }\,{\rho }_{11;00}^{(R;-2)}=\tfrac{1}{2}.$
W(2) via Hamiltonian evolution
First of all, we remark that the construction of W(2) in equation (A1) is not unique, and we look for W(2) as in equation (53). Only the second row of W(2) is of interest. Thus, using N = 4 in equation (51) and t0 = 40 in equation (53), and requiring the following structure for the second row of W(2),
We fix the parameters ${\hat{e}}_{10}={\hat{e}}_{20}=\tfrac{1}{2}$ and $\sigma =\tfrac{1}{2\sqrt{2}}$, so that the pure initial states of two senders are
The unitary transformation W is block-diagonal, W = diag(1, W(1), W(2), ⋯ ,W(6)), and we choose the two-excitation block W(2) in the basis {permut ∣0⨂41⨂2〉} = {∣000011〉, ∣000101〉, ⋯ ,∣110000〉} of the whole system as:
Also, the unitary transformation V = diag(1, V(1), V(2), V(3), V(4)) is block-diagonal where V2 in the basis {permut ∣0⨂21⨂2〉} = {∣0011〉, ∣0101〉, ⋯ ,∣1100〉} of the receiver is
With unitary transformations W and V, we can find the solution of a linear system of equations. [We can choose W(i) = I and V(j) = I for i, j ≠ 2.].
Table B1. Summary of the quantum algorithms for matrix operations.
SimonD1994 On the power of quantum computation Proceedings, XXXV Annual Symposium on Foundations of Computer Science. Los Alamitos, CA IEEE Press https://ieeexplore.ieee.org/document/365701
4
SimonD R1997 On the power of quantum computation SIAM J. Comput.26 5
ShorP W1994 Algorithms for quantum computation: discrete logarithms and factoring Proceedings, XXXV Annual Symposium on Foundations of Computer Science. Los Alamitos, CA IEEE Press https://ieeexplore.ieee.org/document/365700
6
ShorP W1997 Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer SIAM J. Comput.26 5
NielsenM AChuangI L2000Quantum Computation and Quantum Information Cambridge Cambridge University Press
10
GroverL1996 A fast quantum mechanical algorithm for database search Proceedings, XXVIII Annual ACM Symposium on the Theory of Computation New York ACM Press
BarzSKassalIRingbauerMLippY ODakićBAspuru-GuzikAWaltherP2014 A two-qubit photonic quantum processor and its application to solving systems of linear equations Sci. Rep.4 6115
DoroninS IFeldmanE BZenchukA I2020 Solving systems of linear algebraic equations via unitary transformations on quantum processor of IBM quantum experience Quantum Inf. Proc.19 68
BerryD WChildsA MOstranderAWangG M2017 Quantum algorithm for linear differential equations with exponentially improved dependence on precision Commun. Math. Phys.356 1057
Ta-ShmaA2013 Inverting well conditioned matrices in quantum logspace, STOC '13 Proceedings of the forty-fifth annual ACM symposium on Theory of Computing 881-890
TongYAnDWiebeNLinL2021 Fast inversion, preconditioned quantum linear system solvers, fast Green's-function computation, and fast evaluation of matrix functions Phys. Rev. A104 032422