KernelForge

World-class random-feature kernel approximation — a high-precision system for scaling Gaussian/RBF kernels to millions of points.

Author: 晨星 (CJX0712) · License: MIT · Version: 0.1.0

KernelForge provides a family of kernel-feature maps that turn the expensive O(N²) kernel Gram matrix into an O(N·m) explicit feature inner product — without sacrificing accuracy. Its flagship method, KernelFuse, consistently beats the classical Nyström baselines at equal feature budget on every benchmark cell.

---

Why KernelForge

Scaling a kernel method (KRR, SVM, GPR, spectral clustering, MMD, …) to large N is bottlenecked by the N×N Gram matrix. The standard fixes are:

• Random Fourier / ORF / QMC features — Monte-Carlo quadrature of Bochner's integral. Easy, fast, but high-variance at small budgets.
• Nyström — sample landmark columns. Accurate, but a slow O(N·m²) projection and well-known rank-deficiency tails.

KernelFuse takes a third route: a residual-aware pivoted Cholesky (RPCholesky) over-sampling pass builds a wide, well-conditioned basis, then an Eckart–Young SVD compression distills it to exactly m features that capture the top-m spectral energy. The result is a feature map that, at budget m, reproduces the kernel more faithfully than Nyström — verified by independent, machine-precision invariants.

---

Methods

Method	Family	Idea
rff	Monte-Carlo	Random Fourier Features (Rahimi & Recht, 2007)
orf	Monte-Carlo	Orthogonal Random Features (Yu et al., 2016)
qmc	Monte-Carlo	Quasi-Monte-Carlo (Halton/Sobol) spectral sampling
nystrom_uniform	Landmark	Uniform-column Nyström
nystrom_leverage	Landmark	Leverage-score (ridge) Nyström
fuse	Hybrid	RPCholesky over-sample + optimal SVD compression (flagship)

All Monte-Carlo families share one convention: a cos/sin pair with weight w_j produces z_j(x) = √(2w_j/m)·[cos(ω_jᵀx), sin(ω_jᵀx)], so ⟨z(x),z(y)⟩ = Σⱼ (2wⱼ/m) cos(ω_jᵀ(x−y)) and ‖z(x)‖² is exactly constant (independent of x) — a property the invariant suite checks to 1e-16.

---

Benchmarks

Measured on 216 runs (6 datasets × 2 budgets × 3 seeds × 6 methods), metric = relative Frobenius error of the test Gram matrix ‖ẐẐᵀ − K‖ / ‖K‖ (lower is better). KernelFuse wins 12/12 cells (win rate = 100%).

Mean relative Gram error by method (all cells):

method	mean rel Gram err
fuse (flagship)	0.0176
nystrom_leverage	0.0264
nystrom_uniform	0.0271
orf	0.1140
qmc	0.1338
rff	0.1640

KernelFuse vs Nyström-leverage at equal budget (speedup = nys_err / fuse_err):

dataset	fuse @64	nys @64	speedup	fuse @256	nys @256	speedup
mixture	0.0003	0.0015	4.80×	0.0000	0.0000	1.32×
lowrank	0.0068	0.0151	2.23×	0.0016	0.0027	1.68×
aniso_gaussian	0.0117	0.0256	2.19×	0.0040	0.0058	1.43×
iso_gaussian	0.0424	0.0628	1.48×	0.0163	0.0227	1.39×
correlated_gaussian	0.0313	0.0592	1.89×	0.0145	0.0190	1.31×
heavytail	0.0514	0.0666	1.30×	0.0313	0.0353	1.13×

(Full reproducible results: benchmark.json. Quick smoke results: benchmark_quick.json.)

---

Installation

pip install -e .
# or, without an editable install:
pip install numpy scipy scikit-learn   # then use the package from this directory

Requires Python ≥ 3.11 and NumPy ≥ 2.1, SciPy ≥ 1.14, scikit-learn ≥ 1.5.

---

Quick start

import numpy as np
from kernelforge import (
    make_dataset, KernelSpec, rbf_kernel,
    build_approximator, FLAGSHIP, relative_gram_error,
)

# 1. A data-generating process + an exact RBF kernel spec.
data = make_dataset("iso_gaussian", seed=11, n_train=600, n_test=400, dim=16)
ls = np.asarray(data.notes["lengthscale"], dtype=float)
spec = KernelSpec(lengthscale=ls)

# 2. Build the flagship feature map at a small budget (m = 64 features).
app = build_approximator(FLAGSHIP, n_features=64, seed=7).fit(data, spec)

# 3. Transform train/test; inner products approximate the exact kernel.
Ztr, Zte = app.transform_train_test(data)
K_te_true = rbf_kernel(data.X_test, lengthscale=spec.lengthscale)
err = relative_gram_error(Zte, K_te_true)   # ~0.01–0.05
print(f"{FLAGSHIP} test Gram error: {err:.4f}")

# 4. Downstream kernel ridge regression vs the exact-kernel oracle.
from kernelforge import KernelRidgeRegression, exact_kernel_oracle

Command-line benchmark

python -m kernelforge --out benchmark.json      # full benchmark
python -m kernelforge --quick                   # fast smoke run
python -m kernelforge --method fuse --dataset lowrank --budget 64
python run_demo.py                              # 18-run visual demo

---

Verifiable invariants (the hard gold standard)

Every approximation has at least one invariant that can be cross-checked independently of the implementation that produced it. The suite (kernelforge.invariants) is a CI gate:

Invariant	What it proves
exact_symmetry / exact_diagonal_in_unit	The analytic RBF kernel is symmetric and k(x,x)=1.
pair_diagonal_constant	‖z(x)‖² is identical for every x (pair convention).
rff_unbiasedness	E[⟨z,z'⟩] = k — MC average converges to the analytic Bochner form.
fuse_compression_optimal	ΦΦᵀ equals the top-m Eckart–Young truncation of the wide basis.
gram_symmetry / gram_psd	Every produced Gram is symmetric positive-semidefinite.
determinism	Same seed ⇒ bit-identical features.

from kernelforge import run_invariants, build_approximator
rec = run_invariants("fuse", data, spec, n_features=64)
assert all(b["passed"] for b in rec.values())

---

Architecture

kernelforge/
├── core/            # types (Dataset, KernelSpec), seed spawning, config, errors
├── kernels/
│   └── exact.py     # analytic RBF (Bochner) kernel — Tier-0 ground truth
├── data/
│   └── dgps.py      # data-generating processes + isotropic median bandwidth
├── features/
│   ├── base.py      # cos/sin pair convention + QMC/random primitives
│   ├── rff.py orf.py qmc.py nystrom.py
│   ├── fuse.py      # ★ KernelFuse: RPCholesky over-sample + SVD compression
│   └── approximator.py  # registry, build_approximator, FLAGSHIP
├── metrics.py       # relative_gram_error, nmse, r2, ...
├── invariants.py    # the verifiable battery (CI gate)
├── krr.py           # dual Kernel Ridge Regression + exact oracle
├── pipeline.py      # run_one / run / aggregate (benchmark engine)
└── benchmark.py     # CLI + summary table

KernelFuse in one picture: RPCholesky over-samples W = oversample·m columns by always pivoting on the largest residual diagonal (so each new column captures the most remaining spectral energy), then SVD-compresses the wide basis X_sampled @ C_wide to exactly m features via C = Vₘ (the top-m right singular vectors). This is provably the optimal rank-m approximation of the wide Gram by the Eckart–Young theorem.

---

Development

pip install -e ".[dev]"    # ruff + pytest
python -m ruff check .     # lint (must be clean)
python -m pytest -q        # 26 tests, all green
python -m kernelforge --quick

---

License

MIT © 晨星 (CJX0712). See LICENSE.