REVIEW 4 major objections 4 minor 1 cited by
Fast kernel methods: Sobolev, physics-informed, and additive models
T0 review · 4 major / 4 minor · reviewed 2026-08-05 · deepseek-v4-flash
Pith's one-line read Kernel ridge regression on a truncated Fourier basis can be computed exactly in O(n log n) time and memory, and achieves the Sobolev minimax rate n^{-2s/(2s+d)}.
desk verdict Useful algorithmic framework, but the theory leans on an unproved Fourier norm-equivalence lemma and an unstated CG iteration assumption. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing object is the truncated Fourier feature map ϕ(x) = (exp(iπ⟨k,x⟩/(2L)))_{‖k‖∞≤m} combined with the non-uniform fast Fourier transform (NUFFT), which evaluates sums Σ_j Y_j exp(iπ⟨k,X_j⟩/(2L)) in O(n log n). The covariance Φ*Φ/n is Hermitian and d-level block Toeplitz, so it is stored as one row and applied via FFT; with a diagonal (or similarly structured) regularization matrix M, the linear system is solved by conjugate gradient. Proposition 2.1 supplies the norm equivalence and tail-decay estimate that tie Fourier truncation to the Sobolev minimax rate.
What would settle it
Take a bounded Lipschitz domain that is not a cube (for example a disk in 2D), pick s = 1, draw n points from a density on it, fit Equation 2 with m = n^{1/(2s+d)} and λ = n^{-2s/(2s+d)}, and measure the test L2 error against a known f⋆ with rough boundary behavior; the claimed O(n^{-2s/(2s+d)}) decay would be contradicted if the measured exponent is slower, because that decay relies on the Fourier norm equivalence and tail bound of Proposition 2.1 for arbitrary Lipschitz domains.
Extended reading notes
Core claim
The paper's central claim is that Equation 2—the ridge estimator in a finite-dimensional feature space—can be computed exactly with O(n log n) complexity whenever the feature map is a truncated Fourier basis and the regularization matrix is diagonal or has similar fast structure. The covariance matrix Φ*Φ/n is d-level block Toeplitz and Hermitian, so it is determined by a single row that the NUFFT computes in O(n log n); matrix-vector products run in O(m^d log m) by FFT, and conjugate gradient solves the system. Propositions 3.1 and 3.2 show that with m = n^{1/(2s+d)} and λ in [n^{-1}, n^{-2s/(2s+d)}], the L2(PX) error is O(n^{-2s/(2s+d)}), the Sobolev minimax rate; the additive model in Pro
Load-bearing premise
The rates rest on the imported claim from the authors' earlier work that every function in the Sobolev space on a bounded Lipschitz domain has Fourier coefficients that faithfully represent the Sobolev norm and whose truncated tail decays as m^{-2s}; if that fails on non-rectangular domains, the minimax-rate guarantees do not follow.
Editorial extensions
If this is right
- Sobolev kernel regression with m = n^{1/(2s+d)} and λ ≈ n^{-2s/(2s+d)} reaches the minimax L2(PX) rate n^{-2s/(2s+d)} while running in near-linear time.
- The low-bias variant (M = I) achieves the same minimax rate with a smaller regularization bias, and it outperforms the standard Sobolev kernel in the paper's experiments, especially at large n.
- Physics-informed constraints—constant-coefficient linear differential operators—can be added as a diagonal-plus-Toeplitz penalty, so enforcing a PDE does not change the O(n log n) complexity.
- Additive models built from univariate Fourier components recover the univariate rate n^{-2s/(2s+1)} in dimension d, with O(d^2 n log n) total cost, avoiding the curse of dimensionality.
- The covariance assembly dominates cost, so many values of the regularization parameter can be tested in a cheap grid search; the paper demonstrates a 300-point λ grid at n = 10^8 in under 30 seconds on a GPU.
Reading between the lines
- Editorial inference: if the Fourier norm equivalence holds, the same O(n log n) machinery should transfer to other translation-invariant kernels with explicit Fourier symbols (e.g., Matérn or periodic kernels) by changing only the diagonal regularization weights.
- Editorial inference: the near-linear cost and cheap λ-grid suggest that model selection for kernel methods could be re-run routinely at data scales where full kernels are infeasible—something the paper demonstrates but does not promote as a general workflow.
- Editorial inference: because the method is exact on the truncated space, it provides a clean benchmark for approximate large-scale kernel algorithms; comparing random-feature or Nyström error against this estimator on the same data would quantify the cost of approximation.
- Editorial inference: the additive construction separates univariate components but does not model interactions; a natural testable extension is to add a small number of pairwise Fourier blocks to capture interactions while keeping a tractable Toeplitz structure.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a truncated-Fourier / NUFFT framework for kernel ridge regression that, for Sobolev, physics-informed, and additive models, computes the ridge estimator with O(n log n) time and memory complexity on GPUs. The theoretical part proves upper bounds for the L2(PX) risk of the truncated Fourier estimator; with m = n^{1/(2s+d)} and a suitable lambda, the Sobolev and low-bias Sobolev estimators are claimed to reach the minimax rate n^{-2s/(2s+d)}, and the additive estimator the univariate-like rate n^{-2s/(2s+1)}. The empirical section reports experiments up to 10^8–10^10 samples, including comparisons with PyGAM for additive models.
Significance. If the results hold, this is a practically significant contribution: it gives a GPU-friendly, near-linear-time exact ridge regression solver for three important model classes and provides rate guarantees consistent with classical kernel theory. The paper ships code, the theoretical upper-bound proofs follow the standard ridge-regression framework, and the large-scale experiments are a clear strength. The main technical risk is that all rate theorems depend on an imported Fourier representation result for arbitrary bounded Lipschitz domains (Proposition 2.1), which is not proved and is not a routine textbook fact. The additive-model and minimax claims also need sharper qualifications. These issues are local and fixable, so the manuscript is a candidate for major revision rather than rejection.
major comments (4)
- [§2.2, Prop. 2.1 and Appendix A.1] Proposition 2.1 is the load-bearing result: every rate theorem (Props. 3.1, 3.2, 5.1) passes through the norm equivalence and the m^{-2s} Fourier approximation error. The proposition is imported from Doumèche et al. (2024) without a proof or a theorem number, and it is not a standard textbook statement for arbitrary bounded Lipschitz domains: it requires constructing a periodic extension to [−2L,2L]^d whose Fourier coefficients satisfy (ii). This is not merely cosmetic: the proof of Prop. 2.2 integrates f over the full cube although f is defined only on Ω, so the extension is genuinely needed. Please include a self-contained proof, or an exact statement with the precise assumptions, and state whether the pointwise identity in (i) holds on Ω only or on the extension.
- [§3.1, after Prop. 3.1; Assumption (v)] The text says the rate n^{-2s/(2s+d)} is the Sobolev minimax rate under Assumptions (i)–(v). But Assumption (v) only bounds the design density from above; it does not bound it away from zero. Standard minimax lower bounds for L2(PX) risk require a lower bound on the density over the support. Without such a lower bound, designs effectively supported on lower-dimensional or small-volume subsets can enjoy faster rates, so the lower-bound claim is not established under the stated assumptions. Either add a lower-density assumption (e.g., 0 < c ≤ dPX/dx ≤ κ) or present the result as an upper bound with the minimax statement restricted to the appropriate design class.
- [§5, Prop. 5.1 and Appendix A.3] There is an inconsistency between the statement and the proof: the proposition sets m = n^{1/(2s+1)}/d, while the proof uses m = n^{1/(2s+1)}. Since the univariate truncation error is O(m^{-2s}), the extra /d is harmless for fixed d, but the two statements should match. More importantly, the complexity paragraph says the linear solve costs O(d^3 m^3) = O(n) 'provided s≥1'. With the stated m, this cost is n^{3/(2s+1)}, which is not O(n) when 1/2 < s < 1. The assumptions only require s > 1/2, so for this regime the claimed O(d^2 n log n) complexity is not established. Please either use an iterative solver with a rigorous iteration bound for the low-smoothness regime or state the complexity result only for s ≥ 1.
- [Abstract and §2.3] The expression 'exact, fast, and memory-efficient computations' is problematic because the NUFFT is an approximate fast transform, while the statistical rates in Props. 3.1, 3.2, and 5.1 are proved for the exact Fourier sums. To make the theoretical and algorithmic claims consistent, the paper should state the NUFFT accuracy/tolerance needed so that the transform error is negligible relative to the n^{-2s/(2s+d)} (or n^{-2s/(2s+1)}) statistical error, and similarly clarify the conjugate-gradient stopping criterion. As written, the central 'exact O(n log n)' claim is stronger than what the algorithm and analysis establish.
minor comments (4)
- [Prop. 3.1, notation] The proposition defines α = sum_k (1+||k||^{2s})^{-1}, but the bound uses α^2. In the Appendix, α^2 is the sum itself. This is a factor/notation slip; it does not change the rates because α is finite, but the display should be consistent with the proof.
- [Throughout] Typos and small presentation issues: 'RKHS approximation if finite-dimensional' should be 'is finite-dimensional'; 'Untractable domains' should be 'Intractable domains'; 'bonding' in Appendix A.2 should be 'bounding'; '2d-NUFFT' in Section 5 should be '2D-NUFFT'.
- [Figures 1 and 4] The risk curves are averaged over resamples but no error bars or confidence intervals are shown. Given the noisiness of the displayed MSE at small n, error bars would help the reader judge the fit to the n^{-2/3} line.
- [References] The citation 'Golub & Loan (Theorem 4.8.2, 2013)' should be 'Golub & Van Loan'. Also, the reference to Doumèche et al. (2024) in Prop. 2.1 should give a theorem/lemma number so the reader can verify the imported statement.
Circularity Check
No significant circularity: the rate claims are derived from an explicit ridge objective plus parameter-free Fourier lemmas; self-citations count as independent support, not circular inputs.
full rationale
The paper's central derivation is not circular. Equation (2) defines an explicit ridge estimator in a truncated Fourier basis; no parameter is fitted to the target quantity and then relabeled as a prediction. The convergence-rate propositions (3.1, 3.2, 5.1) are proved in Appendices A.2 and A.3 via a standard bias-variance-approximation decomposition, with the approximation error controlled by Proposition 2.2. Proposition 2.2 is proved in Appendix A.1 using Proposition 2.1 and a counting lemma from Doumèche et al. (2024). These are parameter-free mathematical statements with stated assumptions (s ∈ N*, Ω a bounded Lipschitz domain, PX density bounded) and do not include the target minimax rate. Under the provided rubric, such a self-citation is independent support, even though it is load-bearing. The main scientific risk is correctness/verification of Proposition 2.1 for arbitrary bounded Lipschitz domains — if that lemma failed, the bias terms would not decay as claimed — but that is a correctness concern, not circularity. No equation reduces to itself by construction, no fitted parameter is renamed as a prediction, no uniqueness theorem is imported to forbid alternatives, and no known empirical pattern is merely renamed. The O(n log n) claim is an algorithmic complexity analysis based on NUFFT and Toeplitz structure, not a fitted prediction. Empirical comparisons use external baselines (PyGAM), so there is no circular benchmark. Overall, the derivation chain is self-contained in the circularity sense.
Assumptions & free parameters
free parameters (2)
- regularization parameter lambda =
n^{-2s/(2s+d)} (theory), n^{-2s/(2s+1)} (additive)
- PDE penalty weight mu =
1 (used in physics-informed experiments)
assumptions (5)
- domain assumption Fourier representation of H^s(Omega) (Proposition 2.1, Doumeche et al. 2024): for any f in H^s(Omega) there is a coefficient sequence with norm equivalence and the truncation error decays as m^{-2s}.
- domain assumption Assumptions (i)-(v): bounded Lipschitz domain Omega subset [-L,L]^d, f* in H^s(Omega) with s > d/2, noise with zero conditional mean and bounded variance, and marginal density bounded by kappa.
- domain assumption Assumption (vi): known linear constant-coefficient differential operator D with D f* = 0 on Omega.
- domain assumption Assumption (vii): additive decomposition f*(x) = g1(x1) + ... + gd(xd) with each g_l in H^s.
- ad hoc to paper Conjugate-gradient convergence for the system (Sigma_hat + lambda R): the number of iterations is O(n^{d/(2s+d)}) so the total solve cost is O(n^{2d/(2s+d)} log n).
Cite this review
Pith. "Pith review of Fast kernel methods: Sobolev, physics-informed, and additive models." pith.science (2026). https://pith.science/paper/N7Y5QTS7
@misc{pith2026250902649,
author = {Pith},
title = {Pith review of: Fast kernel methods: Sobolev, physics-informed, and additive models},
year = {2026},
howpublished = {\url{https://pith.science/paper/N7Y5QTS7}},
note = {Machine review of arXiv:2509.02649}
}
read the original abstract
Kernel methods are powerful tools in statistical learning, but their cubic complexity in the sample size n limits their use on large-scale datasets. In this work, we introduce a scalable framework for kernel regression with O(n log n) complexity, fully leveraging GPU acceleration. The approach is based on a Fourier representation of kernels combined with non-uniform fast Fourier transforms (NUFFT), enabling exact, fast, and memory-efficient computations. We instantiate our framework in three settings: Sobolev kernel regression, physics-informed regression, and additive models. When known, the proposed estimators are shown to achieve minimax convergence rates, consistent with classical kernel theory. Empirical results demonstrate that our methods can process up to tens of billions of samples within minutes, providing both statistical accuracy and computational scalability. These contributions establish a flexible approach, paving the way for the routine application of kernel methods in large-scale learning tasks.
Figures
Figures from the paper (3 more)
Forward citations
Cited by 1 Pith paper
-
PIKS: Universal Physics-Informed Kernel Methods
PIKS is universally consistent for linear differential constraints: with universal kernels it recovers both the target and the physics residual in the misspecified regime, with rates under source conditions.
Reference graph
Works this paper leans on
-
[2]
+ 3E(∥ √ Σ( ˆΣ + λR)−1Φ∗ϵ/n∥2 2) + 3E(∥ √ Σ( ˆΣ + λR)−1Φ∗δY/n∥2 2), where we used the fact that (x1 + x2 + x3)2 ≤ 3x2 1 + 3x2 2 + 3x2
-
[3]
The first term E(∥ √ Σ( ˆΣ + λR)−1λRθ∥2
-
[4]
The second term E(∥ √ Σ( ˆΣ + λR)−1Φ∗ϵ/n∥2
is a bias term depending on the regularization λR. The second term E(∥ √ Σ( ˆΣ + λR)−1Φ∗ϵ/n∥2
-
[5]
The last term E(∥ √ Σ( ˆΣ + λR)−1Φ∗δY/n∥2
is a variance term depending on the noise ϵ. The last term E(∥ √ Σ( ˆΣ + λR)−1Φ∗δY/n∥2
-
[6]
In what follows, we will successively bound these three errors terms
is another approximation error term, which measures the impact of the inter- action between the regularization R and the approximation errors δY. In what follows, we will successively bound these three errors terms. Bias term. According to Bach (Lemma 7.1, 2024), one has E(∥ √ Σ( ˆΣ + λR)−1λRθ∥2
work page 2024
-
[7]
Observe that R1/2 = S for the Sobolev regression and R1/2 = I for the low-bias Sobolev regression
= λ2E(∥ √ ΣR−1/2(R−1/2 ˆΣR−1/2 + λ)−1R1/2θ∥2 2) ≤ λ(1 + α2 λn )2⟨R1/2θ, R−1/2ΣR−1/2(R−1/2ΣR−1/2 + λ)−1R1/2θ⟩ ≤ λ(1 + α2 λn )2∥R1/2θ∥2 2, where α = max x∈Rd ∥R−1/2ϕ(x)∥2. Observe that R1/2 = S for the Sobolev regression and R1/2 = I for the low-bias Sobolev regression. Thus, α2 = P k∈Zd 1 1+∥k∥2s 2 < ∞ for the Sobolev regression, while α2 = (2m + 1)d for t...
work page 2024
-
[8]
≤ σ2n−1(1 + α2 λn )tr(E((R−1/2ΣR−1/2 + λ)−1R−1/2ΣR−1/2)) ≤ σ2n−1(1 + α2 λn )(2m + 1)d. Approximation term. Let (Xn+1, Yn+1) be a new sample, drawn independently from (X1, Y1), . . . ,(Xn, Yn). Let C = n+1X j=1 R−1/2ϕ(Xj)ϕ(Xj)⋆R−1/2, so that C = nR−1/2 ˆΣR−1/2 + R−1/2ϕ(Xn+1)ϕ(Xn+1)⋆R−1/2. The rationale behind this ap- proach is that (Xn+1, Yn+1) relates to...
work page 2024
-
[9]
≤ (1 + α2 λn )2 E(∥f ⋆ − fθ∥2 L2(PX )). Conclusion. Putting everything together, we deduce from inequality 3 that E(∥fˆθ − f ⋆∥2 L2(PX )) ≤ inf θ∈Hm (2 + 6(1 + α2 λn )2)E(∥f ⋆ − fθ∥2 L2(PX )) + λ(1 + α2 λn )2∥R1/2θ∥2 2 + σ2n−1(1 + α2 λn )(2m + 1)d, where • α2 = P k∈Zd 1 1+∥k∥2s 2 < ∞ for the Sobolev regression, while α2 = (2 m + 1)d for the low-bias Sobol...
Show all 14 references
-
[10]
Therefore, in the Sobolev regression, setting λ = O(n−2s/(2s+d)) such that λ ≥ n−1, and taking m = n1/(2s+d) leads to Sobolev mini- max rate
≤P k∈CZd (1 + ∥k∥2s 2 )|θ(f ⋆)k|2 ≤ CΩ∥f ⋆∥2 H s(Ω). Therefore, in the Sobolev regression, setting λ = O(n−2s/(2s+d)) such that λ ≥ n−1, and taking m = n1/(2s+d) leads to Sobolev mini- max rate. Moreover, in the low-bias Sobolev regression, setting λ = Θ( n−2s/(2s+d)) and taki...
-
[11]
Similarly, the bias term is bounded by E(∥ √ Σ( ˆΣ + λ)−1λθ∥2
+ 3E(∥ √ Σ( ˆΣ + λ)−1Φ∗ϵ/n∥2 2) + 3E(∥ √ Σ( ˆΣ + λ)−1Φ∗δY/n∥2 2), where Σ = E(ϕ(X)ϕ(X)⋆). Similarly, the bias term is bounded by E(∥ √ Σ( ˆΣ + λ)−1λθ∥2
-
[12]
As for the variance term, it is bounded by E(∥ √ Σ( ˆΣ + λ)−1Φ∗ε/n∥2
≤ λ(1 + α2 λn )2∥θ∥2 2, where α = maxx∈Rd ∥ϕ(x)∥2 = p d(2m + 1). As for the variance term, it is bounded by E(∥ √ Σ( ˆΣ + λ)−1Φ∗ε/n∥2
-
[13]
Finally, the approximation term is bounded by E(∥ √ Σ( ˆΣ + λR)−1Φ∗δY/n∥2
≤ σ2n−1(1 + α2 λn )tr(E((Σ + λ)−1Σ)) ≤ σ2 d(2m + 1) n (1 + α2 λn ). Finally, the approximation term is bounded by E(∥ √ Σ( ˆΣ + λR)−1Φ∗δY/n∥2
-
[14]
Putting everything together, we obtain E(∥fˆθ − f ⋆∥2 L2(PX )) ≤ inf θ∈Cd(2m+1) (2 + 6(1 + α2 λn )2)E(∥f ⋆ − fθ∥2 L2(PX )) + λ(1 + α2 λn )2∥θ∥2 2 + σ2 α2 n (1 + α2 λn )
≤ (1 + α2 λn )2 E(∥f ⋆ − fθ∥2 L2(PX )). Putting everything together, we obtain E(∥fˆθ − f ⋆∥2 L2(PX )) ≤ inf θ∈Cd(2m+1) (2 + 6(1 + α2 λn )2)E(∥f ⋆ − fθ∥2 L2(PX )) + λ(1 + α2 λn )2∥θ∥2 2 + σ2 α2 n (1 + α2 λn ). Convergence rate. The convexity of the function x 7→ x2 leads to E(...
-
[2024]
Nathan Doum`eche, Francis Bach, ´Eloi Bedek, G´erard Biau, Claire Boyer, and Yannig Goude
URL https://arxiv.org/abs/2409.13786. Nathan Doum`eche, Francis Bach, ´Eloi Bedek, G´erard Biau, Claire Boyer, and Yannig Goude. Fore- casting time series with constraints, 2025. URL https://arxiv.org/abs/2502.10485. Gene Howard Golub and Charles F. Van Loan.Matrix computation...
2025
Reviewed August 5, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.