Kernel Gram Matrices: Analysis & Algorithms

Updated 8 February 2026

Kernel Gram matrices are defined by positive-definite kernels that capture pairwise similarities, forming the basis for interpolation and spectral analysis in various applications.
They enable efficient computation in high-dimensional settings through structural compression, randomized algorithms, and sparse algebraic operations.
Advances such as subquadratic algorithms, deep kernel processes, and asymptotic spectral analyses enhance scalability and practical performance on large datasets.

A kernel Gram matrix is a fundamental object in mathematical analysis, statistics, and machine learning, encoding pairwise similarities between a discrete set of points according to a symmetric, positive (semi-)definite function known as a kernel. The matrix’s spectral, algebraic, and structural properties govern the efficacy and scalability of numerous algorithms across numerical linear algebra, high-dimensional statistics, and spatial statistics. Modern research into kernel Gram matrices focuses both on their asymptotic spectral properties and on enabling efficient computations through structural compression, randomized algorithms, or algebraic factorization.

1. Formal Definition and Fundamental Properties

Given a collection of points $\{x_1,\dots,x_N\}\subset \Omega\subset\mathbb{R}^d$ and a symmetric, positive-definite kernel $k:\Omega\times\Omega\to\mathbb{R}$ , the kernel Gram matrix $K\in\mathbb{R}^{N\times N}$ is defined by

$K_{ij} = k(x_i, x_j).$

Key properties include:

Positive semidefiniteness: For all $v\in\mathbb{R}^N$ , $v^T K v\geq 0$ .
Mercer’s theorem: For continuous, positive-definite $k$ , there exist an orthonormal basis $\{\phi_\ell\}$ for $L^2(\Omega)$ and non-negative $\{\lambda_\ell\}$ such that $k(x,y) = \sum_{\ell=1}^\infty \lambda_\ell \phi_\ell(x)\phi_\ell(y)$ .
Spectral decomposition: $K=V\Lambda V^T$ with $\Lambda = \operatorname{diag}(\lambda_1, ..., \lambda_N)$ and orthonormal $V$ , yielding a spectrum $\{\lambda_j\}_{j=1}^N$ (Bakshi et al., 2022).
Representer theorem: In kernel-based learning, solutions to regularized empirical risk minimization problems belong to the span of the kernel functions evaluated at the data points.

A significant class of Gram matrices arises from random or structured point sets and various smoothness classes of $k$ , including radial and inner-product kernels. In reproducing kernel Hilbert spaces (RKHS), the Gram matrix records inner products between evaluation functionals, making it central in interpolation, approximation, and variance analysis (Spitzer, 2024).

2. Multiresolution and S-Format Sparse Algebra

High-dimensional and large-sample settings lead to kernel matrices that are prohibitively costly ( $O(N^2)$ memory, $O(N^3)$ arithmetic) in direct computation. The multiresolution samplet framework achieves optimal sparse representations and nearly linear scaling (Harbrecht et al., 2022).

Construction:

Build a balanced binary cluster tree $\mathcal{T}$ over $X$ .
At each level $j$ , define scaling functions $\phi_{j,k}$ (supported on clusters) and detail (samplet) functions $\sigma_{j,k}$ (orthogonal complements with vanishing moments).
The orthonormal samplet basis $\Sigma$ allows a change of basis $K^\Sigma = \Psi K \Psi^T$ .
Partitioned into blocks by cluster association, blocks between "far" clusters (well-separated) are set to zero (S-compression).

Algebraic operations in S-format:

Addition: Entrywise on the sparsity pattern; $O(N\log N)$ .
Multiplication: Only on nonzero pattern; $O(N\log^2 N)$ .
Inversion: Via selected inversion on the sparsified Cholesky or LDL $^T$ factor, leveraging the closure of S-format under inversion.
Matrix functions: Holomorphic functional calculus (e.g., $K^\alpha$ , $\exp(K)$ ) via contour integration and repeated sparse inversion or multiplication.

Theoretical guarantees:

For kernels of smoothness $q+1$ , S-format achieves $\|K^\Sigma-K^\eta\|_F \leq C\eta^{-2(q+1)}N\sqrt{\log N}$ and $O(N\log N)$ storage.
Compression, algebraic operations, and inverses maintain this computational scale.
Application to Gaussian process learning reduces computational complexity from $O(N^3)$ to $\widetilde{O}(N\log^{1.5}N)$ , with empirical errors $\lesssim 10^{-5}$ for F-norm on datasets up to $N\sim 10^6$ (Harbrecht et al., 2022).

3. Subquadratic Algorithms, KDE Reductions, and Randomized Methods

Recent techniques leverage kernel density estimation (KDE) data structures to avoid explicit Gram matrix formation, enabling subquadratic or even sublinear algorithms for core tasks (Bakshi et al., 2022, Backurs et al., 2021).

Kernel Density Estimation framework:

After $O(Nd/(\epsilon^2\tau^p))$ pre-processing, KDE queries allow fast estimation of row/column sums.
Reductions enable:
- Weighted vertex/edge sampling (graph sparsification, random walks).
- Importance sampling for low-rank approximation (using norm of rows).
- Approximation of spectral quantities (top eigenvalue/vector) via noisy power method.

Complexity and applications:

Spectral sparsification: Construction of an $m=O(N\log N/(\epsilon^2\tau^3))$ -edge graph approximating the Laplacian in near-linear time.
Low-rank approximation: Additive-error sketching in $O(N)$ KDE queries plus nearly linear algebraic overhead.
Empirical studies demonstrate $9\times$ to $41\times$ reductions in evaluations or storage over standard baselines (Bakshi et al., 2022).
Sublinear or subquadratic routines for trace, top eigenvalues/vectors, and dense-matrix operations are possible for kernels admitting fast KDE with $p<1$ (Backurs et al., 2021).

4. Spectral and Limit Theorems for Kernel Gram Matrices

The spectral properties of kernel Gram matrices are central in understanding their limiting behavior, algorithmic stability, and statistical concentration.

Asymptotic expansions and perturbation theory:

In settings where points are i.i.d. from a probability law $P$ , the empirical Gram matrix $K_n$ approximates the spectrum of the associated integral operator $T$ on $L^2(P)$ (Bae et al., 1 Feb 2026).
First-order expansions for eigenvalues and projections:

$\hat\lambda_{i,n} = \lambda_i + \langle(H_n-T)\phi_i,\phi_i\rangle + O(\|H_n-H\|_{\mathrm{op}}^2),$

where $H_n$ is the empirical kernel covariance operator; $\|H_n-H\|_{\mathrm{op}} = O(n^{-1/2})$ with high probability.

Weak convergence to Gaussian limits for eigenvalue fluctuations and projections holds under only Mercer-type conditions, with minimal regularity requirements.

Flat limit and eigenstructures:

As kernel width shrinks ( $\varepsilon\to 0$ ), spectra concentrate: the largest eigenvalue scales as $O(1)$ , the rest decay as powers of $\varepsilon^{2(s-1)}$ (infinitely smooth case), or $\varepsilon^{2r-1}$ for finite smoothness of order $r$ (Barthelmé et al., 2019).
Limiting eigenvectors correspond to discrete orthogonal polynomials (analytic kernel) or to the eigenspace of projected distance matrices for finite smoothness.
Determinant asymptotics and singularity structure are characterized via Vandermonde and Wronskian determinants, relevant for interpolation and preconditioning.

Random kernel matrix spectrum:

For Gram matrices from i.i.d. Gaussian vectors and inner-product kernels $f(\langle X_i, X_j\rangle)$ , the empirical spectral distribution converges to a canonical law governed by a cubic equation for the Stieltjes transform.
Smooth kernels generically yield the Marčenko–Pastur law; non-smooth kernels admit distinct (non-MP) limiting spectral distributions (Cheng et al., 2012).

5. Matrix Factorizations, Gram Matrices in Algebraic and Combinatorial Settings

Gram matrices associated with specific algebraic or combinatorial structures display factorization properties leading to explicit minors, inverses, and extremal entries.

Orthogonal polynomial bases:

In classical orthogonal polynomial systems, the Gram matrix in the polynomial basis is diagonal; under basis change (e.g., to monomials), Davis’s lemma provides an explicit factorization and inversion (Spitzer, 2024).
Exact rational forms for the inverse are obtainable for the Legendre, Hermite, and Laguerre cases, with interpolation kernels constructed as quadratic forms in $G^{-1}$ .

Number-theoretic and graph-theoretic kernels:

The GCD Gram matrix $S_{ij} = \gcd(x_i, x_j)$ is positive-definite, with total nonnegativity characterized by the monotonicity of exponents in the prime factorizations of $\{x_i\}$ . For such TN cases, $S$ is a Green’s matrix and admits tridiagonal inverses and explicit minor formulas (Guillot et al., 2019).
Over graphs, the Gram matrix of a RKHS induced by the Laplacian admits a decomposition into a sum involving the inverse of the non-zero Laplacian spectrum. Spectral and entrywise operators are bounded by combinatorial extremal graphs (e.g., path, star, complete), controlling pointwise self-similarity and smoothness (Seto et al., 2012).

6. Advanced Constructions: Deep Kernel Processes and Stochastic Gram Matrix Dynamics

Beyond fixed kernel matrices, modern stochastic process models operate directly on sequences of Gram matrices via deterministic kernel transforms and stochastic sampling.

Deep kernel processes (DKPs):

DKPs iterate kernel nonlinearities and Wishart/inverse-Wishart sampling layers, yielding a Markov chain on positive definite matrices.
Standard deep Gaussian processes (DGPs), infinite Bayesian neural networks, and deep inverse Wishart processes are all special cases of DKPs constructed as chains of Gram matrices (Aitchison et al., 2020).
Variational inference can be carried out directly on the space of Gram matrices (rather than feature vectors), with tractable evidence lower bound calculations via the conjugacy of the inverse-Wishart posterior.
The chain-of-Gram-matrix approach enables fully kernel-based, permutation- and rotation-invariant models, and empirical work reports improved predictive performance and evidence lower bounds over standard DGPs on MNIST and related datasets.

7. Applications, Limitations, and Frontiers

Kernel Gram matrices are core to Gaussian process regression, manifold learning, spectral clustering, semi-supervised learning on graphs, and more. Compressed representations such as the S-format, fast KDE-based reductions, and algebraic factorization methods overcome the $O(N^2)$ barrier, enabling work with $N\sim 10^6$ in practice (Harbrecht et al., 2022, Bakshi et al., 2022). Asymptotic expansions guide regularization and subspace estimation in PCA, while extremal and factorization results support explicit solution-formulas in arithmetic and combinatorial settings.

Limitations include the need for kernel-dependent data structures (for subquadratic algorithms), the dependence of errors and runtime on kernel parameter lower bounds (e.g., $\tau$ ), and the exponential growth of condition numbers in certain Gram matrices (notably in monomial bases), mandating careful numerical handling (Spitzer, 2024). Extending these results to adaptive, streaming, or higher-dimensional kernels remains a subject of ongoing research.