Signed random Fourier features for fast density estimation with indefinite kernels

Xie Wang; Nicolas Langrené; Wen Chen

Source-linked AI summary

Signed random Fourier features for fast density estimation with indefinite kernels

Xie Wang, Nicolas Langrené, Wen Chen

arXiv:2608.29265v1Aug 29, 2026stat.COcs.LGmath.PRstat.ML

TL;DR

Large-scale KDE is costly because direct pairwise evaluation is quadratic, while standard random Fourier features exclude many non-positive-definite kernels. The paper introduces signed random Fourier features, which use signed spectral sampling to approximate such kernels, and reports accurate, substantially faster KDE on one million points.

Problem
Direct KDE requires quadratic computation, while standard random Fourier features are limited to positive definite kernels and therefore exclude many kernels used in KDE.
Method
SRFF generalizes RFF to kernels with absolutely integrable inverse Fourier transforms by sampling from normalized absolute spectral density and attaching spectral signs.
Results
On one million points, SRFF computed KDE with M = 1000 in about 2 seconds, more than a hundred times faster than direct KDE evaluation.
Takeaways & Limitations
SRFF provides a practical fast-KDE approach for indefinite compact kernels, including the multivariate Kuttner-Golubov class and several common KDE kernels.
Takeaways & Limitations
The signed Monte Carlo framework assumes absolute integrability, and the fast KDE implementation has an author-identified limitation for certain evaluation settings that do not benefit from the cosine decomposition.

Abstract

from arXiv · show

Kernel density estimation (KDE) is one of the most fundamental statistical estimators of density functions. Its direct implementation on a dataset of $N$ points incurs an $\mathcal{O}(N^{2})$ computational cost, which is prohibitive for large-scale datasets. Kernel approximation techniques can be applied to bring the computational cost down to $\mathcal{O}(N)$. The random Fourier features (RFF) technique, based on sampling from the spectral density of the kernel function, has become popular to speed up kernel estimators for machine learning applications. Unfortunately, it is restricted to positive definite kernels, while the majority of kernel functions popular in KDE, such as the parabolic kernel, do not satisfy this property. To overcome this limitation, this article introduces the signed random Fourier features (SRFF) technique. It is a generalization of RFF compatible with indefinite kernels whose inverse Fourier transform is absolutely integrable. The motivation for introducing this method is to speed up KDE in the case of multivariate compact kernels, which are generally not positive definite. We detail how to implement SRFF for both product kernels and isotropic kernels. For the class of Kuttner-Golubov kernels $K(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=(1-\left\Vert \boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\Vert ^α)^β\mathbf{1}_{\{\left\Vert \boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\Vert \leq1\}}$ where $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$, $\boldsymbol{x}_{j}\in\mathbb{R}^{d}$, $α>0$, $β>0$, which includes the triangular, parabolic, biweight, triweight, and other kernel functions of interest for KDE as particular examples, we provide an explicit acceptance-rejection algorithm to sample from its signed spectral density. Our numerical tests on a dataset of one million points confirm the computational efficiency and accuracy of SRFF for large-scale KDE.

1. Introduction

Kernel estimators are accurate but difficult to scale because direct pairwise computation is quadratic, while standard kernel approximations generally require positive definiteness. The paper introduces SRFF to extend Fourier-feature acceleration to indefinite KDE kernels.

Direct evaluation of all pairwise kernel values requires O(N^2) operations and O(N^2) memory, limiting large-scale kernel density estimation.The quadratic cost becomes prohibitive as dataset size grows.
Nyström and standard random Fourier features reduce computational and memory costs but require the kernel to be positive definite.This requirement restricts their use for many KDE kernels.
Most KDE kernels of interest, including the parabolic kernel, are not positive definite, making common kernel decomposition techniques inapplicable.The Gaussian kernel is identified as one of the few exceptions.
SRFF samples from the normalized absolute spectral density, attaches each frequency the spectrum's sign, and avoids computing spectral-mass weights.An acceptance-rejection algorithm is used to simulate frequencies.
The paper develops SRFF for compact multivariate kernels, including triangular, parabolic, biweight, triweight, and tricube kernels, and evaluates speed and accuracy numerically.The paper is organized around signed Monte Carlo, SRFF, fast KDE, numerical evaluation, and conclusions.

2. Signed Monte Carlo

This section generalizes Monte Carlo integration from probability densities to absolutely integrable signed functions. The resulting estimator samples from a normalized absolute density and weights each evaluation by the sampled function's sign.

2.1. Signed Monte Carlo formula: Signed Monte Carlo relaxes the usual nonnegativity requirement and applies to absolutely integrable functions f over R^d.The function f need not be a probability density.
2.1. Signed Monte Carlo formula: For an absolutely integrable nonzero f, f can be represented using samples from |f| normalized by its L1 norm and the sign of f.The sign function contributes +1, 0, or −1 according to the value of f.
2.1. Signed Monte Carlo formula: The preferred representation uses an explicitly known integral because the L1 norm generally lacks a closed-form expression.An odd number of Monte Carlo samples can avoid a zero denominator event in the estimator.
2.2. Signed Monte Carlo estimator: The signed Monte Carlo estimator uses M independent vector samples and multiplies each g(X_m) by sign(f(X_m)).The samples are drawn from the distribution with density |f| normalized by its L1 norm.
2.2. Signed Monte Carlo estimator: As M tends to infinity, the estimator converges almost surely under the stated sampling and integrability conditions.The convergence follows from the strong law of large numbers and the signed representation.

3. Signed random Fourier features

This section extends random Fourier features to continuous, symmetric indefinite kernels whose inverse Fourier transforms are absolutely integrable. It also describes simplified constructions for product and isotropic kernels.

SRFF extends the spectral Monte Carlo foundation of RFF to kernel functions that are not positive definite.The extension uses signed Monte Carlo theory to construct random Fourier features for indefinite kernels.
Standard RFF samples random projections from a nonnegative spectral density, so its kernel representation requires positive definiteness.Under symmetry, the inverse Fourier transform is the density of the random projection used in the Monte Carlo approximation.
Theorem 1 requires continuity, symmetry, K(0) ≠ 0, and absolute integrability of the inverse Fourier transform, but not its nonnegativity.This yields a probabilistic representation for a broader class of kernels than standard RFF.
SRFF samples frequencies from |f| / ∥f∥1 and attaches the spectrum's sign to each simulated frequency.Acceptance-rejection simulation avoids computing the generally non-analytical scaling constant ∥f∥1.
3.3.1. Product kernels: For product kernels, the inverse Fourier transform factorizes across dimensions, allowing independent signed-frequency sampling from each univariate factor.The multivariate construction applies when each factor's inverse Fourier transform is absolutely integrable.
3.3.2. Isotropic kernels: For isotropic kernels, the spectral distribution admits a scale-mixture representation using a uniformly distributed unit-sphere direction and an independent nonnegative radius.This representation supports simplified treatment of radially symmetric kernels.

4. Fast kernel density estimation

SRFF reduces kernel MVM and KDE computation from direct quadratic cost to O(MN) by combining signed spectral sampling with fast feature decompositions. The method covers compact multivariate kernels, while adaptive balloon bandwidths remain outside the direct decomposition.

4.1. Fast KDE computation: Direct kernel MVM and KDE evaluation require O(N^2) operations, making large datasets computationally expensive.The kernel density estimator is a special case of the kernel MVM.
4.1. Fast KDE computation: Using M random projections, SRFF evaluates the kernel MVM at all dataset points in O(MN) operations.Algorithm 1 implements this fast computation using the SRFF approximation.
4.1. Fast KDE computation: The fast SRFF algorithm samples random projections from the absolute spectral density and attaches each frequency the sign of the inverse Fourier transform.For isotropic kernels, the projection is represented through an independently sampled spectral radius and a uniform random direction.
4.1. Fast KDE computation: Adaptive balloon bandwidths do not benefit from the cosine decomposition because source and target variables cannot be disentangled, so the direct fast algorithm is inapplicable.A rank-based uniformization followed by constant-bandwidth SRFF is proposed as an alternative approach.
4.2. Choice of multivariate kernel function: The Kuttner-Golubov compact-kernel class includes triangular, parabolic, biweight, and other kernels commonly used in density estimation.The class is relevant because many commonly used KDE kernels are not positive definite or are positive definite only in low dimensions.
4.3. Spectral rejection sampling: For isotropic compact kernels, acceptance-rejection sampling uses an envelope whose tail matches the inverse-power decay of the spectral-radius density.The normalized envelope combines beta, uniform, and Pareto components and can be sampled by composition.

5. Numerical experiments

The experiments evaluate SRFF accuracy and computational efficiency against direct KDE on one million points in dimensions 1 and 2. Increasing M reduces approximation error while increasing SRFF runtime approximately linearly, enabling an accuracy–cost trade-off.

Experimental design: The experiments compare SRFF with direct summation for computational efficiency and approximation precision.Accuracy is measured by mean absolute error between the kernel matrix-vector multiplication and its SRFF approximation.
Accuracy evaluation: For N = 1,000,000, the accuracy study tests 28 values of M from 17 to 196,609 in d = 1 and d = 2, averaging five random trials.Within each trial, frequencies are nested by using prefixes of one master sequence.
Accuracy evaluation: MAE decreases overall as M increases, with log–log slopes close to −1/2; at M = 196,609, average MAE is approximately 4.85 × 10−4 for d = 1 and 3.0 × 10−4 for d = 2.The reported trend matches the expected Monte Carlo approximation behavior.
Runtime evaluation: SRFF runtime has complexity O(NM) and increases approximately linearly with M in both dimensions, unlike direct KDE’s O(N2) complexity.Runtime measurements average seven repeated experiments using the parabolic kernel.
Runtime evaluation: At M = 1000, SRFF computes KDE on N = 1,000,000 points in about 2 seconds, more than a hundred times faster than direct KDE evaluation.The accuracy and runtime results can be combined to select M for a desired trade-off between approximation accuracy and computational cost.

6. Conclusion

The conclusion presents SRFF as an extension of RFF to non-positive-definite kernels and summarizes algorithms for signed spectral sampling and fast kernel matrix-vector multiplication. Numerical tests confirm its speed and accuracy for large-scale KDE.

Contribution: SRFF extends classical RFF to kernels that are not positive definite, including most compact kernels used in KDE.The method targets kernels whose inverse Fourier transforms can have negative values.
Contribution: For Kuttner-Golubov kernels, the paper provides an explicit Fourier-transform formula and an algorithm for simulating their signed spectral distribution.The sampling algorithm uses acceptance-rejection for the absolute value of the inverse Fourier transform.
Contribution: Sampling from |f| / ∥f∥1 by acceptance-rejection avoids analytically computing ∥f+∥1, ∥f−∥1, or ∥f∥1.The conclusion states this approach is more efficient than independently sampling from the positive and negative parts.
Conclusion: The paper constructs an SRFF algorithm for fast kernel matrix-vector multiplication, with KDE as a particular case, and reports numerical confirmation of speed and accuracy.

A. Fourier transforms of multivariate compact kernels

For continuous isotropic kernels, the Fourier transform is radial and can be expressed using a Bessel function. Bochner’s theorem links this transform being a density to positive definiteness of the kernel.

Fourier-transform structure: A continuous isotropic kernel K(u) = k(∥u∥) has a radial Fourier transform.The transform is represented through a Bessel function of the first kind.
Positive definiteness: Bochner’s theorem states that the Fourier transform is a spectral density if and only if the kernel K is positive definite.The appendix uses this characterization to analyze kernels relevant to density estimation.

A.1. Kuttner-Golubov kernels

This section develops normalization and Fourier-transform formulas for the compact Kuttner-Golubov kernel family. The formulas use radial integration, Bessel-function series, and the Fox–Wright generalized hypergeometric function.

Kernel definition: The Kuttner-Golubov kernel is parameterized by α > 0 and β > 0.The section recalls the compact isotropic form before deriving its integral and Fourier transform.
Normalization: The kernel integral is computed explicitly to provide the normalizing constant needed when a kernel estimator scales K as a density.The proof uses a polar change of variables.
Relation to prior formulas: The appendix notes that the formulas can also be deduced from previously reported formulas after a change of notation and simplification.The paper nevertheless supplies more detailed proofs.
Fourier transform: The Fourier transform of K(u) = (1 − ∥u∥α)β 1{∥u∥≤1} is derived from the radial Fourier-transform representation.The derivation substitutes a series representation of the Bessel function and then changes variables.
Fourier transform: The resulting Fourier-transform formulas are expressed using the Fox–Wright generalized hypergeometric function.This function is defined by its numbers of numerator and denominator parameters.

A.2. Symmetric beta kernels

Setting α = 2 yields multivariate symmetric beta kernels, including parabolic, biweight, and triweight kernels. Their Fourier transform is expressed using the confluent hypergeometric function 0F1 and rising factorials.

Setting α = 2 in equation (21) produces the class of multivariate symmetric beta kernels.
This class includes the multivariate parabolic, biweight, and triweight kernels.
The Fourier transform of K(u) = (1 − ∥u∥^2)^β 1_{∥u∥≤1} is given in terms of 0F1.Here, 0F1 is the confluent hypergeometric function with zero numerator and one denominator parameter.
The rising factorial used in the transform is defined by (b)_n = Γ(b + n)/Γ(b).

A.3. Askey kernels

Setting α = 1 yields the Askey functions, including the multivariate triangular kernel. Their Fourier transform is represented with 1F2, and duplication formulas generalize these expressions to α = 2/n.

Setting α = 1 in equation (21) gives the Askey functions.
The Askey functions include the multivariate triangular kernel as a particular case.
The Fourier transform of K(u) = (1 − ∥u∥)^β 1_{∥u∥≤1} is expressed using 1F2.The notation 1F2 denotes a confluent hypergeometric function with one numerator and two denominator parameters.
Applying the Legendre duplication formula simplifies the α = 1 Fourier-transform expression and its multiplicative constant.
The Gauss duplication formula allows equation (28) to be written with pFq when α = 2/n, generalizing the α = 2 and α = 1 cases.Equation (36) generalizes equation (31) for α = 2, n = 1 and equation (34) for α = 1, n = 2.