Implementation of the Principal Component Analysis onto High-Performance Computer Facilities for Hyperspectral Dimensionality Reduction: Results and Comparisons

4.1. Preliminary Issues

4.2. CUDA Implementation of PCA Jacobi

• Data Initialization. The CUDA implementation of the PCA Jacobi algorithm requires the transfer of the hyperspectral image from the CPU (host) memory to the global memory of the GPU (device). Because some operations of the PCA Jacobi algorithm are implemented by means of cuBLAS library, which uses column-major storage, the hyperspectral image in column-major order has been transferred from the host to the GPU memory. Moreover, in order to simplify the transfer operation, the Thrust vector containers have been used for both, host and device. With this approach, we could completely avoid the need to manually manage the memory on the host and device using a Thrust vector for storing our data.
  In Code 4.1 the hyperspectral image is loaded into a Thrust host vector, h_image_in, by means of an iterator (lines 1.1, 1.2, 1.3, 1.4): in line 1.1, the filename with the hyperspectral image is passed as an argument of the command line; in lines 1.2, 1.3, and 1.4, the original image is loaded into a Thrust host vector container with BANDS × PIXELS elements, h_image_in. Finally, this image is transferred to the device memory, d_image_in, by means of thrust::copy (lines 1.5, 1.6). The result is that the hyperspectral image is loaded into the device as a sequence of bands. This is, the image matrix is row-major ordered. This fact must be considered whenever a cuBLAS function is called because cuBLAS uses column-major storage. This issue is achieved by using the transpose parameter whenever a cuBLAS function is called with the hyperspectral image.

  Code 4.1. Data initialization
  1.1	ifstream ifile(argv [1]);
  1.2	thrust::host_vector<float> h_image_in(BANDS * PIXELS);
  1.3	istream_iterator<float> beg(ifile), end;
  1.4	thrust::copy(beg, end, h_image_in.begin());
  1.5	thrust::device_vector<float> d_image_in(BANDS * PIXELS);
  1.6	thrust::copy(h_image_in.begin(), h_image_in.end(), d_image_in.begin());

• Image preprocessing. In this stage, the original hyperspectral image is going to be mean-centered. In order to parallelize this step, we have divided it into two sequential calculations: the average per every band of the hyperspectral image and the subtraction of this average from the original hyperspectral image.

  2.1.

  Average per band calculation. For each band of the original image (BANDS × PIXELS matrix), its average is calculated and stored into a BANDS array. In order to calculate the average per band in parallel with the cuBLAS library, the BANDS × PIXELS matrix is multiplied with ones vector (vector with PIXELS elements) by means of the cublasSgemv function. According to the cuBLAS documentation, cublasSgemv performs the following matrix-vector multiplication:

  $$y=\alpha op\left(A\right)x+\beta y$$

  (11)

  where A is a m × n matrix that is stored in column-major format located in the GPU (device), x and y are vectors, which must be also located in the GPU, and α and β are scalars. The operand op(…) denotes the use or not of the transpose operation. In our average calculation:

  • A is the BANDS×PIXELS matrix which contains the original hyperspectral image in column-major format. According to Code 4.1, this matrix is already in the device, d_image_in (see lines 1.5 and 1.6 in Code 4.1). Moreover, because this matrix is in row-major format, it needs to be transposed in the cuBLAS function.
  • x is a vector with PIXELS elements filled to one. This vector is called d_ones and it has been filled to 1 by means of Thrust (see line 2.1.2 of Code 4.2.1).
  • α is set to 1/ PIXELS, and β is set to 0 (line 2.1.3 of Code 4.2.1).
  • y is the array with BANDS elements that stores the average per band. In Code 4.2.1 this vector is called d_means and it is filled with the result of cublasSgemv function (see lines 2.1.1 and the next to last parameter of cuBLAS function in line 2.1.4).

  Code 4.2.1. Average per band calculation
  2.1.1	thrust::device_vector<float> d_means(BANDS); // to store average per band
  2.1.2	thrust::device_vector<float> d_ones(PIXELS, 1.f); //ones vector
  2.1.3	float alpha = 1.f / (float) PIXELS; float beta = 0.f;
  2.1.4	cublasSgemv(handle,         //handle for use cuBLAS       CUBLAS_OP_T,       //the d_image_in is transposed       PIXELS, BANDS,      //rows and columns of d_image_in       &alpha, thrust::raw_pointer_cast 1(d_image_in.data()),       PIXELS,         //leading dimension of d_image_in       thrust::raw_pointer_cast(d_ones.data()),       1,          //leading dimension of d_ones       &beta,        //β scalar       thrust::raw_pointer_cast(d_means.data()),       1 );          //leading dimension of d_means array

  ¹ thrust::raw_pointer_cast converts the device vector container to a raw pointer.

  2.2.

  Subtraction of the average per band from the original hyperspectral image. For each band of the original image, the average of this band must be removed. The subtraction of the average of a band from every pixel in such band is carried out by means of a Thrust transformation. Transformations are algorithms that apply an operation to each element in a set of (zero or more) input ranges and then store the result in a destination range. The transformation thrust::transform in line 2.2.2 of Code 4.2.2, uses the function object thrust::minus to subtract a device vector, d_image_in, from another device vector, d_means. As d_means is a vector with BANDS elements and d_image_in contains PIXELS columns with BANDS elements, the d_means array must be subtracted PIXEL times for all of the pixels in the original hyperspectral image, d_image_in. This purpose is achieved by combining three Thrust fancy iterators: permutation_iterator, transform_iterator and counting_iterator. Each of these iterators is created by a function called make_XXX, where XXX denotes the name of the Thrust fancy iterator. The permutation_iterator permutes the range of values in d_means by the index scheme obtained by the transform_iterator. In general, permutation_iterator allows for operating on a specific set of values in a sequence instead of the entire sequence. The transform_iterator changes the range generated by the counting_iterator (linear index beginning in 0) to its associated row index by calling the linear_index_to_row_index functor (A C++ functor is an object that acts just like a function, but it has the advantage of being stateful, meaning that it can keep data reflecting its state between calls) (see line 2.2.1 in Code 4.2.2). This last functor transforms the linear index that is generated by the counting_iterator to the row index of a matrix.

  Code 4.2.2. Subtraction the average per band from the original hyperspectral image
  2.2.1	template<typename T> struct linear_index_to_row_index: public thrust::unary_function<T, T> {T Ncols; // --- Number of columns__host__ __device__ linear_index_to_row_index(T Ncols):Ncols(Ncols) {}__host__ __device__ T operator()(T i) {return i / Ncols;}};
  2.2.2	thrust::transform(d_image_in.begin(), d_image_in.end(),          thrust::make_permutation_iterator             (d_means.begin(),             thrust::make_transform_iterator(                thrust::make_counting_iterator(0),                linear_index_to_row_index<int>(PIXELS))),          d_image_in.begin(),thrust::minus<float>());

• Covariance computation. This stage computes the covariance matrix associated with the original image in parallel in the GPU (device) by means of Thrust and the cuBLAS library. The process to calculate this covariance matrix is carried out by the cublasSgemm function, which multiplies the preprocessed matrix by its transpose. According to the cuBLAS documentation, cublasSgemm performs the matrix-matrix multiplication, following the next operation:

  $$C=\alpha op\left(A\right)op\left(B\right)+\beta C$$

  (12)

  where A, B, and C are matrices in column-major format located in the GPU (device), and α and β are scalars. A is an m × k matrix, B is a k × n matrix, and C is an m × n matrix. As mentioned before, the op(…) denotes the use or not of the transpose operation. In our covariance matrix calculation:
  • A and B are BANDS × PIXELS matrices that contain the preprocessed hyperspectral image. According to our Data Initialization code, these matrices are already in the device as d_image_in (see lines 1.5 and 1.6 of Data Initialization code). Moreover, because d_image_in is in row-major format, it needs to be transposed in cuBLAS function. This is the reason why op(A) must be CUBLAS_OP_T (transposed d_image_in ) and op(B) must be CUBLAS_OP_N (original d_image_in, not transposed) (see call to cublasSgemm function in line 3.3).
  • m is equivalent to BANDS, n to BANDS and k to PIXELS.
  • α is set to 1, and β is set to 0 (line 3.2 of the Covariance computation code).
  • C is the result matrix that corresponds with the covariance matrix, d_cov, a BANDS × BANDS matrix (see lines 3.1 and the next to last parameter in the cuBLAS function).
  The last step for the covariance computation stage is to normalize the result of cublasSgemm using a Thrust transformation along with a Thrust iterator, as it is shown in line 3.4 of code 4.3. In this normalization, the elements of the d_cov device vector are divided by PIXELS-1.

  Code 4.3. Covariance computation
  3.1	thrust::device_vector<float> d_cov(BANDS * BANDS);
  3.2	alpha = 1.f; beta = 0.f;
  3.3	cublasSgemm (handle, CUBLAS_OP_T, //the d_image_in is transposed to get column-major ordered format op(A) CUBLAS_OP_N, //op(B) BANDS,   //rows of d_image_in (A) and d_cov (C) BANDS,   //columns of d_image_in (B) and d_cov (C) PIXELS,   //columns of d_image_in (A) and rows of d_image_in (B) &alpha,  //scalar set to 1 thrust::raw_pointer_cast(d_image_in.data()), //matrix d_image_in (A) PIXELS,  //leading dimension of d_image_in (A) thrust::raw_pointer_cast(d_image_in.data()), //matrix d_image_in (B) PIXELS,  //leading dimension of d_image_in (B) &beta,  //scalar set to 0 thrust::raw_pointer_cast(d_cov.data()),  //matrix d_image_in (C) BANDS  //leading dimension of d_image_in (C));
  3.4	thrust::transform (d_cov.begin(), d_cov.end(), thrust::make_constant_iterator((float) (PIXELS - 1)), d_cov.begin(), thrust::divides<float>());

• Eigenvector decomposition. This stage gets the eigenvectors associated to the covariance matrix computed in the previous stage by means of the cyclic Jacobi, which is divided into three parts:

  4.1.

  The first part of this stage deals with a sequential code on the host side, which is in charge of finding the maximum element of the covariance matrix that is not greater than a provided stop factor (ε). The calculation of the maximum value and its position ($ro{w}_{max}$, $co{l}_{max}$) are carried out over the upper triangular matrix. After obtaining the maximum value position, the Equations (2)–(5) are calculated, where Equation (2) in this implementation is as the Equation (7), and i and j corresponds with the position of the maximum element in the covariance matrix ($ro{w}_{max}$, $co{l}_{max}$). Moreover, the set of operations in (13) is computed over the identity matrix.

  $$\begin{gathered} identity\left[ro{w}_{max},ro{w}_{max}\right]=cosA \\ identity\left[co{l}_{max},co{l}_{max}\right]=cosA \\ identity\left[ro{w}_{max},co{l}_{max}\right]=sinA \\ identity\left[co{l}_{max},ro{w}_{max}\right]=-sinA \end{gathered}$$

  (13)

  where l = 1. BANDS (columns of the covariance matrix).

  4.2.

  The second part of this stage is executed on the GPU (device). A CUDA kernel is launched with the maximum number of threads per block and so many blocks as the relation between the number of BANDS and the $GP{U}_{BLOCKSIZE}$ (1024 threads):

  $$\begin{array}{l}threadsPerBlock=GP{U}_{BLOCKSIZE}\\ numberBlocks=\frac{BANDS+threadsPerBlock-1}{threadsPerBlock}\\ kernel\ll \langle numberBlocks,threadsPerBlock\gg \rangle (\dots )\end{array}$$

  (14)

  This kernel computes in parallel the Jacobi method with the equations in (15):

  $$\begin{array}{l}{d}_{cov}\left[ro{w}_{max},l\right]=cosA\cdot d_{cov}\left[l,ro{w}_{max}\right]-sinA\cdot d_{cov}\left[l,co{l}_{max}\right]\\ d_{cov}\left[co{l}_{max},l\right]=cosA\cdot d_{cov}\left[l,co{l}_{max}\right]+sinA\cdot d_{cov}\left[l,ro{w}_{max}\right]\\ d_{cov}\left[l,ro{w}_{max}\right]=d_{cov}\left[ro{w}_{max},l\right]\\ d_{cov}\left[l,co{l}_{max}\right]=d_{cov}\left[co{l}_{max},l\right]\end{array}$$

  (15)

  As the goal of this stage is zeroing the off-diagonal elements, both parts of this stage are repeated until a stop condition is fulfilled. Because each iteration of this stage may undo the zeroes that are obtained in previous iterations, these iterations are repeated until the maximum value of the covariance matrix becomes smaller than a provided stop factor (ε). For each iteration, the covariance matrix becomes a diagonal matrix whose main diagonal represents the eigenvalues. Moreover, the eigenvector matrix is computed in parallel on the host side, as the equation shown in (16) using the cuBLAS library.

  $$P=identit{y}_1\cdot identit{y}_2\cdot identit{y}_3\cdot \cdot \cdot \cdot \cdot identit{y}_k$$

  (16)

  where k represents the last iteration of the eigenvector decomposition step.

  4.3.

  The last part of this stage consists of extracting and sorting both eigenvalues and eigenvectors. The eigenvalues are located in the main diagonal of the covariance matrix obtained in parallel in the last stage. A kernel is in charge of extracting these eigenvalues in parallel from the covariance matrix (see line 4.3.2 of Code 4.4.3). This kernel is launched with one block of BANDS threads and each thread extracts its value of the main diagonal of the covariance matrix and stores it in a Thrust device vector (see lines 4.3.1.), d_diagonal.

  Code 4.4.3. Extraction and sort of eigenvalues and eigenvectors
  4.3.1	thrust::device_vector<float> d_diagonal(BANDS);
  4.3.2	diagonalMain<<<1,BANDS>>> (…, …); //parameters covariance matrix and d_diagonal
  4.3.3	thrust::device_vector<int> d_pos_data(BANDS*BANDS);
  4.3.4	thrust::sequence(thrust::device, d_pos_data.begin(), d_pos_data.end(), 0);
  4.3.5	thrust::sort_by_key (d_diagonal.begin(),d_diagonal.end(),d_pos_data.begin(),thrust::greater<float>());
  4.3.6	sortEigenvectors <<<1,BANDS>>> (…,…,…) //parameters: eigenvector matrix, the order of eigenvalues and the result of the sort

  The main diagonal is sorted in descending order by means of two Thrust transformations: thrust::sequence and thrust::sort_by_key. The first transformation, thrust::sequence, fills a range with a sequence of numbers using the thrust::device execution policy for parallelization; the second transformation, thrust::sort_by_key, performs a key value sort, where the key value is provided by the thrust::sequence transformation. The result is a sorted device vector, d_diagonal, in descending order while using the thrust::greater comparison operator.

  Moreover, the eigenvectors must be sorted according to the order that is followed by the eigenvalues. As each eigenvector contains BANDS elements, the sort operation is carried out changing the order of each eigenvector in matrix P, according to the order followed for the eigenvalues, which is kept in the d_pos_data device vector after sorting the eigenvalues. This sort is computed by a kernel, sortEigenvectors, which is launched with a block with BANDS threads. Each thread is in charge of an eigenvector.

• Projection and reduction. In this stage, the original hyperspectral image is projected onto the first principal component, that is, onto the first eigenvector. Although the original hyperspectral image must be projected onto the set of eigenvectors, for simplicity this projection is performed only onto the first eigenvector. The projection is computed by multiplying the original hyperspectral image by a vector, the first eigenvector stored in matrix P. This operation is performed by a call to the cublasSgemm function (see Code 4.5).

  Code 4.5. Projection and reduction

  thrust::device_vector<float> d_projection (PIXELS,0); thrust::device_vector<float> d_eigenvector (BANDS); //copy the first eigenvector into d_eigenvector: thrust::copy_n (thrust::device, d_G_ordered.begin(), BANDS, d_eigenvector.begin()); alpha = 1.0f; beta = 0.0f; //scalars for cublasSgemv //d_projection = d_image_in x d_eigenvector (the first eigenvector): cublasSgemm(handle,  //handle for use cuBLAS CUBLAS_OP_N, CUBLAS_OP_N, PIXELS, BANDS,1, &alpha, //α scalar, thrust::raw_pointer_cast(d_image_in.data()), thrust::raw_pointer_cast(d_eigenvector.data()), &beta,  //β scalar thrust::raw_pointer_cast(d_projection.data()));

7.1. Performance of the Sequential PCA Jacobi Algorithm

7.2. Performance of the CUDA PCA Jacobi Algorithm

7.3. Performance of the MPPA PCA Jacobi Algorithm

7.4. Comparisons