Nyström Approximation Overview

• Nyström approximation is a method for discretizing linear operators by selecting landmark nodes to construct low-rank approximations with controlled error.
• It leverages matrix properties like positive semidefiniteness and smoothness, incorporating QR-based refinements and randomized numerical linear algebra techniques.
• This framework underpins scalable kernel learning, integral equation discretization, and operator approximations in both finite and infinite-dimensional settings.

The Nyström approximation is a unifying framework for discretizing and efficiently approximating linear operators, most notably large kernel matrices and integral operators of the second kind. It is central to scalable kernel learning, boundary- and integral-equation methods, and numerical analysis of both finite- and infinite-dimensional linear and nonlinear problems. The method exploits structure—typically positive semidefiniteness or smoothness—and judiciously selected "landmarks" or quadrature nodes to construct low-rank compressions or interpolants with quantifiable error. Its modern theory is shaped by matrix sketching, randomized numerical linear algebra, and operator theory, with algorithms now supporting generalized, highly-accurate, mixed-precision, and even infinite-dimensional settings.

1. Core Principles of the Nyström Approximation

The classical finite-dimensional Nyström approximation considers a symmetric positive semidefinite (SPSD) matrix $K\in\mathbb{R}^{n\times n}$, often a kernel (Gram) matrix in data analysis or scientific computing. Given a set of $l$ sampled "landmarks" or indices, the method forms:

• $C = K(:, \{i_j\}_{j=1}^l)\in\mathbb{R}^{n\times l}$: columns corresponding to the landmarks,
• $W=K(\{i_j\},\{i_j\})\in\mathbb{R}^{l\times l}$: the landmark submatrix.

The full-rank Nyström approximation is

$G = C\,W^{\dagger}\,C^T$

where $W^{\dagger}$ is the Moore–Penrose pseudoinverse, yielding at most rank $l$. For a desired fixed rank $k < l$, classical approaches truncate $W$ by its top $k$ spectral components, but this is suboptimal: the truncation should occur after reconstructing the projected matrix, not before. QR-based modifications optimally achieve this by (i) orthogonalizing $C$, (ii) projecting the core $W$, and (iii) performing an eigen-decomposition at the most informative stage (Pourkamali-Anaraki et al., 2017).

Beyond matrices, the Nyström idea underpins collocation methods for integral equations, where a continuous operator $K$ is approximated by discretizing integrals (often through quadrature rules), leading to systems of the form

$f(y) + \int_D k(x, y)f(x) d\mu(x) = g(y)$

with discretization replaced by weighted sums at interpolation nodes. The method generalizes further to infinite-dimensional settings, e.g., approximating Hilbert–Schmidt operators via random processes (Persson et al., 2024).

2. Mathematical Formulations and Algorithmic Variants

2.1 Standard and QR-Based Nyström Methods

Let $K$ be SPSD and choose $l$ landmarks ($l \geq k$). Define

$C = K(:, P),\quad W = P^T K P$

with $P$ the sampling matrix. The optimal fixed-rank Nyström via QR is constructed as: $C = QR,\quad R W^{\dagger} R^T = V \Sigma V^T$ Truncate $V, \Sigma$ to top $k$ and reconstruct

$$G_{(k)}^{\mathrm{opt}} = (QV_k)\,\Sigma_k\,(QV_k)^T$$

which minimizes the approximation error in trace norm among all Nyström approximants of rank $k$ (Pourkamali-Anaraki et al., 2017). Standard methods (truncating $W$ before projecting back) are generally suboptimal, both theoretically and empirically.

2.2 Generalized and Infinite-Dimensional Nyström

For general (possibly rectangular or indefinite) matrices $A\in\mathbb{R}^{m\times n}$, the "generalized Nyström" uses left and right sketches $S,T$: $\hat{A} = (A S) (T^* A S)^{\dagger} (T^* A)$ This encompasses CUR decompositions and admits extensions to cases where $A$ is only accessible through matrix-vector products (single-pass, mixed-precision, or streaming scenarios) (Lazzarino et al., 16 Jan 2026, Carson et al., 2022).

For integral operators and infinite-dimensional trace class or Hilbert–Schmidt operators, the randomized Nyström approach employs random function-valued sketches drawn from Gaussian processes and constructs operator-valued Nyström approximations with rigorous expectation/tail error bounds (Persson et al., 2024).

3. Error Analysis and Theoretical Guarantees

3.1 Error Bounds in the Finite Setting

Traditional Nyström error analysis (Frobenius norm) yields

$$\|K - \hat{K}_r\|_F \leq \|K-K_r\|_F + O\left(\frac{N}{m^{1/4}}\right)$$

for uniform random selection of columns (landmarks), where $N$ is matrix size and $m$ is the number of samples. When there is a large eigengap at rank $r$ (i.e., $\lambda_r - \lambda_{r+1}$ is large), the error improves to $O(N/m^{1/2})$ (Mahdavi et al., 2012). Under spectral incoherence and a $p$-power law decay of eigenvalues, the spectral-norm error further improves to $O(N/m^{p-1})$ (Jin et al., 2011). These results depend on concentration inequalities (e.g., Hilbert–Schmidt) and perturbation theory.

When sampling is i.i.d. from a probability measure, fine-grained spectral tail bounds become available. The expectation of the diagonal (pointwise) kernel error satisfies

$$\mathbb{E}\int_X (k(x,x) - k_{Z,s}(x,x)) d\mu(x) \leq \sum_{i>s} \lambda_i + O\left(\frac{(\log\ell)^{2d+1}}{\ell}\right)$$

for landmark set size $\ell$ and dimension $d$, improving on older rates $O(\ell^{-1/4})$ (Hayakawa et al., 2023).

3.2 Relative Error and Stability for Indefinite and Infinite-Dimensional Operators

For symmetric indefinite matrices, classical Nyström approximations can fail catastrophically due to cancellation in the core $W$ matrix. The remedy is to truncate the core to its top $r$ singular values prior to pseudoinversion, yielding relative-error nuclear norm guarantees if singular values decay rapidly (Park et al., 2022). For infinite-dimensional non-negative self-adjoint trace-class operators, the expected Hilbert–Schmidt, trace, and operator norm errors are controlled by the tail of the operator spectrum, with explicit structural constants (Persson et al., 2024).

4. Numerical Schemes for Integral and Nonlinear Equations

The Nyström method is also a workhorse for discretizing integral equations, both linear and nonlinear. In the solution of linear Fredholm equations of the second kind

$f(y) + (Kf)(y) = g(y)$

on intervals or domains, nodes and weights are chosen by Gauss or Lobatto quadrature, with collocation at quadrature nodes (Fermo et al., 2023, Fermo et al., 2013). The convergence rate is governed by the smoothness of $k$ and $g$.

For nonlinear integral equations (e.g., Hammerstein equations)

$$u(x) = g(x) + \int_{-1}^1 K_1(x,t)u(t)\,dt + \int_{-1}^1 K_2(x,t)h(t,u(t))\,dt$$

the discrete Nyström scheme applies quadrature and converts the continuous problem into a nonlinear system, solvable by Newton or fixed-point iterations. For weakly singular kernels, specialized product rules (Legendre, Mastroianni–Milovanović) are used for stability and higher-order convergence (Fermo et al., 2024).

Fredholm equations defined with weighted measures benefit from averaged or weighted–averaged quadrature rules, which provide improved stability and intrinsic a posteriori error estimation (Fermo et al., 2023).

5. High-Accuracy, Adaptive, and Modern Algorithmic Developments

Standard Nyström methods often stagnate at moderate accuracy, particularly for non-PSD, rectangular, or highly ill-conditioned problems. High-accuracy schemes employ alternating row and column skeletonization using strong rank-revealing (SRR) factorizations, pivoted subset expansion, Schur complement refinement, and progressive error estimation (Xia, 2023). These adaptive strategies can match SVD accuracy down to machine precision with modest sample sizes, achieving essentially optimal singular value decay and robustness to poorly conditioned matrices.

Single-pass and mixed-precision schemes, critical for streaming or ultra-large-scale matrices, perform the dominant $A$-sketch in lower precision (e.g., fp16) and the rest in higher precision, with careful shifts and stability analysis tailored to the approximation rank and matrix spectral properties (Carson et al., 2022). Error and conditioning bounds guide the safe selection of precisions and shifts.

For generalized Nyström approximations (including CUR-type and oblique extensions), rigorous and computationally efficient leave-one-out error estimators can be derived—such as LPO, LTO, and LRO—using only the already-computed sketches, obviating the need for extra matrix accesses (Lazzarino et al., 16 Jan 2026).

6. Practical Implementation: Algorithms and Guidelines

Implementation of Nyström approximations is context-dependent. For kernel matrices, the choice of landmarks is central: uniform random, leverage-score, determinantal point processes, and clustering all have distinct statistical properties and computational costs. Regularization (adding a ridge term to $W$ or core matrices) is essential in ill-conditioned or empirical-Mercer settings, especially with rapidly decaying spectra (Hayakawa et al., 2023). The landmark subset size typically scales as $O(s \log s)$ for target rank $s$.

When discretizing integral operators, careful node placement (e.g., Gauss, Lobatto, sigmoidal transformations) and weight construction (e.g., product quadrature for singular kernels) protect against instability and maximize convergence. Error estimation—using averaged/weighted rules or built-in residual/leave-one-out estimators—enables adaptive, high-accuracy solvers (Fermo et al., 2023, Fermo et al., 2024).

Pseudocode for core algorithms is standardized:

Given: Data matrix X, landmarks Z, kernel κ, target rank k 1. Form C_{ij} = κ(x_i, z_j), W_{ij} = κ(z_i, z_j) 2. Perform thin-QR: [Q, R] = qr(C) 3. Eigendecompose R * pinv(W) * R^T = V * Σ * V^T 4. Truncate to k: V_k = V(:,1:k), Σ_k = Σ(1:k,1:k) 5. Return U_k = Q*V_k, Λ_k = Σ_k

as in (Pourkamali-Anaraki et al., 2017), with extensions for rectangular, mixed-precision, or operator-valued settings as appropriate.

7. Applications and Empirical Performance

Nyström methods are foundational for scalable kernel learning, Gaussian process regression, kernel PCA, spectral clustering, and preconditioning in conjugate gradient methods. Adaptive and QR-based variants consistently outperform classical truncation schemes in both accuracy and robustness (Pourkamali-Anaraki et al., 2017, Park et al., 2022, Xia, 2023). In numerical studies, trace-norm, Frobenius-norm, and max-norm errors are uniformly smaller for QR/refined or high-accuracy Nyström, with clustering accuracy and matrix conditioning matching or closely tracking those of the full SVD in many cases. In machine learning, spectral and generalization error rates can be explicitly controlled as functions of landmark size, spectrum, and data coherence, with sublinear numbers of support vectors attainable for rapidly decaying spectra (Jin et al., 2011).

In integral and boundary-value contexts, adaptive Nyström discretizations yield highly accurate solutions for linear and nonlinear equations, including those with singularities or on irregular domains, and with rigorous control on inversion and error propagation (Fermo et al., 2013, Fermo et al., 2023, Fermo et al., 2024, Anand et al., 2017). For infinite-dimensional operator approximations, empirical errors decay in lockstep with the truncated singular value spectrum, as predicted by theory (Persson et al., 2024).

---

References (arXiv IDs)

• (Pourkamali-Anaraki et al., 2017) Improved Fixed-Rank Nyström Approximation via QR Decomposition: Practical and Theoretical Aspects
• (Mahdavi et al., 2012) An Improved Bound for the Nystrom Method for Large Eigengap
• (Jin et al., 2011) Improved Bound for the Nystrom's Method and its Application to Kernel Classification
• (Park et al., 2022) Randomized low-rank approximation for symmetric indefinite matrices
• (Xia, 2023) Making the Nyström method highly accurate for low-rank approximations
• (Persson et al., 2024) Randomized Nyström approximation of non-negative self-adjoint operators
• (Fermo et al., 2013) A Nystrom method for a boundary integral equation related to the Dirichlet problem on domains with corners
• (Lazzarino et al., 16 Jan 2026) Efficient error estimators for Generalized Nyström
• (Carson et al., 2022) Single-pass Nyström approximation in mixed precision
• (Fermo et al., 2024) A global approximation method for second-kind nonlinear integral equations
• (Fermo et al., 2023) Averaged Nyström interpolants for the solution of Fredholm integral equations of the second kind
• (Hayakawa et al., 2023) Sampling-based Nyström Approximation and Kernel Quadrature
• (Anand et al., 2017) A Nyström-based finite element method on polygonal elements