Title: LEARNING EFFICIENT TENSOR REPRESENTATIONS WITH RING STRUCTURE NETWORKS

Abstract: Tensor train (TT) decomposition is a powerful representation for high-order tensors, which has been successfully applied to various machine learning tasks in recent years. In this paper, we propose a more generalized tensor decomposition with ring structure network by employing circular multilinear products over a sequence of lower-order core tensors, which is termed as TR representation. Several learning algorithms including blockwise ALS with adaptive tensor ranks and SGD with high scalability are presented. Furthermore, the mathematical properties are investigated, which enables us to perform basic algebra operations in a computationally efficiently way by using TR representations. Experimental results on synthetic signals and real-world datasets demonstrate the effectiveness of TR model and the learning algorithms. In particular, we show that the structure information and high-order correlations within a 2D image can be captured efficiently by employing an appropriate tensorization and TR decomposition.

1. INTRODUCTION:

Tensor decompositions aim to represent a higher-order (or multi-dimensional) data as a multilinear product of several latent factors, which attracted considerable attentions in machine learning (Yu & Liu, 2016; Anandkumar et al., 2014; Romera-Paredes et al., 2013; Kanagawa et al., 2016; Yang et al., 2017) and signal processing (Zhou et al., 2016) in recent years. For a dth-order tensor with "square" core tensor of size r, standard tensor decompositions are the canonical polyadic (CP) decomposition (Goulart et al., 2015) which represents data as a sum of rank-one tensors by O(dnr) parameters and Tucker decomposition (De Lathauwer et al., 2000; Xu et al., 2012; Wu et al., 2014; Zhe et al., 2016) which represents data as a core tensor and several factor matrices by O(dnr + r d ) parameters. In general, CP decomposition provides a compact representation but with difficulties in finding the optimal solution, while Tucker decomposition is stable and flexible but its number of parameters scales exponentially to the tensor order.

Recently, tensor networks have emerged as a powerful tool for analyzing very high-order tensors (Cichocki et al., 2016) . A powerful tensor network is tensor train / matrix product states (TT/MPS) representation (Oseledets, 2011) , which requires O(dnr 2 ) parameters and avoid the curse of dimensionality through a particular geometry of low-order contracted tensors. TT representation has been applied to model weight parameters in deep neural network and nonlinear kernel learning (Novikov et al., 2015; Stoudenmire & Schwab, 2016; Tsai et al., 2016) , achieving a significant compression factor and scalability. It also has been successfully used for feature learning and classification (Bengua et al., 2015) . To fully explore the advantages of tensor algebra, the key step is to efficiently represent the real-world dataset by tensor networks, which is not well studied. In addition, there are some limitations of TT including that i) the constraint on TT-ranks, i.e., r 1 = r d+1 = 1, leads to the limited representation ability and flexibility; ii) TT-ranks are bounded by the rank of k-unfolding matricization, which might not be optimal; iii) the permutation of data tensor will yield an inconsistent solution, i.e., TT representations and TT-ranks are sensitive to the order of tensor dimensions. Hence, finding the optimal permutation remains a challenging problem.

In this paper, we introduce a new structure of tensor networks, which can be considered as a generalization of TT representations. First of all, we relax the condition over TT-ranks, i.e., r 1 = r d+1 = 1, leading to an enhanced representation ability. Secondly, the strict ordering of multilinear products between cores should be alleviated. Third, the cores should be treated equivalently by making the model symmetric. To this end, we add a new connection between the first and the last core tensors, yielding a circular tensor products of a set of cores (see Fig. 2 ). More specifically, we consider that each tensor element is approximated by performing a trace operation over the sequential multilinear products of cores. Since the trace operation ensures a scalar output, r 1 = r d+1 = 1 is not necessary. In addition, the cores can be circularly shifted and treated equivalently due to the properties of the trace operation. We call this model tensor ring (TR) decomposition and its cores tensor ring (TR) representations. To learn TR representations, we firstly develop a non-iterative TR-SVD algorithm that is similar to TT-SVD algorithm (Oseledets, 2011) . To find the optimal lower TR-ranks, a block-wise ALS algorithms is presented. Finally, we also propose a scalable algorithm by using stochastic gradient descend, which can be applied to handling large-scale datasets.

Another interesting contribution is that we show the intrinsic structure or high order correlations within a 2D image can be captured more efficiently than SVD by converting 2D matrix to a higher order tensor. For example, given an image of size I × J, we can apply an appropriate tensorization operation (see details in Sec. 5.2) to obtain a fourth order tensor, of which each mode controls one specific scale of resolution. To demonstrate this, Fig. 1 shows the effects caused by noise corruption of specific tensor cores. As we can see, the first mode corresponds to the small-scale patches, while the 4th-mode corresponds to the large-scale partitions. We have shown in Sec. 5.2 that TR model can represent the image more efficiently than the standard SVD.

2. TENSOR RING DECOMPOSITION:

The TR decomposition aims to represent a high-order (or multi-dimensional) tensor by a sequence of 3rd-order tensors that are multiplied circularly. Specifically, let T be a dth-order tensor of size ×n d , TR representation is to decompose it into a sequence of latent tensors Z k ∈ R r k ×n k ×r k+1 , k = 1, 2, . . . , d, which can be expressed in an element-wise form given by

n 1 ×n 2 ×• • •×n d , denoted by T ∈ R

T (i 1 , i 2 , . . . , i d ) = Tr {Z 1 (i 1 )Z 2 (i 2 ) • • • Z d (i d )} = Tr d k=1 Z k (i k ) .

(1) T (i 1 , i 2 , . . . , i d ) denotes the (i 1 , i 2 , . . . , i d )th element of the tensor. Z k (i k ) denotes the i k th lateral slice matrix of the latent tensor Z k , which is of size r k × r k+1 . Note that any two adjacent latent tensors, Z k and Z k+1 , have a common dimension r k+1 on their corresponding modes. The last latent tensor Z d is of size r d × n d × r 1 , i.e., r d+1 = r 1 , which ensures the product of these matrices is a square matrix. These prerequisites play key roles in TR decomposition, resulting in some important numerical properties. For simplicity, the latent tensor Z k can also be called the kth-core (or node).

The size of cores, r k , k = 1, 2, . . . , d, collected and denoted by a vector r = [r 1 , r 2 , . . . , r d ] T , are called TR-ranks. From (1), we can observe that T (i 1 , i 2 , . . . , i d ) is equivalent to the trace of a sequential product of matrices {Z k (i k )}. Based on (1), we can also express TR decomposition in the tensor form, given by

T = r1,...,r d α1,...,α d =1 z 1 (α 1 , α 2 ) • z 2 (α 2 , α 3 ) • • • • • z d (α d , α 1 ),

where the symbol '•' denotes the outer product of vectors and z k (α k , α k+1 ) ∈ R n k denotes the (α k , α k+1 )th mode-2 fiber of tensor Z k . The number of parameters in TR representation is O(dnr 2 ), which is linear to the tensor order d as in TT representation.

The TR representation can also be illustrated graphically by a linear tensor network as shown in Fig. 2 . A node represents a tensor (including a matrix and a vector) whose order is denoted by the number of edges. The number by an edge specifies the size of each mode (or dimension). The connection between two nodes denotes a multilinear product operator between two tensors on a specific mode. This is also called tensor contraction, which corresponds to the summation over the indices of that mode. It should be noted that Z d is connected to Z 1 by the summation over the index α 1 , which is equivalent to the trace operation. For simplicity, we denote TR decomposition by

T = (Z 1 , Z 2 , . . . , Z d ).

Theorem 1 (Circular dimensional permutation invariance). Let T ∈ R n1×n2×...×n d be a dthorder tensor and its TR decomposition is given by T

= (Z 1 , Z 2 , . . . , Z d ). If we define ← - T k ∈ R n k+1 ×•••×n d ×n1×•••×n k

as the circularly shifted version along the dimensions of T by k, then we have ← -

T k = (Z k+1 , . . . , Z d , Z 1 , . . . Z k ).

A proof of Theorem 1 is provided in Appendix B.1.

It should be noted that circular dimensional permutation invariance is an essential feature that distinguishes TR decomposition from TT decomposition. For TT decomposition, the product of matrices must keep a strictly sequential order, yielding that the tensor with a circular dimension shifting does not correspond to the shifting of tensor cores.

3. LEARNING ALGORITHMS:



3.1. SEQUENTIAL SVDS:

We propose the first algorithm for computing the TR decomposition using d sequential SVDs. This algorithm will be called the TR-SVD algorithm. Theorem 2. Let us assume T can be represented by a TR decomposition. If the k-unfolding matrix T k has Rank(T k ) = R k+1 , then there exists a TR decomposition with TR-ranks r which satisfies that ∃k, r 1 r k+1 ≤ R k+1 .

Proof. We can express TR decomposition in the form of k-unfolding matrix,

T k (i 1 • • • i k , i k+1 • • • i d ) = Tr    k j=1 Z j (i j ) d j=k+1 Z j (i j )    = vec   k j=1 Z j (i j )   , vec   k+1 j=d Z T j (i j )   .

(2)

It can also be rewritten as

EQUATION

where we defined the subchain by merging multiple linked cores as

Z <k (i 1 • • • i k-1 ) = k-1 j=1 Z j (i j ) and Z >k (i k+1 • • • i d ) = d j=k+1 Z j (i j ). Hence, we can obtain T k = Z ≤k (2) (Z >k [2] ) T , where the subchain Z ≤k (2) is of size k j=1 n j × r 1 r k+1 , and Z >k [2] is of size d j=k+1 n j × r 1 r k+1 . Since the rank of T k is R k+1 , we can obtain r 1 r k+1 ≤ R k+1 .

According to (2) and (3), TR decomposition can be written as

T 1 (i 1 , i 2 • • • i d ) = α1,α2 Z ≤1 (i 1 , α 1 α 2 )Z >1 (α 1 α 2 , i 2 • • • i d ).

Since the low-rank approximation of T 1 can be obtained by the truncated SVD, which is T 1 = UΣV T + E 1 , the first core Z 1 (i.e., Z ≤1 ) of size r 1 × n 1 × r 2 can be obtained by the proper reshaping and permutation of U and the subchain Z >1 of size r 2 × d j=2 n j × r 1 is obtained by the proper reshaping and permutation of ΣV T , which corresponds to the remaining d -1 dimensions of T . Note that this algorithm use the similar strategy with TT-SVD (Oseledets, 2011), but the reshaping and permutations are totally different between them. Subsequently, we can further reshape the subchain Z >1 as a matrix Z >1 ∈ R r2n2× d j=3 nj r1 which thus can be written as

Z >1 (α 2 i 2 , i 3 • • • i d α 1 ) = α3 Z 2 (α 2 i 2 , α 3 )Z >2 (α 3 , i 3 • • • i d α 1 ).

By applying truncated SVD, i.e., Z >1 = UΣV T + E 2 , we can obtain the second core Z 2 of size (r 2 × n 2 × r 3 ) by appropriately reshaping U and the subchain Z >2 by proper reshaping of ΣV T . This procedure can be performed sequentially to obtain all d cores Z k , k = 1, . . . , d.

As proved in (Oseledets, 2011), the approximation error by using such sequential SVDs is given by

T -(Z 1 , Z 2 , . . . , Z d ) F ≤ d-1 k=1 E k 2 F .

Hence, given a prescribed relative error p , the truncation threshold δ can be set to p √ d-1 T F . However, considering that E 1 F corresponds to two ranks including both r 1 and r 2 , while E k F , ∀k > 1 correspond to only one rank r k+1 . Therefore, we modify the truncation threshold as

EQUATION

A pseudocode of the TR-SVD algorithm is summarized in Alg. 1. Note that the cores obtained by the TR-SVD algorithm are left-orthogonal, which is

Z T k 2 Z k 2 = I for k = 2, . . . , d -1.

3.2. BLOCK-WISE ALTERNATING LEAST-SQUARES (ALS):

The ALS algorithm has been widely applied to various tensor decomposition models such as CP and Tucker decompositions (Kolda & Bader, 2009; Holtz et al., 2012) . The main concept of ALS is optimizing one core while the other cores are fixed, and this procedure will be repeated until some convergence criterion is satisfied. Given a dth-order tensor T , our goal is optimize the error function as min

EQUATION

According to the TR definition in (1), we have

T (i 1 , i 2 , . . . , i d ) = α1,...,α d Z 1 (α 1 , i 1 , α 2 )Z 2 (α 2 , i 2 , α 3 ) • • • Z d (α d , i d , α 1 ) = α k ,α k+1 Z k (α k , i k , α k+1 )Z =k (α k+1 , i k+1 • • • i d i 1 • • • i k-1 , α k ) , where Z =k (i k+1 • • • i d i 1 . . . i k-1 ) = d j=k+1 Z j (i j )

k-1 j=1 Z j (i j ) denotes a slice matrix of subchain tensor by merging all cores except kth core Z k . Hence, the mode-k unfolding matrix of T can be expressed by

T [k] (i k , i k+1 • • • i d i 1 • • • i k-1 ) = α k α k+1 Z k (i k , α k α k+1 )Z =k (α k α k+1 , i k+1 • • • i d i 1 • • • i k-1 ) .

By applying different mode-k unfolding operations, we can obtain that T

[k] = Z k(2) Z =k [2] T

, where Z =k is a subchain obtained by merging d -1 cores.

The objective function in ( 5) can be optimized by solving d subproblems alternatively. More specifically, having fixed all but one core, the problem reduces to a linear least squares problem, which is min

Z k(2) T [k] -Z k(2) Z =k [2] T F , k = 1, . . . , d.

This optimization procedure must be repeated till the convergence, which is called TT-ALS.

Here, we propose a computationally efficient block-wise ALS (BALS) algorithm by utilizing truncated SVD, which facilitates the self-adaptation of ranks. The main idea is to perform the blockwise optimization followed by the separation of a block into individual cores. To achieve this, we consider merging two linked cores, e.g., Z k , Z k+1 , into a block (or subchain) Z (k,k+1) ∈ R r k ×n k n k+1 ×r k+2 . Thus, the subchain Z (k,k+1) can be optimized while leaving all cores except Z k , Z k+1 fixed. Subsequently, the subchain Z (k,k+1) can be reshaped into Z(k,k+1) ∈ R r k n k ×n k+1 r k+2 and separated into a left-orthonormal core Z k and Z k+1 by a truncated SVD:

EQUATION

where Z k 2 ∈ R r k n k ×r k+1 is the 2-unfolding matrix of core Z k , which can be set to U, while Z k+1 1 ∈ R r k+1 ×n k+1 r k+2 is the 1-unfolding matrix of core Z k+1 , which can be set to ΣV T . This procedure thus moves on to optimize the next block cores Z (k+1,k+2) , . . . , Z (d-1,d) , Z (d,1) successively in the similar way. Note that since the TR model is circular, the dth core can also be merged with the first core yielding the block core Z (d,1) .

The key advantage of our BALS algorithm is the rank adaptation ability which can be achieved simply by separating the block core into two cores via truncated SVD, as shown in ( 6). The truncated rank r k+1 can be chosen such that the approximation error is below a certain threshold. One possible choice is to use the same threshold as in the TR-SVD algorithm, i.e., δ k described in (4). However, the empirical experience shows that this threshold often leads to overfitting and the truncated rank is higher than the optimal rank. This is because the updated block Z (k,k+1) during ALS iterations is not a closed form solution and many iterations are necessary for convergence. To relieve this problem, we choose the truncation threshold based on both the current and the desired approximation errors, which is

δ = max T F / √ d, p T F / √ d .

A pseudo code of the BALS algorithm is described in Alg. 2.

3.3. STOCHASTIC GRADIENT DESCENT:

For large-scale dataset, the ALS algorithm is not scalable due to the cubic time complexity in the target rank, while Stochastic Gradient Descent (SGD) shows high efficiency and scalability for matrix/tensor factorization (Gemulla et al., 2011; Maehara et al., 2016; Wang & Anandkumar, 2016) . In this section, we present a scalable and efficient TR decomposition by using SGD, which is also suitable for online learning and tensor completion problems. To this end, we first provide the element-wise loss function, which is

L(Z 1 , Z 2 , . . . , Z d ) = 1 2 i1,...,i d T (i 1 , i 2 , . . . , i d ) -Tr d k=1 Z k (i k ) 2 + 1 2 λ k Z k (i k ) 2 , (7)

where λ k is the regularization parameters. The core idea of SGD is to randomly select one sample T (i 1 , i 2 , . . . , i d ), then update the corresponding slice matrices Z k (i k ), k = 1, . . . , d from each latent core tensor Z k based on the noisy gradient estimates by scaling up just one of local gradients, i.e. ∀k = 1, . . . , d,

EQUATION

We employ Adaptive Moment Estimation (Adam) method to compute adaptive learning rates for each parameter. Thus, the update rule for each core tensor is given by

EQUATION

where

M t = β 1 M t-1 + (1 -β 1 ) ∂L ∂Z t k (i k )

denotes an exponentially decaying average of past gradients and

V t = β 2 V t-1 + (1 -β 2 )( ∂L ∂Z t k (i k )

) 2 denotes exponentially decaying average of second moment of the gradients.

The SGD algorithm can be naturally applied to tensor completion problem, when the data points are sampled only from a sparse tensor. Furthermore, this also naturally gives an online TR decomposition. The batched versions, in which multiple local losses are averaged, are also feasible but often have inferior performance in practice. For each element T (i 1 , i 2 , . . . , i d ), the computational complexity of SGD is O(d 2 r 3 ). If we define N = d k=1 n k consecutive updates as one epoch of SGD, the computational complexity per SGD epoch is thus O(N d 2 r 3 ), which linearly scales to data size. As compared to ALS, which needs O(N dr 4 + dr 6 ), it is more efficient in terms of computational complexity for one epoch. The memory cost of TR-SVD, ALS, and SGD are O(ndr 2 + k n k ), O(ndr 2 +r 2 k n k +r 4 ), and O(ndr 2 +r 2 ), respectively. Thus, TR-SGD requires much less memory. The convergence condition of SGD algorithm follows other stochastic tensor decompositions (Ge et al., 2015; Maehara et al., 2016) .

In practice, how to choose the TR algorithms depends on the task and data at hand. TR-SVD is a non-iterative algorithm leading to the fast computation. The expected approximation error is required from user to choose TR-ranks automatically. However, the learned ranks is not necessary to be the optimal in terms of compression rate. TR-ALS is an iterative algorithm in which the fixed TR-ranks are required to be given by user. The computation time is highly related to TR-ranks, which is very efficient only for small TR-ranks. TR-BALS is an iterative algorithm that can choose TR-ranks in an adaptive manner given the expected approximation error. TR-SGD is an iterative algorithm that is suitable for a large-scale tensor with given TR-ranks. In particular, the memory cost of TR-SGD depends on only the size of core tensors rather than the size of data tensor.

4. PROPERTIES OF TR REPRESENTATION:

By assuming that tensor data have been already represented as TR decompositions, i.e., a sequence of third-order cores, we justify and demonstrate that the basic operations on tensors, such as the addition, multilinear product, Hadamard product, inner product and Frobenius norm, can be performed efficiently by the appropriate operations on each individual cores. We have the following theorems: Property 1. Let T 1 and T 2 be dth-order tensors of size

n 1 × • • • × n d . If TR decompositions of these two tensors are T 1 = (Z 1 , . . . , Z d ) where Z k ∈ R r k ×n k ×r k+1 and T 2 = (Y 1 , . . . , Y d ) where Y k ∈ R s k ×n k ×s k+1

, then the addition of these two tensors, T 3 = T 1 + T 2 , can also be represented in the TR format given by T 3 = (X 1 , . . . , X d ), where X k ∈ R q k ×n k ×q k+1 and q k = r k + s k . Each core X k can be computed by

EQUATION

A proof of Property 1 is provided in Appendix B.2. Note that the sizes of new cores are increased and not optimal in general. This problem can be solved by the rounding procedure (Oseledets, 2011). 

= T × 1 u T 1 × 2 • • • × d u T d

, can be computed by the multilinear product on each cores, which is

EQUATION

A proof of Property 2 is provided in Appendix B.3. It should be noted that the computational complexity in the original tensor form is O(dn d ), while it reduces to O(dnr 2 + dr 3 ) that is linear to tensor order d by using TR representation. Property 3. Let T 1 and T 2 be dth-order tensors of size

n 1 × • • • × n d . If the TR decompositions of these two tensors are T 1 = (Z 1 , . . . , Z d ) where Z k ∈ R r k ×n k ×r k+1 and T 2 = (Y 1 , . . . , Y d )

where Y k ∈ R s k ×n k ×s k+1 , then the Hadamard product of these two tensors, T 3 = T 1 T 2 , can also be represented in the TR format given by T 3 = (X 1 , . . . , X d ), where X k ∈ R q k ×n k ×q k+1 and q k = r k s k . Each core X k can be computed by

EQUATION

The inner product of two tensors can be computed by TR format

T 1 , T 2 = (V 1 , . . . , V d ) where V k = n k i k =1 Z k (i k ) ⊗ Y k (i k ).

Therefore, the Frobenius norm can be also computed by TR format using T F = T , T .

A proof of Property 3 is provided in Appendix B.4. In contrast to O(n d ) in the original tensor form, the computational complexity is equal to O(dnq 2 + dq 3 ) that is linear to d by using TR representation.

5. EXPERIMENTAL RESULTS:



5.1. NUMERICAL ILLUSTRATION:

We consider highly oscillating functions that can be approximated perfectly by a low-rank TT format (Khoromskij, 2015) , as shown in Fig. 3 . We firstly tensorize the functional vector resulting in a dth-order tensor of size

n 1 × n 2 × • • • × n d

, where isometric size is usually preferred, i.e.,

n 1 = n 2 = • • • = n d = n,

with the total number of elements denoted by N = n d . The error bound (tolerance), denoted by p = 10 -3 , is given as the stopping criterion for all compared algorithms. As shown in Table 1 , TR-SVD and TR-BALS can obtain comparable results with TT-SVD in terms of compression ability. However, when noise is involved, TR model significantly outperforms TT model, indicating its more robustness to noises. It should be noted that TT representation has the property that r 1 = r d+1 = 1 and r k , k = 2, . . . , d-1 are bounded by the rank of k-unfolding matrix of T k , which limits its generalization ability and consistency when the tensor modes have been shifted or permuted. To demonstrate this, we consider shifting the dimensions of T of size

n 1 × • • • × n d by k times leading to ← - T k of size n k+1 × • • • × n d × n 1 × • • • × n k .

As shown in Table 2 , the average TT-ranks are varied dramatically along with the different shifts. In particular, when k = 8, rtt becomes 14.6, resulting in a large number of parameters N p = 10376. In contrast to TT, TR can obtain consistent and compact representation with N p = 1360.

We also tested TR-ALS and TR-SGD algorithms on datasets which are generated by a TR model, in which the core tensors are randomly drawn from N (0, 1). As shown in Table 3 , TR-SGD can achieve similar performance as TR-ALS in all cases. In particular, when data is relatively large-scale (10 8 ), TR-SGD can achieve relative error = 0.01 by using 1% of data points only once. 

= r 2 = • • • = 2.

5.2. IMAGE REPRESENTATION BY HIGHER-ORDER TENSOR DECOMPOSITIONS:

An image is naturally represented by a 2D matrix, on which SVD can provide the best low-rank approximation. However, the intrinsic structure and high-order correlations within the image is not well exploited by SVD. In this section, we show the tensorization of an image, yielding a higher-order tensor, and TR decomposition enable us to represent the image more efficiently than SVD. Given an image (e.g. 'Peppers') denoted by Y of size I × J, we can reshape it as I 1 × I 2 × . . . × I d × J 1 × J 2 × . . . × J d followed by an appropriate permutation to I 1 × J 1 × I 2 × J 2 . . . × I d × J d and thus reshape it again to I 1 J 1 × I 2 J 2 × . . . × I d J d , which is a dth-order tensor. The first mode corresponds to small-scale patches of size I 1 × J 1 , while the dth-mode corresponds to large-scale partition of whole image as I d × J d . Based on this tensorization operations, TR decomposition is able to capture the intrinsic structure information and provides a more compact representation. As shown in Table 4 , for 2D matrix case, SVD, TT and TR give exactly same results. In contrast, for 4th-order tensorization cases, TT needs only half number of parameters (2 times compression rate) while TR achieves 3 times compression rate, given the same approximation error 0.1. It should be noted that TR representation provides significantly high compression ability as compared to TT. In addition, Fig. 4 shows TR-SGD results on 'Lena' image by sampling different fraction of data points.

Table 4 : Image representation by using tensorization and TR decomposition. The number of parameters is compared for SVD, TT and TR given the same approximation errors. 4.8e3 4.3e3 7.8e4 7.8e4 1.1e5 9.8e4 1.3e5 1.0e5 n = 2, d = 16 7.4e3 7.4e3 1.0e5 1.0e5 1.5e5 1.5e5 1.7e5 1.7e5

5.3. CIFAR-10:

The CIFAR-10 dataset consists of 60000 32 × 32 colour images. We randomly pick up 1000 images for testing of TR decomposition algorithms. As shown in Fig. 5 , TR model outperforms TT model in terms of compression rate given the same approximation error, which is caused by strict limitation that the mode-1 rank must be 1 for TT model. In addition, TR is a more generalized model, which contains TT as a special case, thus yielding better low-rank approximation. Moreover, all other TR algorithms can also achieve similar results. The detailed results for = 1e -1 are shown in Table 5 . Note that TR-SGD can achieve the same performance as TR-ALS, which demonstrates its effectiveness on real-world dataset. Due to high computational efficiency of TR-SGD per epoch, it can be potentially applied to very large-scale dataset. For visualization, TR-SGD results after 10 and 100 epoch are shown in Fig. 6 .

5.4. TENSORIZING NEURAL NETWORKS USING TR REPRESENTATION:

TT representations have been successfully applied to deep neural networks (Novikov et al., 2015) , which can significantly reduce the number of model parameters and improve computational efficiency. To investigate the properties of TR representation, we applied TR framework to approximate the weight matrix of a fully-connected layer and compared with TT representation. We run the experiment on the MNIST dataset for the task of handwritten-digit recognition. The same setting of neural network (two fully-connected layers with ReLU activation function) as in (Novikov et al., 2015) was applied for comparisons. The input layer is tensorized to a 4th-order tensor of size 4 × 8 × 8 × 4, the weight matrix of size 1024 × 625 is represented by a TR format of size 4 × 8 × 8 × 4 × 5 × 5 × 5 × 5. Through deriving the gradients over each core tensor, all computations can be performed on small core tensors instead of the dense weight matrix by using properties in Sec. 4, yielding the significant improvements of computational efficiency. The experimental results are shown in Fig. 7 . We observe that the TR-layer provides much better flexibility than TT-layer, leading to much lower training and testing errors under the same compression level (i.e., TT/TR ranks). In addition, TR can achieve much better compression rate under the same level of test error. When r 1 = . . . = r 4 = 2, the compression rate of dense weight matrix is up to 1300 times.

We tested the tensorizing neural networks with the same architecture on SVHN dataset (http://ufldl.stanford.edu/housenumbers/). By setting all the TT-ranks in the network to 4, we achieved the test error of 0.13 with compression rate of 444 times, while we can achieve the same test error by setting all the TR-ranks to 3 with compression rate of 592 times. We can conclude that the TR representation can obtain significantly higher compression rate under the same level of test error.

6. CONCLUSION:

We have proposed a novel tensor decomposition model, which provides an efficient representation for a very high-order tensor by a sequence of low-dimensional cores. The number of parameters in our model scales only linearly to the tensor order. To optimize the latent cores, we have presented several different algorithms: TR-SVD is a non-recursive algorithm that is stable and efficient, while TR-BALS can learn a more compact representation with adaptive TR-ranks, TR-SGD is a scalable algorithm which can be also used for tensor completion and online learning. Furthermore, we have investigated the properties on how the basic multilinear algebra can be performed efficiently by operations over TR representations (i.e., cores), which provides a powerful framework for processing large-scale data. The experimental results verified the effectiveness of our proposed algorithms.

A RELATION TO OTHER MODELS:

In this section, we discuss the relations between TR model and the classical tensor decompositions including CPD, Tucker and TT models. All these tensor decompositions can be viewed as the transformed representation of a given tensor. The number of parameters in CPD is O(dnr) that is linear to tensor order, however, its optimization problem is difficult and convergence is slow. The Tucker model is stable and can approximate an arbitrary tensor as close as possible, however, its number of parameters is O(dnr + r d ) that is exponential to tensor order. In contrast, TT and TR decompositions have similar representation power to Tucker model, while their number of paramters is O(dnr 2 ) that is linear to tensor order.

A.1 CP DECOMPOSITION:

The cannonical polyadic decomposition (CPD) aims to represent a dth-order tensor T by a sum of rank-one tensors, given by

EQUATION

where each rank-one tensor is represented by an outer product of d vectors. It can be also written in the element-wise form given by

EQUATION

where •, . . . , • denotes an inner product of a set of vectors, i.e., u

EQUATION

i k ) which is a diagonal matrix for each fixed i k and k, where k = 1, . . . , d, i k = 1, . . . , n k , we can rewrite (14) as

EQUATION

Hence, CPD can be viewed as a special case of TR decomposition T = (V 1 , . . . , V d ) where the cores V k , k = 1, . . . , d are of size r × n k × r and each lateral slice matrix V k (i k ) is a diagonal matrix of size r × r.

A.2 TUCKER DECOMPOSITION:

The Tucker decomposition aims to represent a dth-order tensor T by a multilinear product between a core tensor G ∈ R r1×•••×r d and factor matrices

EQUATION

By assuming the core tensor G can be represented by a TR decomposition G = (V 1 , . . . , V d ), the Tucker decomposition ( 16) in the element-wise form can be rewritten as

EQUATION

where the second step is derived by applying Theorem 2. Hence, Tucker model can be represented as a TR decomposition T = (Z 1 , . . . , Z d ) where the cores are computed by the multilinear products between TR cores representing G and the factor matrices, respectively, which is

EQUATION

A.3 TT DECOMPOSITION

The tensor train decomposition aims to represent a dth-order tensor T by a sequence of cores G k , k = 1, . . . , d, where the first core G 1 ∈ R n1×r2 and the last core G d ∈ R r d ×n d are matrices while the other cores G k ∈ R r k ×n k ×r k+1 , k = 2, . . . , d -1 are 3rd-order tensors. Specifically, TT decomposition in the element-wise form is expressed as

EQUATION

) where g 1 (i 1 ) is the i 1 th row vector of G 1 , g d (i d ) is the i d th column vector of G d , and G k (i k ), k = 2, . . . , d -1 are the i k th lateral slice matrices of G k .

According to the definition of TR decomposition in (1), it is obvious that TT decomposition is a special case of TR decomposition where the first and the last cores are matrices, i.e., r 1 = r d+1 = 1. On the other hand, TR decomposition can be also rewritten as

EQUATION

where z 1 (α 1 , i 1 , :) ∈ R r2 is the α 1 th row vector of the matrix Z 1 (i 1 ) and z d (:, i d , α 1 ) is the α 1 th column vector of the matrix Z d (i d ). Therefore, TR decomposition can be interpreted as a sum of TT representations. The number of TT representations is r 1 and these TT representations have the common cores Z k , for k = 2, . . . , d -1. In general, TR outperforms TT in terms of representation power due to the fact of linear combinations of a group of TT representations. Furthermore, given a specific approximation level, TR representation requires smaller ranks than TT representation.

B PROOFS B.1 PROOF OF THEOREM 1:

Proof. It is obvious that (1) can be rewritten as

T (i 1 , i 2 , . . . , i d ) = Tr(Z 2 (i 2 ), Z 3 (i 3 ), . . . , Z d (i d ), Z 1 (i 1 )) = • • • = Tr(Z d (i d ), Z 1 (i 1 ), . . . , Z d-1 (i d-1 )).

Therefore, we have ← -T k = (Z k+1 , . . . , Z d , Z 1 , . . . , Z k ).

B.2 PROOF OF PROPERTY 1:

Proof. According to the definition of TR decomposition, and the cores shown in (10), the (i 1 , . . . , i d )th element of tensor T 3 can be written as

T 3 (i 1 , . . . , i d ) =Tr d k=1 Z k (i k ) 0 0 d k=1 Y k (i k ) = Tr d k=1 Z k (i k ) + Tr d k=1 Y k (i k ) .

Hence, the addition of tensors in the TR format can be performed by merging of their cores.

B.3 PROOF OF PROPERTY 2:

Proof. The multilinear product between a tensor and vectors can be expressed by

c = i1,...,i d T (i 1 , . . . , i d )u 1 (i 1 ) • • • u d (i d ) = i1,...,i d Tr d k=1 Z k (i k ) u 1 (i 1 ) • • • u d (i d ) =Tr d k=1 n k i k =1 Z k (i k )u k (i k )

.

Thus, it can be written as a TR decomposition shown in (11) where each core X k ∈ R r k ×r k+1 becomes a matrix. The computational complexity is equal to O(dnr 2 ).

Algorithm 1 TR-SVD:

Input: A dth-order tensor T of size (n 1 × • • • × n d ) and the prescribed relative error p .

Output: Cores Z k , k = 1, . . . , d of TR decomposition and the TR-ranks r.

1: Compute truncation threshold δ k for k = 1 and k > 1.

2: Choose one mode as the start point (e.g., the first mode) and obtain the 1-unfolding matrix T 1 .

3: Low-rank approximation by applying δ 1 -truncated SVD:

T 1 = UΣV T + E 1 . 4: Split ranks r 1 , r 2 by min r1,r2 r 1 -r 2 , s. t. r 1 r 2 = rank δ1 (T 1 ). 5: Z 1 ← permute(reshape(U, [n 1 , r 1 , r 2 ]), [2, 1, 3]). 6: Z >1 ← permute(reshape(ΣV T , [r 1 , r 2 , d j=2 n j ]), [2, 3, 1]). 7: for k = 2 to d -1 do 8: Z >k-1 = reshape(Z >k-1 , [r k n k , n k+1 • • • n d r 1 ]). 9: Compute δ k -truncated SVD: Z >k-1 = UΣV T + E k . 10: r k+1 ← rank δ k (Z >k-1 ). 11: Z k ← reshape(U, [r k , n k , r k+1 ]).

12:

Z >k ← reshape(ΣV T , [r k+1 , d j=k+1 n j , r 1 ]). 13: end for B.4 PROOF OF PROPERTY 3

Proof. Each element in tensor T 3 can be written as

T 3 (i 1 , . . . , i d ) =Tr d k=1 Z k (i k ) Tr d k=1 Y k (i k ) = Tr d k=1 Z k (i k ) ⊗ d k=1 Y k (i k ) =Tr d k=1 Z k (i k ) ⊗ Y k (i k )

.

Hence, T 3 can be also represented as TR format with its cores computed by ( 12), which costs O(dnq 2 ).

Furthermore, one can compute the inner product of two tensors in TR representations. For two tensors T 1 and T 2 , it is defined as T 1 , T 2 = i1,...,i d T 3 (i 1 , . . . , i d ), where T 3 = T 1 T 2 . Thus, the inner product can be computed by applying the Hadamard product and then computing the multilinear product between T 3 and vectors of all ones, i.e., T 3 ×

1 u T 1 × 2 • • • × d u T d

, where u k = 1, k = 1, . . . , d, which can be computed efficiently by using Property 2.

C ALGORITHMS:

Pseudo codes of TR-SVD and TR-BALS are provided in Alg. 1 and Alg. 2, respectively.

D ADDITIONAL EXPERIMENTAL RESULTS:



D.1 COIL-100 DATASET:

The Columbia Object Image Libraries (COIL)-100 dataset (Nayar et al., 1996) contains 7200 color images of 100 objects (72 images per object) with different reflectance and complex geometric characteristics. Each image can be represented by a 3rd-order tensor of size 128 × 128 × 3 and then is downsampled to 32 × 32 × 3. Hence, the dataset can be finally organized as a 4th-order tensor of size 32 × 32 × 3 × 7200. The number of features is determined by r 4 × r 1 , while the flexibility of subspace bases is determined by r 2 , r 3 . Subsequently, we apply the K-nearest neighbor (KNN) classifier with K=1 for classification. For detailed comparisons, we randomly select a certain ratio ρ = 50% or ρ = 10% samples as the training set and the rest as the test set. The classification performance is averaged over 10 times of random splitting. In Table 6 , r max of TR decompositions is much smaller than that of TT-SVD. It should be noted that TR representation, as compared to TT, can obtain more compact and discriminant representations. Fig. 8 shows the reconstructed images under different approximation levels. Compute the subchain Z =(k,k+1) .

5:

Obtain the mode-2 unfolding matrix Z

=(k,k+1) [2] of size d j=1 n j /(n k n k+1 ) × r k r k+2 . 6: Z (k,k+1) (2) ← arg min T [k] -Z (k,k+1) (2) Z =(k,k+1) [2] T F . 7:

Tensorization of mode-2 unfolding matrix

EQUATION

).

8::

Reshape the block core by Z(k,k+1) ← reshape(Z (k,k+1) , [r k n k × n k+1 r k+2 ]).

9::

Low-rank approximation by δ-truncated SVD Z(k,k+1) = UΣV T .

10:

Z k ← reshape(U, [r k , n k , r k+1 ]).

11:

Z k+1 ← reshape(ΣV T , [r k+1 , n k+1 , r k+2 ]).

12:

r k+1 ← rank δ ( Z(k,k+1) ).

13:

k ← k + 1. 14: until The desired approximation accuracy is achieved, i.e., ≤ p . We test the TR representation for KTH video database (Laptev & Lindeberg, 2006) containing six types of human actions (walking, jogging, running, boxing, hand waving and hand clapping) performed several times by 25 subjects in four different scenarios: outdoors, outdoors with scale variation, outdoors with different clothes and indoors as illustrated in Fig. 9 . There are 600 video sequences for each combination of 25 subjects, 6 actions and 4 scenarios. Each video sequence was downsampled to 20×20×32. Finally, we can organize the dataset as a tensor of size 20×20×32×600. For extensive comparisons, we choose different error bound p ∈ {0.2, 0.3, 0.4}. In Table 7 , we can see that TR representations achieve better compression ratio reflected by smaller r max , r than that of TT-SVD, while TT-SVD achieves better compression ratio than CP-ALS. For instance, when ≈ 0.2, CP-ALS requires r max = 300, r = 300; TT-SVD requires r max = 139, r = 78, while TR-SVD only requires r max = 99, r = 34.2. For classification performance, we observe that the best accuracy (5 × 5-fold cross validation) achieved by CP-ALS, TT-SVD, TR-SVD are 80.8%, 84.8%, 87.7%, respectively. Note that these classification performances might not be the state-of-the-art on this dataset, we mainly focus on the comparisons of representation ability among CP, TT, and TR decomposition frameworks. To obtain the best performance, we may apply the powerful feature extraction methods to TT or TR representations of dataset. It should be noted that TR decompositions achieve the best classification accuracy when = 0.29, while TT-SVD and CP-ALS achieve their best classification accuracy when = 0.2. This indicates that TR decomposition can preserve more discriminant information even when the approximation error is relatively high. This experiment demonstrates that TR decompositions are effective for unsupervised feature representation due to their flexibility of TR-ranks and high compression ability. 

:



Figures and Tables:

Figure 1: The effects of noise corrupted tensor cores. From left to right, each figure shows noise corruption by adding noise to one specific tensor core.

Figure 3: Highly oscillated functions. The left panel isf 1 (x) = (x + 1) sin(100(x + 1) 2). The middle panel is Airy function:f 2 (x) = x -1 4 sin( 2 3 x32 ). The right panel is Chirp function f 3 (x) = sin x 4 cos(x 2 ).

Figure 4: TR-SGD decomposition with TR-ranks of 12 on the 8th-order tensorization of an image. Iter: 50% indicates that only 50% elements are sampled for learning its TR representation. RSE indicates root relative square error Ŷ -Y F / Y F .

n = 10, d = 4 ( = 0.01, Epoch = 19) ( = 0.01, Epoch = 10 ) n = 10, d = 6 ( = 0.01, Epoch = 10) ( = 0.01, Epoch = 0.4 ) n = 10, d = 8 ( = 0.05, Epoch = 9) ( = 0.01, Epoch = 0.01 )

Figure5: The comparisons of compression rate and approximation error (RSE) on CIFAR-10 dataset by using TR and TT models.

, d = 4 5.1e3 3.8e3 6.8e4 6.4e4 1.0e5 7.3e4 1.3e5 7.4e4 n = 4, d = 8

Figure 6: The reconstructed images by using TR-SGD after 10 and 100 epochs.

Figure 7: The classification performances of tensorizing neural networks by using TR representation.

TR-BALSInput: A d-dimensional tensor T of size (n 1 × • • • × n d )and the prescribed relative error p . Output: Cores Z k and TR-ranks r k , k = 1, . . . , d.1: Initialize r k = 1 for k = 1, . . . , d.2: Initialize Z k ∈ R r k ×n k ×r k+1 for k = 1, . . . , d.3: repeat k ∈ circular{1, 2, . . . , d};4:

Figure 8: The reconstruction of Coil-100 dataset by using TR-SVD. The top row shows the original images, while the reconstructed images are shown from the second to sixth rows corresponding to =0.1, 0.2, 0.3, 0.4, 0.5, respectively.

Figure 9: Video dataset consists of six types of human actions performed by 25 subjects in four different scenarios. From the top to bottom, six video examples corresponding to each type of actions are shown.

n1×

Table content: <table/>

The functional data f 1 (x), f 2 (x), f 3 (x) is tensorized to 10th-order tensor (4×4×. . .×4). In the table, , r, N p denote relative error, average rank, and the total number of parameters, respectively.

Table content: <table><tr><td/><td>r</td><td>f1(x) Np</td><td>Time (s)</td><td/><td>r</td><td>f2(x) Np</td><td>Time (s)</td><td>r</td><td colspan="2">f3(x) Np Time (s)</td><td colspan="2">f1(x) + N (0, σ), SN R = 60dB r Np Time (s)</td></tr><tr><td>TT-SVD</td><td colspan="2">3e-4 4.4 1032</td><td>0.17</td><td>3e-4</td><td>5</td><td>1360</td><td>0.16</td><td colspan="2">3e-4 3.7 680</td><td>0.16</td><td colspan="2">1e-3 16.6 13064</td><td>0.5</td></tr><tr><td colspan="3">TR-SVD 3e-4 4.4 1032</td><td>0.17</td><td>3e-4</td><td>5</td><td>1360</td><td>0.28</td><td colspan="2">5e-4 3.6 668</td><td>0.15</td><td>1e-3 9.7</td><td>4644</td><td>0.4</td></tr><tr><td>TR-ALS</td><td colspan="2">3e-4 4.4 1032</td><td>13.2</td><td>3e-4</td><td>5</td><td>1360</td><td>18.6</td><td colspan="2">8e-4 3.6 668</td><td>4.0</td><td>1e-3 4.4</td><td>1032</td><td>11.8</td></tr><tr><td colspan="3">TR-BALS 9e-4 4.3 1052</td><td>4.6</td><td colspan="3">8e-4 4.9 1324</td><td>5.7</td><td colspan="2">5e-4 3.7 728</td><td>3.4</td><td>1e-3 4.2</td><td>1000</td><td>6.1</td></tr></table>

The results under different shifts of dimensions on functional data f 2 (x) with error bound set at 10 -3 . For the 10th-order tensor, all 9 dimension shifts were considered and the average rank r is compared.

Table content: <table><tr><td/><td/><td/><td/><td/><td>r</td><td/><td/><td/></tr><tr><td/><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td></tr><tr><td colspan="5">TT 5.2 5.8 6 6.2</td><td>7</td><td colspan="4">7 8.5 14.6 8.4</td></tr><tr><td>TR</td><td>5</td><td colspan="5">4.9 5 4.9 4.9 5</td><td>5</td><td colspan="2">4.8 4.9</td></tr><tr><td colspan="2">Iter: 50%,</td><td/><td>Iter: 250%,</td><td/><td/><td>Iter: 500%,</td><td/><td colspan="2">Iter: 5000%,</td></tr><tr><td colspan="2">RSE=0.2058</td><td/><td colspan="2">RSE=0.1430</td><td/><td>RSE=0.1270</td><td/><td colspan="2">RSE=0.1116</td></tr></table>

Results on synthetic data with fixed ranks r 1

Table content: <table/>

The comparisons of different algorithms on Coil-100 dataset. , r max ,r denote relative error, the maximum rank and the average rank, respectively.

Table content: <table><tr><td/><td/><td>r max</td><td>r</td><td>Acc. (%) (ρ = 50%)</td><td>Acc. (%) (ρ = 10%)</td></tr><tr><td/><td>0.19</td><td>67</td><td>47.3</td><td>99.05</td><td>89.11</td></tr><tr><td/><td>0.28</td><td>23</td><td>16.3</td><td>98.99</td><td>88.45</td></tr><tr><td>TT-SVD</td><td>0.37</td><td>8</td><td>6.3</td><td>96.29</td><td>86.02</td></tr><tr><td/><td>0.46</td><td>3</td><td>2.7</td><td>47.78</td><td>44.00</td></tr><tr><td/><td>0.19</td><td>23</td><td>12.0</td><td>99.14</td><td>89.29</td></tr><tr><td/><td>0.28</td><td>10</td><td>6.0</td><td>99.19</td><td>89.89</td></tr><tr><td>TR-SVD</td><td>0.36</td><td>5</td><td>3.5</td><td>98.51</td><td>88.10</td></tr><tr><td/><td>0.43</td><td>3</td><td>2.3</td><td>83.43</td><td>73.20</td></tr></table>

The comparisons of different algorithms on KTH dataset. denotes the obtained relative error; r max denotes maximum rank; r denotes the average rank; and Acc. is the classification accuracy.

Table content: <table><tr><td/><td colspan="2">r max 0.20 300</td><td>r 300</td><td>Acc. (5 × 5-fold) 80.8 %</td></tr><tr><td>CP-ALS</td><td>0.30</td><td>40</td><td>40</td><td>79.3 %</td></tr><tr><td/><td>0.40</td><td>10</td><td>10</td><td>66.8 %</td></tr><tr><td/><td colspan="3">0.20 139 78.0</td><td>84.8 %</td></tr><tr><td>TT-SVD</td><td>0.29 0.38</td><td>38 14</td><td>27.3 9.3</td><td>83.5 % 67.8 %</td></tr><tr><td/><td>0.20</td><td>99</td><td>34.2</td><td>78.8 %</td></tr><tr><td>TR-SVD</td><td>0.29</td><td>27</td><td>12.0</td><td>87.7 %</td></tr><tr><td/><td>0.37</td><td>10</td><td>5.8</td><td>72.4 %</td></tr></table>

Peer review (iclr2018):
 The paper addresses the problem of tensor decomposition which is relevant and interesting. The paper proposes Tensor Ring (TR) decomposition which improves over and bases on the Tensor Train (TT) decomposition method. TT decomposes a tensor in to a sequences of latent tensors where the first and last tensors are a 2D matrices. The proposed TR method generalizes TT in that the first and last tensors are also 3rd-order tensors instead of 2nd-order. I think such generalization is interesting but the innovation seems to be very limited. The paper develops three different kinds of solvers for TR decomposition, i.e., SVD, ALS and SGD. All of these are well known methods. Finally, the paper provides experimental results on synthetic data (3 oscillated functions) and image data (few sampled images). I think the paper could be greatly improved by providing more experiments and ablations to validate the benefits of the proposed methods. Please refer to below for more comments and questions. -- The rating has been updated. Pros: 1. The topic is interesting. 2. The generalization over TT makes sense. Cons: 1. The writing of the paper could be improved and more clear: the conclusions on inner product and F-norm can be integrated into "Theorem 5". And those "theorems" in section 4 are just some properties from previous definitions; they are not theorems. 2. The property of TR decomposition is that the tensors can be shifted (circular invariance). This is an interesting property and it seems to be the major strength of TR over TT. I think the paper could be significantly improved by providing more applications of this property in both theory and experiments. 3. As the number of latent tensors increase, the ALS method becomes much worse approximation of the original optimization. Any insights or results on the optimization performance vs. the number of latent tensors? 4. Also, the paper mentions Eq. 5 (ALS) is optimized by solving d subproblems alternatively. I think this only contains a single round of optimization. Should ALS be applied repeated (each round solves d problems) until convergence? 5. What is the memory consumption for different solvers? 6. SGD also needs to update at least d times for all d latent tensors. Why is the complexity O(r^3) independent of the parameter d? 7. The ALS is so slow (if looking at the results in section 5.1), which becomes not practical. The experimental part could be improved by providing more results and description about a guidance on how to choose from different solvers. 8. What does "iteration" mean in experimental results such as table 2? Different algorithms have different cost for "each iteration" so comparing that seems not fair. The results could make more sense by providing total time consumptions and time cost per iteration. also applies to table 4. 9. Why is the .epsion in table 3 not consistent? Why not choose .epsion = 9e-4 and .epsilon=2e-15 for tensorization? 10. Also, table 3 could be greatly improved by providing more ablations such as results for (n=16, d=8), (n=4, d=4), etc. That could help readers to better understand the effect of TR. 11. Section 5.3 could be improved by providing a curve (compression vs. error) instead of just providing a table of sampled operating points. 12. The paper mentions the application of image representation but only experiment on 32x32 images. How does the proposed method handle large images? Otherwise, it does not seem to be a practical application. 13. Figure 5: Are the RSE measures computed over the whole CIFAR-10 dataset or the displayed images? Minor: - Typo: Page 4 Line 7 "Note that this algorithm use the similar strategy": use -> uses