Parallel GPU Implementation of Iterative PCA Algorithms

Principal component analysis (PCA) is a key statistical technique for multivariate data analysis. For large data sets the common approach to PCA computation is based on the standard NIPALS-PCA algorithm, which unfortunately suffers from loss of orthogonality, and therefore its applicability is usually limited to the estimation of the first few components. Here we present an algorithm based on Gram-Schmidt orthogonalization (called GS-PCA), which eliminates this shortcoming of NIPALS-PCA. Also, we discuss the GPU (Graphics Processing Unit) parallel implementation of both NIPALS-PCA and GS-PCA algorithms. The numerical results show that the GPU parallel optimized versions, based on CUBLAS (NVIDIA) are substantially faster (up to 12 times) than the CPU optimized versions based on CBLAS (GNU Scientific Library).

💡 Research Summary

The paper addresses two fundamental challenges in applying Principal Component Analysis (PCA) to large‑scale, high‑dimensional data: computational speed and the loss of orthogonality inherent in the widely used iterative NIPALS‑PCA algorithm. NIPALS computes principal components one at a time by alternating regressions of the data matrix on scores and loadings, but rounding errors accumulated during the many matrix‑vector multiplications gradually destroy the orthogonal relationship between successive components. Consequently, NIPALS is reliable only for the first few components, especially when the data columns are highly collinear.

To overcome this limitation, the authors propose a new iterative algorithm, GS‑PCA, which incorporates a Gram‑Schmidt (GS) re‑orthogonalization step at each iteration. The algorithm proceeds as follows: for each component k, it computes provisional scores V ← R U and loadings U ← Rᵀ V, then, if k > 0, it orthogonalizes the new vectors against all previously computed ones by subtracting their projections (U ← U − U (UᵀU) and V ← V − V (VᵀV)). After orthogonalization, the vectors are normalized, the eigenvalue λ is updated, and convergence is tested against a tolerance ε. This process guarantees perfect orthogonality of both scores and loadings, and when the number of components K equals the data rank N, GS‑PCA yields a full Singular Value Decomposition (SVD) of the original matrix with the same numerical accuracy as a conventional SVD routine.

Implementation-wise, the authors develop both CPU and GPU versions of NIPALS‑PCA and GS‑PCA. The CPU version uses the GNU Scientific Library (GSL) CBLAS interface, while the GPU version relies on NVIDIA’s CUDA‑based CUBLAS library. Both implementations exploit the same BLAS operations: Level‑1 (daxpy, dnrm2), Level‑2 (dgemv), and Level‑3 (dger). The GPU code allocates matrices on device memory, transfers data from host to device, performs the BLAS kernels in parallel across thousands of CUDA threads, and finally copies the results back to host memory. The authors note that double‑precision support requires a GPU of at least the GTX 280/Tesla C1060 class.

Performance experiments were conducted on an AMD Phenom 99 50 2.6 GHz CPU paired with an NVIDIA GTX 280 GPU. Random matrices of size M × N (M ranging from 5 × 10² to 1.5 × 10⁴, N = M/2) were used, with K = 10 components and convergence tolerance ε = 10⁻⁷. The results show that the GPU implementation scales dramatically better than the CPU version: for the largest problem (M = 1.5 × 10⁴) the GPU is roughly 12× faster. The GS‑PCA algorithm incurs only a modest overhead (5–7 % slower) compared with NIPALS, due to the additional orthogonalization steps, but this cost is negligible relative to the overall speedup obtained from GPU parallelism. Moreover, the authors observe that the GPU advantage becomes pronounced only when the problem size is large enough to keep thousands of threads busy, confirming the need for sufficiently big data to exploit massive parallelism.

In summary, the paper makes two key contributions: (1) it introduces GS‑PCA, an iterative PCA method that eliminates the orthogonality loss of NIPALS by integrating a stable Gram‑Schmidt re‑orthogonalization, and (2) it demonstrates that both NIPALS‑PCA and the new GS‑PCA can be efficiently accelerated on modern GPUs using CUBLAS, achieving up to a twelve‑fold speedup over optimized CPU implementations based on CBLAS. These findings are highly relevant for practitioners who need fast, numerically reliable PCA on large datasets in fields such as chemometrics, image analysis, and machine‑learning preprocessing.