Low-Rank Kernel Approximation

• Low-Rank Kernel Approximation is a framework that compresses dense kernel matrices using analytic, algebraic, and hybrid methods to reduce computational complexity.
• The approach leverages explicit separable expansions and error guarantees to balance rank, accuracy, and dimensionality in high-dimensional RBF kernels.
• Practical strategies such as block clustering and combinatorial singular value analysis enable efficient algorithm implementations in large-scale scientific computing.

Low-rank kernel approximation is a broad framework encompassing analytic, algebraic, and hybrid methodologies for efficiently compressing large, dense kernel matrices and operators. The primary goal is to reduce computational and storage complexity while retaining sufficient accuracy for downstream scientific computing or statistical tasks. This article gives a comprehensive technical overview of key principles, explicit constructions, error guarantees, and modern algorithmic strategies, with an emphasis on radial basis function (RBF) and analytic kernels in high dimension, as detailed in "On the numerical rank of radial basis function kernels in high dimension" (Wang et al., 2017) and complementary research.

1. Separable Low-rank Expansions of Kernel Functions

Let $K(x,y) = f(\|x-y\|_2^2)$ be a radial basis function kernel with $x, y \in \mathbb{R}^d$. A low-rank kernel approximation seeks an expansion

$K(x, y) \approx \sum_{i=1}^R \phi_i(x)\, \psi_i(y)$

with functions $\phi_i, \psi_i: \mathbb{R}^d \to \mathbb{R}$ and minimal rank $R$ for a prescribed error.

When $f$ is analytic on $[0, D^2]$ and extends to a Bernstein ellipse in the complex domain, such expansions admit explicit analytic controls. In particular, a degree-$n$ expansion yields a rank

$R(n, d) = \binom{n + d + 2}{d + 2}$

and a uniform error bound in $L_\infty$ norm: $$\|K - \widetilde{K}\|_\infty \le \varepsilon_n(D) := \frac{2 C_D}{\rho_D^2 - 1}\, \rho_D^{-2n}$$ where $C_D$ bounds $|f|$ and $\rho_D^2 > 1$ is the Bernstein parameter.

For kernels $f$ of finite smoothness (i.e., $q$ derivatives with total variation $V_q$ on $[0,D^2]$), one still obtains algebraic convergence in $n$: $$\|K - \widetilde{K}\|_\infty \le \frac{2 V_q D^{2q}}{\pi q [2(n-q)]^q}, \quad n > q$$ with nearly the same combinatorial growth for $R(n,d)$.

2. Error-Rank-Dimension Trade-offs

The relation between rank, accuracy, smoothness, and dimension is central. For analytic $f$ and fixed $\varepsilon$,

$$R(\varepsilon, d) = \binom{d + 2 + \frac{-\ln(c_1\varepsilon)}{\ln\rho_D^2}}{d + 2}$$

which is polynomial in $d$ for fixed $\varepsilon$. More succinctly, for large $d$,

$$R(\varepsilon, d) = O\Bigl(d^{\,\frac{-\ln(c_1\varepsilon)}{\ln\rho_D^2}}\Bigr)$$

Thus, RBF kernels admit polynomial (rather than exponential) growth in rank as a function of $d$. In settings where $f$ has only finite smoothness $q$, the error satisfies

$\|K-K_R\|_\infty = O\Bigl(D^{2q} R^{-q/d}\Bigr)$

due to $R \sim n^d$, $n \sim R^{1/d}$, and error $\sim n^{-q}$. Consequently, higher smoothness and smaller domain diameter $D$ enable more rapid rank reduction.

For block-partitioned domains (Fourier–Taylor expansions), the error reflects both the smoothness and the geometry of the clusters: $$\|K-K_R\|_\infty \le f_\infty\left(\frac{D_x D_y}{D^2}\right)^{M_t+1} + \frac{V_q}{\pi q} \left(\frac{2D^2}{\pi M_f}\right)^q$$ where $D_x, D_y$ are source and target cluster diameters. Hence, reducing cluster diameters directly improves low-rank approximation error for a fixed rank.

3. Singular Value Plateaux and Group Structure

RBF kernel matrices exhibit distinctive spectral decay patterns. Rather than a simple exponential tail, the singular values $\sigma_i$ present plateaux separated by sharp drops at indices

$I_k = \sum_{\ell=0}^k n_\ell = \binom{k + d}{d}$

where $n_k = \binom{k + d - 1}{d - 1}$ is the count of separable $x$-$y$ terms in the $k$th order Taylor expansion. These plateaux reflect the grouping of separable basis functions arising from polynomial or Fourier–Taylor expansions and align with the practical behavior of SVD, Nyström, randomized SVD, and other decomposition methods. Larger drops represent the exhaustion of all combinations of a given polynomial degree, and the empirical singular value spectrum closely matches this combinatorial structure.

4. Practical Guidance: Algorithmic Strategies and Block Partitioning

Rank selection: For specified $\varepsilon$, select $n$ so that $\rho_D^{-2n}\lesssim \varepsilon$, then set $R\approx \binom{n + d + 2}{d + 2}$.

Block clustering: Clustering data into spatially localized blocks of small diameter $D_x$ or $D_y$ significantly reduces necessary rank, as guided by Fourier–Taylor error bounds. This approach is foundational in fast multipole methods, hierarchical matrices, and modern scalable kernel learning.

Algorithms and implementations:

• For analytic $f$, construct an explicit polynomial or Fourier–Taylor expansion truncated to degree $n$.
• For smooth non-analytic $f$, employ Taylor expansion up to order $n$ with $n$ chosen for the desired algebraic decay.
• For block-wise settings, exploit the geometric dimensions of the clusters to reduce storage and computational cost. Numerical experiments confirm that Monte Carlo–based algorithms (randomized SVD, Nyström) display drops in reconstruction error at thresholds $R = \binom{k+d}{d}$, while implementation simplicity and speed benefit from leveraging small-diameter clustering.

5. Implications and Empirical Verification in High Dimension

Despite the apparent curse of dimensionality (i.e., the $O(r^{2d})$ scaling for tensor-product bases), analytic RBF kernels allow for accurate low-rank separation with $R$ polynomial in $d$ at fixed $\varepsilon$. Empirical tests with up to $10^4$ points and $d=3\dots 20$ demonstrate:

• For a fixed $\varepsilon$, rank grows as $d^n$ when $\varepsilon\approx \rho_D^{-2n}$.
• Thresholds in numerical error decay occur at $R = \binom{k+d}{d}$, matching combinatorial predictions.
• Block clustering directly reduces observed required rank, matching the Fourier–Taylor theoretical bound.

The "group" structure in the spectrum justifies block partitioning and provides an explicit roadmap for rank allocation, hybrid expansions, and blockwise approximation in data-driven and scientific computing applications.

6. Recommendations for High-dimensional Kernel Approximation

• For RBF kernels on $d$-dimensional domains, always assess analyticity or smoothness of $f$ to set feasible rank-accuracy trade-offs.
• Use small-diameter clustering whenever possible to exploit geometric decay in blockwise low-rank error.
• Select $n$ (Taylor or Chebyshev order) based on the desired uniform (or operator) norm tolerance, leveraging $\rho_D^{-2n}\lesssim\varepsilon$ or $n^{-q}\lesssim\varepsilon$ depending on kernel regularity.
• Allocate rank in blocks according to combinatorial singular value plateau structure for maximal efficiency in hierarchical or matrix-free solvers.

These strategies lead to practical, efficient algorithms for large-scale numerical linear algebra, machine learning, Gaussian process regression, and PDE-control problems involving RBF and analytic kernel matrices.

---

In summary, the theory and practice of low-rank kernel approximation for analytic and RBF kernels is now sharply quantified: for fixed error, the function rank grows only polynomially in $d$. The plateaux and group patterns in the singular spectrum, explained by expansion combinatorics, directly inform blockwise algorithms and rank selection, with empirical results corroborating the theoretical predictions on high-dimensional and large-$n$ data (Wang et al., 2017).