Dense High-Dimensional Compressed Kernels

• The paper introduces a near-linear time and space compression method via incomplete Cholesky factorization for dense high-dimensional kernel matrices.
• It demonstrates that truncated Cholesky factorization achieves controlled operator-norm error, enabling accurate low-rank approximations and robust PCA.
• Applications include scalable elliptic PDE solvers, Gaussian process inference, and compressed subspace matching in signal processing.

A dense continuum of high-dimensional compressed kernels refers to both the structure and algorithmic manipulation of families of kernels—often arising as dense matrices from elliptic partial differential equations, Gaussian process covariance functions, or parametric signal subspaces—such that their essential properties and computational utility can be compressed, matched, or approximated efficiently in high dimensions. This concept is pivotal in large-scale numerical analysis, machine learning, and signal processing, especially where direct manipulation of $N \times N$ dense matrices is computationally infeasible. Recent developments provide rigorous frameworks for near-linear time and space compression via incomplete Cholesky factorization, as well as methods for compressed matching in families of high-dimensional subspaces (Schäfer et al., 2017, Mantzel et al., 2014).

1. Dense Kernel Matrices in High Dimensions

Dense kernel matrices $\Theta \in \mathbb{R}^{N \times N}$ arise from point evaluations of a symmetric positive-definite kernel $G(\cdot, \cdot)$ at locations $\{ x_{i} \}_{i = 1}^N \subset \mathbb{R}^d$:

$\Theta_{ij} = G(x_i, x_j).$

For $G$ as the Green’s function of an elliptic operator $L$ over a domain $\Omega \subset \mathbb{R}^d$, such matrices also represent the discretized covariance of spatially indexed Gaussian processes. The growth of $N$ quickly makes direct storage, inversion, or eigendecomposition intractable for both dense PDE solvers and statistical inference (Schäfer et al., 2017).

Complementarily, one considers a continuum of parametric kernel-induced subspaces $$\{S_\theta \subset \mathbb{R}^N: \theta \in \Theta \subset \mathbb{R}^D\}$$, where the index space is continuous (e.g., shift, frequency, or scale parameters). Problems such as template matching and source localization are naturally posed in such kernel families (Mantzel et al., 2014).

2. Sparse Compression via Incomplete Cholesky Factorization

A breakthrough in matrix compression is provided by the zero-fill incomplete Cholesky factorization (ICHOL(0)) tailored to a sparsity pattern $S_p \subset \{1,\dots,N\}^2$. The index set $S_p$ is defined in terms of a “maximin ordering” of the $x_i$ (prioritizing points furthest from the boundary or previous points) and a fill distance $[k]$ at stage $k$:

$$S_p = \left\{ (i,j): \operatorname{dist}(x_i, x_j) \le p \max\{[i], [j]\} \right\}$$

for a parameter $p \sim O(\log(N/\epsilon))$ tuned to the desired accuracy $\epsilon$.

The algorithm proceeds by iteratively computing the Cholesky decomposition of $\Theta$ but fills only positions in $S_p$, discarding candidate fill-ins outside the prescribed sparsity. Theoretical results demonstrate that for elliptic Green’s function kernels, the off-diagonal entries of the full Cholesky factor decay exponentially in a hierarchical pseudo-metric, so truncation to $S_p$ incurs only controlled operator-norm error:

$$\|\Theta - L_{S_p} L_{S_p}^T\|_2 \le \epsilon \|\Theta\|_2$$

provided $p \gtrsim C' \log(N/\epsilon)$ (Schäfer et al., 2017).

3. Continuum of Compressed Kernels and Complexity Guarantees

Varying $\epsilon$ yields a dense continuum of compressed kernel approximations $$\{\widetilde{\Theta}(\epsilon) = L(\epsilon)L(\epsilon)^T\}_{\epsilon > 0}$$, where each factorization achieves error at most $\epsilon \|\Theta\|$ and is supported on $\# S_p = O(N \log N \log^d(N/\epsilon))$ entries:

$$\text{space} = O(N \log N \log^d(N/\epsilon)), \qquad \text{time} = O(N \log^2 N \log^{2d}(N/\epsilon)).$$

The sparsity pattern adapts easily to the intrinsic rather than ambient dimension for points lying on a low-dimensional manifold within $\mathbb{R}^d$. Critically, the algorithm requires only the spatial configuration of $\{x_i\}$ and never needs analytic kernel expressions, making it broadly applicable (Schäfer et al., 2017).

4. Approximate PCA and Low-Rank Structure

The truncated Cholesky factorization directly furnishes a low-rank approximation suitable for sparse principal component analysis:

$\Theta \approx L_{:,1:k} L_{:,1:k}^T,$

where $L_{:,1:k}$ denotes the first $k$ columns of $L$. The quality of the rank-$k$ approximation is tightly controlled; if $v_{k}$ is the $k$th pivot radius (local mesh size), the operator norm satisfies

$$\|\Theta - L_{:,1:k} L_{:,1:k}^T\|_2 \le C v_{k+1}^{2-d}$$

for some constant $C$ (Schäfer et al., 2017). This ensures that near-optimal PCA (in operator norm) is achieved at near-linear computational cost.

5. Inversion and Fast Elliptic PDE Solvers

Fast inversion exploits the structure of both $\Theta$ and its Cholesky factors. The reverse-order Cholesky factorization provides a similarly sparse factorization for $\Theta^{-1}$, supporting direct elliptic PDE solvers:

$$\text{Solve:} \quad \Theta u = f \quad \text{via} \quad \Theta^{-1} = L_{\mathrm{rev}} L_{\mathrm{rev}}^T.$$

The space and time complexities remain $O(N \log N \log^d(N/\epsilon))$ and $O(N \log^2 N \log^{2d}(N/\epsilon))$, matching those of the forward analysis, with rigorously controlled error approximation (Schäfer et al., 2017).

6. Matching over High-Dimensional Kernel Continua via Compression

In the parallel context of compressed subspace matching (Mantzel et al., 2014), a collection of $K$-dimensional subspaces $\{S_\theta\}$ over a continuous parameter space $\Theta$ can be matched to high-dimensional observed signals from a small number of random projections. The collection’s geometric complexity, as quantified by covering numbers and a scalar $\Delta$, determines the requisite number of measurements $M$, with

$M = O\left(K (\Delta + \log K) \right)$

guaranteeing uniform preservation of the matching error landscape:

$$\Big|\|\tilde P_\theta y\|_2^2 - \|P_\theta h_0\|_2^2\Big| \le O\left( \sqrt{ \frac{K(\Delta+\log K)}{M} } \right)$$

hold with high probability for all $\theta$ (Mantzel et al., 2014). This framework applies directly to compressed-domain template matching, time-of-arrival estimation, and matched field processing, even when the template dictionary is uncountably infinite.

7. Applications and Implications

The described frameworks yield efficient, scalable tools for manipulating dense kernel matrices and high-dimensional subspace families:

• Elliptic PDE Solvers: Near-linear time direct solvers for large-scale elliptic boundary value problems.
• Gaussian Process Inference: Scalable PCA, eigendecomposition, and inversion for spatial GP models with Green's function covariance.
• Signal Processing and Source Localization: Accurate template matching and subspace identification from compressed random sketches, robust to the dictionary cardinality and ambient dimension.
• Manifold Adaptivity: Computational costs scale with the intrinsic (not ambient) dimension when data is concentrated according to low-dimensional manifolds.

These results rigorously establish that a dense continuum of compressed kernel approximations can support essential analysis primitives—compression, inversion, PCA, and matching—at computational costs near-linear in $N$ and robust to infinite subspace families, provided the underlying geometric complexity is controlled (Schäfer et al., 2017, Mantzel et al., 2014).