-In this paper, we discuss several possible implementations of GPU-targeted gradient descent algorithm for dimensionality reduction. Four realizations of the gradient descent algorithm are created using HIP, a new framework for GPGPU programming. We get six times performance improvement over a multithreaded CPU version using AMD Radeon RX Vega 56.
At first, we should define a cost function. In the case of Sammon’s mapping the cast function is as follows:
In the present days, data mining becomes more and more popular research area. Data mining techniques process enormous amounts of data in order to detect significant dependencies in data. In many cases, such data can be interpreted as points in some multidimensional space.
Processing multidimensional data causes problems. The first problem is that multidimensional data requires more storage capacity. The second problem is that the time needed to process these data grows fast with the number of dimensions. But in many cases, multidimensional data contains substantial redundancy in terms of information. Therefore, one can map the data into the space of lower dimensionality without a major information loss. This process can be referred to as dimensionality reduction.
Dimensionality reduction techniques are often used in
different problems of image analysis (see [
There is a variety of dimensionality reduction methods.
The most famous methods are Principal component analysis
(PCA) and Independent Component Analysis (ICA), which
are linear methods. The major disadvantage of linear
methods is that they can find only linear dependencies within
data. The nonlinear dimensionality reduction methods are
Curvilinear Component Analysis (CCA) [
While there are GPU implementations for linear
dimensionality reduction techniques [
II. DIMENSIONALITY REDUCTION USING GRADIENT DESCENT
Let N denote the number of points, n denote the dimensionality of a high-dimensional base space and m denote the dimensionality of a low-dimensional target space. Our goal is to map points from the n-dimensional space into where dnij is the distance between points i and j in the base space, dmij is the distance between points i and j in the target space. In this case, the Euclidean distance was used. Gradient descent is a common iterative optimization algorithm. Its common definition is:
x(t+1) = x(t) - f, where the x = {x0, …, xk} denotes the parameters vector, t denotes the iteration number, denotes the speed constant and f denotes the gradient of the cost function f(x) in the current point.
If we take the gradient of the Sammon’s error with the coordinates of all the points in the lower-dimensional space as a vector of parameters, we will get the following equation, describing the Sammon’s mapping procedure: (1) (2) (3) (4) where yi = {y0i, …, ymi-1} denotes the coordinates vector of the i-th point , t denotes the iteration number, dnij denotes the distance between points i and j in the base space, dmij denotes the distance between points i and j in the target space, and denotes the descent speed constant.
The computational complexity of this algorithm is O(N2(n+m)). But it is worth noting that all points are updated independently during an iteration and therefore can be processed in parallel. This idea opens the possibility to use the graphics processing units (GPU).
III. HIP
HIP is a new framework for GPU programming. It is
developed by the GPUOpen initiative [
The terminology and programming model of HIP are inherited from CUDA. The interaction with a GPU is performed using API calls that perform data copy, synchronization, and execution launching. Each command is issued to some command stream. Commands in different command streams are executed in parallel if possible.
A function that was dedicated to perform on a GPU is called a kernel. Each kernel is executed using many GPU threads. GPU threads are grouped into thread blocks, whose size and count are specified on a kernel call. Available memory regions are global memory and shared memory. The global memory is an equivalent of the RAM. The shared memory has very high throughput but also a very small volume, so it was not used in the current work.
Equation 3 shows that points are updated independently from each other at one iteration. Therefore, several threads can be used to process points in parallel. The more threads used, the better performance may be obtained in theory. Since GPUs can process thousands of threads simultaneously, they may be used to accelerate the processing.
To implement the iterative gradient descent, 3 buffers were used:
1. the buffer to store the coordinates of points in the original space;
2. the buffer to store the coordinates of points in the target space obtained at the previous iteration;
3. the buffer to store the coordinates of points in the target space calculated at the current iteration.
We discuss two approaches for the gradient descent implementation. The first approach (denoted below as a standard approach) is to process data according to (3), i.e. for each point all distances to other point was computed, and then the point was updated. The second approach (denoted below as a modified approach) uses a fact that dnij= dnji. Therefore we may compute each distance only once and then perform an atomic update of points i and j. This approach reduces the number of inner cycle iterations from N(N-1) to N(N-1)/2, but also adds an overhead caused by the atomic operations. In this case, the linear indexing of the virtual distances matrix was used to exclude the conditional operators in the program code, which should be avoided in GPU programs.
We also discuss two possible memory layouts of data points. The points set can be stored as a 2D data matrix. The question is how to place the coordinates of one point. The first layout (denoted below as a flat layout) places the point coordinates in rows. This layout allows us to compute the distances between points faster since in this case we are able to use the vector load and store instructions. The second layout (denoted below as a transposed layout) places the coordinates of points in columns of the data matrix. This allows the neighboring threads to read the neighboring memory cells when accessing the data. This facilitates memory caching and allows the memory controller to merge many small memory access requests to one (memory coalescing).
Two implementation approaches and two memory layouts give 4 implementations: 1. 2. 3. 4.
a standard approach with a flat layout; a modified approach with a flat layout; a standard approach with a transposed layout; a modified approach with a transposed layout;
Proposed implementations were tested using a dataset with random points with coordinates from the range [0, 255]. We measured the average execution time for various thread block sizes and datasets sizes. The time of the data transfer from RAM to GPU was not accounted, only the kernel time was measured.
We compared the experimental results with the results of CPU implementation. This implementation was created using an OpenMP as a parallel framework and was compiled with the AVX cycle vectorization allowed. We measured execution time for both single-threaded and 16-threaded versions of this realization.
The experiments configuration: were run on the following
RAM: 8 GB DDR4-3200;
Number of Points found in Table II. The comparison of CPU and GPU final error values is shown in Table III.
Let us analyze the results. It is clear that performance improvement grows along with the task size. The most effective implementation is implementation 2 (a modified approach with a flat memory layout). In our case, the ability to compute distances faster was more important than advanced caching and memory coalescing. The data given by the profiler shows that the main performance limiting factor is the HBM memory performance. operations adds an additional load to the memory system. This partially explains why the thread block size of 16 threads was the best. Since threads of one thread block are executed in groups of 64 threads on the AMD GPU, the thread block of only 16 threads causes only 25% occupancy but also reduces the number of simultaneous memory requests. The results from GPUs with different memory types, such as GDDR5/6, probably will be different.
In this paper, we have discussed different approaches to implement the dimensionality reduction method based on gradient descent for GPUs. Several experiments were taken using an AMD Radeon Vega 56 GPU. The experimental results show approximately 6 times performance improvement in comparison to the multithreaded version executed on the 8-core CPU, and almost 60 times improvement over single-threaded implementation.
Future work will be devoted to the comparison of the proposed implementations using NVIDIA and AMD GPUs.
The work was partly funded by RFBR according to the research project 18-07-01312-a in parts of «2. Method» - «3. Experiments» and by the Russian Federation Ministry of Science and Higher Education within a state contract with the «Crystallography and Photonics» Research Center of the RAS in parts «1. Introduction» and «4. Conclusion».
The HBM (High Bandwidth Memory) features the very wide data bus along with relatively low clocks. Since in our case threads access memory cells located relatively far from each other, the caching is ineffective, and the memory requests often cause delays in execution. The use of atomic