cuPSO: GPU Parallelization for Particle Swarm Optimization Algorithms

The parameters of the Standard Particle Swarm Optimization (SPSO) algorithm are defined in Table 1. In general, the parameters $w$, $c1$, and $c2$ are changed in different practical situations, controlling the behavior and efficiency of the SPSO algorithm. $max\_pos$, $min\_pos$, $max\_v$, and $min\_v$ restrict the lower and upper boundaries of each dimension in the search space for the position and velocity respectively. $max\_iter$ performs the termination criterion of the number of iterations. $particle\_cnt$ is the total number of the particles. Moreover, $pos_{i}$, $v_{i}$, $fit_{i}$, $pbest\_pos_{i}$, and $pbest\_fit_{i}$ are the necessary information of each particle_i, which represent as the position, speed, fitness, best-known position, and best-known fitness respectively. The subscript $i$ is the index of each particle. Finally, $gbest\_pos$ is the best-fitness point ever found by the whole swarm, and $gbest\_fit$ is the best-fitness value obtained by corresponding best-fitness point $gbest\_pos$.

Algorithm 1 shows the full workflow of the SPSO algorithm, and it could be simply decomposed into 5 steps. The algorithm would continue operating from step 2 to step 5 until satisfying the termination condition as follows: (1) initialize the information of each particle, including position, velocity, fitness, the best-known position of each particle, the best-known fitness of each particle, global best-known position, and global best-known fitness. The position and velocity are initialized with a uniformly distributed random vector, constricted in the fixed range to ensure the convergence. (2) update the velocity of each particle as Equation 1 and the position of each particle as Equation 2 within the domain of the fitness function at time $t$ to explore the whole search-space. r1 and r2 are uniformly distributed on random value in [0, 1]. (3) calculate each particle fitness by its position and fitness function. (4) update the best-known position and the fitness, local best for short, while the new fitness is better than the previous one. (5) update the global best-known position and fitness, global best for short, while the fitness of the particle is better than the global best-known fitness.

In order to take advantage of the computing power of GPUs (graphics processing units), a parallel PSO (PPSO) has been developed based on the SPSO introduced in Section 2.1 to improve the execution speed of Steps 2-5 in Algorithm 1. The PPSO program is executed by CPU and GPU. Specifically, the particle information initialization (Step 1) and finalization (to obtain the final results) are done by the CPU, whereas the compute-intensive part (Steps 2-5) is performed on the GPU.

In PPSO (CHEN Feng, 2014), a thread is used to handle the computations for calculating the information of a particle, e.g., velocity, position, and fitness, as defined in Steps 2 and 3, where the workflow is illustrated in the ”1st kernel” (the first parallel code section) in Figure 1. The threads are grouped into blocks to update the best-known position and fitness for a group of particles, where the best data within each block is computed at the end of the ”1st kernel” in Figure 1. It is important to note that in order to speed up the process of obtaining the best data within each block, the parallel reduction operations are performed within each block (i.e., each thread group). Subsequently, the second parallel code (the ”2nd kernel”) is used to further derive the best position and fitness for the particle data across the thread groups (blocks), which is also performed by the parallel reductions. It is worthy to note that the first kernel code is performed by multiple thread blocks and each block performs the reductions internally to find the local best data points. On the contrary, the second kernel code is handed by a thread block to search for the global best data point from the output data derived by the thread blocks for the first kernel.

During the parallel reduction, it incurs extensive arithmetic operations for calculating the proper addresses of the reduction operands and the branching instructions in the reduction loops. A common approach is to unroll the reduction loops to reduce the above overheads since the addresses are computed offline and the sequential unrolled code avoids the branches. Nevertheless, unrolling the loops cannot remove the implicit synchronization between the first and second kernels, and during the iterative search process, the two kernels are invoked repeatedly, which introduces a significant amount of overheads. In this work, we develop the novel algorithm to tackle the problem.

In the practical test, we find that it is seldom the case that the conditions of Steps 4 and 5 in Algorithm 1 are satisfied; the possibility of the satisfied condition may be less than 0.1%. Specifically taking Step 4 as an example, the $fit_{i}$ calculated by the $i$-th thread is usually less than the $pbest\_fit_{i}$ listed in line 14 of the algorithm. Based on the observation, we develop the shared memory based solution, which is referred to as the queue algorithm, to let each thread to update the queue (the in-order accessed arrays and the index to the arrays guarded by the atomic addition operation) if its calculated $fit$ value is larger than that of the global value $gbest\_fit$. The major concept is the sequential accesses (enforced by the atomic add operation) to the shared memory arrays and is listed in line 1-4 in Algorithm 2. To further identify the best configuration of each block, the thread with its thread identifier equal to zero (i.e., $threadIdx.x==0$) is selected to scan the entire queue, finds the local best configuration, and updates the local best information to the auxiliary arrays, as listed in line 7-20 in Algorithm 2. In particular, the local best information for a block (stored in $bestFitQueue[0]$) is updated to the auxiliary arrays indexed by the block identifier ($auxFit[blockIdx.x]$ and $auxPos[blockIdx.x]$), as specified in line 17-18 Algorithm 2. Regarding the complexity of the step for updating the intra-block best, we reduce the complexity from $\mathcal{O}(\log{}n)$ of using the reduction operation to theoretical $\mathcal{O}(1)$ of using the queue method²²2Thanks to the queue method that maintains the best-known fitness of each particle and performs the reduction operation if necessary, the resultant complexity is close to $\mathcal{O}(1)$ at most of the time.. Note that the similar approach is also adopted when using the second kernel to update the global best information.

As described in Section 3.1, the local best information of each block is stored in the corresponding position of the auxiliary arrays (line 17-18 of Algorithm 2). This design allows the first threads of each block to store their best data to the auxiliary arrays concurrently. Later, the second kernel is used to further search for the global best information from within the auxiliary arrays. In such a design, it requires the two kernels to compute the global best configuration, as illustrated in Figure 1. In order to alleviate the synchronization overhead (between the two kernels) and shorten the time of launching the second kernel, we fuse the two kernels together via computing the best configuration (originally done in the second kernel) at the end of the first kernel. That is, Algorithm 3 is proposed to replace the line 17-18 of Algorithm 2.

Algorithm 3 reduces the execution time because each thread block compares its local best configuration against the global best one as soon as its local best configuration is available. In the fused kernel, each of the thread blocks finds its local best configuration and updates to the global configuration if necessary. Such a design eliminates the need to write to the auxiliary arrays and relaxes the synchronization of the thread blocks (i.e., the second kernel should be started after the threads in the first kernel are terminated), where a lock is introduced to ensure the global best information ($gbest\_fit$ and $gbest\_pos$) is accessed sequentially by utilizing the atomic lock. Therefore, it would deliver better performance, especially for the one dimensional data for $gbest\_fit$, which has no bad side effort.

In general, the developers often use an array of structures (AoS) layout to design a program, since the AoS layout is easily understood and supported directly by most programming languages. However, AoS layout is almost the worst case scenarios to develop a parallel program as CUDA. On the other hand, structure of arrays (SoA) layout can reduce memory access time because each thread can do the coalesced access for continuous memory address. Accordingly, we change the data structure as the SoA layout in our implementation.

The data structure how we implement the SoA layout of the PSO algorithm. With the SoA adjustment, every element is seen as to be saved in an array with a fixed value, and it guarantees all threads of the same warp that could access consecutive values in global memory and obtain a better performance than AoS layout. In the high dimension case of the PSO algorithm, we follow the SoA rule in our implementation. As Figure 2, all threads accessing at the same dimension, which are still satisfied with the condition of the coalesced memory access, making the memory access each element efficiently.

NVIDIA GPUs provide 64KB of constant memory that is treated differently from standard global memory. The constant memory has its own cache. When all threads access the same memory address simultaneously and threads in the same block re-access that address continuously, the constant memory can reach the most effective performance. Thus, using constant memory instead of global memory, which can conserve the memory bandwidth. Although the constant memory can complete the memory accessing task with the less latency, it constrains the data usage to be read-only. In our study case, we can place the parameters $w$, $c1$, $c2$, $max\_pos$, $min\_pos$, $max\_v$, $min\_v$, $particle\_cnt$ in constant memory, making all threads accessed to the constant memory in the same time, with an eye to reducing the memory latency time.

The shared memory buffers reside physically on the GPU as opposed to residing in off-chip DRAM, which is much faster than the global memory. In fact, if there are no bank conflicts, the shared memory latency is about 100 times lower than the uncached global memory latency. Because the shared memory is shared by threads in a thread block, it has a mechanism for threads to cooperate. As subsection 2.2, using reduction would access the global memory frequently so that we can utilize shared memory to reduce the memory bandwidth latency. Similarly, shared memory can also be used on the GPU Queue method proposed at subsection 3.1. Owing to shared memory providing 48KB for each stream multiprocessors (SMs) only, we should store the index of the particle number rather than all axes positions on high dimension cases. After getting the best fitness of the particle, we can use the index to fill all axes position back, which not only saves the latency time but also reduces the shared memory usage to ensure maximum parallelization for SMs.

The PSO algorithm relies on the sequence of random numbers to find the best solution within the designated space. To generate the high-quality pseudorandom and quasirandom numbers on the GPU side, it is important to use the cuRAND (Corporation, 2010) library offered by the CUDA toolkit. While the users can provide the random number generator function on the GPU side by porting that function seen in the CPU side, it is time-consuming and inefficient for GPU execution since the custom-made library should be a thread-safe implementation and it may not take advantage from the low-level primitives offered by the CUDA toolkit. Based on our experiments, the cuRAND library outperforms the custom-made implementation on the PPSO program by a factor of 1.1. Therefore,in our experimental results, we adopt the cuRAND library (i.e., the $curand\_uniform\_double()$ function) to run the experiments by default.

The hardware and software platforms used for the experiments are listed in Table 2. The five PSO algorithms have been implemented and tested in the experiments, 1) CPU for the serial version defined in Section 2.1, 2) Reduction for the parallel reduction method introduced in Section 2.2, 3) Loop Unrolling for the loop unrolling method introduced in Section 2.2, 4) Queue for the proposed queue algorithm defined in Section 3.1, and 5) Queue Lock for the proposed queue lock algorithm defined in Section 3.2, where all the algorithms are running on the GPU, except for the first serial version. The double precision floating point numbers are adopted in the PSO implementations. Each of the reported data is the average numbers of the execution time for 10 runs, removing the maximum and minimum numbers. In order to observe the impacts of the improved algorithms, the running time at the GPU side is measured and reported. The above implementations are used to solve the problems with both the low and high dimension space, denoted as 1D problem and 120D problem, respectively.

As listed in Table 3, the parallel algorithms running on the GPU outperform the serial version on the CPU. The execution time of the CPU version increases almost linearly with the growth of the computation load (related to the number of particles and iteration). On the contrary, the parallel versions have relatively stable execution times across the different workloads since the problem size is relatively small which can be covered by GPU. Specifically, each particle is assigned with a GPU thread to handle the required computation. The maximum number of particles is 2,048 in this experiment, and the Nvidia GPU has 3,584 CUDA threads. The ranking, in terms of the execution time, of the five algorithms are plotted in Figure 3. Obviously, the Queue Lock algorithm achieves the best performance than all other parallel versions, and the Queue algorithm is the second best choice. Particularly, the Queue Lock algorithm is 2.2 times faster than the Reduction algorithm, which is the state-of-the-art work. Table 4 shows the speedup ratio of the serial version running on the CPU to the Queue Lock algorithm on the GPU. It shows a maximum of 195x speedup when handling 65,536 particles, and the speedup drops for handling 131,072 particles. Based on our analysis, the GPU is overloaded in such a case, where the GPU capacity is full and it cannot handle the workload efficiently (context switching among thread groups do not help the performance).

We configure the PSO algorithms to handle the problem defined by the 120 dimensional space. In the experiment setup, we select the Queue algorithm to compare against the serial version since the time saved by the Queue Lock algorithm is negligent in the high dimension space since the execution of the first kernel dominate the total execution time and the synchronization overhead saved by the Queue Lock is relatively small. Hence, we think the Queue algorithm is a better solution for the high dimension problem. The speedups achieved by the Queue algorithm are listed in Table 5. Because of a higher dimension space, the peak speedup is achieved (225x) when the number of the particles is 32,768, which is only the half of the number (65,535) reported in Table 4.