## Quantum Higher Order Singular Value Decomposition

Lejia Gu, ∗ Xiaoqiang Wang, † H. W. Joseph Lee, ‡ and Guofeng Zhang § Department of Applied Mathematics, The Hong Kong Polytechnic University, Hong Kong

(Dated: April 7, 2020)

Higher order singular value decomposition (HOSVD) is an important tool for analyzing big data in multilinear algebra and machine learning. In this paper, we present two quantum algorithms for HOSVD. Our methods allow one to decompose a tensor into a core tensor containing tensor singular values and some unitary matrices by quantum computers. Compared to the classical HOSVD algorithm, our quantum algorithms provide an exponential speedup. Furthermore, we introduce a hybrid quantum-classical algorithm of HOSVD model applied in recommendation systems.

## I. INTRODUCTION

Matrix computations are vital to many optimization and machine learning problems. Nowadays, due to the rise of neural networks in machine learning methods, the elements of a network are usually described by tensors which can have more than two indices. Tensors (or hypermatrices), as a higher order generalization of matrices, have found widespread applications in scientific and engineering fields [23, 27, 30, 31]. Tensor decomposition describes a tensor as a sequence of elementary operations acting on other, often simpler tensors. Usually, key information can be extracted from the decomposed tensor, and less space is required to store the original tensor to some accuracy. Tensor network, as a countable collection of tensors connected by contractions, has been widely employed in training machine learning models. A quantum state has a tensor representation. Hence, a quantum network, namely a multipartite system, can be represented by a tensor network. Indeed, quantum circuits are a special class of tensor networks, where the arrangement of the tensors and their types are restricted [5, 11, 18]. Moreover, tensor analysis has been applied to the quantum entanglement problem and classicality problem of spin states [17, 32, 33, 47].

Quantum computers are devices that perform calculations by utilizing quantum mechanical features including superposition and entanglement. Although largescale quantum computers are not built yet, theoretical research on quantum algorithms has been conducted for a few decades. In 1994, Shor's algorithm [38], a hybrid quantum-classical algorithm, is proved to be able to solve integer-factorization problem with polynomial time, while it is widely regarded as an NP hard problem in classical computing. In 1996, Grover's search algorithm [13] finds an entry from an unstructured database quadratically faster than classical algorithms. In 2009,

∗ le-jia.gu@connect.polyu.hk

† xiaoqiang.wang@connect.polyu.hk

‡ joseph.lee@polyu.edu.hk

§ guofeng.zhang@polyu.edu.hk

Harrow, Hassidim and Lloyd put forward a quantum algorithm for solving linear systems of equations, which is famous as the HHL Algorithm [15]. Physical demonstrations of the HHL algorithm can be found in [8, 28, 48]. Base on this algorithm, many quantum versions of classical machine learning methods are designed, such as quantum least-squares linear regression [45] and support vector machines [34]. Also, the quantum singular value decomposition of nonsparse low-rank matrices is introduced in [35]. A different method to implement singular value decomposition, named as quantum singular value estimation (QSVE), is proposed in [20]. The runtimes of the mentioned quantum algorithms are both polylogarithmic in the dimensions of the matrix, so they provide exponential speedups over their classical counterparts.

There are several types of tensor decompositions, such as CANDECOMP/PARAFAC (CP) decomposition [9, 16], tensor-train (TT) decomposition [27], Tucker decomposition [42], and etc. However, currently there are no quantum tensor decomposition algorithms. In this paper, we propose a quantum higher order singular value decomposition (Q-HOSVD). HOSVD is a specific orthogonal Tucker decomposition, and can be considered as an extension of SVD from matrices to tensors.

Classical HOSVD has been well studied, see, e.g., De Lathauwer, De Moor, and Vandewalle in 2000 [10], and it has been successfully applied to signal processing [26] and pattern recognition [43] problems. Furthermore, HOSVD has shown its strong power in quantum chemistry, especially in the second order Møller Plesset perturbation theory calculations [2]. In addition, HOSVD is used in [46] to derive the output m photon state of a quantum linear passive system which is driven by an m photon input state; more specifically, the wave function of the output is expressed in terms of the HOSVD of the input wave function.

Since HOSVD deals with high dimensional data, it has been put into practice in some machine learning methods. For example, it has been successfully applied in recommendation systems [19, 40]. In [36], HOSVD representation for neural networks is proposed. By applying HOSVD the parameter-varying system can be expressed in a tensor product form by locally tuned neural network models. Additionally in [21], HOSVD is applied for com-

pressing convolutional neural networks (CNN).

Our Q-HOSVD algorithms are based upon the quantum matrix singular value decomposition algorithm [35], quantum singular value estimation [20], and several other quantum computing techniques. The input can be a tensor of any order and dimension. By our Q-HOSVD algorithms, it is possible to perform singular value decomposition on tensors exponentially faster than classical algorithms. It can be directly applied to quantum machine learning algorithms, and may help solve computationally challenging problems arising in quantum mechanics and chemistry.

Compared to our preliminary research [14], the current paper provides two different algorithms to implement quantum HOSVD. Algorithm 1 is similar to that in [14]; however, here we present more details on the derivations, and in particular build the SWAP operator for quantum principal component analysis (qPCA) in a direct way. An alternative quantum HOSVD algorithm, Algorithm 2, is proposed in the current paper which is based on quantum singular value estimation [20]. Moreover, a thorough explanation for the computational complexity of these two algorithms is presented. Finally, we propose a hybrid quantum-classical algorithm of HOSVD model applied in recommendation systems.

The remainder of this paper is organized as follows. Some preliminaries are given in Section II. Two quantum higher order singular value decomposition algorithms are presented in Sections III and IV respectively. The computational complexity is discussed in Section V. In Section VI, we give an application of HOSVD model on quantum recommendation systems. In the last section, we summarize the results and compare the quantum HOSVD algorithm with the classical counterpart.

## II. PRELIMINARIES

First, we would like to add a comment on the notation that is used. Different symbols are used to facilitate the distinction among scalars, vectors, matrices, and tensors. Scalars are denoted by both lower-case letters ( a, b, . . . ; α, β, . . . ) and capital letters ( A,B,... ). Boldface lower-case letters ( a b , , . . . ) represent vectors. Since the algorithms we present are quantum algorithms, vectors are represented as quantum states after this section, ket |·〉 denotes a column vector, and bra 〈·| denotes a row vector. Bold-face capitals ( A B , , . . . ) correspond to matrices or operators, and tensors are written as calligraphic letters ( A B , , . . . ).

For a matrix M ∈ C m n × , its singular value decomposition (SVD) reads

$$M = U \Sigma V ^ { \dagger }, \quad \ \ ( 1 )$$

where U is an m × m unitary matrix, Σ is a diagonal m × n matrix with non-negative real values on the diagonal, V is an n × n unitary matrix and V † is the conjugate transpose of V . The decomposition can also be written as

$$M = \sum _ { i = 1 } ^ { r } \sigma _ { i } \mathbf u _ { i } \mathbf v _ { i } ^ { \dagger },$$

where r is the rank of M , σ i is the i -th largest singular value, and u i and v i are the corresponding left and right singular vectors respectively.

Denote [ m ] ≡ { 1 2 , , · · · , m } . An m th-order tensor A = ( a i 1 ··· i m ) is a multi-array of Π m j =1 I j entries, where i j ∈ [ I j ] for j ∈ [ m ]. ( I 1 , I 2 , . . . , I m ) is the dimension of A . When I 1 = I 2 = · · · = I m = n , A is called an m th-order n -dimensional tensor.

Given a tensor A ∈ C I 1 × I 2 ×···× I m and a matrix B ∈ C J k × I k , their k -mode tensor-matrix multiplication

$$& ( \mathcal { A } \times _ { k } \mathbf B ) _ { i _ { 1 } i _ { 2 } \cdots i _ { k - 1 } j _ { k } i _ { k + 1 } \cdots i _ { m } } \\ = & \sum _ { i _ { k } = 1 } ^ { I _ { k } } a _ { i _ { 1 } i _ { 2 } \cdots i _ { k - 1 } i _ { k } i _ { k + 1 } \cdots i _ { m } } b _ { j _ { k } i _ { k } }, \quad \quad ( 3 )$$

produces an I 1 × I 2 ×···× I k -1 × J k × I k +1 ×···× I m tensor. The inner product of two tensors A B ∈ , C I 1 × I 2 ×···× I m , denoted as A· B , is defined as

$$\mathcal { A } \cdot \mathcal { B } = \sum _ { i _ { 1 } = 1 } ^ { I _ { 1 } } \sum _ { i _ { 2 } = 1 } ^ { I _ { 2 } } \cdots \sum _ { i _ { m } = 1 } ^ { I _ { m } } a _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } ^ { * } b _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }. \quad ( 4 )$$

Similar to the matrix case, the induced norm √ A· A is called the Frobenius norm of A , denoted as ‖A‖ F . The l 1 -norm of the tensor A is defined as ‖A‖ 1 = ∑ I 1 i 1 =1 ∑ I 2 i 2 =1 · · · ∑ I m i m =1 | a i 1 2 i ··· i m | .

For tensor A ∈ C I 1 × I 2 ×···× I m , if there exist matrices X ( k ) = [ x ( k ) 1 x ( k ) 2 · · · x ( k ) I k ] ∈ C I k × I k with ‖ x ( k ) i k ‖ = 1 for k ∈ [ m ] and i k ∈ [ I k ] such that

$$\mathcal { A } = \mathcal { S } \times _ { 1 } \mathbf X ^ { ( 1 ) } \times _ { 2 } \mathbf X ^ { ( 2 ) } \times _ { 3 } \cdots \times _ { m } \mathbf X ^ { ( m ) }, \quad ( 5 )$$

then (5) is said to be a Tucker decomposition of A , and S = ( s i 1 2 i ··· i m ) is called the core tensor of A . Higher order singular value decomposition is a specific orthogonal Tucker decomposition. Specifically, for A ∈ C I 1 × I 2 ×···× I m , the HOSVD [10] is defined as

$$\mathcal { A } = \mathcal { S } \times _ { 1 } \mathbf U ^ { ( 1 ) } \times _ { 2 } \mathbf U ^ { ( 2 ) } \times _ { 3 } \cdots \times _ { m } \mathbf U ^ { ( m ) }, \quad ( 6 )$$

where the k -mode singular matrix U ( k ) = [ u ( k ) 1 u ( k ) 2 · · · u ( k ) I k ] is a complex unitary I k × I k matrix, the core tensor S ∈ C I 1 × I 2 ×···× I m and its subtensors S i k = α , of which the k th index is fixed to α ∈ [ I k ], have the following properties.

(i) all-orthogonality:

Two subtensors S i k = α and S i k = β are orthogonal for k = 1 2 , , · · · , m :

glyph[negationslash]

$$\mathcal { S } _ { i _ { k } = \alpha } \cdot \mathcal { S } _ { i _ { k } = \beta } = 0 \quad \text{when} \ \alpha \neq \beta, \quad \text{(7)}$$

FIG. 1. Block diagram of the HOSVD for a third-order tensor. The full lines indicate the full HOSVD in (6). The dashed lines and the small block in S indicate the truncated HOSVD.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00719v2-figure-69403414-8cc6-4ea1-a27e-d898f4047323.jpg)

FIG. 2. The tensor network notation of equation (6).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00719v2-figure-8d3ab2ce-2abd-43de-a1f1-a7061cca4506.jpg)

(ii) ordering:

Similar to the matrix case, the tensor singular values are defined as the Frobenius norms of the ( m -1)th-order subtensors of the core tensor S :

$$\sigma _ { \alpha } ^ { ( k ) } = \| \mathcal { S } _ { i _ { k } = \alpha } \| _ { F }, \quad \ \ \ \ ( 8 ) \quad \ \ \ \ \ r i t ^ { \cdot }$$

for k ∈ [ m ] and α ∈ [ I k ]. Furthermore, these tensor singular values have the following ordering property

$$\sigma _ { 1 } ^ { ( k ) } \geq \sigma _ { 2 } ^ { ( k ) } \geq \cdots \geq \sigma _ { I _ { k } } ^ { ( k ) } \geq 0 \quad \quad ( 9 ) \quad \frac { \infty } { \text{Im} }$$

for k ∈ [ m ]. The block diagram of the HOSVD for a third-order tensor A ∈ C I 1 × I 2 × I 3 is described in Fig. 1. When m = 2, i.e., A is a matrix, the HOSVD is degenerated to the well-known matrix SVD.

HOSVD performs orthogonal coordinate transformations for a higher-order tensor. Here, the unitary matrix U ( k ) is also called the k -mode factor matrix and considered as the principal components in k th mode. Moreover, the entries of the core tensor S show the level of interaction among different components.

The tensor network notation of HOSVD is depicted in Fig. 2. It is known that a tensor corresponds to a multipartite quantum state. Any local unitary transformation on the original tensor can be considered as the local unitary transformation on the corresponding singular matrices, which is vital and useful in quantum computation. The core tensor and singular matrices can also be considered as the two layers in the neural network, with local operations in the first layer and global operations in the second layer.

glyph[negationslash]

For an m th-order tensor A ∈ C I 1 × I 2 × ... × I m , the modek matrix unfolding A ( k ) ∈ C I k × (Π j = k I j ) contains the element a i 1 ··· i m at the position with row number i k and column number

$$& ( i _ { k + 1 } - 1 ) I _ { k + 2 } I _ { k + 3 } \cdots I _ { m } I _ { 1 } I _ { 2 } \cdots I _ { k - 1 } \\ + & ( i _ { k + 2 } - 1 ) I _ { k + 3 } I _ { k + 4 } \cdots I _ { m } I _ { 1 } I _ { 2 } \cdots I _ { k - 1 } + \cdots \\ + & ( i _ { m } - 1 ) I _ { 1 } I _ { 2 } \cdots I _ { k - 1 } + ( i _ { 1 } - 1 ) I _ { 2 } I _ { 3 } \cdots I _ { k - 1 } \\ + & ( i _ { 2 } - 1 ) I _ { 3 } I _ { 4 } \cdots I _ { k - 1 } + \cdots + i _ { k - 1 }.$$

By the above construction, the rank of A ( k ) is at most I k . Clearly, the elements of tensor A and unfolding matrix A ( k ) have a one-to-one correspondence to each other.

In HOSVD, the columns of U ( k ) have been sorted such that the j th column u ( k ) j corresponds to the j th largest nonzero singular value of A ( k ) . Then, we can similarly define the truncated (or compact) HOSVD [41]. For k ∈ [ m ], we remain the first r k columns of U ( k ) , then U ( k ) ∈ C I k × r k . Finally, the core tensor S is of size r 1 × r 2 ×·· ·× r m , and the tuple of numbers ( r , r 1 2 , · · · , r m ) is called a multilinear rank. The block diagram of the truncated HOSVD for a third-order tensor is depicted in Fig. 1. This truncation is widely used in big data problems. Since the data may be sparse or low-rank, we can take the value of r k such that r k glyph[lessmuch] I k . Denote r = max k ∈ [ m ] r k , and I = max k ∈ [ m ] I k . The total number of entries reduces from I m to r m + mIr .

## III. Q-HOSVD ALGORITHM 1

In this section, we present our first Q-HOSVD algorithm.

Algorithm 1 Quantum Higher Order Singular Value Decomposition (Q-HOSVD)

Input : A ∈ C I 1 ×···× I m , glyph[epsilon1], | b 〉

:

Output

S

,

U

(1)

,

U

(2)

, . . . ,

U

1. Load A into qRAM.

for

k

= 1

, . . . , m do

2. Implement qPCA (18) by the SWAP operator S ( k ) ˜ A in (17).

3. Initialize | 0 〉 | glyph[vector] 1 〉 | b 〉 , then apply phase estimation to obtain

$$\mathfrak { l } _ { \psi } \rangle = \sum _ { j = 0 } ^ { r - 1 } \beta _ { j } \left | \tilde { \lambda } _ { j } / N \right \rangle \left | \tilde { u } _ { j } \right \rangle. \quad \quad ( 1 0 )$$

4. Perform measurement on | ˜ λ /N j 〉 and extract | ˜ u j 〉 to compose U ( k ) .

$$\text{end for}$$

$$\int _ {. } ^ { \infty } \int _ { 5. \ \mathcal { S } \leftarrow \mathcal { A } \times _ { 1 } \mathbf U ^ { ( 1 ) } \times _ { 2 } \mathbf U ^ { ( 2 ) } \times _ { 3 } \cdots \times _ { m } \mathbf U ^ { ( m ) }. }$$

Several techniques and subroutines are applied in Algorithm 1. First, tensor A to be decomposed is loaded into the quantum register by qRAM. For a fixed k , we design a SWAP operator S ( k ) ˜ A based on matrix unfolding and Hermitian extension, and harness it to apply qPCA. After that, we initialize the state | 0 〉 | glyph[vector] 1 〉 | b 〉 , where | b 〉 could be any state and considered as a superposition of eigenstates of ˜ A ( k ) , then apply phase estimation on it to obtain

(

m

)

the state | ψ 〉 which is a superposition state composed of eigenvalues and eigenstates. Then, if we hope to obtain the core tensor S , quantum measurement is performed to reconstruct the singular matrices. Finally, S is calculated by the quantum tensor-matrix multiplication among tensor A and the singular matrices.

In the following subsections, we explain the implementation of Algorithm 1 step by step. Without loss of generality, we assume ‖A‖ F = 1, and I 1 = I 2 = · · · = I m = n .

## A. Step 1

A vector x = ( x , x 1 2 , . . . , x n ) T ∈ R n with unit norm can be loaded into a quantum register by an oracle named quantum random access memory (qRAM) [12]:

↦

$$\mathbf x \mapsto | x \rangle = \sum _ { i = 0 } ^ { n - 1 } x _ { i + 1 } | i \rangle \quad \quad ( 1 1 ) \quad \ \ t r e _ { \substack { \mathbf t o _ { F } \\ \infty \varepsilon } }$$

with preparation time O (log n ) by a unitary operator U x . Note that in quantum computing the indices usually count from 0.

Similarly, a tensor A = ( a i 1 2 i ··· i m ) can be accessed by the following multipartite state

$$| \Psi \rangle = \sum _ { i _ { 1 }, i _ { 2 }, \dots, i _ { m } = 0 } ^ { n - 1 } a _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } | i _ { 1 } i _ { 2 } \cdots i _ { m } \rangle \,, \quad ( 1 2 ) \quad \text{ due } c \quad \text{ due } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{\quad } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text$$

where i k = 0 , . . . , n -1 for k ∈ [ m ]. This procedure can be achieved using O m ( log n ) operations by U A .

A tree data structure is presented to implement the operation of qRAM [20, 29], which is a classical data structure with quantum access. In this paper, we extend this data structure for real matrices to both real and complex m th-order tensors and summarize it in the following theorem.

Theorem 1. For a tensor A ∈ R C ( ) I 1 × I 2 ×···× I m , there exists a data structure for storing this tensor with quantum access in time O (polylog( I 1 I 2 · · · I m )) .

Proof. We prepare a series of unitary operators and append an ancilla | 0 〉 at the front when applying every op- erator

$$\underset { n - 1 } { \text{to} } \underset { \dots } { \mathbf U } _ { m } \colon | 0 \rangle \rightarrow | \phi _ { m } \rangle = \frac { 1 } { \| \mathcal { A } \| _ { F } } \sum _ { i _ { m } = 0 } ^ { I _ { m } - 1 } \| \mathcal { A } ( \colon, \dots, \colon, i _ { m } ) \| _ { F } \, | i _ { m } \rangle$$

U

$$\text{of} \, \text{arator} \, \text{in} \, \text{to} \, \text{$U_{m} \, \colon | 0 \rangle \rightarrow | \phi _ { m } \rangle = \frac { 1 } { \| \mathcal { A } \| _ { F } } \sum _ { i _ { m } = 0 } ^ { I _ { m } - 1 } \| \mathcal { A } ( \cdot, \dots, \cdot ; i _ { m } ) \| _ { F } | i _ { m } \rangle \\ \text{in} \, \text{$U_{m-1} \, \colon | 0 \rangle \, | \phi _ { m } \rangle \rightarrow | \phi _ { m - 1 } \rangle \\ \text{in} \, \text{$=} \, \text{$1 \| \mathcal { A } \| _ { F } } \sum _ { i _ { m - 1 } = 0 } ^ { I _ { m - 1 } - 1 } \sum _ { i _ { m } = 0 } ^ { I _ { m } - 1 } \frac { \| \mathcal { A } ( \cdot, \dots, \cdot ; i _ { m - 1 }, i _ { m } ) \| _ { F } } { \| \mathcal { A } ( \cdot, \dots, \cdot ; i _ { m } ) \| _ { F } } \\ \text{in} \, \text{$=} \, \text{$1 \| \mathcal { A } \| _ { F } } \sum _ { i _ { m - 1 } = 0 } ^ { I _ { m - 1 } - 1 } \sum _ { i _ { m } = 0 } ^ { I _ { m } - 1 } \| \mathcal { A } ( \cdot, \dots, \cdot ; i _ { m - 1 }, i _ { m } ) \| _ { F } \\ \text{in} \, \text{$=} \, \text{$1 \| \mathcal { A } \| _ { F } } \sum _ { i _ { m - 1 } = 0 } ^ { I _ { m - 1 } - 1 } \sum _ { i _ { m } = 0 } ^ { I _ { m } - 1 } \| \mathcal { A } ( \cdot, \dots, \cdot ; i _ { m - 1 }, i _ { m } ) \| _ { F } \\ \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{out} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \, \text{in} \,
\end{cases}$$

$$= & \frac { 1 } { \| \mathcal { A } \| _ { F } } \sum _ { i _ { 1 } = 0 } ^ { I _ { 1 } - 1 } \cdots \sum _ { i _ { m } = 0 } ^ { I _ { m } - 1 } a _ { i _ { 1 } \cdots i _ { m } } \, | i _ { 1 } \rangle \cdots | i _ { m } \rangle \,. \quad ($$

$$( 1 3 )$$

Basically, the data structure consists of several binary trees named as B ( k ) i , i = 0 , . . . , I k -1 , k ∈ [ m ]. The top root stores the Frobenius norm of A . The root of each B ( k ) i stores the value ‖A (: , . . . , : , i k +1 , . . . , i m ) ‖ 2 F for k ∈ [ m -1] and the weight of an interior node is just the sum of the weights of its children. The leaf nodes at the bottom store the weights A ( i 1 , · · · , i m ) 2 and its sign sgn( A ( i 1 , · · · , i m )) for real tensors. For complex tensors, there is one more layer of binary trees for storing the real and imaginary part of a complex entry A ( i 1 , · · · , i m ) and their signs. Therefore, one more qubit is required for storing a complex tensor, and the extra time is O (1) so we can omit it. If a quantum algorithm has access to this data structure, a series of controlled rotations are applied on the initial state in time O ( glyph[ceilingleft] log I 1 glyph[ceilingright] · · · glyph[ceilingleft] log I m glyph[ceilingright] ). We consider it as O (polylog( I 1 I 2 · · · I m )).

For demonstration, the graph illustrations of the data structure for a real and complex 2 × × 2 2 tensor are given in Fig. 3 and 4 respectively.

## B. Step 2

The quantum unfolding matrix A ( k ) = ( a ′ ( k ) i k j k ) can be processed by a SWAP operator U ( k ) SP :

$$& \text{ processed by a SWAP operator } \mathbf U _ { \text{SP} } ^ { ( k ) } \colon \\ & \sum _ { i _ { 1 }, i _ { 2 }, \dots, i _ { m } = 0 } ^ { n - 1 } a _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } | i _ { 1 } i _ { 2 } \cdots i _ { m } \rangle \\ \quad \overset { \mathbf U _ { \text{SP} } ^ { ( k ) } } { \longrightarrow } \sum _ { i _ { 1 }, i _ { 2 }, \dots, i _ { m } = 0 } ^ { n - 1 } a _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } | i _ { k } i _ { k + 1 } \cdots i _ { m } i _ { 1 } \cdots i _ { k - 1 } \rangle \\ & = \sum _ { i _ { k } = 0 } ^ { n - 1 } \sum _ { j _ { k } = 0 } ^ { n ^ { m - 1 } - 1 } a _ { i _ { k } j _ { k } } ^ { \prime ( k ) } | i _ { k } j _ { k } \rangle \,,$$

FIG. 3. The data structure for a real 2 × 2 × 2 tensor.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00719v2-figure-a24222ce-a594-48a1-99e2-7088c33530d9.jpg)

FIG. 4. The data structure for a complex 2 × 2 × 2 tensor.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00719v2-figure-0c953673-0330-44fc-84fd-e882d2203e14.jpg)

where | j k 〉 = | i k +1 · · · i m i 1 · · · i k -1 〉 . For example, for a 2 × 2 × 2 tensor A , the entries correspond to those of mode-3 unfolding matrix A (3) by

$$a _ { 0 0 0 } \left | 0 0 0 \right \rangle & \rightarrow a _ { 0 0 } ^ { \prime ( 3 ) } \left | 0 0 \right \rangle \\ a _ { 0 1 0 } \left | 0 1 0 \right \rangle & \rightarrow a _ { 0 1 } ^ { \prime ( 3 ) } \left | 0 1 \right \rangle \\ & \vdots \quad \quad \quad$$

The corresponding SWAP operator is U (3) SP = (SWAP ⊗ I )( I ⊗ SWAP), where I is a 2 × 2 identity matrix, and SWAP is the well-known 4 × 4 SWAP operator SWAP = ∑ 1 j,glyph[lscript] =0 | glyph[lscript] 〉〈 j | ⊗ | j 〉〈 glyph[lscript] | . The circuit of the above operations and tensor input is given in Fig. 5. Except for mode1 unfolding which does not require SWAP operations, other modek unfoldings require m -1 SWAP operations. Combined with the complexity of input in Step 1, the total complexity is still O m ( log n ).

$$\text{Denote } N \ = \ n + n ^ { m - 1 }. \ \text{Since } \mathbf A ^ { ( k ) } \ = \ \left ( a _ { i j } ^ { \prime ( k ) } \right ) \ \in \\ \mathbb { C } ^ { n \times n ^ { m - 1 } } \text{ is not a Hermitian matrix, we consider the fol- } \ \text{ where}$$

is not a Hermitian matrix, we consider the fol-

FIG. 5. The quantum circuit to perform mode-3 matrix unfolding of a 2 × 2 × 2 tensor. The unitary operator U A is the one used in Step 1.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00719v2-figure-38cad557-41ab-4031-8e72-7e009253746e.jpg)

lowing extended matrix

$$\mathbb { P } \otimes \text{ extended matrix} \\ \mathbb { P } & \otimes \left [ \begin{matrix} 0 & \mathbb { A } ^ { ( k ) } \\ \mathbb { A } ^ { ( k ) } \end{matrix} \right ] \\ \text{ions} \\ \text{ode-} \\ \text{ions}, \text{ the} \end{matrix} \right ] \$$

where | i k 〉 ∈ C N is the i k -th computational basis. Note

that i k runs from 0 to n -1. Then ˜ A ( k ) is an N × N Hermitian matrix. For Hermitian matrices, the singular values are the absolute value of eigenvalues, so phase estimation [25] can be used to apply the singular value decomposition. Let r = rank( ˜ A ( k ) ). Since rank( A ( k ) ) ≤ n , r ≤ 2 n .

For the Hermitian matrix ˜ A ( k ) , we define a SWAP-like operator S ( k ) ˜ A ∈ C N 2 × N 2 based on the entries of ˜ A ( k ) :

$$( 1 7 )$$

$$\text{operator $S_{A}^{(k)} \in \mathbb{C}^{N^{2}\times N^{2}}$ based on the entries of $\tilde{A}^{(k)$} \quad = t \\ S _ { A } ^ { ( k ) } \coloneqq \sum _ { \ell, j = 0 } ^ { N - 1 } \tilde { A } _ { \ell j } ^ { ( k ) } | j$$

This operator is one-sparse in a quadratically bigger space, i.e., there is no more than one non-zero entry in every row and column, and its entries are efficiently computable. Therefore, the matrix exponentiation e -i S ( k ) ˜ A ∆ t is efficiently implemented [4].

The SWAP-like operator (17) combines the modek unfolding (14) and Hermitian extension (16), since they are all related to SWAP operations. It only requires access to the entries of original tensor A .

For simplicity, in the following we use A to represent A ( k ) when k is fixed. Let ρ 1 and ρ 2 be two distinct density matrices, where ρ 1 = 1 | glyph[vector] 〉 〈 glyph[vector] 1 | with | glyph[vector] 1 〉 = 1 √ N ∑ N -1 glyph[lscript] =0 | glyph[lscript] 〉 .

Lemma 1. [24] By quantum principal component analysis (qPCA), the unitary e -i ˜ A N ∆ t is simulated using S ˜ A through

$$\text{tr} _ { 1 } \{ e ^ { - \text{i} \mathbb { S } _ { \lambda } \Delta t } \rho _ { 1 } \otimes \rho _ { 2 } e ^ { \text{i} \mathbb { S } _ { \lambda } \Delta t } \} \approx e ^ { - \text{i} \mathbb { A } \Delta t } \rho _ { 2 } e ^ { \text{i} \mathbb { A } \Delta t }. \ \ ( 1 8 ) \quad \text{with}$$

Let glyph[epsilon1] 0 be the trace norm of the error term O (∆ t 2 ) in (18). For s steps, the resulting error is glyph[epsilon1] 1 = sglyph[epsilon1] 0 ≤ 2 s ‖A‖ 2 max ∆ , t 2 where ‖A‖ max = max i 1 ,...,i m | a i 1 ··· i m | . The simulated time is t = ∆ . Then, s t

$$\frac { \epsilon _ { 1 } } { s } \leq 2 \| \mathcal { A } \| _ { \max } ^ { 2 } \left ( \frac { t } { s } \right ) ^ { 2 }. \quad \quad \ \ ( 1 9 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Thus,

$$s = O \left ( \frac { t ^ { 2 } } { \epsilon _ { 1 } } \| \mathcal { A } \| _ { \max } ^ { 2 } \right ) \text{ \quad \ \ } ( 2 0 )$$

steps are required to simulate e -i ˜ A N ∆ t if glyph[epsilon1] 1 and t are fixed.

Since we have assumed ‖A‖ F = 1, then ‖A‖ max = O (1). Hence, s = O t /glyph[epsilon1] ( 2 1 ). Applying the output in equation (18) again in the second register, we obtain

$$\ u e c - \quad & \text{equation (18) again in the second register, we obtain} \\ ) \leq n, \\ \Pr \text{-like} ^ { ( k ) } \colon = & \text{tr} _ { 1 } \left \{ e ^ { - i \mathbf S _ { \tilde { \lambda } } \Delta t } \rho _ { 1 } \otimes \left ( \rho _ { 2 } - i \frac { \Delta t } { N } [ \tilde { \alpha }, \rho _ { 2 } ] + O ( \Delta t ^ { 2 } ) \right ) e ^ { i \mathbf S _ { \tilde { \lambda } } \Delta t } \right \} \\ = & \text{tr} _ { 1 } \{ e ^ { - i \mathbf S _ { \tilde { \lambda } } \Delta t } \rho _ { 1 } \otimes \rho _ { 2 } e ^ { i \mathbf S _ { \tilde { \lambda } } \Delta t } \} \\ & - i \frac { \Delta t } { N } \text{tr} _ { 1 } \{ e ^ { - i \mathbf S _ { \tilde { \lambda } } \Delta t } ( \rho _ { 1 } \otimes [ \tilde { \alpha }, \rho _ { 2 } ] ) e ^ { i \mathbf S _ { \tilde { \lambda } } \Delta t } \} + O ( \Delta t ^ { 2 } ) \\ = & \rho _ { 2 } - i \frac { \Delta t } { N } [ \tilde { \alpha }, \rho _ { 2 } ] - i \frac { \Delta t } { N } \text{tr} _ { 1 } \{ \rho _ { 1 } \otimes [ \tilde { \alpha }, \rho _ { 2 } ] \} + O ( \Delta t ^ { 2 } ) \\ = & \rho _ { 2 } - i \frac { 2 \Delta t } { N } [ \tilde { \alpha }, \rho _ { 2 } ] + O ( \Delta t ^ { 2 } ). \\ & \text{ Thus, by-continuously using $k$ copies of $\rho_{1}$ we can simu-}$$

Thus, by continuously using k copies of ρ 1 we can simulate e -i( ˜ A /N k ) ∆ t .

## C. Step 3

Next, we use the quantum phase estimation algorithm to estimate the eigenvalues of e -i( ˜ A /N )∆ t .

Lemma 2. [22] Let U be an n × n unitary operator, with eigenvectors | v j 〉 and eigenvalues e i θ j for θ j ∈ [ -π, π , ] i.e. we have U | v j 〉 = e i θ j | v j 〉 for j = 0 1 , , . . . , n -1 . For a precision parameter glyph[epsilon1] &gt; 0 , there exists a quantum phase estimation algorithm that runs in time O T ( ( U ) log n/glyph[epsilon1] ) and with probability 1 -1 / poly( n ) maps a state ∑ n -1 j =0 α j | v j 〉 to the state ∑ n -1 j =0 α j | v j 〉 ∣ ∣ ¯ θ j 〉 such that ¯ θ j ∈ θ j ± glyph[epsilon1] for all j = 0 1 , , . . . , n -1 .

Theorem 2. For the input | 0 〉 ⊗ d | glyph[vector] 1 〉 | b 〉 , apply (18) in Lemma 1 to the input, and apply phase estimation, then the superposition state (10), where | ˜ λ /N j 〉 is the estimated eigenvalue of ˜ A /N encoded in basis qubits.

Proof. Given an initial quantum state

$$\left | 0 \right \rangle ^ { \otimes d } \left | \vec { 1 } \right \rangle \left | b \right \rangle = \left | \underbrace { 0 0 \cdots 0 } _ { d } \right \rangle \left | \vec { 1 } \right \rangle \left | b \right \rangle \text{ \quad \quad ( 2 2 ) }$$

with d = O ( glyph[ceilingleft] log(1 /glyph[epsilon1] 2 ) glyph[ceilingright] ) control qubits, where | b 〉 is the superposition of eigenvectors | ˜ u j 〉 corresponding to ˜ λ j :

$$| b \rangle = \sum _ { j = 0 } ^ { N - 1 } \beta _ { j } \, | \tilde { u } _ { j } \rangle \,, \quad \sum _ { j = 0 } ^ { N - 1 } | \beta _ { j } | ^ { 2 } = 1, \quad ( 2 3 )$$

glyph[epsilon1] 2 is the accuracy for approximating the eigenvalues. Let ρ 2 = | b 〉 〈 b | . We first apply Hadamard operators to the first register, then the state (22) becomes

$$\frac { 1 } { \sqrt { 2 ^ { d } } } \sum _ { \ell = 0 } ^ { 2 ^ { d } - 1 } \left | \ell \right \rangle \left | \vec { 1 } \right \rangle \left | b \right \rangle,$$

whose density matrix has the following form

$$\frac { 1 } { 2 ^ { d } } \sum _ { \ell = 0 } ^ { 2 ^ { d } - 1 } | \ell \rangle \left \langle \ell \right | \otimes \rho _ { 1 } \otimes \rho _ { 2 }. \quad \ \ ( 2 5 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Then we multiply ∑ 2 d -1 glyph[lscript] =0 | glyph[lscript] 〉 〈 glyph[lscript] | ⊗ ( e -i S A ∆ t ) glyph[lscript] and ∑ 2 d -1 glyph[lscript] =0 | glyph[lscript] 〉 〈 glyph[lscript] | ⊗ ( e i S A ∆ t ) glyph[lscript] to both sides of (25) to obtain

$$\sum _ { \ell = 0 } ^ { 2 ^ { d } - 1 } | \ell \rangle \, \langle \ell | \otimes \left ( ( e ^ { - i \mathbf S _ { A } \Delta t } ) ^ { \ell } \rho _ { 1 } \otimes \rho _ { 2 } ( e ^ { i \mathbf S _ { A } \Delta t } ) ^ { \ell } \right ). \quad ( 2 6 ) \quad \text{Afi}$$

Next, we perform a partial trace to the second register using (18) resulting in

$$\sum _ { \ell = 0 } ^ { 2 ^ { d } - 1 } | \ell \rangle \, \langle \ell | \otimes \left ( ( e ^ { - i \tilde { \frac { n } { N } } \Delta t } ) ^ { \ell } \rho _ { 2 } ( e ^ { i \tilde { \frac { n } { N } } \Delta t } ) ^ { \ell } \right ). \quad ( 2 7 ) \, \begin{array} { c } \text{$tens} \\ \text{$qual} \\ [ 3 7 ]. \end{array}$$

After that, we apply the phase estimation algorithm in Lemma 2 to obtain the estimated eigenvalues of ˜ A /N , since

$$\intertext { e } e ^ { - i \tilde { \mathfrak { A } } \Delta t } \left | b \right \rangle & = \sum _ { j = 0 } ^ { N - 1 } \beta _ { j } e ^ { - i \tilde { \mathfrak { A } } \Delta t } \left | \tilde { u } _ { j } \right \rangle \\ & = \sum _ { j = 0 } ^ { N - 1 } \beta _ { j } e ^ { - i \lambda _ { j } ( \tilde { \mathfrak { A } } ) \Delta t } \left | \tilde { u } _ { j } \right \rangle. \quad \text{(28)} \quad \text{where} \\. \int _ { j = 0 } ^ { N - 1 } \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \end{cases}$$

At last, we implement the inverse quantum Fourier transform [25] and remove the first register, the final state (10)

$$\left | \psi \right \rangle = \sum _ { j = 0 } ^ { r - 1 } \beta _ { j } \left | \tilde { \lambda } _ { j } / N \right \rangle \left | \tilde { u } _ { j } \right \rangle$$

is obtained, where | ˜ λ /N j 〉 is the estimated eigenvalue of ˜ A /N encoded in basis qubits. The corresponding eigenvector | ˜ u j 〉 is proportional to ( u j ; ± v j ) ∈ C N , where u j and v j are the left and right singular vectors of ˜ , A and the norm of each subvector u j and v j is 1 / √ 2, independent of their respective lengths n and n m -1 .

## D. Step 4

Since A is of size n × n m -1 , A has at most n singular values { σ j } . As a result, ˜ A has at most 2 n nonzero eigenvalues ˜ λ j ∈ {± σ j } . Next, we measure the first register of state (10) in the computational basis {| 0 〉 , · · · , | 2 d - 〉} 1 , all eigenpairs | ˜ λ /N j 〉 | ˜ u j 〉 are obtained with probability | β j | 2 . Discarding the first register and projecting | ˜ u j 〉 onto the u j part by using projection operators P u = ∑ n -1 i =0 | i 〉 〈 i | and P v = ∑ n m -1 + n -1 i = 1 n | i 〉 〈 i | result in | u j 〉 with probability 〈 ˜ u j | u j , 0 〉 = 2 . Then, the singular matrix U is calculated by

$$\mathbf U = \sum _ { j = 1 } ^ { n } | \mathbf u _ { j } \rangle \left \langle j \right |.$$

Repeating measurements with the initial state | b 〉 = | 0 〉 , | 1 〉 , · · · , | n - 〉 1 and applying amplitude amplification [1], we can obtain all the singular vectors in T U = O n ( 3 / 2 ) times with probability close to 1. Thus, the singular matrix U ( k ) is reconstructed.

## E. Step 5

After we get all U ( k ) for k = 1 2 , , . . . , m , in this step we calculate the core tensor S :

$$\mathcal { S } = \mathcal { A } \times _ { 1 } \mathbf U ^ { ( 1 ) ^ { \dagger } } \times _ { 2 } \mathbf U ^ { ( 2 ) ^ { \dagger } } \times _ { 3 } \cdots \times _ { m } \mathbf U ^ { ( m ) ^ { \dagger } }. \quad ( 3 0 )$$

Here, the calculation is accelerated by the quantum tensor-matrix multiplication, which is similar to the quantum matrix multiplication algorithm by swap test [37]. We may calculate A× k U ( k ) † through the following state

$$\ n & \text{ state} \\ \n & \frac { 1 } { \| \mathcal { A } \| _ { F } \| \mathbf U ^ { ( k ) } \| _ { F } } \sum _ { i _ { 1 }, \dots, i _ { k - 1 }, j _ { k }, i _ { k + 1 }, \dots, i _ { m } = 0 } ^ { n - 1 } \| \mathbf U ^ { ( k ) } _ { \bullet } j _ { k } \| _ { 2 } \\ \| \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet } i _ { k + 1 } \cdots i _ { m } \| _ { 2 } \langle \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet } i _ { k + 1 } \cdots i _ { m } | \mathbf U ^ { ( k ) } _ { \bullet } j _ { k } \rangle \\ & \quad | i _ { 1 }, \cdots, i _ { k - 1 }, j _ { k }, i _ { k + 1 }, \cdots, i _ { m } \rangle \, | 0 \rangle + | 0 \rangle ^ { \perp } \,, \quad ( 3 1 ) \\ \text{where} \, + \, A. \quad. \quad. \ \, \text{$\,$i$ on $n$ local comment $contents $to $len$}.$$

where |A i 1 ··· i k -1 · i k +1 ··· i m 〉 is an n -level quantum state ( n -entry vector) if i 1 , . . . , i k -1 , i k +1 , . . . , i m are all fixed.

By (31), the success probability is

$$& \text{By (31), the success probability is} \\ \underset { \begin{array} { c } { n } { 1 }, \dots, i _ { k - 1 }, j _ { k }, i _ { k + 1 }, \dots, i _ { m } = 0 } { \sum } \\ { i _ { 1 }, \dots, i _ { k - 1 }, j _ { k }, i _ { k + 1 }, \dots, i _ { m } = 0 } \end{array} } { \sum } \| \mathsf U _ { \bullet j _ { k } } ^ { ( k ) } \| _ { 2 } ^ { 2 } \| \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet i _ { k + 1 } \cdots i _ { m } } \| _ { 2 } ^ { 2 } \\ & \langle \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet i _ { k + 1 } \cdots i _ { m } } | \mathsf U _ { \bullet j _ { k } } ^ { ( k ) } \rangle _ { ^ { 2 } } \Big / \Big ( \| \mathcal { A } \| _ { F } ^ { 2 } \| \mathsf U ^ { ( k ) } \| _ { F } ^ { 2 } \Big ) \\ e \text{ of } \quad = & \frac { \| \mathcal { A } \times _ { k } \mathsf U ^ { ( k ) ^ { \dagger } } \| _ { F } ^ { 2 } } { \| \mathcal { A } \| _ { F } ^ { 2 } \| \mathsf U ^ { ( k ) } \| _ { F } ^ { 2 } }.$$

Note that unitary matrices preserve norms, and we have assumed that ‖A‖ F = 1, therefore

$$\| \mathcal { A } \times _ { k } \mathbf U ^ { ( k ) ^ { \dagger } } \| _ { F } = \| \mathcal { A } \| _ { F } = 1. \quad \quad ( 3 3 )$$

Thus, after post-selecting | 0 , the state (31) becomes 〉

$$\sin _ { 1 } ^ { 2 } \sum _ { \substack { | \psi ^ { ( k ) } \rangle \coloneqq \\ \sum _ { i _ { 1 }, \dots, i _ { k - 1 }, j _ { k }, i _ { k + 1 }, \dots, i _ { m } = 0 } } } ^ { n - 1 } \| \mathbf U _ { \bullet j _ { k } } ^ { ( k ) } \| _ { 2 } \\ \sum _ { \substack { 2 n \\ \beta = 1 } } ^ { 2 } \sum _ { \substack { | \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet i _ { k + 1 } \cdots i _ { m } } \| _ { 2 } \langle \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet i _ { k + 1 } \cdots i _ { m } } | \mathbf U _ { \bullet j _ { k } } ^ { ( k ) } \rangle \\ \sum _ { \substack { | i _ { 1 }, \cdots, i _ { k - 1 }, j _ { k }, i _ { k + 1 }, \cdots, i _ { m } \rangle \,. \\ \dots / \nwarrow } } ^ { n - 1 } \sum _ { \substack { | \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet i _ { k + 1 } \cdots i _ { m } } \| _ { 2 } \langle \mathcal { A } _ { i _ { 1 } \cdots i _ { k - 1 } \bullet i _ { k + 1 } \cdots i _ { m } } \rangle \,. \\ \dots / \nwarrow } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrows } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrower } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack \\ \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow \ } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { n - 1 } } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \ dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum_{ \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1. } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarRow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack ${ \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \\ \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum _ { \substack { \dots \nwarrow } } ^ { n - 1 } \sum$$

| ψ ( k ) 〉 corresponds to the tensor A× k U ( k ) † , whose amplitudes are exactly the entries of tensor A× k U ( k ) † .

After applying amplitude amplification [1], the computational complexity is

$$\dots \dots \dots \dots _ { F } \dots _ { \nu } \dots \\ T _ { M } & = \tilde { O } \left ( \frac { \| \mathcal { A } \| _ { F } \| \mathbf U ^ { ( k ) } \| _ { F } } { \epsilon _ { 3 } \| \mathcal { A } \times _ { k } \mathbf U ^ { ( k ) ^ { \dagger } } \| _ { F } } \right ) \\ & = \tilde { O } \left ( \frac { \| \mathbf U ^ { ( k ) } \| _ { F } } { \epsilon _ { 3 } } \right ) = \tilde { O } \left ( \frac { \sqrt { n } } { \epsilon _ { 3 } } \right ) \quad \quad ( 3 5 )$$

to accuracy glyph[epsilon1] 3 .

Repeating the multiplication between A and U ( k ) † for k = 1 2 , , . . . , m , we obtain the final state

$$| \psi _ { f } \rangle = \sum _ { j _ { 1 }, \dots, j _ { m } = 0 } ^ { n - 1 } s _ { j _ { 1 } \cdots j _ { m } } \, | j _ { 1 }, \cdots, j _ { m } \rangle \,, \quad ( 3 6 )$$

corresponding to the core tensor S . The total complexity is ˜ O ( m n √ glyph[epsilon1] 3 ) .

Without loss of generality, we can regard the accuracy of matrix exponentiation, phase estimation and tensormatrix multiplication as the same value, i.e. glyph[epsilon1] 1 = glyph[epsilon1] 2 = glyph[epsilon1] 3 =: glyph[epsilon1] .

## IV. Q-HOSVD ALGORITHM 2

For Algorithm 1 to be efficiently implemented, the unfolding matrices are required to be low-rank. This result is analyzed in Section V for complexity analysis. However, in some cases the input may not have such good structure. In this section we propose an alternative quantum HOSVD algorithm which is based on quantum singular value estimation (QSVE) [20]. In this method, the input is a general matrix, not required to be sparse, lowrank or well-conditioned.

In [20, 29], the authors introduce a classical tree data structure with quantum access, stated in Lemma 3, where the quantum states are efficiently prepared corresponding to the rows and columns of matrices. Based on this data structure, a fast quantum algorithm to perform singular value estimation stated in 4 is designed.

Lemma 3. [20] Consider a matrix A = ( a ij ) ∈ R N 1 × N 2 with ω nonzero entries. Let A i be its i -th row of A , and ˆ A = 1 ‖ A ‖ F [ ‖ A 0 ‖ 2 , ‖ A 1 ‖ 2 , · · · , ‖ A N 1 -1 ‖ 2 ] glyph[latticetop] . There exists a data structure storing the matrix A in O ω ( log 2 ( N N 1 2 )) space such that a quantum algorithm having access to the data structure can perform the mapping

$$\ a a t a \ s t r u c t u r e \, c a n \, p e r \rho n t u e \ m a p p i n g & \, \stackrel { \alpha \colon } { \text{data} } \\ \mathbf U _ { P } \colon | i \rangle \, | 0 \rangle \rightarrow | i \rangle \, | \mathbf A _ { i } \rangle = \frac { 1 } { \| \mathbf A _ { i } \| _ { 2 } } \sum _ { j = 0 } ^ { N _ { 2 } - 1 } a _ { i j } \, | i \rangle \, | j \rangle & \sum _ { \substack { \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \pi \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \binom { 1 } { \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta\beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta\binom { 1 } { \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \pi \beta \beta \pi \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \pi \beta \beta \binom { 1 } { \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta\beta \beta \beta \beta \beta \pi \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta\beta \beta \beta \beta \beta \Pi \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \Pi \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta\beta \beta \beta \beta \beta \binom { 1 } { \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \beta \binom { 1 } { \beta \beta \beta \beta \beta } }$$

in time polylog( N N 1 2 ) .

In simple terms, there exists a unitary operator U A that prepares A by

$$\mathbf U _ { A } ( | 0 \rangle ^ { \log N _ { 1 } } \left | 0 \right \rangle ^ { \log N _ { 2 } } ) = \frac { 1 } { \left \| \mathbf A \right \| _ { F } } \sum _ { i, j } a _ { i j } \left | i \right \rangle \left | j \right \rangle \quad ( 3 8 ) \quad \begin{pmatrix} \text{ros } \right | \\ \text{conju} \\ \text{on tb} \end{pmatrix}$$

in O (polylog( N N 1 2 )) time.

Define two isometries P ∈ C N N 1 2 × N 1 and Q ∈ C N N 1 2 × N 2 related to U P and U Q :

$$\mathbf P = \sum _ { i = 0 } ^ { N _ { 1 } - 1 } \left | i \right \rangle \left | \mathbf A _ { i } \right \rangle \left \langle i \right |, \ \mathbf Q = \sum _ { j = 0 } ^ { N _ { 2 } - 1 } \left | \hat { \mathbf A } \right \rangle \left | j \right \rangle \left \langle j \right |. \quad ( 3 9 )$$

It can be proved that P P † = I N 1 , Q Q † = I N 2 , 2 PP † -I N N 1 2 is unitary and it can be efficiently implemented in time O (polylog( N N 1 2 )) in the form of U P and U Q . Actually,

$$2 \text{PP} ^ { \dagger } - I _ { N _ { 1 } N _ { 2 } } & = 2 \sum _ { i } | i \rangle \left | \mathbf A _ { i } \right \rangle \left \langle i | \left \langle \mathbf A _ { i } | - \mathbf I _ { N _ { 1 } N _ { 2 } } \\ & = \text{U} _ { P } \text{G} _ { P } \text{U} _ { P } ^ { \dagger },$$

where G P := 2 ∑ i | i 〉 | 0 〉 〈 i | 〈 0 | -I N N 1 2 is a reflection. Similarly, 2 QQ † -I N N 1 2 = U G U Q Q † Q , where G Q := 2 ∑ j | 0 〉 | j 〉 〈 0 | 〈 j | -I N N 1 2 .

Now denote

$$W & = ( 2 \text{PP} ^ { \dagger } - \mathbf I _ { N _ { 1 } N _ { 2 } } ) ( 2 \text{QQ} Q ^ { \dagger } - \mathbf I _ { N _ { 1 } N _ { 2 } } ) \\ & = U _ { P } G _ { P } U _ { P } ^ { \dagger } U _ { Q } G _ { Q } U _ { Q } ^ { \dagger }.$$

All the factors are unitary operators. Then, W is used for phase estimation to obtain the singular values. The eigenvalue of W is e i θ i such that cos( θ / i 2) = σ / i ‖ A ‖ F , where σ i is the singular value of A . The result of QSVE is summarized in the following lemma:

Lemma 4. [20] Let A ∈ R N 1 × N 2 stored in the data structure stated in Lemma 3, and the singular value decomposition of A can be written as A = ∑ r -1 glyph[lscript] =0 σ glyph[lscript] | u glyph[lscript] 〉 〈 v glyph[lscript] | , where r = min( N ,N 1 2 ) . Let glyph[epsilon1] &gt; 0 be the precision parameter. Then the quantum singular value estimation runs in O (polylog( N N 1 2 ) /glyph[epsilon1] ) and achieves ∑ i β i | v i 〉 | 0 〉 ↦→ ∑ i β i | v i 〉 | σ i 〉 , and | σ i -σ i | ≤ glyph[epsilon1] ‖ A ‖ F for all i with probability at least 1 -1 / poly( N 2 ) .

The circuit of QSVE on a single matrix is shown in Fig. 6. Note that QSVE obtains the superposition state of estimated singular values and the corresponding right singular vectors. Since in HOSVD we aim to obtain the left singular vectors, we can obtain the left singular vectors | u ( k ) i 〉 of A ( k ) from the right singular vectors of its conjugate transpose A ( k ) † . Thus, we perform the QSVE on the unfolding matrix A ( k ) † . The Q-HOSVD algorithm based on QSVE is given in Algorithm 2.

FIG. 6. Circuit of QSVE on matrix A . The second register has t qubits, indicating that the accuracy is 2 -t , and j takes the value of 2 0 , 2 1 , . . . , 2 t -1 respectively. U f maps the singular values of W to those of A .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00719v2-figure-9dfa1435-6a21-43c0-a34d-976f0acb5af8.jpg)

## Algorithm 2 Q-HOSVD by QSVE

$$\overline { \text{Input} \colon \mathcal { A } \in \mathbb { C } ^ { I _ { 1 } \times \dots \times I _ { m } }, \epsilon, | b \rangle } _ { \text{Output} \colon \mathcal { S }, U ^ { ( 1 ) }, U ^ { ( 2 ) }, \dots, U ^ { ( m ) } } _ { \text{Output} \colon \quad \dots }$$

- 1. Prepare the initial state 1 √ m ∑ m -1 k =0 | b 〉 | 0 〉 | k 〉 .
- 2. Implement the controlledk QSVE by the last register to obtain the superposition state

$$| \psi \rangle = \frac { 1 } { \sqrt { m } } \sum _ { k = 0 } ^ { m - 1 } \sum _ { i = 0 } ^ { r - 1 } \beta _ { i } ^ { ( k ) } \, | u _ { i } ^ { ( k ) } \rangle \, | \overline { \sigma } _ { i } ^ { ( k ) } \rangle \, | k \rangle \,.$$

3. Post-select k and perform measurement on | σ ( k ) i 〉 extract | u ( k ) 〉 to compose U ( k ) .

$$4. \, \mathcal { S } \leftarrow \stackrel { \cdot } { \mathcal { A } } \times _ { 1 } ^ { \cdot } \mathbf U ^ { ( 1 ) ^ { \dagger } } \stackrel { \cdot } { \times _ { 2 } } \mathbf U ^ { ( 2 ) ^ { \dagger } } \times _ { 3 } \cdots \times _ { m } \mathbf U ^ { ( m ) ^ { \dagger } }.$$

and i .

Algorithm 2 is similar to Algorithm 1, but we do not need to apply phase estimation on the extended Hermitian matrices.

For Step 1, we prepare the initial state

$$\left | \psi _ { 0 } \right \rangle = \frac { 1 } { \sqrt { m } } \sum _ { k = 0 } ^ { m - 1 } \left | b \right \rangle \left | 0 \right \rangle \left | k \right \rangle, \quad \quad \left ( 4 2 \right ) \quad \text{ we} \\ \text{ is } r \\ \text{ set to}.$$

where the first register could be any state and always expressed as a superposition of eigenstates in (23), the second register stores the eigenvalues after phase estimation, and the last register is the index for modek unfolding.

For Step 2, assume that the tensor A is m th-order n -dimensional for simplicity. Recall that for fixed k , | j k 〉 = | i k +1 · · · i m i 1 · · · i k -1 〉 as the same in (14), where j k = 0 , . . . , n m -1 -1. Denote | A j k 〉 be the tube |A i 1 ··· i k -1 · i k +1 ··· i m 〉 . Different from Algorithm 1, we directly prepare the modek unfolding matrix through the unitary operators U ( k ) P and U ( k ) Q as in Lemma 3 according to the mode of unfolding:

$$\mathbf U _ { P } ^ { ( k ) } \colon | j _ { k } \rangle | 0 \rangle \to | j _ { k } \rangle | \overline { \mathbf A } _ { j _ { k } } \rangle = \frac { 1 } { \| \mathbf A _ { j _ { k } } \| _ { 2 } } \sum _ { i _ { k } = 0 } ^ { n - 1 } \overline { a } _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } | j _ { k } \rangle | i _ { k } \rangle \quad & \text{omul\ } \\ & \text{for } \ j _ { k } = 0, \cdots, n ^ { m - 1 } - 1 ; \quad & \text{sort exp} \\ \underset { m - 1 } { \text{measure} }.$$

$$\mathbf U _ { Q } ^ { ( k ) } \colon | 0 \rangle | i _ { k } \rangle & \to | \hat { \mathbf A } _ { k } \rangle | i _ { k } \rangle = \frac { 1 } { \| \mathcal { A } \| _ { F } } \sum _ { j _ { k } = 0 } ^ { n ^ { m - 1 } - 1 } \| \mathbf A _ { j _ { k } } \| _ { 2 } | j _ { k } \rangle | i _ { k } \rangle \quad \text{measu} \\ & \text{for } \quad i _ { k } = 0, \cdots, n - 1. \quad \text{(43)} \quad \text{specific} \quad \text{$$

By this way, the above two operators are prepared in time O m ( polylog n ), corresponding to A ( k ) † . Then we implement the controlledk singular value estimation in parallel by applying the operator ∑ m -1 k =0 W ( k ) ⊗ | k 〉 〈 k | on the initial state | ψ 0 〉 . Finally, we undo the phase estimation and apply the inverse of operator U ( k ) P , obtaining the state

$$| \psi \rangle = \frac { 1 } { \sqrt { m } } \sum _ { k = 0 } ^ { m - 1 } \sum _ { i = 0 } ^ { r - 1 } \beta _ { i } ^ { ( k ) } \, | u _ { i } ^ { ( k ) } \rangle \, | \overline { \sigma } _ { i } ^ { ( k ) } \rangle \, | k \rangle \,. \quad ( 4 4 )$$

After that, similar to the Steps 4 and 5 in Algorithm 1, we can make measurements and obtain the singular matrices and core tensor.

## V. COMPLEXITY ANALYSIS

For simplicity, we assume the input is an m th-order n -dimensional tensor.

For the first Q-HOSVD algorithm, the computational complexity mainly comes from data access, matrix exponential simulation and phase estimation. The data input time is O m ( polylog n ). At a simulation time t , only the eigenvalues of ˜ A ( k ) /N with | ˜ λ j | /N = Ω(1 /t ) matter [35], and the eigenvalues smaller than glyph[epsilon1] are omitted. Note that A ( k ) is an n × n m -1 matrix. For a fixed k , let the number of these eigenvalues be r ≤ 2 n , then by

$$\begin{cases} \ t r \{ \tilde { \mathbf A } _ { r } ^ { 2 } / N ^ { 2 } \} = \sum _ { j = 0 } ^ { r - 1 } \tilde { \lambda } _ { j } ^ { 2 } / N ^ { 2 } = \Omega ( r / t ^ { 2 } ), \\ \ t r \{ \tilde { \mathbf A } _ { r } ^ { 2 } / N ^ { 2 } \} \leq \tilde { \tau } \{ \tilde { \mathbf A } ^ { 2 } / N ^ { 2 } \} = \| \tilde { \mathbf A } \| _ { F } ^ { 2 } / N ^ { 2 } \leq \| \mathcal { A } \| _ { \max } ^ { 2 }, \\ \quad \text{$\max$ and that $two$ and $two$ to $two$ to $multi$ and $multi$} \end{cases}$$

we find that the rank of the effectively simulated matrix is r = O ( ‖A‖ 2 max t 2 ). By (20), there are O t ( 2 ‖A‖ 2 max /glyph[epsilon1] ) steps required to simulate e -i ˜ A N t , where ‖A‖ max = O (1), and 1 /glyph[epsilon1] can be chosen as O (polylog n ). To make this algorithm efficient, t = O (polylog n ), then the rank r = O (polylog n ), i.e. the matrices have to be low-rank. Thus, the time to implement phase estimation is O (polylog n · log( n /glyph[epsilon1] ) ). Therefore, the total computational complexity of obtaining (10) is O m ( 2 polylog n ).

For the second Q-HOSVD algorithm, for each unfolding the time to access the data structure is O m ( polylog n ), and the time to implement QSVE is O m ( polylog n/glyph[epsilon1] ), so the total complexity of obtaining (44) is O m ( 3 polylog n ). Furthermore, there are no requirement for the structure of the input tensor.

Usually, in practical problems, m glyph[lessmuch] n . Thus, we can omit the order m , and the algorithm runs polylogarithmically in the dimensions.

If we want to obtain the singular matrices and core tensor explicitly in the quantum register, we need to make measurements on the states (10) and (44), and reconstruct singular matrices and finally calculate the core tensor by quantum tensor matrix multiplication, the complexities are O m n ( 3 2 polylog n ) and O m n ( 4 2 polylog n ) respectively.

## VI. APPLICATION ON RECOMMENDATION SYSTEMS

In [40], the authors make use of HOSVD for tag recommendations. Given an initial third-order tensor with usage data triplets (user, item, tag), they implement HOSVD and do truncations to obtain the core tensor and reconstructed tensor with smaller dimensions. Then, based on the entries of the reconstructed tensor, the tags are recommended to users. We have carried out the similar SVD and truncation operations in [44] by another tensor decomposition called t-svd.

In this section, we introduce a hybrid quantum-classical recommendation method for context-aware collaborative filtering (CF) based on tensor factorization (TF), named as multiverse recommendation [19]. TF is an extension of matrix factorization (MF) to multiple dimensions. HOSVD is chosen as our TF approach to analyze the recommendation systems, due to its relevance among the different categories. Given the known preference tensor, we use HOSVD model to find out the missing information. This problem is well-known as the completion problem in recommendation systems [39]. Our contribution is designing a hybrid quantum-classical recommendation algorithm to accelerate this process.

Context has been universally acknowledged as an important factor for analyzing recommendation systems. A pair (user, item) is extended to a triplet (user, item, context) or even larger multiplets, where context denotes the factor that may influence a user's preference on a specified item, e.g. time, location, and we consider the interactions among them. Generally, any number of contexts can be added to this recommendation system, and the correlation is described by a relevance score function f rs as follows:

$$f _ { r s } \colon & \text{User} \times \text{Item} \times \text{Context} _ { 1 } \times \dots \times \text{Context} _ { \ell } \quad \frac { \text{Algo} } { \text{Input} } \\ & \longrightarrow \text{Relevance Score} \quad ( 4 6 ) \quad \text{Input} _ { F }$$

with glyph[lscript] different contexts, so that the total number of dimensions is glyph[lscript] +2. Thus, we can use a tensor with order m = glyph[lscript] + 2 to express the set of relevance scores, and utilize HOSVD model to describe it. Such methods are widely applied in recommendation systems like Netflix prize problems [3] and so on.

Denote the given preference tensor Y ∈ { 0 , . . . , 5 } I 1 × I 2 ×···× I m containing the observed ratings ranging from 1 to 5, and value 0 indicates that the item has not been rated yet. The aim is to find out such missing values and give a good recommendation to users. Denote the factor matrices U ∈ R I 1 × d 1 , M ∈ R I 2 × d 2 , C ( k ) ∈ R I k +2 × d k +2 , for k = 1 , . . . , glyph[lscript] . Then, S ∈ R d 1 × d 2 ×···× d m . Let d = [ d , d 1 2 , · · · , d m ] glyph[latticetop] , and d = max j ∈ [ m ] d j .

To obtain the recommendations based on HOSVD, we design a loss function and optimize over it. The loss function is characterized as L ( T ( θ , ) Y ), where θ is the model parameter, i.e., θ := ( S , U M C , , (1) , C (2) , . . . , C ( glyph[lscript] ) ). De- note a set D := { ( i 1 , i 2 , . . . , i m ) | y i 1 2 i ··· i m &gt; 0 } an observation history, and K := | D | the number of observed ratings. The total loss function is defined as

$$\text{$;rec$} \\ \text{with} \quad L ( \mathcal { T } ( \theta ), \mathcal { Y } ) \coloneqq \frac { 1 } { \| \mathcal { S } \| _ { 1 } } \sum _ { ( i _ { 1 }, i _ { 2 }, \dots, i _ { m } ) \in D } l \left ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } \right ), \\ \text{ensor} \\ \text{hen}. \quad \text{which only online on the overhead while in $\mathcal{V}$                                                                                                                                                                                                         } (47)$$

which only applies on the observed values in Y . Function l ( t, y ) is a pointwise loss function, that can be based on l 2 norm, e.g., l ( t, y ) = ( t -y ) 2 / 2, or other types of distance measure. By adding a regularization term to avoid overfitting, we establish the objective function

$$J ( \theta ) \coloneqq L ( \mathcal { T } ( \theta ), \mathcal { Y } ) + \Omega ( \theta )$$

with trivial regularizers

$$\text{sions.} _ { \substack { \ z e \text{ the} \\ \arg \text{ the} \\ \text{ensor}, } } \Omega ( \theta ) = \lambda _ { 1 } \| \mathbf U \| _ { F } ^ { 2 } + \lambda _ { 2 } \| \mathbf M \| _ { F } ^ { 2 } + \sum _ { k = 1 } ^ { \ell } \lambda _ { k + 2 } \| \mathbf C ^ { ( k ) } \| _ { F } ^ { 2 } + \lambda _ { \mathcal { S } } \| \mathcal { S } \| _ { F } ^ { 2 }. \\ \text{ensor}, } \text{form} _ { \substack { \text{the normalize} \text{ of matrix} \text{sion} \text{ can be chosen as the} } } ( 4 9 )$$

Usually, the parameters of matrices can be chosen as the same value, i.e., λ 1 = λ 2 = · · · = λ m =: λ .

We optimize these matrices and the core tensor by stochastic gradient descent (SGD) method [6]. SGD randomly picks a sample at one time and perform gradient descent. It usually compares with batch gradient descent (BGD) which runs over all the samples each iteration. BGD converges globally in every step but it is computationally prohibitive for our problem. The cost of SGD is low, but it usually converges in a local minimum. For big data problems, SGD often converges without running over all the samples. The whole tensor completion algorithm based on the HOSVD model is given in Algorithm 3.

## Algorithm 3 Tensor Completion by HOSVD

Input:

Y

,

d

, η, λ, λ

S

Output:

S

,

U M C

,

,

(1)

,

· · ·

,

C

(

glyph[lscript]

)

Initialize U M C , , (1) , · · · , C ( glyph[lscript] ) , S with small values, T with all zeros.

$$\text{ Initalize } U, M, C ^ { ( i ) }, \cdots, C ^ { ( e ) }, \mathcal { S } \text{ with small values, } \mathcal { T } \text{ with } \\ \text{ of } \\ \text{$der} \\ \text{ and } \\ \text{$i_{1}_{2} \cdots i_{m}=\mathcal{ }S\times_{1}U_{i}$} \times_{2}M_{i}$} \times_{3}C_{i}^{(1)}\times_{4} \cdots \times_{m}C_{i}^{(e)}\times_{5}
are \text{ $U_{i}$} \mapsto \mathcal{ }i_{i} \mapsto - \eta \lambda U_{i}$} \mapsto - \eta \partial _ { U _ { i } } \mapsto l ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } ) \\ \text{$i_{1}$} \mapsto \leftarrow C _ { i _ { 3 } } ^ { ( 1 ) } \mapsto - \eta \lambda C _ { i _ { 3 } } ^ { ( 1 ) } \mapsto - \eta \partial _ { C ^ { ( 1 ) } _ { i _ { 3 } } } l ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } ) \\ \in \text{$C_{i}^{(1)}\mapsto \leftarrow C_{i _ { 3 } } ^ { ( 1 ) } - \eta \lambda C _ { i _ { 3 } } ^ { ( 1 ) } - \eta \partial _ { C ^ { ( 1 ) } _ { i _ { 3 } } } l ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } ) \\ \text{$rot} \\ \text{the} \\ \text{out} \\ \text{$C_{i}^{(e)}\mapsto \leftarrow C_{i _ { m } } ^ { ( e ) } - \eta \lambda C _ { i _ { m } } ^ { ( e ) } - \eta \partial _ { C ^ { ( e ) } _ { i _ { m } } } l ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } ) \\ \text{$ion} \\ \kappa _ { 1 }, \\ \text{$\mathcal{ }S\times_{1}$} \mapsto \mathcal{ }S\times_{1} \mathcal{ }S\times_{1} \mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$} \text{$\mathcal{ }S\times_{1}$}$$

Algorithm 3 can be considered as a training method by SGD. After we obtain the factor matrices and core tensor, T is computed explicitly by

$$\mathcal { T } = \mathcal { S } \times _ { 1 } \mathbf U \times _ { 2 } \mathbf M \times _ { 3 } \mathbf C ^ { ( 1 ) } \times _ { 4 } \cdots \times _ { m } \mathbf C ^ { ( \ell ) } \quad ( 5 0 )$$

as an approximation of the preference tensor Y and we give recommendations to users according to the entries of T .

This algorithm is a hybrid quantum-classical algorithm. The computation of gradients is accelerated by some quantum subroutines, and the rest procedures are performed by classical computers. Denote ◦ the outer product between tensors. The gradients are, e.g.,

$$& \partial _ { \mathbf U _ { i _ { 1 } } \bullet } l \left ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } \right ) & \text{1o calc} \\ = & \partial _ { t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } } l \left ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } \right ) \mathcal { S } \times _ { 2 } \mathbf M _ { i _ { 2 } \bullet } \times _ { 3 } \cdots \times _ { m } \mathbf C _ { i _ { m } \bullet } ^ { ( \ell ) }, \\ ( 5 1 ) & \text{1o calc}$$

$$& \partial _ { \mathcal { S } } l \left ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } \right ) \\ = & \partial _ { t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } } l \left ( t _ { i _ { 1 } i _ { 2 } \cdots i _ { m } }, y _ { i _ { 1 } i _ { 2 } \cdots i _ { m } } \right ) \mathbf U _ { i _ { 1 } \bullet } \circ \mathbf M _ { i _ { 2 } \bullet } \circ \cdots \circ \mathbf C _ { i _ { m } \bullet } ^ { ( \ell ) } \colon \quad \text{order} \quad \\ & \text{the in} \quad \text{$the in}$} \quad \text{$\text{$th}$}$$

For gradient (51), ∂ t i 1 2 i ··· i m l ( t i 1 2 i ··· i m , y i 1 2 i ··· i m ) is a simple function, if the loss function takes l 2 norm, then ∂ l t ( t, y ) = t -y . Define g U ( i 1 · ) = S × 2 M i 2 · × ··· × 3 m C ( glyph[lscript] ) i m · . The entry g j ( U i 1 · ) is equivalent to

$$\mathcal { S } _ { i _ { 1 } = j } \cdot ( \mathbf M _ { i _ { 2 } \bullet } \circ \dots \circ \mathbf C _ { i _ { m } \bullet } ^ { ( \ell ) } ) = \colon \mathcal { S } _ { i _ { 1 } = j } \cdot \mathcal { Z }, \quad ( 5 3 ) \quad \text{is de}$$

the inner product of two ( m -1)th order tensors. By quantum matrix multiplication algorithm [37], the outer product of two vectors | a 〉 and | b 〉 can be performed in time O ( ‖ a ‖ 3 F ‖ ‖ b 3 F polylog n/glyph[epsilon1] 4 ‖ a ◦ b ‖ 3 F ) = O (polylog n/glyph[epsilon1] 4 ) to accuracy glyph[epsilon1] 4 , which holds for classical vectors. Thus, we can perform the quantum outer product M i 2 · ◦· · · ◦ C ( glyph[lscript] ) i m · in time O (polylog n/glyph[epsilon1] 4 ). The resulting state is

$$| z \rangle = \frac { 1 } { \| \mathcal { Z } \| _ { F } } \sum _ { i _ { 2 }, i _ { 3 }, \dots, i _ { m } = 0 } ^ { n - 1 } z _ { i _ { 2 } i _ { 3 } \cdots i _ { m } } \, | i _ { 2 } \rangle \, | i _ { 3 } \rangle \cdots | i _ { m } \rangle \,. \ \ ( 5 4 ) \quad \text{$corq1$}$$

Next, we load the subtensor S i 1 = j into the quantum register as

$$| s \rangle = \frac { 1 } { \| \mathcal { S } _ { i _ { 1 } = j } \| _ { F } } \sum _ { i _ { 2 }, i _ { 3 }, \dots, i _ { m } = 0 } ^ { n - 1 } s _ { j i _ { 2 } i _ { 3 } \cdots i _ { m } } \, | i _ { 2 } \rangle \, | i _ { 3 } \rangle \cdots | i _ { m } \rangle \,. \quad \text{ is not } \text{comp} \\ \text{Then } \text{ was can content to follow any anonymous action to} \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{  \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{$$

Then, we can construct the following superposition state:

$$| \phi _ { 1 } \rangle & = \frac { 1 } { \sqrt { 2 } } ( | + \rangle \, | s \rangle + | - \rangle \, | z \rangle ) \quad \text{D} \\ & = \sin \theta \, | 0 \rangle \, | u \rangle + \cos \theta \, | 1 \rangle \, | v \rangle \quad \text{(56)} \quad \text{dat} \quad \text{d} \quad \text{d}$$

with cos θ = √ (1 -〈 | 〉 s z ) / 2. Then, by applying quantum amplitude estimation algorithm [7], we can obtain h such that

$$\left | \frac { 1 - \langle s | z \rangle } { 2 } - h \right | \leq \epsilon _ { 5 } \, \quad \ \ ( 5 7 ) \, \quad \ ]$$

in time O m ( log d/glyph[epsilon1] 5 ). Therefore, ‖S i 1 = j ‖ F ‖Z‖ F (1 -2 h ) gives a 2 glyph[epsilon1] 5 ‖S i 1 = j ‖ F ‖Z‖ F -approximate of g U ( i 1 · ). Then, if we take glyph[epsilon1] = 2 glyph[epsilon1] 5 ‖S i 1 = j ‖ F ‖Z‖ F , we can obtain the value of g j ( U i 1 · ) in time O m ( log d ‖S i 1 = j ‖ F ‖Z‖ F /glyph[epsilon1] ) to accuracy glyph[epsilon1] . Since the gradient g has d entries, and we have to repeat the above procedure for all the singular matrices, then the total complexity of matrix optimization is O Km d ( 2 polylog d ). Compared to the corresponding classical algorithm which takes O Kmd ( m ) classical calculations, our quantum algorithm is exponentially faster. To calculate the core tensor S , we can directly use the classical computation, the complexity is O Kd ( m ).

## VII. SUMMARY AND DISCUSSION

Wehave introduced two quantum algorithms for higher order singular value decomposition. For the first method, the input has to be low-rank. For the second one, the input has no constraints on its structure. The output is a core tensor including tensor singular values and singular matrices stored in the quantum register. Some subroutines of the first method has been realized by quantum physicists, such as sparse Hamiltonian simulation and qPCA. For the second method, the QSVE algorithm is derived strictly by using quantum operators.

For an m th-order n -dimensional tensor, the complexity of the classical HOSVD is O mn ( m +1 ), while the complexities of our two quantum HOSVD algorithms for obtaining the superposition state of eigenvalues and eigenstates are O m ( 2 polylog n ) and O m ( 3 polylog n ) respectively. Moreover, the complexities of obtaining singular matrices and core tensor are O m n ( 3 2 polylog n ) and O m n ( 4 2 polylog n ) respectively. Usually, the order m is much smaller than the dimension n . In this sense, our quantum HOSVD algorithms provide an exponential speedup over the classical counterpart.

Furthermore, we have provided an application of HOSVD on recommendation systems by quantum computers. Based on the given preference tensor, we may find out the missing values and give recommendations to users by this HOSVD model. One day the size of datasets is not enough for classical computers to handle, quantum computers will be able to solve these problems efficiently.

## VIII. DATA AVAILABILITY

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.

## IX. AUTHOR CONTRIBUTIONS

L.G., X.W., H.W.J.L and G.Z. developed the two quantum HOSVD algorithms and analyzed the computational complexity. L.G. and G.Z designed the hybrid quantum-classical recommendation method. L.G. sum-

marized the results and wrote the manuscript. X.W., H.W.J.L and G.Z. reviewed and modified the manuscript.

## X. COMPETING INTERESTS

The authors declare that there are no competing interests.

- [1] Andris Ambainis. Variable time amplitude amplification and quantum algorithms for linear algebra problems. In STACS'12 (29th Symposium on Theoretical Aspects of Computer Science) , volume 14, pages 636-647. LIPIcs, 2012.
- [2] Franziska Bell, Daniel S Lambrecht, and Martin HeadGordon. Higher order singular value decomposition in quantum chemistry. Molecular Physics , 108(19-20):27592773, 2010.
- [3] James Bennett and Stan Lanning. The Netflix Prize. In Proceedings of KDD cup and workshop , volume 2007, page 35. New York, NY, USA., 2007.
- [4] Dominic W Berry, Graeme Ahokas, Richard Cleve, and Barry C Sanders. Efficient quantum algorithms for simulating sparse Hamiltonians. Communications in Mathematical Physics , 270(2):359-371, 2007.
- [5] Jacob Biamonte and Ville Bergholm. Tensor networks in a nutshell. arXiv preprint arXiv:1708.00006 , 2017.
- [6] L´ eon Bottou. Large-scale machine learning with stochastic gradient descent. In Proceedings of COMPSTAT'2010 , pages 177-186. Springer, 2010.
- [7] Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp. Quantum amplitude amplification and estimation. Contemporary Mathematics , 305:53-74, 2002.
- [8] X-D Cai, Christian Weedbrook, Z-E Su, M-C Chen, Mile Gu, M-J Zhu, Li Li, Nai-Le Liu, Chao-Yang Lu, and Jian-Wei Pan. Experimental quantum computing to solve systems of linear equations. Physical review letters , 110(23):230501, 2013.
- [9] J Douglas Carroll and Jih-Jie Chang. Analysis of individual differences in multidimensional scaling via an N-way generalization of 'Eckart-Young' decomposition. Psychometrika , 35(3):283-319, 1970.
- [10] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. A multilinear singular value decomposition. SIAM journal on Matrix Analysis and Applications , 21(4):12531278, 2000.
- [11] Vedran Dunjko and Hans J Briegel. Machine learning &amp; artificial intelligence in the quantum domain: a review of recent progress. Reports on Progress in Physics , 81(7):074001, 2018.
- [12] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone. Quantum random access memory. Physical review letters , 100(16):160501, 2008.
- [13] Lov K Grover. A fast quantum mechanical algorithm for database search. In Proceedings of the Twenty-eighth Annual ACM Symposium on Theory of Computing , pages 212-219. New York, NY, USA., 1996.
- [14] Lejia Gu, Xiaoqiang Wang, and Guofeng Zhang. Quantum higher order singular value decomposition. In 2019
- IEEE International Conference on Systems, Man and Cybernetics (SMC) , pages 1166-1171. IEEE, 2019.
- [15] Aram W Harrow, Avinatan Hassidim, and Seth Lloyd. Quantum algorithm for linear systems of equations. Physical review letters , 103(15):150502, 2009.
- [16] Richard A Harshman et al. Foundations of the PARAFAC procedure: Models and conditions for an 'explanatory' multimodal factor analysis. 1970.
- [17] Shenglong Hu, Liqun Qi, and Guofeng Zhang. Computing the geometric measure of entanglement of multipartite pure states by means of non-negative tensors. Physical Review A , 93(1):012304, 2016.
- [18] William Huggins, Piyush Patil, Bradley Mitchell, K Birgitta Whaley, and E Miles Stoudenmire. Towards quantum machine learning with tensor networks. Quantum Science and technology , 4(2):024001, 2019.
- [19] Alexandros Karatzoglou, Xavier Amatriain, Linas Baltrunas, and Nuria Oliver. Multiverse recommendation: n-dimensional tensor factorization for context-aware collaborative filtering. In Proceedings of the fourth ACM conference on Recommender systems , pages 79-86. ACM, 2010.
- [20] Iordanis Kerenidis and Anupam Prakash. Quantum recommendation systems. In Christos H. Papadimitriou, editor, 8th Innovations in Theoretical Computer Science Conference (ITCS 2017) , volume 67 of Leibniz International Proceedings in Informatics (LIPIcs) , pages 49:1-49:21, Dagstuhl, Germany, 2017. Schloss DagstuhlLeibniz-Zentrum fuer Informatik.
- [21] Maksym Kholiavchenko. Iterative low-rank approximation for CNN compression. arXiv preprint arXiv:1803.08995 , 2018.
- [22] A Yu Kitaev. Quantum measurements and the Abelian stabilizer problem. arXiv preprint quant-ph/9511026 , 1995.
- [23] Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review , 51(3):455-500, 2009.
- [24] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost. Quantum principal component analysis. Nature Physics , 10(9):631, 2014.
- [25] Michael A Nielsen and Isaac Chuang. Quantum computation and quantum information, 2002.
- [26] Larsson Omberg, Gene H Golub, and Orly Alter. A tensor higher-order singular value decomposition for integrative analysis of DNA microarray data from different studies. Proceedings of the National Academy of Sciences , 104(47):18371-18376, 2007.
- [27] Ivan V Oseledets. Tensor-train decomposition. SIAM Journal on Scientific Computing , 33(5):2295-2317, 2011.

- [28] Jian Pan, Yudong Cao, Xiwei Yao, Zhaokai Li, Chenyong Ju, Hongwei Chen, Xinhua Peng, Sabre Kais, and Jiangfeng Du. Experimental realization of quantum algorithm for solving linear systems of equations. Physical Review A , 89(2):022313, 2014.
- [29] Anupam Prakash. Quantum algorithms for linear algebra and machine learning . PhD thesis, UC Berkeley, 2014.
- [30] Liqun Qi, Haibin Chen, and Yannan Chen. Tensor eigenvalues and their applications , volume 39. Springer, 2018.
- [31] Liqun Qi and Ziyan Luo. Tensor analysis: spectral theory and special tensors , volume 151. SIAM, 2017.
- [32] Liqun Qi, Guofeng Zhang, Daniel Braun, Fabian BohnetWaldraff, and Olivier Giraud. Regularly decomposable tensors and classical spin states. Communications in Mathematical Sciences , 15(6):1651-1665, 2017.
- [33] Liqun Qi, Guofeng Zhang, and Guyan Ni. How entangled can a multi-party system possibly be? Physics Letters A , 382(22):1465-1471, 2018.
- [34] Patrick Rebentrost, Masoud Mohseni, and Seth Lloyd. Quantum support vector machine for big data classification. Physical review letters , 113(13):130503, 2014.
- [35] Patrick Rebentrost, Adrian Steffens, Iman Marvian, and Seth Lloyd. Quantum singular-value decomposition of nonsparse low-rank matrices. Physical review A , 97(1):012327, 2018.
- [36] Andr´ as R¨ ovid, L´ aszl´ o Szeidl, and P´ eter V´ arlaki. On tensor-product model based representation of neural networks. In 2011 15th IEEE International Conference on Intelligent Engineering Systems , pages 69-72. IEEE, 2011.
- [37] Changpeng Shao. Quantum algorithms to matrix multiplication. arXiv preprint arXiv:1803.01601 , 2018.
- [38] Peter W Shor. Algorithms for quantum computation: discrete logarithms and factoring. In Proceedings 35th annual symposium on foundations of computer science , pages 124-134. IEEE, 1994.
- [39] Qingquan Song, Hancheng Ge, James Caverlee, and Xia Hu. Tensor completion algorithms in big data analytics. ACM Transactions on Knowledge Discovery from Data (TKDD) , 13(1):1-48, 2019.
- [40] Panagiotis Symeonidis, Alexandros Nanopoulos, and Yannis Manolopoulos. Tag recommendations based on tensor dimensionality reduction. In Proceedings of the 2008 ACM conference on Recommender systems , pages 43-50. ACM, 2008.
- [41] L´ aszl´ o Szeidl, Peter Baranyi, Zolt´n Petres, a and Peter Varlaki. Numerical reconstruction of the HOSVD based canonical form of polytopic dynamic models. In 2007 International Symposium on Computational Intelligence and Intelligent Informatics , pages 111-116. IEEE, 2007.
- [42] Ledyard R Tucker. Some mathematical notes on threemode factor analysis. Psychometrika , 31(3):279-311, 1966.
- [43] M Alex O Vasilescu and Demetri Terzopoulos. Multilinear analysis of image ensembles: Tensorfaces. In European Conference on Computer Vision , pages 447-460. Springer, 2002.
- [44] Xiaoqiang Wang, Lejia Gu, Heung-wing Joseph Lee, and Guofeng Zhang. Quantum tensor singular value decomposition with applications to recommendation systems. arXiv preprint arXiv:1910.01262 , 2019.
- [45] Nathan Wiebe, Daniel Braun, and Seth Lloyd. Quantum algorithm for data fitting. Physical review letters , 109(5):050505, 2012.
- [46] Guofeng Zhang. Dynamical analysis of quantum linear systems driven by multi-channel multi-photon states. Automatica , 83:186-198, 2017.
- [47] Mengshi Zhang, Guyan Ni, and Guofeng Zhang. Iterative methods for computing U-eigenvalues of non-symmetric complex tensors with application in quantum entanglement. Computational Optimization and Applications , pages 1-20, 2019.
- [48] Zhikuan Zhao, Alejandro Pozas-Kerstjens, Patrick Rebentrost, and Peter Wittek. Bayesian deep learning on a quantum computer. Quantum Machine Intelligence , 1(1-2):41-51, 2019.
Lejia Gu Xiaoqiang Wang H. W. Joseph Lee Guofeng Zhang
Quantum Higher Order Singular Value Decomposition
Higher order singular value decomposition (HOSVD) is an important tool for
analyzing big data in multilinear algebra and machine learning. In this paper,
we present two quantum algorithms for HOSVD. Our methods allow one to decompose
a tensor into a core tensor containing tensor singular values and some unitary
matrices by quantum computers. Compared to the classical HOSVD algorithm, our
quantum algorithms provide an exponential speedup. Furthermore, we introduce a
hybrid quantum-classical algorithm of HOSVD model applied in recommendation
systems.