[Paper Review] Randomized LU Decomposition

Gil Shabat, Yaniv Shmueli+2 more|arXiv (Cornell University)|Oct 27, 2013

Sparse and Compressive Sensing Techniques9 references14 citations

TL;DR

This paper introduces a randomized LU decomposition algorithm that efficiently computes low-rank matrix approximations using random projections, leveraging subgaussian random matrix theory to derive tight error bounds. The method achieves superior performance on GPUs with minimal CPU-GPU data transfer, outperforming randomized SVD and ID in speed while maintaining high approximation accuracy close to the best possible SVD results.

ABSTRACT

We present a fast randomized algorithm that computes a low rank LU decomposition. Our algorithm uses random projections type techniques to efficiently compute a low rank approximation of large matrices. The randomized LU algorithm can be parallelized and further accelerated by using sparse random matrices in its projection step. Several different error bounds are proven for the algorithm approximations. To prove these bounds, recent results from random matrix theory related to subgaussian matrices are used. As an application, we also show how the algorithm can be utilized to solve problems such as the rank-deficient least squares problem. Numerical examples, which illustrate the performance of the algorithm and compare it to other decomposition methods, are presented.

Motivation & Objective

• To develop a fast, parallelizable randomized LU decomposition algorithm for large-scale low-rank matrix approximation.
• To extend the randomized SVD framework to LU factorization, enabling efficient computation on modern hardware like GPUs.
• To derive rigorous error bounds based on random matrix theory, particularly for subgaussian matrices.
• To minimize computational cost and data transfer by fully parallelizing the algorithm for GPU execution.
• To empirically validate the algorithm’s accuracy and performance against existing methods like randomized SVD, ID, and Lanczos SVD.

Proposed method

• The algorithm uses random projection via a Gaussian matrix $ G $ to compress the input matrix $ A $, reducing its dimension while preserving its range.
• It applies a pivoted LU decomposition on the projected matrix $ AG $, ensuring numerical stability and rank-revealing properties.
• The method constructs low-rank approximations $ LU $ such that $ \|LU - PAQ\|_2 \leq C(m,n,k)\sigma_{k+1}(A) $, where $ \sigma_{k+1} $ is the $ (k+1) $-th singular value.
• Error bounds are derived using recent results on the extreme singular values of subgaussian random matrices, improving upon existing randomized SVD bounds.
• The algorithm is fully parallelized using standard GPU-accelerated BLAS operations, avoiding CPU-GPU data transfers for performance gains.
• Two variants are proposed: standard randomized LU (Algorithm 4.1) and a faster, sparser variant (Algorithm 4.4) using structured random projections.

Experimental results

Research questions

• RQ1 Can a randomized LU decomposition be designed to achieve low-rank approximation with error bounds comparable to randomized SVD?
• RQ2 How do the error bounds of the randomized LU method scale with matrix dimensions and rank, based on subgaussian random matrix theory?
• RQ3 To what extent can the algorithm be parallelized and accelerated on GPU architectures without CPU-GPU data transfer?
• RQ4 How does the performance (accuracy and speed) of the randomized LU compare to randomized SVD, ID, and Lanczos SVD on real-world and synthetic matrices?
• RQ5 What is the trade-off between approximation accuracy and computational efficiency in the fast randomized LU variant?

Key findings

• The randomized LU algorithm achieves a PSNR of approximately 32.5 dB on a 2124×7225 image at rank $ k=200 $, outperforming randomized ID and approaching the performance of the Lanczos SVD.
• The algorithm runs significantly faster than randomized SVD, ID, and Lanczos SVD on GPU, with execution times reduced by up to 30% compared to the next fastest method.
• The error of the fast randomized LU (Algorithm 4.4) is higher than that of the standard version (Algorithm 4.1), but both decay at the same rate with increasing $ k $.
• Theoretical error bounds are proportional to $ \sigma_{k+1}(A) $, confirming the method's reliability for low-rank approximation.
• The algorithm’s performance is stable across matrices with both rapidly and slowly decaying singular values, demonstrating robustness.
• The use of structured random projections in Algorithm 4.4 reduces memory and computation cost but requires a larger $ l $ to maintain accuracy, as shown in Lemma 4.13.