Matrix Encoding Method in Variational Quantum Singular Value Decomposition

Alexander I. Zenchuk; Wentao Qi; Junde Wu

doi:10.2478/qic-2025-0020

Introduction

1.

Quantum algorithms has become an attractive field of practical realization of quantum physics in everyday life. Among other quantum algorithms we concentrate on the family of algorithms manipulating with matrices and thus representing the quantum analogues of classical counterparts. Quantum Fourier transform [1–3] and phase estimation [3,4] must be distinguished as most popular and well recognized quantum algorithms used as subroutines included in many algorithms for data processing.

Set of algorithms reaches the goal using exclusively quantum approach. Apart from Fourier transform and Quantum phase estimation quoted above we refer to the algorithms for matrix manipulations (addition, multiplication, Kronecker sum, tensor product, Hadamard product) based on the Trotterization method [5–8] and the Baker-Champbell-Hausdorff [9] approximation for exponentiating matrices [10]. In [11], the matrix operations (addition, multiplication, inner product of vectors) are realized via action of special unitary operations on the matrix encoded into the either mixed or pure superposition state of some quantum system. Later, such unitary transformations where realized in terms of simple one- and two-qubit operations [12]. The matrix-encoding approach was implemented in the quantum algorithms for determinant calculation, matrix inversion and solving linear systems [13].

Another large family of algorithms includes algorithms combining both quantum and classical subroutines. To such algorithms one can refer the original version of the well-known Harrow-Hassidim-Lloyd (HHL) algorithm for solving systems of linear algebraic equations [14–21]. In this algorithm the classical subroutine is required for inversion of eigenvalues λ_j of the matrix A in equation A_x = b. However, later the method of approximate calculation of the inverse eigenvalues 1/λ_j was developed [22,23]. This method can be incorporated into HHL algorithm removing necessity of classical operations. Variational algorithms represent widely acknowledged class of hybrid algorithms for solving problems based on optimization methods. In particular, such algorithm was developed for the singular value decomposition (SVD) [24–27], where the loss (or objective) function at fixed optimization parameters was calculated by the quantum algorithm, while the iteration of parameters was performed using the classical gradient method. In Ref. [26], all simulations were implemented via Paddle Quantum [28] on the PaddlePaddle Deep Learning Platform [29,30].

We have to note that some principal issues on quantum algorithms for SVD are referred to Refs. [31–34]. Thus, the principle of realization of the quantum gradient descent algorithm is described in [31], but no explicit representation for unitary transformation performing iteration steps in this algorithm is given, this problem requires further study. In [32], the problem of fixing the phase of the singular vector associated with the appropriate singular value is explored. The data analysis involving the quantum SVD-based data representation is proposed in [33]. The quantum singular value estimation algorithm (i.e., estimation of the singular value associated with each singular vector) is described in [31,34]. However, the particular realizations of quantum algorithm are not discussed there. This fact motivates research on development of quantum SVD algorithms which can be validated on near-term quantum processors. The importance of the algorithms for SVD is determined by the wide applicability of SVD as a subroutine in various algorithms including some variants of matrix inversion [22,23], solving systems of linear equations [35], quantum recommendation systems [31,36,37].

In our paper, we modify the algorithm developed in [26,27] implementing the encoding the elements of the N × N matrix A into the probability amplitudes of the pure state of some quantum system with subsequent application of two parametrized unitary transformations and matrix multiplications using the algorithm proposed in [12]. Implementing this encoding we avoid representation of A as a linear combination of unitary transformations utilized in [26,27]. In comparison with the above references, we reduce the number of measurements required to get the complete information for calculating the objective function. In our case, each run of the algorithm is supplemented with O(1) measurements of the uncillae states (subsystems K and B below), while the number of measurements in [26,27] is O(N), here N is the number of diagonal elements in the considered matrix. Of course, because of the probabilistic method of obtaining the objective function, one has to perform series of runs of the algorithm.

The paper is organized as follows. In Section 2 we present the detailed description of our version of the quantum part of the variational SVD and briefly discuss the complete hybrid algorithm. Conclusions are given in Section 3.

Singular Value Decomposition

2.

Preliminaries

2.1.

We consider the singular value decomposition of an arbitrary square N × N matrix M assuming N = 2ⁿ,

1

M = \hat{U} D {\hat{V}}^{^{†}},

where D = diag(d₁, …, d_N) is the diagonal matrix of singular values (some of those values might be zero) and Û and

{\hat{V}}^{†}

are unitary matrices. To find the singular values (entries of the diagonal matrix D) we, first of all, introduce the following objective function [26]:

2

L (α, β) = \sum_{j = 0}^{N - 1} q_{j} \times Re (〈 ψ_{j} | U^{T} (α) M U (β) | ψ_{j} 〉), U = {b_{i j} : i, j = 0, \dots, N - 1},

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

3

\max_{α, β} L (α, β) = L (α^{*}, β^{*}) .

Then

4

\begin{array}{l} d_{j} = 〈 ψ_{j} | U^{T} (α^{*}) M U (β^{*}) | ψ_{j} 〉, j = 0, \dots, N - 1, \\ {\hat{U}}^{†} = U^{T} (α^{*}), \hat{V} = U (β^{*}) . \end{array}

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} = \sum_{j = 0}^{N - 1} q_{j} | j 〉_{q}, \sum_{j = 0}^{N - 1} q_{j}^{2} = 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 $\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 $\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(β):

6

{\hat{U}}^{†} = U^{T} (α) = {\hat{V}}^{†} = U^{†} (β), \Rightarrow U (α) = \bar{U} (β),

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.

Quantum Algorithm Preparing Objective Function

2.2.

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}:

7

A = \frac{e^{- i \arg (m_{00})}}{\sqrt{\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {| m_{i j} |}^{2}}} M .

Now we can encode the elements of the matrix A into the superposition state of R and C as follows:

8

\begin{array}{l} | A 〉 = \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} a_{i j} | i 〉_{R} | j 〉_{C}, \\ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {| a_{i j} |}^{2} = 1, \end{array}

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

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_{χ} H_{K}$ :

10

| Φ_{1} 〉 = W_{χ K}^{(0)} | Φ_{0} 〉 = \frac{1}{2^{(n + 1) / 2}} \sum_{k = 0}^{N - 1} | A 〉 | k 〉_{χ} | 0 〉_{ψ} | φ 〉_{q} (| 0 〉_{K} + | 1 〉_{K}),

thus creating the systems of orthonormal states |k〉_χ, k = 0, …, N – 1, and initializing the superposition state of the controlling qubit K.

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

11

P_{χ_{i} K} = | 1 〉_{χ_{i}} | 1 〉_{K χ_{i}} {1 |}_{K} 〈 1 |, i = 1, \dots, n, 〉,

and the controlled operator

12

W_{χ ψ q K}^{(1)} = \prod_{i = 1}^{n} (P_{χ_{i} K} \otimes σ_{ψ_{i}}^{(x)} σ_{q_{i}}^{(x)} + (I_{χ_{i} K} - P_{χ_{i} K}) \otimes I_{ψ_{i} q_{i}}) .

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_{χ ψ q K}^{(1)}$ we obtain

13

| Φ_{2} 〉 = W_{χ ψ q K}^{(1)} | Φ_{1} 〉 = \frac{1}{2^{(n + 1) / 2}} (\sum_{k = 0}^{N - 1} q_{0} | A 〉 | k 〉_{χ} | k 〉_{ψ} | 0 〉_{q} | 0 〉_{K} + \sum_{k = 0}^{N - 1} q_{k} | A 〉 | k 〉_{χ} | k 〉_{ψ} | 0 〉_{q} | 1 〉_{K}) + | g_{2} 〉,

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:

14

W_{χ ψ K}^{(2)} = | 1 〉_{K K} 〈_{1 | \otimes U_{χ} (α) U_{ψ} (β) + | 0 〉 K K} 〈 0 | \otimes I_{χ ψ} .

We can represent the action of the operator U on the vectors |k〉_χ and |k〉_ψ in terms of its elements as follows:

15

U_{χ} (α) | k 〉_{χ} = \sum_{l = 0}^{N - 1} b_{l k} (α) | l 〉_{χ}, U_{ψ} (β) | k 〉_{ψ} = \sum_{l = 0}^{N - 1} b_{l k} (β) | l 〉_{ψ} .

Then, applying $W_{χ ψ K}^{(2)}$ to |Φ₂〉 we obtain

16

\begin{array}{l} | Φ_{3} 〉 & = W_{χ ψ K}^{(2)} | Φ_{2} 〉 \\ = \frac{1}{2^{(n + 1) / 2}} (\sum_{k = 0}^{N - 1} q_{0} | A 〉 | k 〉_{χ} | k 〉_{ψ} | 0 〉_{q} | 0 〉_{K} + \sum_{k = 0}^{N - 1} \sum_{l_{1} = 0}^{N - 1} \sum_{l_{2} = 0}^{N - 1} q_{k} b_{l_{1} k} (α) b_{l_{2} k} (β) | A 〉 {| l_{1} 〉}_{χ} {| l_{2} 〉}_{ψ} | 0 〉_{q} | 1 〉_{K}) + | g_{3} 〉 . \end{array}

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} (α)$ , 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

17

P_{ψ_{i} K} = | 1 〉_{ψ_{i}} | 1 〉_{K} ψ_{i} {1 |}_{K} 〈 1 |, i = 1, \dots, n, 〉

we introduce the following controlled operators:

18

\begin{array}{l} C_{χ K R}^{(1)} & = \prod_{i = 1}^{n} (P_{χ_{i} K} \otimes σ_{R_{i}}^{(x)} + (I_{χ_{i} K} - P_{χ_{i} K}) \otimes I_{R_{i}}), \\ C_{ψ K C}^{(2)} & = \prod_{i = 1}^{n} (P_{ψ_{i} K} \otimes σ_{C_{i}}^{(x)} + (I_{ψ_{i} K} - P_{ψ_{i} K}) \otimes I_{C_{i}}), \end{array}

Here, the operator $C_{χ K R}^{(1)}$ is required for multiplying U^T(α) and A, the operator $C_{ψ K C}^{(2)}$ serves for multiplying A and U(β), Applying the operator $W_{R C_{χ ψ K}}^{(3)} = C_{ψ K C}^{(2)} C_{χ K R}^{(1)}$ to the state |Φ₃〉 we obtain

19

\begin{array}{l} | Φ_{4} 〉 = W_{R C χ ψ K}^{(3)} | Φ_{3} 〉 \\ = \frac{1}{2^{(n + 1) / 2}} (\sum_{k = 0}^{N - 1} q_{0} a_{00} | 0 〉_{R} | 0 〉_{C} | k 〉_{χ} | k 〉_{ψ} | 0 〉_{q} | 0 〉_{K} \\ + \sum_{k = 0}^{N - 1} \sum_{l_{1} = 0}^{N - 1} \sum_{l_{2} = 0}^{N - 1} q_{k} b_{l_{1} k} (α) b_{l_{2} k} (β) a_{l_{1} l_{2}} | 0 〉_{R} | 0 〉_{C} | 〉 {| l_{1} 〉}_{χ} {| l_{2} 〉}_{ψ} | 0 〉_{q} | 1 〉_{K}) + | g_{4} 〉 \end{array}

where the first part in the rhs collects the terms needed for further calculations (these terms will be labelled later on by the operator

W_{R C χ ψ q B \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

20

W_{χ ψ}^{(4)} = H_{χ} H_{ψ},

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_{χ \tilde{χ} ψ}^{(4)}

to the state |φ₄〉, selecting only the needed terms and moving others to the garbage |g₅〉, we obtain

21

\begin{array}{l} | Φ_{5} 〉 & = W_{χ \tilde{χ} ψ}^{(4)} | Φ_{4} 〉 \\ = \frac{1}{2^{(3 n + 1) / 2}} ({\tilde{a}}_{00} | 0 〉_{K} + \sum_{l = 0}^{N - 1} A_{l} | 1 〉_{K}) | 0 〉_{R} | 0 〉_{C} | 0 〉_{χ} | 0 〉_{ψ} | 0 〉_{q} + | g_{5} 〉, \end{array}

where,

22

{\tilde{a}}_{00} = 2^{n} q_{0} a_{00},

23

A_{l} (α, β) = \sum_{k = 0}^{N - 1} \sum_{m = 0}^{N - 1} q_{l} a_{k m} b_{k l} (α) b_{m l} (β) = q_{l} {(U^{T} (α) M U (β))}_{l l} .

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 $\tilde{B}$ in the ground state and the controlled operator

24

\begin{array}{l} W_{R C χ ψ q B \tilde{B}}^{(5)} = P \otimes σ_{B}^{(x)} σ_{\tilde{B}}^{(x)} + (I_{R C χ ψ q B \tilde{B}} - P) \otimes I_{B \tilde{B}} \\ P = {| 0 〉}_{R} {| 0 〉}_{C} {| 0 〉}_{χ} {| 0 〉}_{ψ} {| 0 〉}_{q} R 〈 0 |_{C} 〈 0 |_{χ} 〈 0 |_{ψ} 〈 {0 |}_{q} 〈 0 | . \end{array}

Then, applying $W_{R C χ ψ q B \tilde{B}}^{(5)}$ to the state $| Φ_{5} 〉 | 0 〉_{B} | 0 〉_{\tilde{B}}$ we obtain

25

\begin{array}{l} | Φ_{6} 〉 = W_{R C χ ψ q B \tilde{B}}^{(5)} | Φ_{5} 〉 | 0 〉_{B} | 0 〉_{\tilde{B}} = \\ \frac{1}{2^{(3 n + 1) / 2}} {({\tilde{a}}_{00} | 0 〉}_{K} + \tilde{L} (α, β) | 1 〉_{K}) | 0 〉_{R} | 0 〉_{C} | 0 〉_{χ} | 0 〉_{ψ} | 0 〉_{q} | 1 〉_{B} | 1 〉_{\tilde{B}} \\ + | g_{5} 〉 | 0 〉_{B} | 0 〉_{\tilde{B}}, \tilde{L} (α, β) = \sum_{l = 0}^{N - 1} A_{l} (α, β) . \end{array}

Now we can remove the garbage by introducing the controlled measurement [38],

26

W_{B \tilde{B}}^{(6)} = | 1 〉_{B}_{B} {〈 1 | \otimes M_{\tilde{B}} + | 0 〉}_{B}_{B} 〈 0 | \otimes I_{\tilde{B}},

where

M_{\tilde{B}}

is the measurement operator applied to

\tilde{B}

. Applying

W_{B \tilde{B}}^{(6)}

to |Φ₆〉 we obtain

27

\begin{array}{l} | Φ_{7} 〉 = W_{B \tilde{B}}^{(6)} | Φ_{6} 〉 = {| Ψ_{o u t} 〉}_{K B} | 0 〉_{R} | 0 〉_{C} | 0 〉_{χ} | 0 〉_{ψ} | 0 〉_{q}, \\ {| Ψ_{o u t} 〉}_{K B} = G^{- 1} ({\tilde{a}}_{00} | 0 〉_{K} + \tilde{L} (α, β) | 1 〉_{K}) | 1 〉_{B}, \end{array}

where

G = \sqrt{{\tilde{a}}_{00}^{2} + | \tilde{L} (α, β) |^{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.

Remark 1.

To obtain the output state |ψ_out〉_KB we involve so-called controlled measurement. Of course, we could obtain this state just using multiple running of the algorithm each time measuring the state of the ancilla $\tilde{B}$ till obtaining the state $| 1 〉_{\tilde{B}}$ as the required result of measurement. In this case we again obtain the state |Φ₇〉 (27) with the output state |ψ_out〉_KB. However, the probability of such output of the measurement is O(2^–(3n+1)) which exponentially decreases with the number of qubits n characterizing the space of the circuit. This is a disadvantage of the algorithm that can be completely overcome referring to the controlled measurement. Although we can not present a particular realization of this operator in terms of well-known quantum and/or classical operations, the controlled measurement reflects the deep relation between the quantum and classical physics. The measurement of state of the qubit $\tilde{B}$ will be performed if only the excited state |1〉_B exists in the considered superposition state. Since our superposition state |Φ₆〉 (25) includes |1〉_B by construction, the measurement must be performed with the predictable result $| 1 〉_{\tilde{B}}$ . In addition, the controlled measurement is the last element in the following scheme. We know the controlled operations, where both controlling and controlled subsystems are quantum (CNOT is the simplest representative). The controlled operations with controlling classical subsystem and controlled quantum subsystem are also acknowledged [3]. The proposed concept of controlled measurement represent the last possible type of control operations in quantum-classical system with the control by the state of a quantum subsystem. All these arguments do not prove realizability of controlled measurement but motivate the research for its realization that would not contradict the basic postulates and theorems of quantum mechanics including the no-cloning theorem [39].

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

28

C_{B K} = | 1 〉_{B}_{B} 〈_{1 | \otimes H_{K} + | 0 〉 B}_{B} 〈 0 | \otimes I_{K} .

Thus, we have

29

\begin{array}{l} | {\tilde{Ψ}}_{out} 〉 = C_{B K} H_{B} {| Ψ_{out} 〉}_{K B} \\ = \frac{1}{\sqrt{2} G} ({\tilde{a}}_{00} | 0 〉_{K} + \tilde{L} (α, β) | 1 〉_{K}) | 0 〉_{B} \\ - \frac{1}{2 G} (({\tilde{a}}_{00} + \tilde{L} (α, β)) | 0 〉_{K} + ({\tilde{a}}_{00} - \tilde{L} (α, β)) | 1 〉_{K}) | 1 〉_{B} . \end{array}

Now we measure the both qubits K and B with probabilities p_ij for fixing the state $| i 〉_{K} | j 〉_{\tilde{B}}$ , i, j = 0,1, thus having

30

\begin{array}{l} p_{00} (α, β) = \frac{{\tilde{a}}_{00}^{2}}{2 G^{2} (α, β)}, \\ p_{10} (α, β) = \frac{| \tilde{L} (α, β) |^{2}}{2 G^{2} (α, β)}, \\ p_{01} (α, β) = \frac{G^{2} (α, β) + 2 {\tilde{a}}_{00} Re (\tilde{L} (α, β))}{4 G^{2} (α, β)}, \\ p_{11} (α, β) = \frac{G^{2} (α, β) - 2 {\tilde{a}}_{00} Re (\tilde{L} (α, β))}{4 G^{2} (α, β)} . \end{array}

From this system we obtain the expressions for the needed quantities:

31

G^{2} (α, β) = \frac{{\tilde{a}}_{00}^{2}}{2 p_{00} {(α, β)}^{'}},

32

L (α, β) = Re (\tilde{L} (α, β)) = \frac{{\tilde{a}}_{00}}{2 p_{00} (α, β)} (p_{01} (α, β) - p_{11} (α, β)) .

The depth of the operator $W_{R C χ ψ q B \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)}$ 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:

33

{\tilde{a}}_{00}^{2} = 2^{2 n} q_{0}^{2} a_{00}^{2} ≫ ε \Leftrightarrow q_{0} a_{00} ≫ 2^{- n} \sqrt{ε} .

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 ${\tilde{a}}_{00}^{2} ~ ε$ . 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} 〉 = \frac{1}{\sqrt{2}} (| 0 〉_{R} | 0 〉_{C} | 0 〉_{χ} | 0 〉_{ψ} | 0 〉_{q} | 1 〉_{K} + | A 〉 | 0 〉_{χ} | 0 〉_{ψ} | φ 〉_{q} | 1 〉_{K}) = | Φ_{1} 〉 .

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)}$ for Real Matrix A

2.3.

The operator $W_{χ ψ K}^{(2)}$ in (14) can be conveniently represented as a product of two operators

35

W_{χ ψ K}^{(2)} = W_{χ K}^{(2)} W_{ψ K}^{(2)},

36

W_{χ K}^{(2)} = | 1 〉_{K K} 〈 1 | \otimes U_{χ} (α) + | 0 〉_{K K} 〈 0 | \otimes I_{χ},

37

W_{ψ K}^{(2)} = | 1 〉_{K}_{K} 〈 1 | \otimes U_{ψ} (β) + | 0 〉_{K}_{K} 〈 0 | \otimes I_{ψ}

as shown in Figure 2a. Notice that two operators

W_{χ K}^{(2)}

and

W_{ψ 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)}

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

38

\begin{array}{l} U (α) = \prod_{k = 1}^{Q} R_{k} (α_{(k - 1) n + 1}, \dots, α_{k n}), \\ R_{k} (α_{(k - 1) n + 1}, \dots, α_{k n}) = \prod_{m = 1}^{n - 1} C_{χ_{m} χ_{m + 1}} \prod_{j = 1}^{n} R_{y χ_{j}} (α_{(k - 1) n + j}), \\ C_{χ_{m} χ_{m + 1}} = | 1 〉_{χ_{m} χ_{m}} {〈 1 | \otimes σ_{χ_{m + 1}}^{(x)} + | 0 〉}_{χ_{m} χ_{m}} 〈 0 | \otimes I_{χ_{m + 1}} . \end{array}

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)}$ given in (36). To this end we substitute U determined in (38) into (36) and transform it to

39

\begin{array}{l} W_{χ K}^{(2)} = \prod_{k = 1}^{Q} \prod_{m = 1}^{n - 1} C_{K χ_{m} χ_{m + 1}} \\ \times \prod_{j = 1}^{n} R_{y χ_{j}} (α_{(k - 1) n + j} / 2) C_{K χ_{j}}^{(z)} R_{y χ_{j}} (α_{(k - 1) n + j} / 2) C_{K χ_{j}}^{(z)} \end{array}

where

40

\begin{array}{l} C_{K χ_{m} χ_{m + 1}} = P_{K χ_{m}} \otimes σ_{χ_{m + 1}}^{(x)} + (I_{K χ_{m}} - P_{K χ_{m}}) \otimes I_{χ_{m + 1}}, \\ P_{K_{χ_{m}}} = | 1 〉_{K} | 1 〉_{χ_{m}} K {1 |}_{χ_{m}} 〈 1 |, \\ C_{K χ_{j}}^{(z)} = | 0 〉_{K K} {〈 0 | \otimes σ_{χ_{j}}^{(z)} + | 1 〉}_{K K} 〈 1 | \otimes I_{χ_{j}} . \end{array}

In (39), the controlled rotations R_yχj are represented by the second product since

41

R_{y χ_{j}} (α_{(k - 1) n + j} / 2) C_{K χ_{j}}^{(z)} R_{y χ_{j}} (α_{(k - 1) n + j} / 2) C_{K χ_{j}}^{(z)} = {\begin{array}{l} I_{χ_{j}} & for | 0 〉_{K} \\ R_{y χ_{j}} (α_{(k - 1) n + j}) & for | 1 〉_{K} \end{array},

The circuit for the operator $W_{χ K}^{(2)}$ (and also for $W_{ψ K}^{(2)}$ ) is shown in Figure 3. The depth of each operator $W_{χ K}^{(2)}$ and $W_{ψ 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).

The circuit for the operators $W_{ψ K}^{(2)} (α)$ (and $W_{ψ K}^{(2)} (β)$ ). Here the set of parameters γ is either α (for $W_{ψ K}^{(2)} (α)$ ) or β (for $W_{ψ K}^{(2)} (β)$ ), Z ≡ σ^(z).

Remark 2.

For the case of complex matrix A, the R_y-rotation should be replaced with the operator R(η, φ, θ) = R_z(η)Ry(φ)R_z(θ).

Derivatives of Objective Function

2.4.

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 $\hat{γ} = {{\hat{γ}}_{1}, \dots, {\hat{γ}}_{2 n Q}}$ be the set of all parameters α and β: $\hat{γ} = {α, β}$ . Since all the parameters $\hat{γ}$ are introduced through the R_y-rotation, it is simple to calculate any-order derivative of L with respect to the parameters ${\hat{γ}}_{j}$ [26]. For instance, for the first- and second-order derivatives we have

42

\begin{array}{l} \frac{\partial L (\hat{γ})}{\partial {\hat{γ}}_{k}} = \frac{1}{2} L ({\hat{γ}}^{(k)}), \\ \frac{\partial^{2} L (\hat{γ})}{\partial {\hat{γ}}_{k} \partial {\hat{γ}}_{m}} = \frac{1}{4} L ({\hat{γ}}^{(k, m)}), k, m = 1, \dots, 2 n Q . \end{array}

Here ${\hat{γ}}^{(k)}$ means the set of parameters $\hat{γ}$ , in which ${\hat{γ}}_{k}$ is replaced with ${\hat{γ}}_{k} + π$ and ${\hat{γ}}^{(k, m)}$ is the set of parameters ${\hat{γ}}^{(k)}$ , in which ${\hat{γ}}_{m}^{(k)}$ is replaced with ${\hat{γ}}_{m}^{(k)} + π$ , i.e.,

43

\begin{array}{l} {\hat{γ}}^{(k)} = {\hat{γ} |}_{{\hat{γ}}_{k} \to {\hat{γ}}_{k} + π}, \\ {\hat{γ}}^{(k, m)} = {{\hat{γ}}^{(k)} |}_{{\hat{γ}}_{m} \to {\hat{γ}}_{m} + π}, k, m = 1, \dots, 2 n Q . \end{array}

In order to probabilistically find $L (\hat{γ})$ 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 ({\hat{γ}}^{(k)})$ , 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 ({\hat{γ}}^{(k, m)})$ , k, m = 1,…, 2nQ, we have to perform $\frac{2 n Q (2 n Q + 1)}{2} = 2 (^{n Q) 2} + n Q$ sets of runs in addition to the above runs (we take into account that $L ({\hat{γ}}^{(k, m)}) = L ({\hat{γ}}^{(m, k)})$ ). 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{γ}$ .

Hybrid Algorithm for SVD: Brief Discussion

2.5.

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} = \sqrt{\sum_{i, j} {| a_{i j} |}^{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{γ}$ . The algorithm can be briefly described as follow. For some initial values of the parameters ${\hat{γ}}^{[0]}$ we calculate the values of the objective function $L ({\hat{γ}}^{[0]})$ and its derivatives with respect to the parameters $\hat{γ}$ 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 ${\hat{γ}}^{[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{γ}$ and so on till we reach the required convergence criterion ε ≪ 1, i.e., till the following condition is satisfied: $Δ L = | L ({\hat{γ}}^{[k + 1]}) - L ({\hat{γ}}^{[k]}) | < ε$ . The scheme for this algorithm is shown in Figure 4.

The hybrid algorithm for calculating SVD. The matrices *Û, D* and $\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.

Conclusions

3.

We present a new version of the quantum part of the variational SVD algorithm based on the matrix encoding approach, when the entries of the matrix M are encoded into the probability amplitudes of the superposition state of a quantum system, in our case, subsystems R, C. The one-qubit subsystem K is an auxiliary subsystem that controls set of unitary operations. Eventually, the resulting state |Ψ_out〉_KB is the superposition state of this auxiliary system K multiplied by the excited state |1〉_B of the ancilla B. This state is obtained as a result of controlled measurement of the ancilla $\tilde{B}$ which removes the problem of small success probability that unavoidably appears in the case of ordinary measurement of the ancilla state because of the Hadamard transformation used in this algorithm. At the moment, we cannot suggest a particular realization of the controlled measurement. However, this operator means the control of the classical operation (measurement of the qubit $\tilde{B}$ ) by the quantum state of another qubit which is the qubit B in our case, and there is no any particular reason to discard possibility of such control. This control would reflect the deep relation between the classical and quantum phenomena without strong border between them. The controlled measurement means that the measurement of $\tilde{B}$ will be performed if only the superposition state includes the excited qubit |1〉_B. Otherwise, the quantum system remains unchanged. Thus, the realization of controlled measurement remains an open problem. It is quite possible that this concept requires some modification to be realizable.

To measure the value of the objective function we use the state |Ψ_out〉_KB and, after the Hadamard and controlled Hadamard transformations, we find the probabilities of the states |i〉_K|j〉_B and then required objective function L(α, β). Since the result is probabilistic we have to run the algorithm many times to obtain the required precision for the objective function. However, this multiple running is a necessary part of any probabilistic algorithm. In a similar way we can calculate all derivatives of the objective function required for running the classical optimization algorithm. The depth of the whole quantum algorithm can be estimated as O(Qlog(N)/ε), it linearly depends on the number of runs N_r which, in turn, is inverse proportional to the precision ε required for the probabilistic measurement of the objective function. The space of the algorithm is O(log N). We also notice that the different type of the matrix encoding is used in [27] yielding certain privileges for that algorithm over the algorithm in [26].

Although our algorithm deals with square matrices, it can be applied to the rectangular matrices as well because the rectangular matrix can be written in a square form by adding appropriate number of zero rows or columns. We also have to note that SVD is also a key for constructing the inverse or pseudoinverse of the matrix [34] because the pseudoinverse matrix for any given matrix A having SVD, $A = \sum_{i = 1}^{r} s_{i} | y_{i} 〉〈 x_{i} |$ , can be written as $A^{+} = \sum_{i = 1}^{r} \frac{1}{s_{i}} | x_{i} 〉〈 y_{i} |$ .

The fact that matrix-encoding approach is applicable to the variational Quantum SVD algorithm confirms the wide applicability of this approach which has already been used in algorithms for matrix manipulations including addition, multiplication, determinant calculation, inverse matrix calculation and solving systems of linear equations [11–13,38].

Acknowledgments

The manuscript is supported by the National Natural Science Foundation of China (Grants No. 12031004, No. 12271474 and No. 61877054). The work was partially funded by a state task of Russian Fundamental Investigations (State Registration No. 124013000760-0).

Notes

[1] Contributed by Author Contributions

All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

[2] Conflicts of interest Conflicts of Interest

The authors claim that they do not have any Conflict of Interest.

[3] Data Availability Statement

All necessary data are included into the manuscript.

Matrix Encoding Method in Variational Quantum Singular Value Decomposition

Full Article

Introduction

1.

Singular Value Decomposition

2.

Preliminaries

2.1.

1

2

3

4

5

Figure 1.

6

Quantum Algorithm Preparing Objective Function

2.2.

7

8

9

10

Figure 2.

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

Remark 1.

28

29

30

31

32

33

34

Realization of Operator WχψK(2) for Real Matrix A

2.3.

35

36

37

38

39

40

41

Figure 3.

Remark 2.

Derivatives of Objective Function

2.4.

42

43

Hybrid Algorithm for SVD: Brief Discussion

2.5.

Figure 4.

Conclusions

3.

Acknowledgments

Notes

Realization of Operator $W_{χ ψ K}^{(2)}$ for Real Matrix A