Low-Rank Approximation by Randomly Pivoted LU

TL;DR

This work presents Randomly Pivoted LU (RPLU), a randomized complete-pivoting variant for low-rank matrix approximation that samples pivots proportional to squared residual entries. It proves geometric convergence in expectation when singular values decay geometrically with rate ρ<1/2, via a one-step Gram-bound that yields a 4^k residual growth factor, and introduces an efficient CUR-form, low-memory implementation suitable for large structured matrices and Cauchy-like matrices. The paper also introduces C2PLU as a competitive deterministic counterpart, demonstrates practical GPU-friendly, memory-light variants, and shows substantial speedups and high-rank capabilities in applications such as preconditioning for large linear systems and fast rational approximation. Collectively, these results enable high-quality, feature-preserving, low-rank approximations on hardware with limited memory and for matrices with exploitable structure, broadening the toolkit beyond traditional SVD-based methods.

Abstract

The low-rank approximation properties of Randomly Pivoted LU (RPLU), a variant of Gaussian elimination where pivots are sampled proportional to the squared entries of the Schur complement, are analyzed. It is shown that the RPLU iterates converge geometrically in expectation for matrices with rapidly decaying singular values. RPLU outperforms existing low-rank approximation algorithms in two settings: first, when memory is limited, RPLU can be implemented with $\mathcal{O}(k^2 + m + n)$ storage and $\mathcal{O}( k(m + n)+ k\mathcal{M}(\mat{A}) + k^3)$ operations, where $\mathcal{M}(\mat{A})$ is the cost of a matvec with $\mat{A}\in\mathbb{C}^{n\times m}$ or its adjoint, for a rank-$k$ approximation. Second, when the matrix and its Schur complements share exploitable structure, such as for Cauchy-like matrices. The efficacy of RPLU is illustrated with several examples, including applications in rational approximation and solving large linear systems on GPUs.

Figures (10)

• Figure 1: Performance of RPLU and three greedy schemes CPLU, C2PLU and IterativeCUR pritchard2025fast ($b=5$) on a Cauchy matrix $A_{ij} = 1 / (x_i - y_j)$ for points $\{x_i\}$ (blue) and $\{y_j\}$ (green) forming a smiley (top) and two spirals (bottom). Left: pivots selected by CPLU. Middle: pivots from one run of RPLU. Right: approximation error vs. rank $k$, with mean over 30 runs and min/max as a shaded region for RPLU and IterativeCUR. In the spiral example, greedy schemes select only "outlier" points, whereas RPLU also samples from the cluster, yielding lower error. Similar behavior is observed for PSD matrices rpchol.
• Figure 2: Comparison of low-rank approximation performance on kernel matrices $\mathbf{A}_{ij} = \exp(-0.5\|x_i-x_j\|^2)$ for two point sets shown on left. The nuclear-norm error is averaged over 30 runs, with standard deviation indicated by shaded regions.
• Figure 3: Illustration of the bounds in \ref{['theorem:rpchol_bound']} and \ref{['lem:rchol_doubling']} for RPCholesky and \ref{['theorem:doubling']} for RPLU on random $200 \times 200$ matrices with different singular value decay rates. Top row: RPCholesky with polynomial decay $\sigma_k \propto 1/k^2$ and geometric decay $\sigma_k = \rho^k$ for $\rho \in \{0.75, 0.4, 0.25\}$. Bottom row: RPLU on the same decay profiles. To compute the bound in \ref{['theorem:rpchol_bound']}, for each $k$ we perform a grid search over $\varepsilon$ and $r$ to find the values that minimize the bound. When the singular values decay geometrically with rate $\rho < \frac{1}{2}$, \ref{['lem:rchol_doubling']} and \ref{['theorem:doubling']} provide a loose bound with geometric decay for RPCholesky and RPLU, respectively, and no guarantees when $\rho \geq \frac{1}{2}$. In contrast, the RPCholesky bound in \ref{['theorem:rpchol_bound']} provides a sub-geometric bound in all cases. In all experiments with geometric decay, the observed residuals (mean over 30 trials with min/max bands) track closely with the optimal rank-$k$ approximation error, significantly outperforming all theoretical bounds.
• Figure 4: Approximation quality and timing of the various algorithms on a multi-level Toeplitz matrix with $128^3$ and $320^3$ grid points on a GPU. "OOM" denotes when the randomized SVD runs out of memory in each case. The $2$-norm is estimated using a Lanczos iteration using matrix-vector products with $\mathbf{A}$ and $\mathbf{\hat{A}}$.
• Figure 5: Convergence and timing of three GPU-based GMRES solvers for an integro-differential equation discretized on a $320^3$ grid ($n \approx 3.2 \times 10^7$). We compare unpreconditioned GMRES (black), and GMRES with a rank-$512$ CUR preconditioner computed by RPLU (blue), C2PLU (orange), CPQR (red) and RPQR (purple).

Theorems & Definitions (7)

• Lemma 1: Error doubling, Lemma 5.5 in rpchol
• Theorem 1: simplified bound, Theorem 2.3 in rpchol
• Theorem 2
• Proposition 1
• Proposition 2
• proof
• proof