[Paper Review] Parallel GPU Implementation of Iterative PCA Algorithms

M. Andrecut|ArXiv.org|Nov 7, 2008

Neural Networks and Applications1 references4 citations

TL;DR

This paper proposes a GPU-accelerated, orthogonal iterative PCA algorithm (GS-PCA) based on Gram-Schmidt re-orthogonalization to overcome the loss of orthogonality in the standard NIPALS-PCA algorithm. Implemented using CUBLAS on NVIDIA GPUs, GS-PCA achieves up to 12× speedup over CPU-optimized CBLAS implementations while ensuring numerical stability for computing full sets of principal components.

ABSTRACT

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).

Motivation & Objective

• To address the loss of orthogonality in iterative PCA algorithms like NIPALS-PCA, which limits their use to only the first few components.
• To develop a numerically stable iterative PCA algorithm that maintains orthogonality across all components through Gram-Schmidt re-orthogonalization.
• To accelerate PCA computation for large-scale datasets using GPU parallelization via CUBLAS.
• To compare the performance and stability of GPU-optimized GS-PCA and NIPALS-PCA against CPU-optimized CBLAS implementations.

Proposed method

• Proposes GS-PCA, an iterative PCA algorithm that applies Gram-Schmidt orthogonalization at each component extraction step to preserve orthogonality.
• Uses a two-step iterative process: first estimating loadings via matrix-vector multiplication, then updating scores using orthogonal projection.
• Applies re-orthogonalization by subtracting projections onto previously computed components during each iteration to maintain orthogonality.
• Implements the algorithm on GPU using CUBLAS for high-performance linear algebra operations, including GEMV, GEMM, NRM2, and SCAL.
• Performs mean centering of input data on GPU using CUBLAS DAXPY and DCOPY operations before PCA computation.
• Transfers results between host and device using cublasGetMatrix and cublasSetMatrix, ensuring data consistency across CPU and GPU memory.

Experimental results

Research questions

• RQ1 Can a GPU-parallelized iterative PCA algorithm maintain orthogonality across all components, unlike standard NIPALS-PCA?
• RQ2 How much faster is a GPU-optimized GS-PCA implementation compared to a CPU-optimized CBLAS version?
• RQ3 Does the use of Gram-Schmidt re-orthogonalization eliminate the numerical instability caused by accumulated errors in iterative PCA?
• RQ4 What is the maximum speedup achievable using CUBLAS-based GPU kernels for PCA on large datasets?
• RQ5 Can GS-PCA reliably compute the full set of principal components without loss of accuracy or orthogonality?

Key findings

• The GPU-optimized GS-PCA implementation achieved up to 12 times speedup over the CPU-optimized CBLAS version for large datasets.
• GS-PCA successfully maintains orthogonality of principal components across all iterations, eliminating the numerical instability of NIPALS-PCA.
• The algorithm correctly computes the full set of principal components, not just the first few, due to the orthogonalization process.
• The use of CUBLAS on NVIDIA GPUs significantly accelerates key linear algebra operations such as matrix-vector multiplication and norm computation.
• The convergence criterion based on eigenvalue difference (|λ′ − λ| < ε) was effective in detecting stable component estimation in both NIPALS-PCA and GS-PCA.
• The implementation correctly reconstructs the data matrix as X = TP^T + R, with residuals R close to zero after convergence, confirming algorithmic correctness.