# High Performance Graph Data Imputation on Multiple GPUs
## 4 Efficient Convolutional Imputation Algorithm on GPU
### Design and Implementation of the Baseline GPU Convolutional Imputation Algorithm
#### 4.1.3 Parallelization of the Algorithm  
In Algorithm 1, the major computations are graph Fourier transform, inverse graph Fourier transform, and singular value soft-threshold in lines 8-10 in each iteration. Since the process of inverse graph Fourier transform is quite similar to the process of graph Fourier transform, we only discuss graph Fourier transform in the following sections for conciseness. For the outermost for-loop, we set the configurable parameter \(C\) as the number of iterations, and \(\lambda^{j}(j=1,2,...,C)\) is the regularization parameter for iteration \(j\). The result for iteration \(j\) with \(\lambda^{j}\) is a start for the next iteration \(j+1\) with \(\lambda^{j+1}\). For the nested while-loop, the program ends in the maximum iterations \(T\) or the convergence condition satisfies the formula \(\|\mathcal{G}_{t}^{\text{est}}-\mathcal{G}_{t-1}^{\text{est}}\|^{2}/\|\mathcal{ G}_{t-1}^{\text{est}}\|^{2}\leq\epsilon\), where \(\epsilon>0\) is an adjustable parameter that can affect the result of computing and \(t\geq 2\). In the while-loop:  
* For the graph Fourier transforms computation in line 8 of Algorithm 1, the algorithm accesses data along the third dimension of the graph-tensor. That is to access the graph-tensor model-3 tube by model-3 tube. Based on the data storage mentioned earlier, it cannot provide coalesced memory access for the GFT computation on GPU. The premise of coalesced memory accesses is that accesses must be sequential and addresses aligned. Therefore, we introduce the step of data reorganization. First,  
Figure 3: Reconstruction of graph-tensor data for computing graph Fourier transforms in the batched scheme. we utilize the eigenvalue solver method provided in the cuSLOVER library [24] to get the graph Fourier transform matrix \(\mathbf{U}\), which is composed of eigenvectors of the graph Laplacian matrix \(L\). Then, we propose a mapping to reconstruct the graph-tensor \(\mathcal{G}^{\text{impute}}\), where the original slice-by-slice layout (the index of \(\mathcal{G}^{\text{impute}}_{ijk}\) is \((i-1)+m*(j-1)+m*(k-1)\)) of data is reorganized into the tube-by-tube layout (the index of \(\mathcal{G}^{\text{impute}}_{ijk}\) is \(N*(i-1)+m*N*(j-1)+(k-1)\)). Finally, we utilize batched matrix-matrix multiplication in the cuBLAS library [24] to calculate the GFT of \(\mathcal{G}^{\text{impute}}\), and then convert the result back to the original data layout in the graph spectral domain to get \(\widetilde{\mathcal{G}}^{\text{impute}}\). Since the benefits obtained from batched matrix multiplication far outweigh the overhead introduced by reorganizing the data, the overall algorithm performance is improved. We design a batched scheme for the GFT computation by reorganizing the graph-tensor data to improve parallelism. Figure 3 shows an example to illustrate the detailed scheme. For a graph-tensor \(\mathcal{G}\) with size \(2\times 3\times 4\) and graph Fourier transform matrix \(\mathbf{U}\) with size \(4\times 4\), the original computation is that each frontal slice in the graph spectral domain is a linear combination of data matrices on the graph, where data access is random. The result of \(\mathcal{G}(1,1,1)\) is the dot product of vector \(\mathbf{U}(1,:)\) and \(\mathcal{G}(1,1,:)\), \(0*0+2*6+1*2+1*8=22\). However, it is more time consuming because random accesses to the data are slower than sequential accesses. Therefore, we reorganize the graph-tensor data to achieve sequential accesses, such as the \((1,1)\)-th model-3 tube is reorganized into the 1 st-column of 1-st frontal slice, the \((2,1)\)-th model-3 tube is reorganized into the 2 nd-column of 1-st frontal slice, and so on. The graph Fourier transform matrix \(\mathbf{U}\) can now multiply each frontal slice matrix in batched by exploiting batched matrix-matrix multiplication. * For the singular value soft-threshold computation in line 9 of Algorithm 1, it includes singular value decomposition and matrix multiplication for each frontal slice, as shown in the the pseudo code in Algorithm 2.