Matrix Encoding Method in Variational Quantum Singular Value Decomposition

We consider the singular value decomposition of an arbitrary square N × N matrix M assuming N = 2ⁿ, 1M=U^DV^†,M = \hat UD{{\hat V}^{^\dag }}, where D = diag(d₁, …, d_N) is the diagonal matrix of singular values (some of those values might be zero) and Û and V^†{{\hat V}^\dag } are unitary matrices. To find the singular values (entries of the diagonal matrix D) we, first of all, introduce the following objective function [26]: 2L(α,β)=∑j=0N−1qj×Re(〈ψj|UT(α)MU(β)|ψj〉),U={ bij:i,j=0,…,N−1 },L(\alpha ,\beta ) = \sum\limits_{j = 0}^{N - 1} {{q_j}} \times {\mathop{\rm Re}\nolimits} \left( {\langle {\psi _j}|{U^T}(\alpha )MU(\beta )\left| {{\psi _j}} \right\rangle } \right),\quad U = \left\{ {{b_{ij}}:i,j = 0, \ldots ,N - 1} \right\}, where |ψ_j〉, j = 0, …, N − 1, is a set of orthogonal vectors, α = {α₀, …, α_nQ} and β = {β₀, …, β_nQ} represent two sets of optimization parameters, q₀ > ⋯ > q_N−1 are real weights, Q is some integer associated with the subroutine used for preparing the unitary transformation U in Section 2.3. We also note that the sum in (2) is over all N singular values including possible zeros unlike Ref. [26], where the sum is truncated keeping only T ≤ N largest singular values. This truncation can be simply realized in our algorithm just equating to zero all q_j with j > T. The reason to take transposition T of the matrix U in (2) will be clarified letter, see eq. (23). Here we emphasize that, although the objective function is a sum of N terms, its expectation value will be found by measuring the states of two qubis, see eqs. (29), (30), unlike Refs. [26,27], where the expectation value of each term must be measured separately. We appeal to the gradient method to find such parameters α* and β* that maximize the objective function, i.e., 3maxα,βL(α,β)=L(α*,β*).\mathop {\max }\limits_{\alpha ,\beta } L(\alpha ,\beta ) = L\left( {{\alpha ^*},{\beta ^*}} \right).

Then 4dj=〈ψj|UT(α*)MU(β*)|ψj〉,j=0,…,N−1,U^†=UT(α*),V^=U(β*).\matrix{ {{d_j} = \langle {\psi _j}|{U^T}\left( {{\alpha ^*}} \right)MU\left( {{\beta ^*}} \right)\left| {{\psi _j}} \right\rangle ,\quad j = 0, \ldots ,N - 1,} \hfill \cr {{{\hat U}^\dag } = {U^T}\left( {{\alpha ^*}} \right),\quad \hat V = U\left( {{\beta ^*}} \right).} \hfill \cr }

As for the set of orthogonal vectors |ψ_j〉, we take the vectors of computational basis |j〉, j = 0,…, N – 1, where the integer j is written in the binary form. We introduce five n-qubit subsystems, see Figure 1: the subsystems R and C serve to enumerate, respectively, the rows and column of M, the subsystems χ and ψ are needed to operate with 〈ψ_j| and |ψ_j〉 in (2), these two subsystems are also used to organize the weighted sum over j in (2), and the subsystem q encodes the normalized vector of weights, 5|φ〉q=∑j=0N−1qj|j〉q,∑j=0N−1qj2=1.|\varphi {\rangle _q} = \sum\limits_{j = 0}^{N - 1} {{q_j}} |j{\rangle _q},\quad \sum\limits_{j = 0}^{N - 1} {q_j^2} = 1.

Figure 1.

Structure of subsystems required for performing the quantum part of variational SVD. It includes five n-qubit subsystems for encoding the matrix M(subsystems R, C), orthonormal eigenvectors 〈ψ_j| and |ψ_j〉 (subsystems χ, ψ) and weights q_i (subsystem q). In addition, the one-qubit subsystem K serves as a controlling qubit in controlled operators of the algorithm and two one-qubit ancillae B and B˜{\tilde B} serve for the controlled measurement.

In addition, the single qubit of the subsystem K serves as a controlling qubit in succeeding controlled operations. At the last steps of the algorithm we will introduce two one-qubit ancillae B and B˜{\tilde B} to properly organize garbage removal and required measurements.

We shall note that the proposed method can be also used to construct the eigensystem for the square positive semidefinite matrix R (for instance, for the density matrix), i.e., we can factorize R as R = UDU^†. For this purpose, we just have to involve the following relation between the matrices U(α) and U(β): 6U^†=UT(α)=V^†=U†(β), ⇒ U(α)=U¯(β),{{\hat U}^\dag } = {U^T}(\alpha ) = {{\hat V}^\dag } = {U^\dag }(\beta ),\; \Rightarrow \;U(\alpha ) = \bar U(\beta ), where the bar means complex conjugate. This condition can be simply satisfied in the case when all the parameters α_j and β_j are introduced via the x-, y-, or z-rotation R_θj = e^{–σ(θ)α_j/2}, θ = x, y, z, σ^(θ) are the Pauli matrices. If θ = y, then α_j = β_j. In the case θ = x, z, we take α_j = –β_j.

First of all, we have to prepare the above mentioned matrix M = {m_ij : i, j = 0, …, N – 1} for encoding into the state of a quantum system. To this end we normalize M and make real the first diagonal element assuming that this element does not equal to zero (m₀₀ ≠ 0), i.e., replace M with the matrix A = {a_ij : i,j = 0, …, N – 1}: 7A=e−iarg(m00)∑i=0N−1∑j=0N−1| mij |2M.A = {{{e^{ - i\arg \left( {{m_{00}}} \right)}}} \over {\sqrt {\sum\nolimits_{i = 0}^{N - 1} {\sum\nolimits_{j = 0}^{N - 1} {{{\left| {{m_{ij}}} \right|}^2}} } } }}M.

Now we can encode the elements of the matrix A into the superposition state of R and C as follows: 8|A〉=∑i=0N−1∑j=0N−1aij|i〉R|j〉C,∑i=0N−1∑j=0N−1| aij |2=1,\matrix{ {|A\rangle = \sum\limits_{i = 0}^{N - 1} {\sum\limits_{j = 0}^{N - 1} {{a_{ij}}} } |i{\rangle _R}|j{\rangle _C},} \hfill \cr {\sum\limits_{i = 0}^{N - 1} {\sum\limits_{j = 0}^{N - 1} {{{\left| {{a_{ij}}} \right|}^2}} } = 1,} \hfill \cr } where the normalization is provided by eq. (7). In other words, we have quantum access to the matrix A [33]. Here we shall note that the matrix encoding problem is a complicated task by itself and, in general, the depth of the algorithm encoding the N × N matrix is at least O(N²). The same holds for encoding the state |φ〉_q in (5). However, both this problems can be referred to the initial state preparation and will not be detailed in our paper. Subsystems χ, ψ and K are in the ground state initially and the state of q is defined in (5), i.e., the initial state of the whole system reads: 9|Φ0〉=|A〉|0〉χ|0〉ψ|φ〉q|0〉K.\left| {{\Phi _0}} \right\rangle = |A\rangle |0{\rangle _\chi }|0{\rangle _\psi }|\varphi {\rangle _q}|0{\rangle _K}.

Hereafter in this paper the subscript means the subsystem to which the operator is applied.

Now we proceed to description of the quantum algorithm, which is also illustrated by the circuit in Fig.2. As the first step, we apply the Hadamard transformation to each qubit of χ and K (we denote this set of transformations as WχK(0)=HχHKW_{\chi K}^{(0)} = {H_\chi }{H_K}: 10|Φ1〉=WχK(0)|Φ0〉=12(n+1)/2∑k=0N−1|A〉|k〉χ|0〉ψ|φ〉q(|0〉K+|1〉K),\left| {{\Phi _1}} \right\rangle = W_{\chi K}^{(0)}\left| {{\Phi _0}} \right\rangle = {1 \over {{2^{(n + 1)/2}}}}\sum\limits_{k = 0}^{N - 1} | A\rangle |k{\rangle _\chi }|0{\rangle _\psi }|\varphi {\rangle _q}\left( {|0{\rangle _K} + |1{\rangle _K}} \right), thus creating the systems of orthonormal states |k〉_χ, k = 0, …, N – 1, and initializing the superposition state of the controlling qubit K.

Figure 2.

The circuit for the quantum part of the variational SVD algorithm. The depth of this circuit can be estimated as O(Qlog(N)/ε). (a) The circuit for creating the state |ψ_out〉_K given in (27); the operatorsW^(j), j = 0, …,6 are presented without subscripts for brevity. (b) The operators applied to the state |ψ_out〉_KB to probabilistically obtain the normalization G and the objective function L(α, β).

Now we double the state |k〉_χ creating the same state |k〉_ψ of the system ψ and also multiply the obtained state by the weight qk resulting in q_k|k〉_χ|k〉_ψ|0〉_q, see eq. (13). In the last case we use the trick proposed in [13] for matrix product. Both operations are controlled by the state of K, i.e., they are applied only if K is in the excited state |1〉_K. To arrange such control we introduce the projectors 11PχiK=|1〉χi|1〉Kχi1|K〈1|,i=1,…,n,〉,{P_{{\chi _i}K}} = |1{\rangle _{{\chi _i}}}|1{\rangle _{K{\chi _i}}}\left. {{{\left. 1 \right|}_K}\langle 1|,\quad i = 1, \ldots ,n,} \right\rangle , and the controlled operator 12WχψqK(1)=∏i=1n(PχiK⊗σψi(x)σqi(x)+(IχiK−PχiK)⊗Iψiqi).W_{\chi \psi qK}^{(1)} = \prod\limits_{i = 1}^n {\left( {{P_{{\chi _i}K}} \otimes \sigma _{{\psi _i}}^{(x)}\sigma _{{q_i}}^{(x)} + \left( {{I_{{\chi _i}K}} - {P_{{\chi _i}K}}} \right) \otimes {I_{{\psi _i}{q_i}}}} \right)} .

Hereafter the subscript attached to the notation of a subsystem indicates the appropriate qubit of this subsystem. Thus, subscript i mean the ith qubit of the appropriate subsystem in (12).

Applying WχψqK(1)W_{\chi \psi qK}^{(1)} we obtain 13|Φ2〉 = WχψqK(1)|Φ1〉 = 12(n+1)/2( ∑k=0N−1q0|A〉|k〉χ|k〉ψ|0〉q|0〉K + ∑k=0N−1qk|A〉|k〉χ|k〉ψ|0〉q|1〉K )+|g2〉,\left| {{\Phi _2}} \right\rangle \; = \;W_{\chi \psi qK}^{(1)}\left| {{\Phi _1}} \right\rangle \; = \;{1 \over {{2^{(n + 1)/2}}}}\left( {\sum\limits_{k = 0}^{N - 1} {{q_0}} |A\rangle |k{\rangle _\chi }|k{\rangle _\psi }|0{\rangle _q}|0{\rangle _K}} \right.\left. {\; + \;\sum\limits_{k = 0}^{N - 1} {{q_k}} |A\rangle |k{\rangle _\chi }|k{\rangle _\psi }|0{\rangle _q}|1{\rangle _K}} \right) + \left| {{g_2}} \right\rangle , where the garbage |g₂〉 collects the terms containing the states |j〉_q with j > 0, which we don’t need hereafter.

Next, we prepare and apply the unitary operators U(α) and U(β) in (2) controlled by the excited state of the one-qubit subsystem K: 14WχψK(2)= |1〉KK〈1|⊗Uχ(α)Uψ(β)+|0〉KK〈0|⊗Iχψ.W_{\chi \psi K}^{(2)} = \;|1{\rangle _{KK}}{\langle 1| \otimes {U_\chi }(\alpha ){U_\psi }(\beta ) + |0\rangle _{KK}}\langle 0| \otimes {I_{\chi \psi }}.

We can represent the action of the operator U on the vectors |k〉_χ and |k〉_ψ in terms of its elements as follows: 15Uχ(α)|k〉χ=∑l=0N−1blk(α)|l〉χ,Uψ(β)|k〉ψ=∑l=0N−1blk(β)|l〉ψ.{U_\chi }(\alpha )|k{\rangle _\chi } = \sum\limits_{l = 0}^{N - 1} {{b_{lk}}} (\alpha )|l{\rangle _\chi },\quad {U_\psi }(\beta )|k{\rangle _\psi } = \sum\limits_{l = 0}^{N - 1} {{b_{lk}}} (\beta )|l{\rangle _\psi }.

Then, applying WχψK(2)W_{\chi \psi K}^{(2)} to |Φ₂〉 we obtain 16|Φ3〉=WχψK(2)|Φ2〉=12(n+1)/2( ∑k=0N−1q0|A〉|k〉χ|k〉ψ|0〉q|0〉K + ∑k=0N−1∑l1=0N−1∑l2=0N−1qkbl1k(α)bl2k(β)|A〉|l1〉χ|l2〉ψ|0〉q|1〉K )+|g3〉.\matrix{ {\left| {{\Phi _3}} \right\rangle } \hfill & { = W_{\chi \psi K}^{(2)}\left| {{\Phi _2}} \right\rangle } \hfill \cr {} \hfill & { = {1 \over {{2^{(n + 1)/2}}}}\left( {\sum\limits_{k = 0}^{N - 1} {{q_0}} |A\rangle |k{\rangle _\chi }|k{\rangle _\psi }|0{\rangle _q}|0{\rangle _K}} \right.\left. {\; + \;\sum\limits_{k = 0}^{N - 1} {\sum\limits_{{l_1} = 0}^{N - 1} {\sum\limits_{{l_2} = 0}^{N - 1} {{q_k}} } } {b_{{l_1}k}}(\alpha ){b_{{l_2}k}}(\beta )|A\rangle {{\left| {{l_1}} \right\rangle }_\chi }{{\left| {{l_2}} \right\rangle }_\psi }|0{\rangle _q}|1{\rangle _K}} \right) + \left| {{g_3}} \right\rangle .} \hfill \cr }

Here, the garbage |g₂〉 from (13) is transformed to |g₃〉, but we don’t describe explicitly this transformation because we are not interested in the particular structure of the garbage. The same holds for the garbage in other states below. To multiply three matrices UχT(α)U_\chi ^T(\alpha ), A and U_ψ(β) and eventually calculate the sum Σ_j q_j 〈j|U^T(α)AU(β)|j〉_ψ in (2), we follow Refs. [12,13]. Using projectors (11) and projectors 17PψiK=|1〉ψi|1〉Kψi1|K〈1|,i=1,…,n,〉{P_{{\psi _i}K}} = |1{\rangle _{{\psi _i}}}|1{\rangle _K}{\psi _i}\left. {{{\left. 1 \right|}_K}\langle 1|,\quad i = 1, \ldots ,n,} \right\rangle we introduce the following controlled operators: 18CχKR(1)=∏i=1n(PχiK⊗σRi(x)+(IχiK−PχiK)⊗IRi),CψKC(2)=∏i=1n(PψiK⊗σCi(x)+(IψiK−PψiK)⊗ICi),\matrix{ {C_{\chi KR}^{(1)}} \hfill & { = \prod\limits_{i = 1}^n {\left( {{P_{{\chi _i}K}} \otimes \sigma _{{R_i}}^{(x)} + \left( {{I_{{\chi _i}K}} - {P_{{\chi _i}K}}} \right) \otimes {I_{{R_i}}}} \right)} ,} \hfill \cr {C_{\psi KC}^{(2)}} \hfill & { = \prod\limits_{i = 1}^n {\left( {{P_{{\psi _i}K}} \otimes \sigma _{{C_i}}^{(x)} + \left( {{I_{{\psi _i}K}} - {P_{{\psi _i}K}}} \right) \otimes {I_{{C_i}}}} \right)} ,} \hfill \cr }

Here, the operator CχKR(1)C_{\chi KR}^{(1)} is required for multiplying U^T(α) and A, the operator CψKC(2)C_{\psi KC}^{(2)} serves for multiplying A and U(β), Applying the operator WRCχψK(3)=CψKC(2)CχKR(1)W_{R{C_{\chi \psi K}}}^{(3)} = C_{\psi KC}^{(2)}C_{\chi KR}^{(1)} to the state |Φ₃〉 we obtain 19|Φ4〉=WRCχψK(3)|Φ3〉=12(n+1)/2( ∑k=0N−1q0a00|0〉R|0〉C|k〉χ|k〉ψ|0〉q|0〉K + ∑k=0N−1∑l1=0N−1∑l2=0N−1qkbl1k(α)bl2k(β)al1l2|0〉R|0〉C|〉|l1〉χ|l2〉ψ|0〉q|1〉K )+|g4〉\matrix{ {\left| {{\Phi _4}} \right\rangle = W_{RC\chi \psi K}^{(3)}\left| {{\Phi _3}} \right\rangle } \hfill \cr { = {1 \over {{2^{(n + 1)/2}}}}\left( {\sum\limits_{k = 0}^{N - 1} {{q_0}} {a_{00}}|0{\rangle _R}|0{\rangle _C}|k{\rangle _\chi }|k{\rangle _\psi }|0{\rangle _q}|0{\rangle _K}} \right.} \hfill \cr {\left. { + \;\sum\limits_{k = 0}^{N - 1} {\sum\limits_{{l_1} = 0}^{N - 1} {\sum\limits_{{l_2} = 0}^{N - 1} {{q_k}} } } {b_{{l_1}k}}(\alpha ){b_{{l_2}k}}(\beta ){a_{{l_1}{l_2}}}|0{\rangle _R}|0{\rangle _C}|\rangle {{\left| {{l_1}} \right\rangle }_\chi }{{\left| {{l_2}} \right\rangle }_\psi }|0{\rangle _q}|1{\rangle _K}} \right) + \left| {{g_4}} \right\rangle } \hfill \cr } where the first part in the rhs collects the terms needed for further calculations (these terms will be labelled later on by the operator WRCχψqBB˜(5)W_{RC\chi \psi qB\tilde B}^{(5)}, see eq. (24)) and |g₄〉 is the garbage to be removed later. Now, according to the multiplication algorithm (see Appendix in Ref. [13]) we introduce the operator 20Wχψ(4)=HχHψ,W_{\chi \psi }^{(4)} = {H_\chi }{H_\psi }, where H_χ and H_ψ are the sets of Hadamard transformations applied to each qubit of the subsystems χ and ψ respectively. These operators complete the multiplications U^T(α)A and AU(β) respectively, and simultaneously calculate the weighted trace Σ_j q_j(U^T(α)AU(β))_jj. Then, applying Wχχ˜ψ(4)W_{\chi \tilde \chi \psi }^{(4)} to the state |φ₄〉, selecting only the needed terms and moving others to the garbage |g₅〉, we obtain 21|Φ5〉=Wχχ˜ψ(4)|Φ4〉=12(3n+1)/2(a˜00|0〉K+∑l=0N−1Al|1〉K)|0〉R|0〉C|0〉χ|0〉ψ|0〉q+|g5〉,\matrix{ {\left| {{\Phi _5}} \right\rangle } \hfill & { = W_{\chi \tilde \chi \psi }^{(4)}\left| {{\Phi _4}} \right\rangle } \hfill \cr {} \hfill & { = {1 \over {{2^{(3n + 1)/2}}}}\left( {{{\tilde a}_{00}}|0{\rangle _K} + \sum\limits_{l = 0}^{N - 1} {{A_l}} |1{\rangle _K}} \right)|0{\rangle _R}|0{\rangle _C}|0{\rangle _\chi }|0{\rangle _\psi }|0{\rangle _q} + \left| {{g_5}} \right\rangle ,} \hfill \cr } where, 22a˜00=2nq0a00,{{\tilde a}_{00}} = {2^n}{q_0}{a_{00}}, 23Al(α,β)=∑k=0N−1∑m=0N−1qlakmbkl(α)bml(β)=ql(UT(α)MU(β))ll.{A_l}(\alpha ,\beta ) = \sum\limits_{k = 0}^{N - 1} {\sum\limits_{m = 0}^{N - 1} {{q_l}} } {a_{km}}{b_{kl}}(\alpha ){b_{ml}}(\beta ) = {q_l}{\left( {{U^T}(\alpha )MU(\beta )} \right)_{ll}}.

Remark that the factor 2ⁿ in the expression for ã₀₀ (22) appears because of the sum over k in (19) which includes n terms.

Now we label and remove the garbage |g₅〉 from the state (21) via the controlled measurement. To this end we introduce two one-qubit ancillae B and B˜{\tilde B} in the ground state and the controlled operator 24WRCχψqBB˜(5)=P⊗σB(x)σB˜(x)+(IRCχψqBB˜−P)⊗IBB˜P=|0〉R|0〉C|0〉χ|0〉ψ|0〉qR 〈0|C 〈0|χ〈0|ψ〈0|q 〈0|.\matrix{ {W_{RC\chi \psi qB\tilde B}^{(5)} = P \otimes \sigma _B^{(x)}\sigma _{\tilde B}^{(x)} + \left( {{I_{RC\chi \psi qB\tilde B}} - P} \right) \otimes {I_{B\tilde B}}} \hfill \cr {P = {{\left. {|0} \right\rangle }_R}{{\left. {|0} \right\rangle }_C}{{\left. {|0} \right\rangle }_\chi }{{\left. {|0} \right\rangle }_\psi }{{\left. {|0} \right\rangle }_q}R\;\left\langle {0{|_C}} \right.\;\left\langle {0{|_\chi }} \right.\left\langle {0{|_\psi }} \right.\left\langle {{{\left. 0 \right|}_q}} \right.\;\left\langle {0|} \right..\;} \hfill \cr }

Then, applying WRCχψqBB˜(5)W_{RC\chi \psi qB\tilde B}^{(5)} to the state |Φ5〉|0〉B|0〉B˜\left| {{\Phi _5}} \right\rangle |0{\rangle _B}|0{\rangle _{\tilde B}} we obtain 25|Φ6〉=WRCχψqBB˜(5)|Φ5〉|0〉B|0〉B˜=12(3n+1)/2(a˜00|0〉K+L˜(α,β)|1〉K)|0〉R|0〉C|0〉χ|0〉ψ|0〉q|1〉B|1〉B˜+|g5〉|0〉B|0〉B˜, L˜(α,β)=∑l=0N−1Al(α,β).\matrix{ {|{\Phi _6}\rangle = W_{RC\chi \psi qB\tilde B}^{(5)}|{\Phi _5}\rangle |0{\rangle _B}|0{\rangle _{\tilde B}} = } \hfill \cr {{1 \over {{2^{(3n + 1)/2}}}}{{({{\tilde a}_{00}}|0\rangle }_K} + \tilde L(\alpha ,\beta )|1{\rangle _K})|0{\rangle _R}|0{\rangle _C}|0{\rangle _\chi }|0{\rangle _\psi }|0{\rangle _q}|1{\rangle _B}|1{\rangle _{\tilde B}}} \hfill \cr { + |{g_5}\rangle |0{\rangle _B}|0{\rangle _{\tilde B}},\;\tilde L(\alpha ,\beta ) = \sum\limits_{l = 0}^{N - 1} {{A_l}} (\alpha ,\beta ).} \hfill \cr }

Now we can remove the garbage by introducing the controlled measurement [38], 26WBB˜(6)=|1〉BB〈1|⊗MB˜+|0〉BB〈0|⊗IB˜,W_{B\tilde B}^{(6)} = |1{\rangle _B}_B{\langle 1| \otimes {M_{\tilde B}} + |0\rangle _B}_B\langle 0| \otimes {I_{\tilde B}}, where MB˜{M_{\tilde B}} is the measurement operator applied to B˜{\tilde B}. Applying WBB˜(6)W_{B\tilde B}^{(6)} to |Φ₆〉 we obtain 27|Φ7〉=WBB˜(6)|Φ6〉=|Ψout〉KB|0〉R|0〉C|0〉χ|0〉ψ|0〉q,|Ψout〉KB=G−1(a˜00|0〉K+L˜(α,β)|1〉K)|1〉B,\matrix{ {\left| {{\Phi _7}} \right\rangle = W_{B\tilde B}^{(6)}\left| {{\Phi _6}} \right\rangle = {{\left| {{\Psi _{{\rm{out }}}}} \right\rangle }_{KB}}|0{\rangle _R}|0{\rangle _C}|0{\rangle _\chi }|0{\rangle _\psi }|0{\rangle _q},} \hfill \cr {{{\left| {{\Psi _{{\rm{out }}}}} \right\rangle }_{KB}} = {G^{ - 1}}\left( {{{\tilde a}_{00}}|0{\rangle _K} + \tilde L(\alpha ,\beta )|1{\rangle _K}} \right)|1{\rangle _B},} \hfill \cr } where G=a˜002+|L˜(α,β)|2G = \sqrt {\tilde a_{00}^2 + |\tilde L(\alpha ,\beta ){|^2}} and we recall that ã₀₀ in (22) is real because a₀₀ is the real element of the matrix A and q₀ is a positive integer.

We are aimed at finding the objective function L and normalization G. To this end we apply the Hadamard transformation H_B to the ancilla B and then apply the Hadamard transformation H_K, controlled by the excited state of B, to the qubit K, i.e., the controlled operator 28CBK=|1〉BB〈1|⊗HK+|0〉BB〈0|⊗IK.{C_{BK}} = |1{\rangle _B}_B{\langle 1| \otimes {H_K} + |0\rangle _B}_B\langle 0| \otimes {I_K}.

Thus, we have 29|Ψ˜out 〉=CBKHB|Ψout 〉KB=12G(a˜00|0〉K+L˜(α,β)|1〉K)|0〉B−12G((a˜00+L˜(α,β))|0〉K+(a˜00−L˜(α,β))|1〉K)|1〉B.\matrix{ {\left| {{{\tilde \Psi }_{{\rm{out }}}}} \right\rangle = {C_{BK}}{H_B}{{\left| {{\Psi _{{\rm{out }}}}} \right\rangle }_{KB}}} \hfill \cr { = {1 \over {\sqrt 2 G}}\left( {{{\tilde a}_{00}}|0{\rangle _K} + \tilde L(\alpha ,\beta )|1{\rangle _K}} \right)|0{\rangle _B}} \hfill \cr { - {1 \over {2G}}\left( {\left( {{{\tilde a}_{00}} + \tilde L(\alpha ,\beta )} \right)|0{\rangle _K} + \left( {{{\tilde a}_{00}} - \tilde L(\alpha ,\beta )} \right)|1{\rangle _K}} \right)|1{\rangle _B}.} \hfill \cr }

Now we measure the both qubits K and B with probabilities p_ij for fixing the state |i〉K|j〉B˜|i{\rangle _K}|j{\rangle _{\tilde B}}, i, j = 0,1, thus having 30p00(α,β)=a˜0022G2(α,β),p10(α,β)=|L˜(α,β)|22G2(α,β),p01(α,β)=G2(α,β)+2a˜00Re(L˜(α,β))4G2(α,β),p11(α,β)=G2(α,β)−2a˜00Re(L˜(α,β))4G2(α,β).\matrix{ {{p_{00}}(\alpha ,\beta ) = {{\tilde a_{00}^2} \over {2{G^2}(\alpha ,\beta )}},} \hfill \cr {{p_{10}}(\alpha ,\beta ) = {{|\tilde L(\alpha ,\beta ){|^2}} \over {2{G^2}(\alpha ,\beta )}},} \hfill \cr {{p_{01}}(\alpha ,\beta ) = {{{G^2}(\alpha ,\beta ) + 2{{\tilde a}_{00}}{\mathop{\rm Re}\nolimits} (\tilde L(\alpha ,\beta ))} \over {4{G^2}(\alpha ,\beta )}},} \hfill \cr {{p_{11}}(\alpha ,\beta ) = {{{G^2}(\alpha ,\beta ) - 2{{\tilde a}_{00}}{\mathop{\rm Re}\nolimits} (\tilde L(\alpha ,\beta ))} \over {4{G^2}(\alpha ,\beta )}}.} \hfill \cr }

From this system we obtain the expressions for the needed quantities: 31G2(α,β)=a˜0022p00(α,β)′,{G^2}(\alpha ,\beta ) = {{\tilde a_{00}^2} \over {2{p_{00}}{{(\alpha ,\beta )}^\prime }}}, 32L(α,β)=Re(L˜(α,β))=a˜002p00(α,β)(p01(α,β)−p11(α,β)).L(\alpha ,\beta ) = {\mathop{\rm Re}\nolimits} (\tilde L(\alpha ,\beta )) = {{{{\tilde a}_{00}}} \over {2{p_{00}}(\alpha ,\beta )}}\left( {{p_{01}}(\alpha ,\beta ) - {p_{11}}(\alpha ,\beta )} \right).

The depth of the operator WRCχψqBB˜(5)W_{RC\chi \psi qB\tilde B}^{(5)} can be estimated as O(log N) which is indicated in Figure 2. But the depth of the whole algorithm is defined by the operator WχψK(2)W_{\chi \psi K}^{(2)} whose depth is also indicated in Figure 2 and equals O(Qlog N). This estimation will be obtained in Section 2.3 after introducing the parameter Q. However, as for the depth of the whole algorithm for calculating the objective function, one can take into account the probabilistic method implemented for calculating the objective function using formula (32) that includes probabilities p_ij (i, j = 0, 1) given in (30). If the number of runs is N_r, then the depth of the algorithm is O(N_rQlogN). In turn, N_r ~ 1/ε, where ε ≪ 1 is the required precision for the probabilistic calculation of p_ij. Therefore, the depth is O(Olog(N)/ε). We have to note that, also expression (32) for the objective function L has p₀₀ in the denominator, there is no singularity at p₀₀ → 0, because (p₀₁ – p₁₁) → 0 as well in this case. However, to provide calculation with required precision ε, we require p₀₀ ≫ ε. This condition can be replaced by the following one: 33a˜002=22nq02a002≫ε ⇔ q0a00≫2−nε.\tilde a_{00}^2 = {2^{2n}}q_0^2a_{00}^2 \gg \varepsilon \; \Leftrightarrow \;{q_0}{a_{00}} \gg {2^{ - n}}\sqrt \varepsilon .

In this case p₀₀ ~ p₀₁ ~ p₁₁ ≫ ε, while the probability p₁₀ does not appear in (32). If the probabilities p₀₀, p₀₁ and p₁₁ are calculated with the precision ε, then formula (32) yields the objective function L with the precision ~ ε. Then, according to eq. (2), we conclude that the singular values can be calculated with the precision ~ ε/N.

Now we turn to the case when the condition (33) is destroyed and consider the case a˜002~ε\tilde a_{00}^2\~\varepsilon . In particular, ã can be zero. This is the case when we cannot use formulae (30)-(32) for calculating the objective function with precision e and we have to modify the algorithm. The simplest way is to change the starting point of the algorithm and replace formulae (9) and (10) with the following one: 34|Φ0〉=12(|0〉R|0〉C|0〉χ|0〉ψ|0〉q|1〉K+|A〉|0〉χ|0〉ψ|φ〉q|1〉K)=|Φ1〉.\left| {{\Phi _0}} \right\rangle = {1 \over {\sqrt 2 }}\left( {|0{\rangle _R}|0{\rangle _C}|0{\rangle _\chi }|0{\rangle _\psi }|0{\rangle _q}|1{\rangle _K} + |A\rangle |0{\rangle _\chi }|0{\rangle _\psi }|\varphi {\rangle _q}|1{\rangle _K}} \right) = \left| {{\Phi _1}} \right\rangle .

After proper modification of the algorithm we result in the formulae similar to (30) in which the parameter ã₀₀ is replaced with another parameter that does not depend on the elements of A and q. Such formulae allow to calculate the objective function L with any desired precision e. Further details of the modified algorithm will not be discussed here. We note that starting equation (34) can be used in case (33) as well. However, initial state (9) is simpler for preparation in comparison with state (34) and therefore it is recommended in case (33).

The space required for realization of this algorithm in both above cases is O(log N) qubits.

Realization of Operator WχψK(2)W_{\chi \psi K}^{(2)} for Real Matrix A

The operator WχψK(2)W_{\chi \psi K}^{(2)} in (14) can be conveniently represented as a product of two operators 35WχψK(2)=WχK(2)WψK(2),W_{\chi \psi K}^{(2)} = W_{\chi K}^{(2)}W_{\psi K}^{(2)}, 36WχK(2)=|1〉KK〈1|⊗Uχ(α)+|0〉KK〈0|⊗Iχ,W_{\chi K}^{(2)} = |1{\rangle _{KK}}{\langle 1| \otimes {U_\chi }(\alpha ) + |0\rangle _{KK}}\langle 0| \otimes {I_\chi }, 37WψK(2)=|1〉KK〈1|⊗Uψ(β)+|0〉KK〈0|⊗IψW_{\psi K}^{(2)} = |1{\rangle _K}_K{\langle 1| \otimes {U_\psi }(\beta ) + |0\rangle _K}_K\langle 0| \otimes {I_\psi } as shown in Figure 2a. Notice that two operators WχK(2)W_{\chi K}^{(2)} and WψK(2)W_{\psi K}^{(2)} are completely equivalent to each other and defer only by the parameters encoded into them. Therefore we describe only one of them, say WχK(2)W_{\chi K}^{(2)}. To realize the operator U for the real matrix A it is enough to use the one-qubit y-rotations R_y(φ) = exp(–σ^(y)φ/2) (σ^(y) is the Pauli matrix) and C-nots [26]: 38U(α)=∏k=1QRk(α(k−1)n+1,…,αkn),Rk(α(k−1)n+1,…,αkn)=∏m=1n−1Cχmχm+1∏j=1nRyχj(α(k−1)n+j),Cχmχm+1=|1〉χmχm〈1|⊗σχm+1(x)+|0〉χmχm〈0|⊗Iχm+1.\matrix{ {U(\alpha ) = \prod\limits_{k = 1}^Q {{R_k}} \left( {{\alpha _{(k - 1)n + 1}}, \ldots ,{\alpha _{kn}}} \right),} \hfill \cr {{R_k}\left( {{\alpha _{(k - 1)n + 1}}, \ldots ,{\alpha _{kn}}} \right) = \prod\limits_{m = 1}^{n - 1} {{C_{{\chi _m}{\chi _{m + 1}}}}} \prod\limits_{j = 1}^n {{R_{y{\chi _j}}}} \left( {{\alpha _{(k - 1)n + j}}} \right),} \hfill \cr {{C_{{\chi _m}{\chi _{m + 1}}}} = |1{\rangle _{{\chi _m}{\chi _m}}}{{\langle 1| \otimes \sigma _{{\chi _{m + 1}}}^{(x)} + |0\rangle }_{{\chi _m}{\chi _m}}}\langle 0| \otimes {I_{{\chi _{m + 1}}}}.} \hfill \cr }

In this formula, R_k represents a single block of transformations encoding n (the number of qubits in the subsystem χ) parameters α_i, i = 1,…, n, R_yχj is the y-rotation applied to the jth qubit of the subsystem χ. Involving Q blocks R_k, k = 1, …, Q, we enlarge the number of parameters to nQ. We note that this number, in general, may be bigger than the number of free real parameters in the N × N unitary transformation, which is N². Such increase in the number of parameters is caused by the non-standard parametrization of the unitary transformation U which, in turn, is chosen for two reasons: (i) simple realization of U in terms of one- and two-qubit operations and (ii) simple realization of derivatives of the objective function with respect to these parameters, see eq. (42). The depth of the operator U is O(log N).

Now we turn to realization of the controlled operator WχK(2)W_{\chi K}^{(2)} given in (36). To this end we substitute U determined in (38) into (36) and transform it to 39WχK(2)=∏k=1Q∏m=1n−1CKχmχm+1×∏j=1nRyχj(α(k−1)n+j/2)CKχj(z)Ryχj(α(k−1)n+j/2)CKχj(z)\matrix{ {W_{\chi K}^{(2)} = \prod\limits_{k = 1}^Q {\prod\limits_{m = 1}^{n - 1} {{C_{K{\chi _m}{\chi _{m + 1}}}}} } } \hfill \cr { \times \prod\limits_{j = 1}^n {{R_{y{\chi _j}}}} \left( {{\alpha _{(k - 1)n + j}}/2} \right)C_{K{\chi _j}}^{(z)}{R_{y{\chi _j}}}\left( {{\alpha _{(k - 1)n + j}}/2} \right)C_{K{\chi _j}}^{(z)}} \hfill \cr } where 40CKχmχm+1=PKχm⊗σχm+1(x)+(IKχm−PKχm)⊗Iχm+1,PKχm=|1〉K|1〉χmK1|χm〈1|,CKχj(z)=|0〉KK〈0|⊗σχj(z)+|1〉KK〈1|⊗Iχj.\matrix{ {{C_{K{\chi _m}{\chi _{m + 1}}}} = {P_{K{\chi _m}}} \otimes \sigma _{{\chi _{m + 1}}}^{(x)} + \left( {{I_{K{\chi _m}}} - {P_{K{\chi _m}}}} \right) \otimes {I_{{\chi _{m + 1}}}},} \hfill \cr {{P_{{K_{{\chi _m}}}}} = |1{\rangle _K}|1{\rangle _{{\chi _m}}}K{{\left. 1 \right|}_{{\chi _m}}}\langle 1|,} \hfill \cr {C_{K{\chi _j}}^{(z)} = |0{\rangle _{KK}}{{\langle 0| \otimes \sigma _{{\chi _j}}^{(z)} + |1\rangle }_{KK}}\langle 1| \otimes {I_{{\chi _j}}}.} \hfill \cr }

In (39), the controlled rotations R_yχj are represented by the second product since 41Ryχj(α(k−1)n+j/2)CKχj(z)Ryχj(α(k−1)n+j/2)CKχj(z)={ Iχj for |0〉KRyχj(α(k−1)n+j) for |1〉K,{R_{y{\chi _j}}}\left( {{\alpha _{(k - 1)n + j}}/2} \right)C_{K{\chi _j}}^{(z)}{R_{y{\chi _j}}}\left( {{\alpha _{(k - 1)n + j}}/2} \right)C_{K{\chi _j}}^{(z)} = \left\{ {\matrix{ {{I_{{\chi _j}}}} \hfill & {{\rm{ for }}|0{\rangle _K}} \hfill \cr {{R_{y{\chi _j}}}\left( {{\alpha _{(k - 1)n + j}}} \right)} \hfill & {{\rm{ for }}|1{\rangle _K}} \hfill \cr } ,} \right.

The circuit for the operator WχK(2)W_{\chi K}^{(2)} (and also for WψK(2)W_{\psi K}^{(2)}) is shown in Figure 3. The depth of each operator WχK(2)W_{\chi K}^{(2)} and WψK(2)W_{\psi K}^{(2)} is O(Qlog N) and therefore the depth of the whole circuit is O(Qlog N), where the parameter Q depends on the required precision of calculating the objective function. The parameter Q is used in the estimation of the depth of the algorithm in the paragraph below eq. (32).

Figure 3.

The circuit for the operatorsWψK(2)(α)W_{\chi K}^{(2)}(\alpha ) (and WψK(2)(β)W_{\psi K}^{(2)}(\beta )). Here the set of parameters γ is either α (for WψK(2)(α)W_{\chi K}^{(2)}(\alpha )) or β (for WψK(2)(β)W_{\psi K}^{(2)}(\beta )), Z ≡ σ^(z).

The input data for the classical optimization algorithm include not only the value of the objective function at the given values of the parameters α and β, but also derivatives of the objective function with respect to these parameters. It can be shown [26] that the required derivatives can be obtained calculating the objective function at certain values of the parameters α and β using the algorithm presented in Section 2.2.

Let γ^={ γ^1,…,γ^2nQ }\hat \gamma = \left\{ {{{\hat \gamma }_1}, \ldots ,{{\hat \gamma }_{2nQ}}} \right\} be the set of all parameters α and β: γ^={α,β}\hat \gamma = \{ \alpha ,\beta \} . Since all the parameters γ^{\hat \gamma } are introduced through the R_y-rotation, it is simple to calculate any-order derivative of L with respect to the parameters γ^j{{\hat \gamma }_j} [26]. For instance, for the first- and second-order derivatives we have 42∂L(γ^)∂γ^k=12L(γ^(k)),∂2L(γ^)∂γ^k∂γ^m=14L(γ^(k,m)),k,m=1,…,2nQ.\matrix{ {{{\partial L(\hat \gamma )} \over {\partial {{\hat \gamma }_k}}} = {1 \over 2}L\left( {{{\hat \gamma }^{(k)}}} \right),} \hfill \cr {{{{\partial ^2}L(\hat \gamma )} \over {\partial {{\hat \gamma }_k}\partial {{\hat \gamma }_m}}} = {1 \over 4}L\left( {{{\hat \gamma }^{(k,m)}}} \right),\quad k,m = 1, \ldots ,2nQ.} \hfill \cr }

Here γ^(k){{\hat \gamma }^{(k)}} means the set of parameters γ^{\hat \gamma }, in which γ^k{{\hat \gamma }_k} is replaced with γ^k+π{{\hat \gamma }_k} + \pi and γ^(k,m){{\hat \gamma }^{(k,m)}} is the set of parameters γ^(k){{\hat \gamma }^{(k)}}, in which γ^m(k)\hat \gamma _m^{(k)} is replaced with γ^m(k)+π\hat \gamma _m^{(k)} + \pi , i.e., 43γ^(k)=γ^|γ^k→γ^k+π,γ^(k,m)=γ^(k)|γ^m→γ^m+π,k,m=1,…,2nQ.\matrix{ {{{\hat \gamma }^{(k)}} = {{\left. {\hat \gamma } \right|}_{{{\hat \gamma }_k} \to {{\hat \gamma }_k} + \pi }},} \hfill \cr {{{\hat \gamma }^{(k,m)}} = {{\left. {{{\hat \gamma }^{(k)}}} \right|}_{{{\hat \gamma }_m} \to {{\hat \gamma }_m} + \pi }},k,m = 1, \ldots ,2nQ.} \hfill \cr }

In order to probabilistically find L(γ^)L(\hat \gamma ) at fixed values of the parameters, we have to perform one set of runs of quantum algorithm, the number of runs N_r in this set is determined by the required precision ε of L: N_r ~ 1/ε. Similarly, to probabilistically find all the first derivatives L(γ^(k))L\left( {{{\hat \gamma }^{(k)}}} \right), k = 1,…, 2nQ, at fixed values of the parameters, we have to perform 2nQ sets of runs (one set for each derivative). If the optimization algorithm requires also the second derivatives L(γ^(k,m))L\left( {{{\hat \gamma }^{(k,m)}}} \right), k, m = 1,…, 2nQ, we have to perform 2nQ(2nQ+1)2=2(nQ)2+nQ{{2nQ(2nQ + 1)} \over 2} = 2{(nQ)^2} + nQ sets of runs in addition to the above runs (we take into account that L(γ^(k,m))=L(γ^(m,k))L\left( {{{\hat \gamma }^{(k,m)}}} \right) = L\left( {{{\hat \gamma }^{(m,k)}}} \right)). The higher order derivatives of the objective function can be treated similarly. In this way we supply the objective function along with all necessary derivatives of this function to the input of the classical optimization algorithm which calculates the successive values of the parameters γ^{\hat \gamma }.

The variational algorithm for calculating the SVD is a hybrid one. It is described in Refs. [26,27] in details including examples of realization of the algorithm via Paddle Quantum [28] on the PaddlePaddle Deep Learning Platform [29,30]. The accuracy of the SVD obtained via the variational algorithm is estimated by comparing it with the original matrix. In [26], the authors give detailed analysis of variational quantum SVD (VQSVD) for the randomly generated 8 × 8 matrix M. The quality of VQSVD was characterized by the matrix distance ∥M – M_SVD∥₂, M_SVD is given in (1), ‖A‖2=∑i,j| aij |2A{_2} = \sqrt {\sum\nolimits_{i,j} {{{\left| {{a_{ij}}} \right|}^2}} } . It was shown that this distance reduces with an increase in the parameter Q. Several realizations of the unitary transformation U (see eq. (4)) where given and some applications of VQSVD (in particular, in Recommendation systems) were proposed. An alternative VQSVD was proposed in [27] with the principal novelty in encoding the elements of the matrix M into the quantum state of some system, at that the matrix elements must be ordered according to the certain scheme. The objective function was also different therein. The quality of VQSVD was characterized by the matrix distance using the Frobenius norm. In comparison with the VQSVD in [26], the modified algorithm demonstrates certain advantages.

The structure of the hybrid SVD algorithm is shown in Figure 4. We use the superscript ^[j] to label the jth iteration values of the parameters γ^{\hat \gamma }. The algorithm can be briefly described as follow. For some initial values of the parameters γ^[0]{{\hat \gamma }^{[0]}} we calculate the values of the objective function L(γ^[0])L\left( {{{\hat \gamma }^{[0]}}} \right) and its derivatives with respect to the parameters γ^{\hat \gamma } using the quantum algorithm. Then we use the found values of the objective function and its derivatives as the input for the classical optimization algorithm (for instance, for the gradient maximization algorithm) to find the succeeding iterated values of the parameters γ^[1]{{\hat \gamma }^{[1]}}. Next, we put them to the input of quantum algorithm, which calculates the objective function and its derivatives for the new values of the parameters γ^{\hat \gamma } and so on till we reach the required convergence criterion ε ≪ 1, i.e., till the following condition is satisfied: ΔL=| L(γ^[k+1])−L(γ^[k]) |<ε\Delta L = \left| {L\left( {{{\hat \gamma }^{[k + 1]}}} \right) - L\left( {{{\hat \gamma }^{[k]}}} \right)} \right| < \varepsilon . The scheme for this algorithm is shown in Figure 4.

Figure 4.

The hybrid algorithm for calculating SVD. The matrices Û, D and V^{\hat V} are defined in terms of U(α*) and U(β*) according to (4), ΔL = |L^[k+1] – L^[k]|. For simplicity, we indicate only the function L and first derivatives ∂_γL to be transferred from the output of the quantum algorithm to the input of the classical algorithm. However, the higher order derivatives might be required as well.