Scalable Random Wavelet Features: Efficient Non-Stationary Kernel Approximation with Convergence Guarantees
Authors: Sawan Kumar
Souvik Chakraborty
Venue: ICLR 2026 Poster (2026)
Keywords: Random features
Non-stationary kernels
Wavelet features
Gaussian Process
Kernel approximation
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: Modeling non-stationary processes, where statistical properties vary across the input domain, is a critical challenge in machine learning; yet most scalable methods rely on a simplifying assumption of stationarity. This forces a difficult trade-off: use expressive but computationally demanding models like Deep Gaussian Processes, or scalable but limited methods like Random Fourier Features (RFF). We close this gap by introducing Random Wavelet Features (RWF), a framework that constructs scalable, non-stationary kernel approximations by sampling from wavelet families. By harnessing the inherent localization and multi-resolution structure of wavelets, RWF generates an explicit feature map that captures complex, input-dependent patterns. Our framework provides a principled way to generalize RFF to the non-stationary setting and comes with a comprehensive theoretical analysis, including positive definiteness, unbiasedness, and uniform convergence guarantees. We demonstrate empirically on a range of challenging synthetic and real-world datasets that RWF outperforms stationary random features and offers a compelling accuracy-efficiency trade-off against more complex models, unlocking scalable and expressive kernel methods for a broad class of real-world non-stationary problems.

Published as a conference paper at ICLR 2026
SCALABLE RANDOM WAVELET FEATURES: EFFICIENT NON-STATIONARI KERNEL APPROXIMATION WITH CONVERGENCE GUARANTEES
Sawan Kumar
Souvik Chakraborty
Department of Applied Mechanics
Department of Applied Mechanics
Indian Institute of Technology Delhi New Delhi I 10016, India
Yardi School of Artificial Intelligence (ScAl)
Indian Institute of Technology Delhi
sawan. kuma r @am. Sitd.ac. ac. in
New Delhi 1 10001, India
souvik iitd.ac.in
ABSTRACT
Modeling non-stationary processes, where statistical properties vary across the in- put domain, is a critical challenge in machine learning; yet most scalable methods rely simplifying assumption of stationarity. This forces difficult trade- off: use expressive but computationally demanding models like Deep Gaussian Processes, scalable but limited methods like Random Fourier Features (RFF). We close this gap by introducing Random Wavelet Features (RWF), a framework that constructs scalable, ota-ssationary kernel approximations by sampling from wavelet families. By harnessing the inherent localization and multi-resolution structure of wavelets, RWF generates an explicit feature map that captures com- plex, input-dependent patterns. Our framework provides a principled way to gen- eralize RFF to the non non-stationary setting and comes with comprehensive theoret- ical analysis, including positive definiteness, unbiasedness, and uniform conver- gence guarantees. We demonstrate empirically on a range challenging synthetic and real-world datasets that RWF outperforms stationary random features and of- fers a compelling eccuracy-efficiens trade-off against more complex models, un- locking scalable and expressive kernel methods for broad class of real-world non-stationary problems.
INTRODUCTION
The ability to model complex, real-world phenomena is one of central challenges in machine learning. Domains such as geospatial modeling, where terrain varies drastically across regions, or speech analysis, where signals exhibit bursts of volatility, are often characterized by pronounced non-stationarity meaning their statistical properties change across the input space. Gaussian Pro- cesses (GPs) offer a principled framework for such problems, providing robust uncertainty estimates and flexible, oonpaarmettri modeling (Williams Rasmussen, (006). Despite these advantages, exact GPs suffer from two major limitations: their expressivity is often constrained by the choice of kernel, and their computational cost scales cubically with the number of training points, rendering them impractical for modern large-scale applications (Liu et al. 2020).
Most of the current approaches force a trade-off between expressivity and efficiency. On one hand, methods like Random Fourier Features (RFF) achieve impressive scalability by approximating the kernel with a linear-in-date feature map (Cutajar al. 2017)] Avron et al. [007 RRaimi & Recht 2007]. Yet, their dependence on Bochner's s theorem (Bochner) 2005) fundamentally restricts them to stationary kernels, which assume uniform behavior across the entire domain. Applying these models to non-stationary data leads to systematic missspecification, resulting in compromised pre- dictive accuracy and uncalibrated uncertainty estimates (Cheema & Rasmussen 2024| Hensman et al. 2013; 2018). On the other hand, expressive models like Deep GPs (Salimbeni et al. 2019), spectral mixtures (Tompkins et al.) 2020), and input- kernels (Rudner et al,] 2020) can capture non-stationarity, but they often reintroduce prohibitive computational costs, complex infer-

---

Published as a conference paper at ICLR 2026
ence schemes, and challenges in optimization and hyperparameter tuning. The gap between scalable stationary models and complex non-stationary ones still remains.
In this work, we close this gap by introducing Random Wavelet Features (RWF), scalable and expressive framework for non-stationary kernel approximation. Instead of relying on globally sup- ported sinusoidal bases like RFF and its variants, construct random features from wavelets fam- ily of functions that are inherently localized in both space and frequency. By sampling wavelets random scales and shifts, RWF generates an explicit feature map that can adapt to local data characteristics. This multi-resolution structure allows the model to capture sharp, localized events with fine-scale wavelets while simultaneously modeling smooth, long-range trends with coarse-scale wavelets. The result is a principled method that generalizes random features to the -tatiiono setting while preserving the linear-time complexity that makes them elegant and efficient. Our main contributions are threefold. First, we provide a comprehensive theoretical analysis of RWF, including positive definiteness of the induced kernels, unbiasedness and variance bounds, and uni- form convergence guarantees with explicit sample complexity. Second, we show that RWF achieves O(ND*) training complexity, retaining the scalability of random feature methods while directly encoding non-stationarity through wavelet localization. Finally, demonstrate empirically syn- thetic, speech, and large-scale regression benchmarks that RWF consistently improves upon sta- tionary random features and offers atleast competitive accuracy-efficiency trade-off against more complex non-stationary models.
1.1 RELATED WORK
Scalable Kernel Approximations. The random features framework was pioneered by (Rahimi & Recht 2007), showing shift-invariant kernels can approximated using random Fourier features with linear-time computations. This framework has since been extended in several directions, in- cluding computationally efficient variants such as Fast kernel learning (Wilson et al. [0014). theoret- ical guarantees approximation error (Sriperumbudur & Szabó Avro et al. 2017 Rudi &
Rosasco, 2017 Li et al. 2021). and structured sampling schemes (Choromanski et al. 2017). There There
works that extend the random Fourier features beyond classical settings using variational approx-
imations (Hensman et al. 2018). adaptive feature learning (Zhen et al. 2020 Shi et al. 22024). and
even a connection to quantum machine learning (Landman al. 2022). Recent progress extends spectral approximations to capture a wider spectrum of kernel families, thereby enhancing the ex- pressivity of scalable feature maps (Langrené al., 2024). While these methods achieve scalability, their reliance on stationary Fourier bases limits their ability to capture non-stationary (Paciorek & Schervish 2003) or localized phenomena, which are crucial in many scientific domains.
Wavelet-motivatee approximations. Wavelets and continuous wavelet frames (Mallat 1999) pro- vide localized multi-scale representations and have previously been used for kernel design through wavelet support vector machines and wavelet kernel learning (Zhang et al. 2004 Yger & Rako- tomamoniy/ 201T), where kernels are derived analytically, or wavelet transforms are used as prepro- cessing. More recently, Guo et al. (2024) proposed Bayesian kernel model based on fixed wavelet bases for high-dimensional Bayesian linear regression. More closely related, Kereta et al.) (2019) studies Monte Carlo discretizations of continuous wavelet frames for signal reconstruction. How- ever, their construction relies on the eigendecomposition of the empirical kernel matrix, leading to cubic computational complexity in the number of samples. While these approaches illustrate the value of wavelets for capturing local structure, they rely on fixed or predefined wavelets dictionaries and do not provide theoretical guarantees such unbiasedness or uniform convergence.
Hybrid and modern kernel learning. Several approaches have been developed to capture non- stationarity in GPs through spectral mixture kernels (Wilson Adams 2013 Remes et al., 2017)
and deep Gaussian processes (Damianou & Lawrence 2013 (allibeem et al. 2019), though both remain costly for large datasets. Scalable variants include KISS-GP (Wilson Nickisch 2015). which exploits structured interpolation, and deep kernel learning (Wilson et al.) (016) combines neural feature extractors with GPs. More recent efforts include deep random features for spatiotem- poral learning (Chen et al. 2024p. graph-based random Fourier features (Zhang et al.) 2025). and adaptive RKHS constructions (Shi et al. 2024). Despite these advances, existing methods often trade off scalability, expressivity, and interpretability. Our work is positioned at this intersection where we aim to design feature maps that inherit the scalability of random features while enabling flexible, non-stationary modeling.
2

---

Published as a conference paper at ICLR 2026
2
PRELIMINARIES AND BACKGROUND
brief review of Gaussian Process regression (GPR), sparse variational GPs, and andom-feature GPs is provided to ground our wavelet construction.
2.1 GAUSSIAN PROCESSES
Given training inputs X [æ1,..., ÆNl E R and targets E TRN, we consider a Gaussian process prior over latent function f. The observations Yn are assumed to be noisy evaluations of this function at the corresponding inputs In:
f ~~ k),
yn f(x)) En where En ~N(,00
(2.1)
We define the covariance matrix Kxx RNXN with [Kxx]i k( m.)). The log marginal likelihood, used for training the hyperparameters GP, is given by the following expression:
logp(y X) = v (Kxx + I - log det(Kxx x + +o"IN) - log((T))
(2.2)
For test input æ., let x.x [k(x., x1), .e................... ..... ... k(x., .N Aiw k( x.), and
(Kxx + o IN ' y.
(2.3)
The predictive posterior moments for test input takes the following form:
He (xx) k.xa,
(2.4a)
dz = Ken k.x [xx + I2IN)" -
(2.4b)
with 1x. = k x_xx The key bottleneck of exact inference is its computational costs O(N3) time, and ONN2) memory. To address these challenges, several approaches have been introduced in the literature; the most common ones are the sparse approximation of GP.
Stochastic Variational GPs (SVGP). SVGP, we introduce My. inducing inputs Z₁ (z1,..., 2M,]] and inducing variables where [Kuljj = k(z,,z) z,). Defining K fuu E TR Ng X with [K sulum k(@n; 2m) and Q5f
f(Z.) equipped with the prior p(u) = N K.n). Ma
K K.
Ky,Kau Knf. the conditional prior becomes
p( u) N(K K_1u, u, Kss Q11))
(2.5)
A Gaussian variational posterior g(u) N(m, S) induces q(f) = N(, E), with μ Am and Σ Kjf Q5 + ASA where A = K Kml Under a Gaussian likelihood p(u.. fa) Nyyn Jn. 2²), ), the ELBO simplifies to
C
By(s,) (orts)llogp(nn p(unl fa] KL(g(u)
(2.6)
n=l
where (ycs))loogpi fa)] = - -2 [(vn 2n) Enm] 2 (og(2mo²))
Using a minibatch B of size b gives the unbiased estimator Ż (N/b) Enes Eg p(y, fa)l (Lg())((())))))( with per-iteration complexity O(bM?) plus a one-time O ( factorization of K
Predictive moments at a test point æ, follow the closed-form GP equations: ps (x(x)) k.u K m and of (x.) k.. k. Ku kns + k..,K K-1K vnd ? where k... [k(#,, 21 ,..,................... kx 2M. }]-
2.2 RANDOM FOURIER FEATURE GPs (RFF-GP)
The random Fourier features approach introduced by Rahimi & Recht (007) approximates station- ary kernels using explicit feature maps. Consider a zero-mean Gaussian process (0. k) with a stationary kernel k(æ, .²)= = k(x x²). By Bochner's theorem (Bochner 2005), the kernel admits the spectral representation
k(x, x²) = / ww (x==²') (ee=a))(w)) p(w) du.
(2.7)
3

---

Published as a conference paper at ICLR 2026
where p(w) is the normalized spectral density of kernel k. Introducing random phase b Unif[0, 2T]], this can be expressed as an expectation over cosine features:
k( 2') Bu.th cos(w x + b) cos(w x' b]]
(2.8)
Approximating the expectation with D Monte Carlo samples {(w,, b;) yields the random fea- ture map X -.
xxx
(D(=)) Φ)(x)) v2 2os(w]] bj),
(2.9)
such that the approximate kernel is k(2z, x') z(x) z(x').
feature space with the dinensionnn features z(-), leading to a Bayesian linear regression model. Placing a Gaus-
From the GP perspective, this corresponds to replacing the
sian prior ur (0, Ip) on the weights, the Gaussian posterior with covariance and mean given by.
Sw = (ID + 2Z
(2.
m S y,
(2.10b)
where Z E RND collects the feature maps of the training inputs. The Gaussian predictive distri- bution for new test point æ. has the mean and covariance defined as,
pa(x,) xxx)) mw:
(2.11a)
Varly, D 2 xxx) = 0z.
(2.11b)
The RFF-GP framework is thus scalable approximation for stationary kernels. However, its liance on globally supported Fourier features limits its ability to model non-stationarity. For further
details on RFF and examples, Appendix A.I
PROPOSED METHODOLOGY
Random Fourier-based kernel approximation methods, which exploit Bochner's theorem (Rahimi & Rechi 2007), yield scalable approximations for stationary kernels but are inherently incapable of modeling non-stationary covariance structures. Sparse variational GPs model non-stationarity with expressive kernels yet rely on indueing sets and cubic costs in M,, per update. We propose Ran- dom Wavelet Features (RWF), which construct non-stationary kernels via multi-resolution, locally supported wavelets. By sampling wavelet scales and shifts, RWF provides an explicit feature map =(-) that: (i) induces positive definite non-stationary kernel; (ii) preserves linear-time training and prediction as in RFF-GPs; and (iii) captures localized, multi-resolution structure that stationary RFF lacks.
3.1 WAVELET-BASED KERNEL. CONSTRUCTION
To model non-stationarity, kernel's properties must adapt across the input domain. Stationary ker- nels, often approximated by Random Fourier Features (RFF), rely on globally supported sinusoidal bases that are inherently spatially invariant. In contrast, wavelets offer a natural alternative by pro- viding basis that is localized in both space and frequency. By randomizing the scale (controlling frequency) and shift (controlling spatial location) of wavelet atoms, we can construct a flexible, non-stationary kernel.
Our construction begins with mother wavelet vi Rdd - R2 a function with zero mean and unit T,2 norm (see Appendix A.2 for details). From v/, we generate family of wavelet atoms via isotropic scaling and translation:
( = -d/2pb 1 x- -)
for scale > 0 and shift E TRd
(3.1)
4

---

Published as a conference paper at ICLR 2026
Each atom Ma.e is a localized "wave packet" centered at t with spatial extent proportional Let θ = (0, oo) x Rd be the parameter space of scales and shifts. We define non-stationary kernel by integrating over this space with respect to non-negative measure p(ds dt):
k< y) = / (a.la)))..e()) p(ds dt).
(3.2)
This construction guarantees positive definiteness, as the integrand is product of scalar features. If (a has a density p(8, 1) 2 O, kernel becomes:
k( v1z,))=AA J.:0///00/
J (.) p(s, f) dt ds.
(3.3)
The density p(s, t) governs the kernel's properties. common choice is factorized form p(s, (s), where Ps (c.g., (og-uniform) spans multiple resolutions and (e.g.. uniform over the data's 's convex hull) provides spatial coverage.
3.2 RANDOM WAVELET FEATURE SAMPLING STRATEGY
The integral in equation 3.3 is typically intractable. We approximate it via Monte Carlo sampling, which forms the basis of our random features.
Definition 3.1 (Random Wavelet Features). Sample (s;, t.) i..d. from distribution with density p(s, t) and define the random feature map Rd + RD as:
1
(a)
DD ( (z)) ..... 1. tn (a:) (x))
(3.4)
The corresponding kernel approximation is k(x y) =<x) =(y))
By construction, k(xx,, y) is an unbiased estimator of k(xx, y). This formulation transforms the kernel method into Bayesian linear model, enabling efficient training and prediction. The full procedure detailed Algorithm[]
Algorithm RWF-GP Training and Prediction
1: Input: Training data (X, y), test inputs X., number of features D, wavelet v/, sampling distri-
bution t)s,t).
2: Hyperparameters: Noise variance o², parameters of p(s,
3: Training:
4: Sample (8) ~ p(s, t) for i- 1 1...., D.
5: Construct feature matrix NxD where Zni vo t'ants (xn).
6: Compute weight posterior: Sw = (lD + -2 2)- and mw 2 y.
7: Optimize hyperparameters (c.g., α², params of p) by maximizing the marginal likelihood of the Bayesian linear model.
8: Prediction:
9: Construct test feature matrix Z. E N.D. where [Z.]ji (onn,,
Compute predictive mean: μ4, 7.
11: Compute predictive variance: o² = diag Z. S Z.T) + +
12: Output: Predictive distribution NV (N,,, )))
3.3 PRACTICAL CONSIDERATIONS
Computational Complexity. RWF efficient because computational cost scales linearly with the dataset size. Constructing D random wavelet features over N inputs of dimension costs O Dd), after which training reduces to the primal form GP regression in D-dimensional feature space. Forming Z Z requires O(ND* and the resulting D x D system is solved in O(D*), so for N >> D the overall training cost is dominated by O(ND2). Predictions require O(D²) per test point. In contrast, Exact GPs scale as O(N3) and SVGP incurs ONNM²) per optimization step due to iterative variational updates. RWF computes its posterior a single closed-form solve, yielding substantial wall-clock speedups for large-scale non-stationary learning.
5

---

Published as a conference paper at ICLR 2026
The key modeling non-stationarity lies in the practical choices for the wavelet family and sampling distribution. The choice of mother wavelet v (c.g., Morlet for time-freqeenc analysis or Daubechies for sharp transitions) and the sampling distribution p(s, t) (e.g.. log- log-uniform for scales, uniform for shifts) (Bergstra & Bengio 2012 Jeftreys 1946) allows the model to adapt to multi-resolution signal structures. For stable training, beneficial to regularize the model by constraining the sampling range for scales and applying weight decay to the linear model.
THEORETICAL ANALYSIS
To analyze the quality our approximation, we establish uniform convergence guarantees. Our analysis relies bounding the complexity of the function class induced by the wavelet features. We define the following key quantities: = sup,,t, (<)) as the uniform bound on the feature magnitude, and = sup= k(æ,a') x²) as the maximum kernel value.
4.1 POSITIVE DEFINITENESS OF WAVELET KERNELS
Theorem 4.1 (Positive Definiteness of Wavelet-Basee Kernels). Let 2/ Rd + R be mother wavelet function, and define the family of wavelets as (x)((x)=
d/2 (x- for scale
s > O and translation E gd. Let p(s, t).): IR x Rd + [0, oo) be a non-negative measure such that the integral is well well-defined and finite for all y. Rd. Then, the function
k(x,y) 3) - I vo, p(s, dt ds
(4.1)
a positive definite kernel on TRd x TRd.
(Proof in Appendix A.44
4.2 UNBIASEDNESS AND VARIANCE BOUNDS
Lemma 4.1 (Unbiasedness). For all æ, y E X, the wavelet random feature approximation is unbi- ased: Exk(x, y)] k(x,y))
(Proof in Appendix A.5
Lemma 4.2 (Variance Bound). For all y E X, the variance of the approximation is bounded:
B²
Var [[ 8) $ D
(4.2)
(Proof in Appendix A.6
4.3 UNIFORM CONVERGENCE GUARANTEES
Theorem 4.2 (Uniform Convergence of Random Wavelet Features). Let C. Rd be compact with diameter di. Let (x. y) be positive definite kernel as in Theorem 4.1 and define the random feature map Rdd R.D by independently sampling (st,t,) for i 1,..., D and setting
1
=(x)=
(W).))()). b.........p......... to)
(4.3)
VD
Assume l and the feature map are Lipschitz continuous with constants La and L,, respectively. Then, for any e > 0,
|x(2) 2(y) - k(xx, 4)1 2 22 4 [4diam(M)I
2d
De²
sup LENEM
E exp
(4.4)
813
(Proof in Appendix A.7

---

Published as a conference paper at ICLR 2026
Doct dP
SVIET
IDF-22
Deep C2
Dm
11
Figure 1: Predictive performance of different GP methods on a step function regression task. Each panel shows the predictive mean (solid line) with =2o confidence intervals (shaded), training data (dots). RWF-GP (ours) captures the discontinuities sharply while maintaining calibrated uncer- tainty. In contrast, Exact GP, Sparse Variational GP, and RFF-GP struggle with sharp transitions, either oversmoothing or miscalibrating the uncertainty.
4.4 SAMPLE COMPLEXITY ANALYSIS
The above theorem provides explicit sample complexity bounds. To achieve approximation with probability at least 5, suffices to choose
8B2 2
diam diiamMM), (M)Lz (M))a
2 5 )
D2
2dd
+ log
(4.5)
This result is derived by inverting the probability bound in Theorem
The constants B and L depend on the choice of wavelet and are discussed in Appendix[A. [A.2 This shows that the number of required features scales logarithmically with the desired accuracy and confidence level.
EXPERIMENTS
This section presents experiments designed to evaluate the performance of the proposed approach. We begin by examining the approximation quality of our approach on non-stationary synthetic data, and then proceed to evaluate it on a highly non-stationary speech signal dataset and benchmark regression tasks, comparing it to various baseline models. Further details about all the experiments can be found in Appendix
Baselines: We compare against scalable and/or expressive variants: SVGP (Hensman et al. 2013 RFF-GP (Rahimi Recht 2007). Deep GPs (Salimbeni et al. 2019). and exact GPs (when feasible) well specialized GP for non-stationary data: Spectral Mixture kernels (Langrené al. 2024
DRF (Chen et al. 2024). IDD-GP (Rudner et al. 2020), and Adaptive RKHS Fourier Feature GPs
(Shi al. 2024)..
EVALUATION ON SYNTHETIC DATA
We first evaluate RWF non-stationary multi-step function, a setting where shallow GPs with stationary kernels fail to capture input-dependent variations (Rudner et al., 2020). Deep GPs, al- though offer more expressiveness, struggle with sharp discontinuities. In contrast, RWF enables shallow GPs to fit accurately: Figure [ shows that RWF-GP captures the non-stationary structure, whereas baselines yield overly smooth or oscillatory fits due limited kernel flexibility. Table illustrates the superior performance of the proposed approach, both in terms of accuracy and train- ing time, over its competitors. Figure summarizes wall-clock time and memory footprints for the compared methods, illustrating the scalability of the proposed approach. Ablation study illustrating convergence proposed approach with feature size shown in Appendix
5.2 TIMIT SPEECH SIGNAL
We evaluate our approach on a regression task derived from the TIMIT corpus, following prior GP- based studies (Shi et al., 2024). TIMIT poses challenge due to strong non-stationarities in the audio signal, such as localized consonant bursts and slowly varying regions. Models relying on stationary kernels struggle capture these variations without either over-smoothing or requiring

---

Published as a conference paper at ICLR 2026
Training
Memory Usage
Training time (min)
Memory (GB)
20
500
1K
10c
20k
50K
500
1K
10x
20k
50c
Number of samples
Number of samples
Figure 2: Scalability on the multi-step function. Time and memory vs. number of training samples on the multi-step function: RWF is most efficient; SVGP and Deep GP incur higher cost.
Table I: Performance comparison of GP baselines on the multi-step function over five runs (mean + std; lower is better). Bold indicates the best result, and underline indicates the second best. Methods:
Exact = Exact GP, SVGP Stochastic Variational GP, RFF = Random Fourier Features,
DRF = Deep-RF GP, DGP = Deep GP, SM Spectral Mixture GP, IDD Inter-domain Deep GP,
A-RKHS = Adaptive RKHS GP. Results for SM, IDD, and A-RKHS are from Shi et al. (2024).
Exact
SVGP
RFF
DRF
DGP
SM
IDD
A-RKHS
RWF (Ours)
0.190
0.231
0.246
0.190
0.162
0.210
0.107
0.095
0.071
RMSE
10.091
=0.144
+0.142
+0.120
#0.110
10.085
+0.00
+0.045
10.011
0.215
0.392
0.238
0.205
0.187
0.201
0.143
0.131
0.112
CRPS
10,030
10.025
10.041
10.032
10.028
10.030
10.020
10.018
10,010
0.042
0.123
0,118
-0.018
-0,268
0.220
-0.820
-1.210
-1.879
NLL
100012
-0.018
1001811
#0.216
+0,211
+0,180
+00.000
10.075
#0.061
Time
12
15
11
18
20
17
17
11
9
large number of features. Unlike RFF, RWF allocates resolution adaptively: small scales capture sharp attacks and large scales capture smooth regions, thus reducing approximation variance for fixed feature size. Results: Table 2 reports RMSE and training time. RWF-GP achieves the lowest error compared to the baselines. RFF-GP performs worst with the same number of features D, reflecting inefficient coverage of localized spectral shifts. Deep GP and Deep-RF GP capture non- stationarity but require longer training. Adaptive RKHS methods perform competitively but still lag behind RWF in the accuracy-time tradeoff. Further details about the experiment are mentioned in Appendix C.2
Table 2: TIMIT regression: RMSE, CRPS, and NLL. (mean + std over 5 runs), and training time. Bold indicates the best result, and underline indicates the second best.
Exact
RFF
SVGP
IDD
A-RKHS
RWF (Ours)
RMSE
2.10:10.008
2.13:100004
2.28:10.005
0.98-10.005
0.54:10.005
0.57:10.15
0.4:100.003
0.42-10.003
CRPS
1.92:00020
1.95=:00118
2.10::0.025
0.85+0.015
0.49::0.100
0.51
0.44+10.099
0.39-10.006
NL.I 3.25110.02
1.31
3.52.L0.14
1.8211.09
1.12:10.66
0.84-10.07
0.70
0.56
133
110
5.3 PERFORMANCE ON UCI DATASET
To evaluate generalization beyond synthetic and domain-specift tasks, we benchmark standard regression datasets from the UCI repository (Dua & Graff 2019), widely used in GP lit- erature. These datasets span range of input dimensions and sample sizes, making them useful benchmark for adaptability and scalability. Following established practice (Salimbeni et al. 2019 Rudner et al. 2020: McDonald & Álvarez 2021). we use a 90/10 train teinttest split, normalize the inputs, and standardize the outputs. Results. Table B reports RMSE and training time. RWF-GP achieves consistently strong predictive performance, yielding the lowest error in five out the seven datasets, and competitive performance on the remaining two datasets. Deep GP and Deep-RF GP
8

---

Published as a conference paper at ICLR 2026
capture non-stationarity but require longer training time. Spectral mixture kemelss provide partial gains on some datasets.
Table Performance on UCI regression benchmarks: RMSE, CRPS, NLL, and training time (min-
utes). Bold indicates the best, and underline indicates the second best.
Data
ENERGY CONCRETE
AIRFOIL 1.5k
STOCK 5k
MOTION
KIN8NM
NAVAL 11k
1k
8k
8k
RFF
0.66.100.33
6.72:10.50
5.34-10.2
1.861:0.03
1.6010002
0.41
0.13:10.002
0.68:10.22
5.92-0117
5.18=0.07
2.13::0.33
1.8740.05
0.10=00.2
0.12::0001
DRF
0.58+00.4
5.01=0001
3.45=0.11
0,95::0.44
0.442:.03
0.12=0.03
0.08+00.011
RMSE
DGP
0.48:10003
4.55
3.66=0.08
0.90::0.03
1.39::0.02
0.09=0.02
0.04:100003
SM
0.67+00.3
5.80+0.19
3,90=0.99
0.92::0.04
1.62+0005
0.11=0.02
0.06+00.001
IDD
0.55:1:0.04
4.20:10008
3.30=0.09
0.88:10004
1.48:0003
0.28
0.07:10.002
1.46.10.03
0.04.10.001
RWF (Ours)
0.42:10.02
4.45
3.20=0.08
0.84::0.03
1.55#0.00
0.09:0.01
0.02:10,001
RFF
0.11+0002
6.21=00.00
4.2==0225
1.72:+0.03
1.51
0.32=0.01
0.11+00.001
SVGP
0.58::0.22
5.20
4.01=0.06
1.70:
1.48:0002
0.11=00.01
0.04
0.06.1:0.00
CRPS
DGP
0.33:10.02
1.25
0.28
0.81110.22
1.32
0.08
0.03
SM
0.63:::003
0,50
3.65=0008
0.84:10.03
1.53
0.09:0.01
0.000100001
IDD
0.49.10.33
3.98
3.008
0.80.10.33
1.42....2
0.2210.01
0.06.10.002
A-RKHS
0:6:::0022
4.10:0.10
2.95-0.09
0.78:10.03
1.11
0.15
0.03:10.001
RWF (Ours)
0.33:00.02
4.2820.12
2.88=0077
0.76:10.03
1.47=00.1
0.07=0.01
0.02:10.001
0.921-0.08
6.85
4.92-0.21
2.02:4:0005
1.724:0.05
0.56-0.03
0.18:100002
SVGP
1.001..07
6.10100.0
4.25.10112
1.98.10.04
1.68
0.32
0.10.10.002
DRF
1.621:005
5.30
3.45 10
1.32,1:0044
0.8211.044
0.42:10.33
0.15:100001
DGP
140+:0.05
5.01
3.68=0.09
0.29::0.33
0.32
0.30=0.02
0.08++0.022
SM
1.98:20008
5.90::020
3.88=0.10
1.36:::.04
1.55:2.00
0.33:00.02
0.12:1:.002
0.52:-0.06
4.85-0.15
3.20=0.08
1.27::0.44
1.88+0.05
0.70=00.3
0.14::0.003
A-RKHS
1.48+00.04
5.01
3.12=0009
1.25+:0.04
1.43
0.33=0002
0.08=000001
RWF (Ours)
0.32
5.10:10.16
2.05 0.08
1.20
1.41
0.28
0.02
0.00:10.001
REF
14
10
10.4
16
14
14
30
DRF
15
16
10.2
15
18
16
33
Time
DGP
15.6
17.8
20
20
35
SM
12.3
21
24
24
IDD
11.1
20
16
19
22
29
A-RKHS
15.3
15
22
35
RWF (Ours)
9.6
10
12
9
24
5.4 PROTEIN DATASET
The Protein dataset has around 45K examples and 9 real-valued input features that originate from biological domain and serve as a practical benchmark for regression tasks. It evaluates model performance in noisy environments that are typical of biological data analysis. Table reports the results. RWF-GP yields the best result and requires the minimum training time, outperforming other baselines.
Table 4: Results the Protein dataset (45K samples). We report RMSE, CRPS, and NLL (mean + std over runs) and training time (minutes). Bold indicates the best result, and underline indicates the second best.
RFF
SVGP
DRF
DGP
SM
IDD
A-RKHS
RWF (Ours)
RMSE CRPS
5.4111001
5.40-1-0001
4.65:10.14
4.35:40.01
4.55+10.02
4.42-10.01
4.32-40.01
4.25-110.02
4.92::0.04
4.88-0003
4.124:.10
3.86:0002
4.01-20.05
3.90-0.03
3.78+0.02
3.65-10.02
NI.I
3.98-10.06
3.92-10.05
3.21:10.08
2.89110.44
3.05110.06
2.95--005
2.82--0.03
2.71110.03
Time (min)
95
120
130
120
133
129
130
90
CONCLUSION
We introduced Random Wavelet Features (RWF), a scalable and principled framework for expres- sive non-stationary kernel approximation. In contrast computationally demanding models like Deep GPs and adaptive convolutional kernels, RWF achieves rare balance of efficiency and ex- pressiveness. By leveraging randomized wavelet families, RWF explicitly encodes the localized,

---

Published as a conference paper at ICLR 2026
multi-resolution patterns inherent in complex real-world processes. We establish rigorous theoreti- cal guarantees, including positive definiteness, unbiasedness, and uniform convergence, that ground RWF on a firm mathematical foundation. Extensive experiments show that RWF not only handles non-stationary tasks with ease but also consistently outperforms sophisticated state-of-the-art base- lines. RWF sets wavelet sampling and integration with deep kernel architectures promising to further expand its reach and impact.
standard for scalable kernel learning, with future directions such adaptive
ETHICS STATEMENT
This work presents methodological advances in scalable random features and kernel approximation using random wavelet features. No human subjects, personally identifiable information, or sensitive data were involved. All experiments use publicly available datasets. The method intended for scientific research; any broader impacts are indirect and depend on the domain-specifie application. We confirm adherence to the ICLR Code of Ethics.
ACKNOWLEDGEMENTS
SK acknowledges the support received from the Ministry of Education (MoE) in the form a Research Fellowship. SC acknowledges the financial support received from Anusandhan National Research Foundation (ANRF) via grant no. CRG/2023/007667.
REPRODUCIBILITY STATEMENT
The proposed Random Wavelet Feature construction is fully specified in Section 3 including the sampling procedure, wavelet parameterization, and kernel approximation formula. All theoretical assumptions are explicitly stated, and complete proofs are provided in the appendix. All code re- quired reproduce the results is publicly available at: /ttps://githus. com/ http//cccccc///////o SRWF .it
REFERENCES
Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Velingker, and Amir Zandich. Random Fourier features for kernel ridge regression: Approximation bounds and statis- tical guarantees. In Doina Precup and Yee Whye Teh (eds.), Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp- 253-262 PMLR, 06-11 Aug 2017.
James Bergstra and Yoshua Bengio. Random search for hyper-paramenn optimization, The journal of machine learning research, 13(1):281-305, 2012.
Salomon Bochner. Harmonic analysis and the theory of probability. Courier Corporation, 2005.
Talay M Cheema and Carl Edward Rasmussen. Integrated variational fourier features for fast spatial modelling with gaussian processes. Transactions on Machine Learning Research, 2024. ISSN 2835-8856.
Weibin Chen, Azhir Mahmood, Michel Tsamados, and So Takao. Deep random features for scalable interpolation of spatiotemporal data. In The Thirteenth International Conference Learning Representations, 2024.
Krzysztof Choromanski et al. Structured evolution with compact architectures for scalable policy optimization. In ICML, 2017.
Kurt Cutajar, Edwin V. Bonilla, Pietro Michiardi, and Maurizio Filippone. Random feature expan- sions for deep Gaussian processes. In Doina Precup and Yee Whye Teh (eds.). Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 884-892 893. PMLR, 2017.
10

---

Published as a conference paper at ICLR 2026
Andreas Damianou and Neil D Lawrence. Deep gaussian processes. In Artificial intelligence and statistics, pp. 207-215. PMLR, 2013.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2019. Irvine, CA: University of California, School Information and Computer Sciences.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrix matrix gaussian process inference with gpu acceleration. Advances in neural information processing systems, 31, 2018.
Wenxing Guo, Xueying Zhang, Bei Jiang, Linglong Kong, and Yaozhong Hu. Wavelet-based bayesian approximate kernel method for high-dimensional data analysis. Computational Statis- fics, (94))3994)2322334, 2024.
James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data. arXiv preprint arXiv:1309.6835, 2013.
James Hensman, Nicolas Durrande, and Arno Solin. Variational fourier features for gaussian pro- cesses. Journal of Machine Learning Research, 18(151):1-52, 2018.
Harold Jeffreys. An invariant form for the prior probability in estimation problems. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 186(1007):453-461, 1946.
Zeljko Kereta, Stefano Vigogna, Valeriya Naumova, Lorenzo Rosasco, and Ernesto De Vito. Monte carlo wavelets: randomized approach frame discretization. In 2019 13th Inter- national conference on Sampling Theory and Applications (SampTA), pp. 1-5, 2019. doi: 10.1 10.1109//ammT744568111 10399903.
Jonas Landman, Slimane Thabet, Constantin Dalyac, Hela Mhiri, and Elham Kashefi. Classically approximating variational quantum machine learning with random fourier features. arXiv preprint 1:Xiv:2210.13200 2022.
Nicolas Langrené, Xavier Warin, and Pierre Gruet. A spectral mixture representation of isotropic kernels to generalize random fourier features. arXiv preprint arXiv:2411.02770 2024.
Zhu Li, Gean-Francois Ton, Dino Oglic, and Dino Sejdinovic. Towards a unified analysis random fourier features. Journal of Machine Learning Research, 22(108):1-51, 2021.
Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfci Cai. When gaussian process meets big data: A review of scalable gps. IEEE transactions on neural networks and learning systems, 31(11); 4405 4423, 2020.
Stephane Mallat. A Wavelet Tour of Signal Processing. Academic Press, 1999.
Thomas McDonald and Mauricio Álvarez. Compositional modeling of nonlinear dynamical systems with ode-based random features. Advances Neural Information Processing Systems, 34:13809 13819, 2021.
Christopher Paciorek and Mark Schervish. Nonstationary covariance functions for gaussian process regression. Advances in neural information processing systems, 16, 2003.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007.
Sami Remes, Markus Heinonen, and Samuel Kaski. Non-stationary spectral kernels. Advances in neural information processing systems, 30, 2017.
Alessandro Rudi and Lorenzo Rosasco. Generalization properties of learning with random features. In I. Guyon, U. Von Luxburg, S. Bengio. H. Wallach, R. Fergus, S. Vishwanathan, and R. Gameet (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
11

---

Published as a conference paper at ICLR 2026
Tim GJ Rudner, Dino Sejdinovie, and Yarin Gal. Inter-domain deep gaussian processes. In Interna- tional Conference on Machine Learning, pp. 8286 8294. PMLR, 2020.
Hugh Salimbeni, Vincent Dutordoir, James Hensman, and Mare Deisenroth. Deep gaussian pro- cesses with importance-weightee variational inference. In International Conference on Machine Learning, pp. 5589-5598. PMLR, 2019.
Xinxing Shi, Thomas Baldwin- Baldwin-MeDonan MoDonll,, and Mauricio A Álvarez. Adaptive rkhs fourier features for compositional gaussian process models. arXiv preprini arXiv:2407.0122 2024.
Bharath Sriperumbudur and Zoltán Szabó. Optimal rates for random fourier features. Advances in neural information processing systems, 28, 2015.
Anthony Tompkins. Rafael Oliveira, and Fabio T Ramos. Sparse spectrum warped input mea- sures for nonstationary kernel learning. Advances in Neural Information Processing Systems, 33: 16153-16164, 2020.
Christopher Williams and Carl Edward Rasmussen. Gaussian processes for machine learning. volume 2. MIT press Cambridge, MA, 2006.
Andrew Wilson and Ryan Adams. Gaussian process kernels for pattern discovery and extrapolation. In International conference on machine learning, pp. 1067-1075. PMLR, 2013.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian pro- cesses (kiss-gp). In International conference on machine learning, pp. 1775-1784. PMLR, 2015.
Andrew G Wilson, Elad Gilboa, Arye Nehorai, and John P Cunningham. Fast kernel learning for multidimensional pattern extrapolation. Advances in neural information processing systems, 27, 2014.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric Xing. Deep kernel learning- In Artificial intelligence and statistics, pp. 370-378. PMLR, 2016.
Florian Yger and Alain Rakotomamoniy. Wavelet kernel learning. Pattern Recognition, 44(10-11); 2614-2629, 2011.
Li Zhang, Weida Zhou, and Licheng Jiao. Wavelet support vector machine. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 34(1):34 39, 2004.
Matthew Zhang, Jihao Andreas Lin, Adrian Weller, Richard E Turner, and Isaac Reid. Graph random features for scalable gaussian processes. arXiv preprint arXiv:2509.03691, 2025.
Xiantong Zhen, Haoliang Sun, Yingjun Du, Jun Xu, Yilong Yin, Ling Shao, and Cees Snoek. Learn- ing to learn kernels with variational random features. In International Conference on Machine Learning, pp. 11409-11419, PMLR, 2020.
A RANDOM FEATURES FOR GAUSSIAN PROCESS
A.1 RANDOM FOURIER FEATURES FOR STATIONARY KERNELS
Let A Rd x 2d + F be stationary kernel, i.c., k(xx, x') k(x x'). By Bochner's theorem, k admits the following representation in terms spectral density p(w):
k(x x²) S ww ((e-x') pmw) dw.
(A.1)
Equivalently,
k(x x') Bwwnp(w) neew x -i T Txx
(A.2)
Expanding the complex exponential into sine and cosine terms gives
kex x') = EEwmp(w) coosww T 2)) cos( T Ta') + sin(w x) sin (wn x¹)]
(A.3)
12

---

Published as a conference paper at ICLR 2026
Introducing auxiliary random phase b Unif 2m), rewrite this
k(x x') E [2 cos (w T "xx b) cos(w> x' + +b]]
(A.4)
Thus, an unbiased Monte Carlo approximation with M samples mm J yields
2
k(x x²)
cos (W,, + bm) cos(wIee + bm).
(A.5)
M
where wm p(w) and bm Unif(O, 2m).
This naturally leads to the random feature mapping
cos(w x + b1) cos(w, x + b2)
((x) /
(A.6)
: cos ( x bm)
so that k(x, x²) x ((x)) ((x').
Example Sqquared-Exponentiae Kernel). The squared- kernel defined as
|
x
k( x') = o² exp
(A.7)
262
where l is the lengthscale and o² the kernel variance. Its Fourier transform (up to normalization) is given by
Ed
p(w)
(2 exp (=p(=1e(((01)))
(A.8)
which corresponds to a Gaussian distribution NO, E L2). Thus, for the -quared-exponential ker- nel, random Fourier features obtained by sampling wm AO, tion.
E (-2Ia) in the above construc-
A.2 WAVELET PRELIMINARIES
Mother wavelet and admissibility. A (real) mother wavelet vb Rd + I satisfies: (i) zero mean Jaa (xx) dæ 0; (ii) square integrability v E I.2(RR) (iii) admissibility constant
f [2( dw co; |w
Cy = /
(A.9)
ensuring invertibility of the continuous wavelet transform (CWT).
Scaled-translated wavelets.
For scale s >> O and translation t E Rdd
(,) -d/2
(A.10)
Energy preserved: I|v),, = | If nb has compact support contained in ball of radius R. then thst has support radius sRt yielding spatial localization.
Continuous wavelet transform. For f (IRd),
{ ()(ss)We,,( Wy(s,)) = [ [())) dæ,, [(x) = CT (/0 Jed
dt ds
(A.11)
sd+s
Vanishing moments. & has M vanishing moments if ())(()))))) dx 0 for all multi-indices [a] M. Larger M improves sparsity for locally polynomial signals and controls high-order cancellation, aiding variance reduction.
13

---

Published as a conference paper at ICLR 2026
Time frequency localization. The Heisenberg-type trade-off bounds the product of spatial vari- ance and spectral variance of tb. Well-localized (e.g.. Morlet, Mexican Hat) wavelets balance this, enabling adaptation to -on-ttationariyy
Bounding feature magnitudes. Suppose scales are sampled in a compact interval E [Sminn Smax] and C with || S C | S c... Then
[.,e(())) S s 4/2 0 I $ min
4/2 C B.
(A.12)
Lipschitzness of wavelets. For any x, x',
((),)()) (=,,((=')l s s- c 1 1 x'1l S 8min 1 || - =x'll =: L x x'1. (A.13)
d/2-11111 (1)
Feature map Lipschitz constant.
Feature map z(xx) = vD
1
Maste (x))) [xm satisfies
D
1
||=(x) - x(x')11 = D (anne(m) Wenne (x) = x' 12,
(A.14)
((a))=((9) =(y) is then jointly Lipschitz with constant 80 b2 S Tyg. Inner product map F( v) =(xr) < 2BLL, under Euclidean metric on Rd x Rd.
Consequences. These bounds verify the assumptions preceding Theorem 4.2 under mild smooth- ness and bounded-scale sampling.
A.3 EXAMPLES OF MOTHER WAVELETS
To ground the proposed framework, we illustrate two specific choices of mother wavelets (x)) used our experiments. Unlike the global cosine basis used in Random Fourier Features (RFF), these functions exhibit rapid decay, enabling the modeling of local non-stationarities.
Mexican Hat Wavelet Defined as the negative normalized second derivative of Gaussian, the Mexican Hat wavelet in d-dimensions is given by:
Mex xz) Ca (1 - 1-|1x|1
(A.15)
where Ca is normalization constant. This wavelet has narrow effective support and exactly zero mean. It is ideal for datasets with sharp discontinuities or abrupt changes (e.g., the Step Function experiment in Section 5.1).
2. Morlet Wavelet. sian window:
The Morlet wavelet consists of a complex plane wave modulated by a Gaus-
MoMor (z) Ca exp
||=|| (-1==11
(os( 2) exp ||اا |]
(A.16)
where wo E Rd is the central frequency. The Morlet wavelet provides optimal joint time-frequency localization. particularly effective for quasi-periodic signals with varying frequencies, such as the TIMIT speech data (Section 5.2).
Comparison with Random Fourier Features. modeling local singularities.
The structural advantage of RWF evident when
RFF (Global Support): Fourier feature (xx) = cos(w x + b) has infinite sup- port. To approximate local step function æ6, RFF requires the superposition of many high-frequency sinusoids to cancel out globally, often leading to oscillations (Gibbs phe- nomenon) in distant regions.
RWF (Local Support): In contrast, wavelet atom (x)(z) is effectively zero outside radius R x s. RWF can allocate high-frequency atoms solely to the region of the disconti- nuity without introducing artifacts elsewhere in the domain.
14

---

Published as a conference paper at ICLR 2026
A.4 PROOF OF THEOREM 4.1
Proof. To show k is positive definite, we must verify that for any finite set of points { C IRd and coefficients {c)i-r IR,
N (cck(2;; 2 0.
(A.17)
=1
Substituting the definition of k(,,,,,))
N N
cock(x; 2z)
ΣΣ (
CiOj
(inle:))a p(s, d) ds
(A.18)
==
Jud
ΣV
(csbs,)()))sss p(s, di ds.
(A.19)
The inner double sum can be rewritten as:
Σ (C) = Σ [avss(os)
N
(A.20)
Thus, the expression simplifies to:
J Jge Σ (ows.)
N
t) dt ds.
(A.21)
Since (EA ca)))))))) 2 0 and p(s,t) 2 0, the integrand is non-negative, proving positive
definiteness.
A.5 PROOF OF LEMMA 4.1
Proof. Define (x, y) = Whases (x) Maane (1) Then
1 D
k(xx, y) D
Σ (x
(A.22a)
2=1 B[ y)1 k(x, y)
(A.22b)
Linearity of expectation yields the result.
A.6 PROOF OF LEMMA 4.2
Proof. Since [Z,x,, y)l = I (x)).... (v)) & Bd almost surely, we have
VarlZ y)] < 14,
(A.23a)
B² D
Varlk(æ, y)] D
Var [Z, <
(A.23b)
However, using the tighter bound Var < E|Z?? < B2, we get the stated result.
A.7 PROOF OF THEOREM 4.2
Outline. (i) Pointwise concentration via Hoeffding; (ii) Cover M x M with an nenet; (iii) Lift bound supremum using Lipschitz continuity; (iv) Optimize n achieve stated constants.
15

---

Published as a conference paper at ICLR 2026
ii) Pointwise concentration.
For fixed (æ, y). define U, 4e,,,, (x) (2),
=(x) xy)
U B[U = k(
|0 S B².
(A.24)
D i=1
Hoeffding yields
De
Pr [1(2)) ((()))==))))))))))))){>)) k( y) 2 2x) S exp
(A.25)
2B*
Noting 2 or tightening via Var [U:] - Ba y) - B and sub-Gaussian refinement) produces equivalent order; re-express constant as S13 in the final statement after net lifting (absorbing improvements from Bernstein-lype refinement).
(iii Covering number. Let (1)) be the minimal cardinality of an 17-net of M in Euclidean norm. Standard volume arguments give
2 diam(MM) d
N(n)= s
(A.26)
Hence M M admits an ne-net with r S (2diam((M) 2d
(iii) Lipschitz lifting. Let (x, y) be arbitrary and choose (æ, 5) E I with | S 77 yy 5 s
7. Write
[=(>)) ==(=))=(() =(y) - k(x, v) < (2(x)) (aar))=(())) - = + z
k(@ 0) (1; k( y)l- (A.27)
By joint Lipschitzness (Section A.2 first and third terms are bounded by
[x(x) (( 2yy) z(²) x - 2BL (1 - कll + | - Ø11 - 4BL
(A.28)
k 3) k( v) < क + 1 e11) 2
(A.29)
Thus, if each net point satisfies
(2(x) T = k( 5) < c/2
(A.30)
and we choose 77 that 4BL, 24477 S e/2, we obtain uniform error < Pick
11 4(2BL2 + Lb)
(using La < 2BL, from Cauchy-Sehwarz).
(A.31)
8BL;
Therefore 2 c/(8B..) suffices; for simplicity we use 17 = e/(4L.) after absorbing constants into exponent.
(iv) Union bound.
With the chosen
D(ee22)3
[ (4dimmMM
De
Pr ssup 1a x - k| e/2) < 2| exp
r
2B4
= 2
71
exp
8B4
(A.32)
Substituting 17 c/(AL:) gives
4 di)
2d
De
Pr sup [x(x) = k(x, w) 2 52 a.y
exp
(A.33)
8B1
Finally, replacing B by B² (tighter variance-barr based constant using Var[U:] S B²b(2, 8) $ B and sub-Gaussian refinement) gives the stated theorem form.
16

---

Published as a conference paper at ICLR 2026
A.8 WAVELET-SPECIES THEORETICAL RESULTS
Lemma A.1 (Stationarity criterion vs. non-stationarity under bounded pc). (s( (t()) with pa independent of t. We define,
Assume p(s,t) =
k<x, y) - (o,1(2) (., (s) pc(t) dt ds,
(x,((x)) s 1 (A.34)
1. If Ps is uniform on a bounded domain D C Rd with nonempty boundary, and db is localized (compactly supported rapidly decaying), then in general k is non-stationary, i.e., there exist (x, y, c) such that
l:(:r + c, c) + k( y-
(A.35)
2. If is trunslation-invariinn on Rd (i.e., Pe(l) = Pe + c) for all shifts c)2 then k is stationary: k(x, y) = k(x - y) which recovers RFF as a special case.
Lemma A.2 (Wavelet localization: explicit feature bounds). Let ab have compact support contained in the ball B(0, R.p)) with || < M and V/V S SGO. Then, for all >> 0 and x. L. E TR...
1 (x) S M 4/2 [1 (|)
(x))) S G --/2--11111 4/2-1
[Im (A.36)
Now for scales s E [smins Smax]. the uniform constants in the concentration bound (Theorem 4.2) might be chosen
d/2 B= Ma min
L G mmn
d/2-1
(A.37)
Corollary A.1 (Wavelet-specific uniform bound with explicit constants). Using Lemma with The- orem 4.2, for S [Smin: Smax] and compactly supported vb,
diam G smin
-d/2-1
2d
D
sup kol( y) - k(x, >5) S a.JEM
exp
8
? Smin
(A.38)
The prefactor and exponential rate depend on (Mg. G Rw. smin) and are therefore wavelet- specific rather than generic constants. This bound quantifies the frequeey trade-off inherent to wavelets but absent in RFF.
Proposition A.2 (Moment cancellation reduces low-scale bias). Assume Pt is locally smooth (CM) around and vi has M vanishing moments. Then
k(x, y) = 5200 Bd
1o (=== (0 (t) gd Px(8) ds
(A.39)
admits Taylor expansion of pe() around I = a where the first M- - terms vanish. Interpretation: Wavelets with higher vanishing moments (e.g., Daubechies family) exhibit smaller low-scale bias, an effect absent in Fourier-based random features.
Corollary A.3 (Comparative constants for specific mother wavelets). For E [8min; Smax). the constants specialize as follows:
Wavelet
Radius (Ry)
Moments (M)
Bound (B) Lipschitz (Lz)
Haar
0.5
dd/ min Ossmin
(
d/2-11
Daubechies-4
Re 1.5
4
d/2
O ( (m
-d/2
-
Mexican Hat
(fast decay)
2
(s
O (smin
Compactly supported wavelets (Haar, Daubechies) yield smaller effective constants, while higher- moment wavelets (e.g., Daubechies) achieve stronger bias reduction of order (M)
B EXPERIMENTAL DETAILS
All the models in the experimental section are implemented using PyTorch and mostly are im- plemented using GPytorch (Gardner et al. 2018), trained by Adam and AdamW Optimizer
17

---

Published as a conference paper at ICLR 2026
NVIDIA A40 GPU. The learning rate for most of the examples is taken to be 0.01 (unless mentioned otherwise) and a batch size of 128. For the Deep-Gr PPpp-PPPo example, we follow the doubly stochastic vari- ational inference as proposed by (Salimbeni et al., 2019) with zero-mean.
Unless specifically stated, we have normalised the imput data for training and inistalized our model with length-scale I 0.1 and o² 0.1kernel variance for TIMIT dataset.
B. EVALUATION METRICS
We evaluate our models using Root Mean Squared Error (RMSE) Let the dataset be denoted as Jn=1- We consider a model f trained on D
D = ((In. yn) ) for training and D {(In: yn )n)] and evaluated using the following criteria. Note here {yn 1 is the ground truth and model predictions f f(X) where X = {X
Root Mean Squared Error (RMSE). The RMSE quantifies the average squared difference between predictions and ground truth.
CRMsE(f:D) V (m.))) (lV D
(B.1)
Empirically estimated as
N Ilyn - 11
FRMSE(f:2
(B.2)
where (Ix,) is the predictive mean at input Zn.-
C ADDITIONAL EXPERIMENTAL DETAILS
C.I SYNTHETIC DATASET
Effect of feature size. Figure reports the convergence behavior of RWF-GP the number of features D increases. As expected. predictive accuracy improves with larger D, but RWF-GP consis- tently attains lower RMSE than RFF-GP across all regimes. Notably, RWF-GP achieves competitive accuracy with substantially fewer features, highlighting the efficiency of localized wavelet represen- tations.
Baseline details. We evaluate all models a dataset consisting of N 4200 training points and Nust = 1800 held-out test points. We utilized the Adam optimizer with a learning rate 0.01. The baseline kernel configurations were chosen follows: the Exact GP. RFF-GP. and SVGP employed a stationary Squared Exponential (RBF) kernel; and the Adaptive-RKHS baseline em- ployed non-stationary convolution kernel. Our proposed RWF-GP utilized Mexican Hat mother wavelet, demonstrating its ability to capture sharp transitions without the stationarity assumptions inherent in the RBF and Matérn baselines.
Comparison with the Non-Stationary Covariance GP. performance of the classical non-stationary covariance model of Paciorek Schervish (2003) on the multi-step function. Results are shown in Table
For completeness, we also report the
Relation prior wavelet-based methods. Prior wavelet kernels (Zhang et al. 2004 ger & Rakotomamonjyy (011) use fixed or learned wavelet dictionaries within deterministic SVM formu- lations, while Kereta al. (Kereta al., 2019) study frame-based discretizations that require spectral decomposition of the empirical kernel matrix. In contrast, RWF constructs an explicit randomized feature map from a continuous wavelet integral representation, yielding an unbiased Monte Carlo estimator with lincar-time primal inference in the number of samples.
C.2 TIMIT SPEECH SIGNAL
Dataset and Preprocessing. We use the TIMIT corpus (630 speakers, 6300 utterances, 16 kHz). For each utterance, 80-dimensional features extracted (25
window, 10 hop, pre-emphasis,
18

---

Published as a conference paper at ICLR 2026
Table 5: Performance comparison of GP baselines the multi-step function five (mean + std; lower is better). Bold indicates the best result and underline indicates the second best. Here, NS-GP Non-stationary Covariance GP (Paciorek Schervish 2003). Kereta et al. (Kereta al. 2019 propose deterministic frame-based discretization that does not yield predictive uncertainty; therefore, CRPS and NLL are not applicable (-).
Method
RMSE
CRPS
NLL
Time
Exact
0.190+0.001
0.215+00.00
0.042++0.12
12
SVGP
0.231+0.014
0.392+00.025
0.123+00018
15
RFF
0.246.10.12
0.238110.041
0.11811.181
11
DRF
0.190.1..12 120
0.20511.002
-0.0181102216
18
NS-GP
0.104110.010
0.16810.007
-1.030-10004
12
DGP
0.162
0.1878
-0.268+02211
20
SM
0,210#0.085
0.201
0.2
17
Kereta
0,822#00009
15
IDD
0.107110.050
0.043
-0.820
17
A-RKHS
0.005-10.045
0.101
-1.210+00.055
11
RWF (Ours)
0.00110.011
0.012
-1.879-0.061
9
4.5
RWF-GP RFF-GP
3.5
3.0
RMS
2.0
0.5
200
Number of features
Figure 3: RMSE vs. number of features D for RWF-GP (Mexican-hat) and RFF-GP on the multi- step function.
CMVN). Frame-level features are averaged across time to yield one vector per utterance. As regres- sion targets, we use either the mean energy a chosen Mel band (mel wwl.bii..knmn. k.nan)) or the mean of PCA component of the spectrogram (mel-pca_k). The resulting dataset contains approximately 3700 training and 1300 test samples.
RWF Configuration. We employ complex Morlet wavelets for time frequency localization. Scales 8 are drawn log-uniformly from [2-1,2] 22) for initialisation, and translations. / are sam- pled uniformly from the input domain. Features are dx(x)) = 1/1 v/s.,, (x). Hyperparameters (Smin; Smax); and noise variance, tuned. Regularization includes (i) clipping extreme scales during warm-up and (ii) ridge penalty Al|w|| 3 with À = 10-4 on Bayesian linear weights. (a) clip- ping extreme scales during warm-up, (b) ridge penalty Al/w|| (with À 10-4) on the Bayesian linear weights' MAP objective surrogate used for hyperparameter inner loops. Wavelet family. Un- less otherwise specified, we employ Morlet and Mexican Hat wavelets as the mother wavelets for constructing random wavelet features.
19


Reviews:
[{"review_id": "QYt5oOgLlc", "reviewer": "Reviewer_Rq9f", "summary": "The paper proposes Random Wavelet Features (RWF) as a Monte-Carlo feature map for Gaussian-process (GP) regression intended to capture non-stationarity while retaining the scalability of random features, and is posed as a generalization of the RFF approach to the case of non-stationary kernels. Concretely, the authors define a kernel by averaging outer products of wavelet atoms over random scale $s>0$ and shift $t \\in \\mathbb{R}^d$ :\n$$\nk(x, y)=\\int_{(s, t)} \\psi_{s, t}(x) \\psi_{s, t}(y) p(s, t) d s d t, \\quad \\psi_{s, t}(x)=s^{-d / 2} \\psi\\left(\\frac{x-t}{s}\\right)\n$$\nand approximate it with $D$ Monte-Carlo samples, yielding an explicit feature map $z(x) \\in \\mathbb{R}^D$ and a linear-in-data GP regression algorithm for large sample sizes. The theory section claims (i) positive definiteness of $k$, (ii) unbiasedness/variance bounds for the Monte-Carlo estimator, and (iii) a uniform approximation bound $\\sup _{x, y \\in \\mathcal{M}}\\left|z(x)^{\\top} z(y)-k(x, y)\\right|$ via covering arguments. Experiments on synthetic data, TIMIT, several UCI datasets, and Protein report lower RMSE and similar or lower training time than RFF-GP, SVGP, deep GPs, spectral mixtures, and adaptive RKHS features.", "questions": "1. Is there anything specific to the wavelet dictionary (other than choices of constants) that can be used to improve the theoeretical results in Sections 4.1 to 4.3?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "The paper derives potentially non-stationary kernels using wavelet based dictionaries together with the random features machinery. Experimental results show slightly better RMSE values for GP regression with slightly reduced training cost compared to baselines on most experiments.", "weaknesses": "Theoretical novelty of the work is extremely limited. The positive-definiteness result (Thm. 4.1) is a direct application of the classic “integral of feature products” construction; nothing wavelet‑specific is used. The unbiasedness and uniform‑convergence proofs follow the standard random‑features and Monte-Carlo estimator analysis tools. Also, the claim about non-stationarity is violated if the weighting probabality distribution $p(s,t)$ in Eqaution 3.3 is taken to be of the form $p_s(s) p_t(t)$ with $p_t$ being translation invariant (For example Lebesgue measure over $\\mathbb{R}^{d}$) The idea of using wavelets in place of Fourier features for scalable kernel approximations also appears in the literature prior to this, including works that explicitly propose random wavelet features for approximate kernels [1]. Closely related Monte‑Carlo wavelet frame constructions exist in learning theory as well [2]. \n\nReferences:\n[1] Guo, W., Zhang, X., Jiang, B. et al. Wavelet-based Bayesian approximate kernel method for high-dimensional data analysis. Comput Stat 39, 2323–2341 (2024). https://doi.org/10.1007/s00180-023-01438-1\n[2] Z. Kereta, S. Vigogna, V. Naumova, L. Rosasco and E. De Vito, \"Monte Carlo wavelets: a randomized approach to frame discretization,\" 2019 13th International conference on Sampling Theory and Applications (SampTA), Bordeaux, France, 2019, pp. 1-5, doi: 10.1109/SampTA45681.2019.9030825.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HfPvksme9L", "reviewer": "Reviewer_QoaL", "summary": "This paper proposes   a scalable method, RWF, for modeling non-stationary processes by sampling from wavelet families. RWF extends RFFs to capture localized, input-dependent patterns, offering theoretical guarantees and strong empirical performance, achieving a superior accuracy-efficiency trade-off for large-scale, non-stationary kernel learning tasks", "questions": "See weaknesses", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1 It provides a rigorous theoretical analysis of Random Wavelet Features (RWF), establishing positive definiteness, unbiasedness, variance bounds, and uniform convergence with explicit sample complexity.\n\n2 RWF achieves O(ND²) training complexity, maintaining the scalability of random feature methods while effectively encoding non-stationarity via wavelet localization.\n\n3 Extensive empirical evaluations on synthetic, speech, and large-scale regression datasets demonstrate that RWF consistently outperforms stationary random features", "weaknesses": "1 While the proposed RWF framework is clearly presented and supported by solid theoretical analysis, I have concerns regarding its novelty. The core formulation of RWF (Eqs. 3.1–3.3) and the sampling procedure (Algorithm 1) appear conceptually similar to existing RWF methods, which also construct kernel approximations using randomly sampled wavelet bases. The authors should clarify what specific differences or innovations distinguish RWF from  earlier approaches such as L. Sun et al., “Wavelet-based Bayesian Approximate Kernel Method for High-Dimensional Data Analysis,” 2023, as well as earlier studies on non-stationary kernel approximations (e.g., Remes et al., “Non-stationary Spectral Kernels,” NeurIPS 2017; Samo & Roberts, “Stochastic Process Regression with Non-stationary Spectral Kernels,” 2015).\n\n2 The manuscript would benefit from better organization. For example, I am very confusing that why so much space is devoted to introducing SVGP in the Preliminaries section, while the discussion on the computational complexity of the proposed method (Section 3.3) is relatively brief and should be expanded.\n\nI am not an expert of GP. I will adjust my score according to the authors' feedback and other reviewer's comments.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vEcgUqVxhu", "reviewer": "Reviewer_daFv", "summary": "This paper presents an approach to build scalable kernel models using wavelet-based decompositions. The method takes a similar approach to random Fourier features (RFFs), but uses wavelet bases, instead of the cosine bases. The result provides a way to build finite-dimensional feature maps whose dot products approximate (possibly non-stationary) positive-definite kernels/covariance functions. Theoretical results are presented on the positive-definiteness of the resulting kernel and on the approximation error due to Monte Carlo approximations. Experiments are presented comparing the proposed RWF against traditional RFFs, sparse variational Gaussian processes, deep GPs, and other baselines.", "questions": "What methods in the literature of wavelet-based kernel machines would be the closest to the approach in this paper?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The paper is well written and organised in a way that makes the presentation easy to follow.\n* The proposed idea is intuitive and the theoretical analysis follows a similar approach to the analysis of RFF methods.\n* The experimental evaluations compare against a wide range of baselines revealing mostly strong performance improvements.\n* Interesting ablation studies are included, comparing for example training time and memory usage, the latter of which is usually rare to find.", "weaknesses": "* My main concern is that the related work discussion fails to cover previous wavelet-based kernel methods and their approximations. A quick literature search reveals wavelet support vector machines (Zhang et al., 2004) and other methods also apparently using wavelet-based kernel decompositions (Guo et al., 2004; Yger, 2011), and potentially others that I'd have missed. Even if these methods are not directly solving the same modelling problem, it'd be important to contrast this paper's approach with them (or at least a few key representatives), to better contextualise and assess the significance of this paper's contribution.\n* In the background on wavelets, a few concrete examples of mother wavelet functions could help to motivate the unfamiliar reader and also demonstrate the capabilities of this modelling approach where, e.g., RFFs would fail.\n\nReferences:\n* Guo, W., Zhang, X., Jiang, B., Kong, L., & Hu, Y. (2024). Wavelet-based Bayesian approximate kernel method for high-dimensional data analysis. Computational Statistics, 39(4), 2323-2341.\n* Yger, F., & Rakotomamonjy, A. (2011). Wavelet kernel learning. Pattern Recognition, 44(10-11), 2614-2629.\n* Zhang, L., Zhou, W., & Jiao, L. (2004). Wavelet support vector machine. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 34(1), 34-39.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HOTFTTdtnA", "reviewer": "Reviewer_28um", "summary": "The manuscript proposes random wavelet features to close the gap between expressibility and computational cost of non-stationary kernels for Gaussian Processes. The objective and approach are laudable, but critical details are missing in the results.", "questions": "(1) It is not immediately clear why \"The RFF-GP framework is thus a scalable approximation for stationary kernels.\" The kernel z(x) z(x') is non-stationary. Please explain the discrepancy.\n(2) \"Random Fourier-based kernel approximation methods exploit Bochner’s theorem (Rahimi & Recht,\n2007) to yield scalable approximations for stationary kernels, but they struggle to capture non-\nstationarity.\" Are they stationary kernels, or do they struggle with non-stationarity? Stationary kernels don't struggle with non-stationarity; they simply cannot model it. Please explain. \n(3) It was difficult for me to decipher the setup of the competing methods. For example, the exact GP in the examples performs pretty poorly, but should not. Given the same setup, it should beat approximate methods. So a stationary kernel in an exact GP should perform better than SVGP with the same kernel. What was the exact setup of the competing methods?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- Broad and impactful scope.\n- Understandable writing.\n- Strong theoretical results.\n- Good diversity of experimental tests.", "weaknesses": "- RMSE only in experiments. At a minimum, CRPS should be added.\n- Unclear setup of the competing methods (see questions). \n- No truly large GP in experiments (> 100,000 data points): There are GPs now that can run on several million points; 45000 max is a little disappointing. \n- Writing has to be polished a bit (missing articles and small typos)\n- Missing comparison to more traditional but powerful ways to encode non-stationarity, for instance, \"Paciorek, Christopher, and Mark Schervish. \"Nonstationary covariance functions for Gaussian process regression.\" Advances in neural information processing systems 16 (2003).\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)