Numerical Methods for Kernel Slicing

Basic Information

• Paper ID: 2510.11478
• Title: Numerical Methods for Kernel Slicing
• Authors: Nicolaj Rux (Chemnitz University of Technology), Johannes Hertrich (Université Paris Dauphine-PSL and Inria Mokaplan), Sebastian Neumayer (Chemnitz University of Technology)
• Classification: math.NA, cs.NA
• Publication Date: October 14, 2025
• Paper Link: https://arxiv.org/abs/2510.11478v1

Abstract

Kernel functions are crucial for modeling interactions in machine learning. However, brute-force computation of kernel sums exhibits quadratic complexity growth with respect to sample size. Recent Fourier slicing methods can reduce this complexity to linear, provided that the kernel can be sliced and its Fourier coefficients are known. To obtain these coefficients, this paper formulates the slicing relationship as an inverse problem and proposes two recovery algorithms. Extensive numerical experiments demonstrate the speed and accuracy of the proposed methods.

Research Background and Motivation

Core Problem

Kernel methods are widely applied in machine learning for density estimation, support vector machine classification, principal component analysis, maximum mean discrepancy (MMD), and other tasks. The computational bottleneck in these applications typically involves evaluating expressions of the form:

$s_m := \sum_{n=1}^N F(\|x_n - y_m\|)w_n, \quad m = 1,\ldots,M$

where $F \in C([0,\infty))$ is a radial basis function, $x_1,\ldots,x_N, y_1,\ldots,y_M \in \mathbb{R}^d$ are sample points, and $w \in \mathbb{R}^N$ are weights.

Computational Complexity Challenges

Direct computation requires $O(NMd)$ operations, which is infeasible for large datasets. Classical methods such as fast Fourier summation and fast multipole methods, while reducing complexity to $O(M+N)$, suffer from exponential dependence on dimension $d > 4$ due to reliance on fast Fourier transforms or spatial partitioning, rendering them impractical.

Advantages of Slicing Algorithms

The fundamental idea of slicing algorithms is to find a function $f \in L^1_{loc}([0,\infty))$ such that:

$F(\|x\|) = \frac{1}{\omega_{d-1}} \int_{S^{d-1}} f(|\langle\xi, x\rangle|)d\xi$

where $\omega_{d-1} = 2\pi^{d/2}/\Gamma(d/2)$ is the surface measure of the $d$-dimensional sphere. By discretizing the integral, kernel summation can be reduced to one-dimensional cases and computed efficiently using fast Fourier summation.

Core Contributions

1. Formalization of slicing function recovery as an inverse problem, establishing a complete theoretical framework
2. Proposal of two numerical algorithms for recovering cosine series coefficients required for fast Fourier summation
3. Provision of rigorous error estimates, including analysis of forward error and slicing error
4. Extensive numerical experiments validating the efficiency and accuracy of the methods on various kernel functions
5. Extension of method applicability to handle unknown slicing functions without analytical knowledge

Methodology Details

Problem Formulation

Given a radial basis function $F: [0,\infty) \to \mathbb{R}$, find a function $f: [0,\infty) \to \mathbb{R}$ such that the slicing relationship $F = S_d[f]$ holds, where $S_d$ is a generalized Riemann-Liouville fractional integral operator:

$S_d[f](s) = \int_0^1 f(ts)\varrho_d(t)dt$

where $\varrho_d(t) := c_d(1-t^2)^{(d-3)/2}$, $c_d := \frac{2\Gamma(d/2)}{\sqrt{\pi}\Gamma((d-1)/2)}$.

Model Architecture

1. Optimization Problem Formulation

The slicing function recovery is transformed into a regularized minimization problem:

$\hat{a} = \arg\min_{a \in \mathbb{R}^K} \|S_d[f_a] - F\|_H^2 + \tau^2\|f_a\|_G^2$

where $f_a = C^{-1}[a]$ is a $K$-term cosine series:

$f_a(t) = a_0 + \sqrt{2}\sum_{k=1}^{K-1} a_k \cos(\pi kt)$

2. Spatial Domain Method (Algorithm 1)

• Matrix Construction: Compute $h_k := S_d[g_k]$, where $g_k$ are cosine basis functions
• Discretization: Approximate integrals using Gauss-Legendre quadrature
• Solution: Solve the least squares problem $\|\hat{H}^T a - \hat{b}\|_2^2 + \tau^2\|Da\|_2^2$

3. Frequency Domain Method (Algorithm 2)

• Operator Representation: Construct matrix representation of operator $S := C \circ S_d \circ C^{-1}$
• Coefficient Calculation: Utilize the relationship $S_{j,k} = S_d[\text{sinc}(\cdot + j) + \text{sinc}(\cdot - j)](k)$
• Optimization Solution: Solve the regularized problem in frequency domain space

Technical Innovations

1. Theoretical Foundation: Establish boundedness theory of slicing operator $S_d$ on different function spaces
2. Numerical Stability: Address ill-posed problems through Tikhonov regularization
3. Error Decomposition: Decompose total error into forward error and slicing error components
4. Convergence Analysis: Prove convergence rates under function smoothness assumptions

Experimental Setup

Datasets

Multiple radial basis functions are tested:

• Gaussian: $F(s) = \exp(-s^2/(2c^2))$
• Laplace: $F(s) = \exp(-c|s|)$
• Inverse Multiquadric (IMQ): $F(s) = (c^2 + s^2)^{-1/2}$
• Thin Plate Spline (TPS): $F(s) = (cs)^2\log(|cs|)$
• Logarithmic Kernel (LOG): $F(s) = \log(|cs|)$
• Bump Function and Multiquadric (MQ)

Evaluation Metrics

• Forward Error: $|F_K(s) - F(s)|$
• Relative L2 Error: $\|s - \hat{s}\|_2/\|s\|_2$
• Runtime Comparison

Comparison Methods

• Direct Method: Truncated Fourier series when analytical solution $f = S_d^{-1}[F]$ is known
• PyKeOps: Highly optimized GPU brute-force computation package
• Three Configurations: S-L2-H1, F-L2-H1, F-H1-H1

Implementation Details

• Use $L = 2^{10}$ quadrature points
• $K = 2^8$ cosine coefficients in domain, $J = 2^{10}$ in range
• Regularization parameter $\tau \in \{10^{-6}, 10^{-7}, 10^{-4}\}$

Experimental Results

Main Results

Forward Error Analysis

For Laplace and Bump functions, forward error $|F_K(s) - F(s)|$ remains below $10^{-2}$ across the entire interval $[0,1]$, with slightly larger errors in irregular regions (e.g., at $s=0$ for Laplace function).

Fast Kernel Summation Accuracy

In tests with $d=1000$ dimensions and $N=M=10^4$ samples:

Runtime Comparison

• Computational Overhead: Coefficient computation time approximately 0.1 seconds (GPU) to 1.3 seconds (CPU)
• Acceleration Effect: Fast summation methods begin to outperform brute-force methods when $N \geq 3 \times 10^3$
• Significant Speedup: Approximately 50-fold acceleration achieved for $N = 5 \times 10^4$ samples

Ablation Studies

The choice of regularization parameter $\tau$ is critical:

• Excessively small $\tau$ leads to numerical instability
• Excessively large $\tau$ results in over-regularization
• Optimal values typically fall within the range $10^{-6}$ to $10^{-4}$

Related Work

Development of Slicing Methods

• Originally appeared in random one-dimensional projections of Wasserstein distance
• Extended to kernel metrics such as MMD
• Closely related to random Fourier features but more general

Fast Kernel Summation Methods

• Traditional Methods: Non-uniform fast Fourier transform, fast multipole methods
• High-Dimensional Challenges: Curse of dimensionality limits applicability of traditional methods
• GPU Implementation: KeOps and similar tools remain competitive at moderate dimensions

Theoretical Foundations

The slicing relationship has multiple names in harmonic analysis and fractional calculus:

• Adjoint Radon transform
• Generalized Riemann-Liouville fractional integral
• Special case of Erdélyi-Kober integral

Conclusions and Discussion

Main Conclusions

1. Theoretical Contribution: Establish complete slicing operator theory, including operator norm estimates and error bounds
2. Numerical Methods: The two proposed algorithms effectively recover coefficients of unknown slicing functions
3. Practical Value: Methods significantly outperform brute-force computation in high dimensions, applicable to large-scale applications

Limitations

1. Dimension Dependence: While complexity is improved, still requires $O(dP)$ computational cost
2. Regularization Sensitivity: Careful tuning of regularization parameters is necessary
3. Smoothness Requirements: Convergence analysis depends on function smoothness assumptions

Future Directions

1. Adaptive Parameter Selection: Develop methods for automatic regularization parameter selection
2. More Efficient Quadrature: Explore specialized quadrature rules to enhance accuracy
3. Application Extensions: Validate method practicality in concrete machine learning tasks

In-Depth Evaluation

Strengths

1. Theoretical Rigor: Provides complete functional analysis framework, including operator boundedness and convergence analysis
2. Method Practicality: Two algorithms each have advantages; spatial domain method is intuitive, frequency domain method is theoretically elegant
3. Comprehensive Experiments: Tests multiple kernel functions from smooth to non-smooth, validating method robustness
4. Excellent Performance: Achieves significant computational acceleration while maintaining accuracy

Weaknesses

1. Parameter Tuning: Regularization parameter selection requires experience, lacking automated methods
2. Memory Requirements: Matrix storage may become a bottleneck in extremely high-dimensional cases
3. Special Case Handling: Method performance is limited for certain ill-conditioned kernels (e.g., LOG)

Impact

1. Academic Value: Provides new theoretical tools and numerical techniques for high-dimensional kernel methods
2. Practical Significance: Holds important value in large-scale machine learning applications
3. Reproducibility: Open-source code provided, facilitating researcher adoption and extension

Applicable Scenarios

• Large-Scale Machine Learning: Particularly suitable for kernel method applications with large sample sizes and high dimensions
• Scientific Computing: Broad application prospects in numerical simulations requiring efficient kernel summation
• Real-Time Systems: Enables fast online inference through pre-computed coefficients

References

The paper cites 52 relevant references spanning kernel methods, fast algorithms, harmonic analysis, and other fields, providing a solid theoretical foundation for the research.