Tensor decomposition is an extension of matrix decomposition in higher dimensions. It is an essential tool in social relation prediction [1], financial encryption [2], and image processing [3]. Where tensor-train (TT) and tensor-ring (TR) decomposition have been widely used in signal processing [4], computer vision [5], and data mining [6]. According to the "Curse of Dimensionality," tensor decomposition's time and space cost increase quickly with the size and dimension of the tensor. It is a critical mission to develop high-performance tensor decompositions. Currently, CPU-based TT and TR decompositions [9,10] do not take full advantage of the algorithms' parallelism, making it difficult to process large amounts of data.
This paper utilizes GPUs to achieve high-performance TT and TR decomposition algorithms. Moreover, we conducted experiments comparing existing CPU algorithms with our optimized GPU algorithms, showing that our optimized TT and TR decomposition implemen-tations on GPUs are efficient.
There are three major contributions to this paper:
• This paper proposes high-performance GPU-based third-order TT and TR decom-positions with the same accuracy as CPUs.
This paper proposes efficient memory access of tensors in GPUs, an efficient diago-nal matrix and matrix multiplication, and utilizes the high-parallel Jacobi-based SVD for replacing the traditional SVD. The optimized algorithms reduce memory footprint and tensor matricization.
This paper conducts experiments to verify the performance of TT and TR decompositions on one Tesla V100 GPU. The optimized TT and TR decomposition implementations get maximum of 6.67 and 6.36 speedups over the GPU basic implementations.
This section describes notations, TT decomposition, and TR decomposition.
This paper utilizes boldface lowercase letters 𝐚 ∈ ℝ n , boldface uppercase letters 𝑨 ∈ ℝ n 1 ×n 2 , and uppercase calligraphic letters 𝒜 ∈ ℝ n 1 ×n 2 ×n 3 to denote vectors, matrices, and tensors, respectively. Tensor contractions is represented with the ∘ symbol. Figure 1 shows that the TT decomposition [9] uses three third-order core tensors to express a thirdorder tensor 𝒜 ∈ ℝ 𝑛 1 ×𝑛 2 ×𝑛 3 by tensor contractions: 𝒜 = 𝒢 (1) ∘ 𝒢 (2) ∘ 𝒢 (3) ,
(1) where 𝒢 (𝑘) ∈ R r k−1 ×n k ×r k expresses the k-th core tensor. [r 0 , r 1 , r 2 , r 3 ] expresses the TT-ranks where r 0 = r 3 = 1. Therefore, 𝒢 (1) and 𝒢 (3) are second-order tensors (matrices).
The tensor-train structure is one of the tensor networks and is represented in Figure 1 by the graphical modeling [7]. The connections between two tensors indicate tensor contractions. The original tensor 𝒜 is obtained by the tensor contractions of all tensors on the left. The steps of third-order TT decomposition [9] are described in Algorithm 1. In the third and tenth lines of Algorithm 1, we convert a tensor 𝒞 (𝑘−1) into a matrix 𝐂 with r k−1 n k rows and ∏ n i 3 i=k+1
columns and a matrix 𝑼 into a tensor 𝒢 (k) with n k columns, r k in the third direction, and r k−1 rows and by reshaping operations, respectively. Figure 2 shows that the TR decomposition [10] uses three third-order core tensors 𝒢 (𝑘) ∈ ℝ r k ×n k ×r k+1 , k = 1,2,3 to express a third-order tensor 𝒜 ∈ ℝ 𝑛 1 ×𝑛 2 ×𝑛 3 by tensor contractions: 𝒜 = 𝒢 (1) ∘ 𝒢 (2) ∘ 𝒢 (3) ,
where tensor contractions are calculated between tensors and also between 𝒢 (1) and 𝒢 (3) . [r 0 , r 1 , r 2 , r 3 ] expresses the TR-ranks. Because of the ring structure of TR tensors, r 0 = r 3 do not need to be forced to equal 1, which is used to distinguish between TR and TT structures. TR structure is another special case of tensor networks and steps of TR decomposition [10] are described in Algorithm 2.
Tensor-Train and Tensor-Ring Decompositions on GPUs
In this section, we propose three parallel optimizations for TT decomposition in Algorithm 1 and TR decomposition in Algorithm 2.
In TT and TR decompositions, the matrix SVD operations take up the most time, reaching 67%. The traditional SVD operation is not matched to GPU structure because of its low parallelism characteristics. As a substitute, Jacobi SVD [8] is adopted to match the GPU's high-parallelism feature. Jacobi SVD needs an iteration number and an accuracy to determine when the algorithm terminates. Under the single precision of data and calculation, this paper set the maximum iteration to 100 and the accuracy to 10e-8 for getting the minimum error in experiments.
The algorithm stops when the number of iterations reaches the maximum number of iterations or the error between the repaired matrix and the original matrix reaches a preset threshold.
Diagonal matrix and matrix multiplication is the operation with the second longest time occupation in the eleventh line of Algorithm 1 and the eleventh and fifth lines of Algorithm 2. The time cost of these operations increases rapidly with the dimension size of data. We find that the traditional processes introduce redundant calcula-tions because values exist only on the diagonal. Therefore, we accelerate the computation using the following parallel computation method:
where s k and 𝑽 k T represent the k-th value and row of matrix 𝑺 on the diagonal and matrix 𝑽 T . This parallel computation method takes advantage of parallelism and reduces redundant computations.
The sixth line to the ninth line of Algorithm 1 and the fifth line to the seventh line of Algorithm 2 are the element-wise products which can be calculated in parallel. We utilize the following parallel method to perform the element-wise product s ⋅ s, with m = #(s): parallel (s m−k+1
Figure 3: A schematic diagram of the tensor's layout in memory.
Figure 3 exhibits a third-order tensor' column-major layout in memory. The tensor data is stored as a front slice of the column master. We adopt the column-major layout in memory to meet the data reading requirements of two libraries: cuSOLVER and cuBLAS. In addition, this kind of memory access method can directly get the mode-1 unfolding of the tensor without the tensor matricization reducing the overhead.
To reduce the memory footprint and the overhead of truncation operations in Algorithm 1 and Algorithm 2, the front truncation sub-sections of matrix 𝑽 T and vector s are calculated directly in the eleventh line of Algorithm 1 and the eleventh and fifth lines of Algorithm 2. We utilize the direct conversion in memory to reduce the overhead of tensor permuting operations in the tenth and eleventh lines of Algorithm 2. Moreover, through this optimized memory access method, the reshape operations are eliminated in Algorithms. These algorithms generate a lot of intermediate variables, which introduces much memory footprint and time overhead. Therefore, we reuse the allocated memory and dynamically delete and allocate intermediate variables in GPU memory to reduce memory consumption.
The input and output data volumes of TT and TR decompositions increase quickly with tensor sizes, which results in high time consumption of data transfer between CPUs and GPUs. To reduce the overhead of access space in transmission, the cores 𝒢 1 , 𝒢 2 , 𝒢 3 are combineded into an array 𝑐 . Meanwhile, [𝑛 1 , 𝑛 2 , 𝑛 3 ] are used to store dimensions of the input tensor and [r 0 , r 1 , r 2 , r 3 ] are used to store TR-ranks or TT-ranks. Therefore, k-th core is acquired through 𝒢 𝓀 = reshape(𝐜(∏ r j n j r j+1 k−1 j=1 , ∏ r j n j r j+1 k j=1
), [r k−1 , n , r k ]) , 𝐜(∏ r j n j r j+1 k−1 j=1
, ∏ r j n j r j+1 k j=1
) denotes the elements from ∏ r j n j r j+1 k−1 j=1
to ∏ r j n j r j+1 k j=1
.
Our experiments run on a server with 80 GB host memory. The server is equipped with two Intel Xeon E5-2640 V4 CPUs. Each CPU has ten cores supporting twenty hardware threads. Moreover, the server is equipped with one Tesla V100 GPU with 32GB device memory and 5,120 CUDA cores @1.53 GHz. We focus on the speedups of our experiment result: speedup = (CPU running time)/(GPU running time). The relative square error (RSE) is utilized to measure the error of data before and after decomposition: RSE = ||𝒜 − 𝒢 (1) ∘ 𝒢 (2) ∘ 𝒢 (3) || F /||𝒜|| F . The experiment tensor data are obtained by tensor contractions of three small tensors. For the Jacobi SVD, under single precision, the accuracy ɛ is set to10e-8, and the max iteration time is set to 100.
The speedups and running time of our optimized third-order TT decomposition are exhibited in Figure 4. The tensor sizes vary from 100 × 100 × 100 to 1,200 × 1,200 × 1,200 . The CPU implemen-tation is referred from MATLAB code [9]. Because of the GPU memory size, the maximum tensor size that can be processed is 1,200 × 1,200 × 1,200 on Tesla V100 GPU. Compared with the CPU implementations, the optimized GPU implement-tation obtains 14.25× on average and up to 24.80× speedups, which are higher than the GPU baseline implementation. The RSE of CPU and GPU are on the 10e-4 level. In our experiment, the speedups of the optimized implementations have a general upward trend. The speedups and running time of our optimized third-order TR decomposition are exhibited in Figure 5. The CPU implementation is referred from MATLAB code [11]. Compared with the CPU implementations, our optimized GPU implementation achieves 11.35 × on average and up to 21.77× speedups, which are higher than GPU baseline implementations. The RSE of CPU and GPU are also on the 10e-4 level. Because of the overhead of iteration and data transfer, the speedup is less than one when the size of tensor is 100 × 100 × 100. The speedups of the optimized TR decomposition keep increasing with the size of tensor.
High-performance third-order tensor-train and tensor-ring decomposition implementations on GPUs are proposed in this paper. To improve the efficiency of the algorithms, three optimization strategies are proposed. First, efficient memory access is proposed to reduce the memory footprint. Second, parallelization strategies are widely adopted in algorithms to match GPUs. Third, we use the highparallel Jacobi SVD to reduce time for critical calculations.
Moreover, we experimentally verify the advantages of our optimized decomposition algorithms. The third-order TT and TR decompositions get maximum 6.67× and 6.36× speedups. Implementing multi-GPU implementations of high-order TT and TR decompositions is our future work. Meanwhile, the optimized third-order TT and TR decomposition algorithms will be combined into the cuTensor library [6].