{"id":"cb8e9a65-6880-4a93-834d-0c514f3715e2","arxiv_id":"2509.02649","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":2,"one_line_summary":"A Fourier-basis plus NUFFT framework solves Sobolev, physics-informed, and additive kernel regressions in O(n log n) time on GPUs while retaining minimax rates.","lead":"This paper shows how to train certain kernel-based regression models on billions of data points in about a minute using a Fourier transform trick on GPUs. It also proves these fast models keep the same statistical accuracy as classical kernel methods for Sobolev and additive-model problems.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Prop. 2.1's imported Fourier norm equivalence for all bounded Lipschitz domains is the load-bearing assumption; if it fails, the m^{-2s} bias bound and the minimax rates in Props. 3.1/3.2/5.1 fail with it.","rationale":"I considered alternatives: the CG iteration count could be repaired by the standard fact that CG terminates in at most D iterations, giving the stated O(n^{2d/(2s+d)} log n) cost; the Toeplitz first-row issue is a matter of ordering/embedding with constants 2^d; the additive-model m inconsistency (m=n^{1/(2s+1)}/d vs m=n^{1/(2s+1)} in the proof) is a typo that does not change the rate for fixed d. None of these threatens the central argument as much as Prop. 2.1, which supports the bias term in every rate theorem. The reader's weakest_assumption already identifies this same point; I agree. If Prop. 2.1 is verified in the cited prior work, the theory is likely sound and the conditional can be lifted; if not, the rates are unsubstantiated. Therefore the reader's CONDITIONAL verdict is appropriate and should remain unchanged.","tokens_in":18698,"tokens_out":26568,"duration_ms":291254,"concrete_test":"Check the cited Doumeche et al. (2024) result: does it prove the two-sided Fourier norm bound for every bounded Lipschitz domain with C_Omega independent of m? If the source covers only rectangular/hypercube domains, the present Assumption (iv) is unsupported. Independent analytic check: for a non-rectangular Omega (e.g., the unit disk) construct the bounded extension operator and verify that for every f in H^s(Omega) the periodized coefficients satisfy sum_k |theta_k|^2 (1+||k||^{2s}) <= C_Omega ||f||^2_{H^s(Omega)}; if such a construction requires extra boundary trace conditions not stated in the paper, the minimax claims in Props. 3.1, 3.2, and 5.1 are not established.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Every theoretical rate in the paper funnels through Prop. 2.1, imported from Doumeche et al. (2024): for f in H^s(Omega) on a bounded Lipschitz Omega, there must exist Fourier coefficients theta(f) with sum_k |theta_k|^2 (1+||k||^{2s}) <= C_Omega ||f||^2_{H^s(Omega)} and with truncation error O(m^{-2s}) in L2(PX). This upper bound is used in the proofs of Props. 3.1/3.2/5.1 to control lambda||S theta||^2 and, via Prop. 2.2, the bias term E(||f* - f_theta||^2). If the equivalence constant depends on m, or the m^{-2s} decay fails for a general Lipschitz domain, the bias term is not O(n^{-2s/(2s+d)}) and the claimed minimax rates do not follow. The proposition is not a standard textbook fact: for Omega=(0,1), the naive periodic Fourier coefficients of f(x)=x have divergent Sobolev-weighted l2 norm, so the result depends on a nontrivial bounded-extension/periodization argument. The paper gives no proof and no theorem number in the citation, while Assumption (iv) claims arbitrary bounded Lipschitz domains. All experiments are on cubes, so they do not exercise the general-domain claim. This is the single most load-bearing gap.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":19128,"tokens_out":20253,"duration_ms":233894,"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":[{"comment":"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.","section":"§2.2, Prop. 2.1 and Appendix A.1"},{"comment":"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.","section":"§3.1, after Prop. 3.1; Assumption (v)"},{"comment":"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.","section":"§5, Prop. 5.1 and Appendix A.3"},{"comment":"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.","section":"Abstract and §2.3"}],"minor_comments":[{"comment":"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.","section":"Prop. 3.1, notation"},{"comment":"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'.","section":"Throughout"},{"comment":"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.","section":"Figures 1 and 4"},{"comment":"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.","section":"References"}],"recommendation":"major_revision","confidential_remarks":"The main risk is Proposition 2.1. Since the authors are also the authors of the cited COLT paper, asking for a self-contained proof or precise theorem reference is reasonable and should be straightforward. The empirical and algorithmic parts are strong, and the upper-bound proofs are otherwise standard. I do not see a reason for rejection, but the minimax claim under the current density assumption and the additive-model complexity statement need to be corrected."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Re: arXiv:2509.02649. My take: the algorithmic package is genuinely useful and the experiments are impressive, but the paper currently leans on two things it doesn't show. First, Proposition 2.1, imported from the authors' own COLT paper, states a Fourier series representation with Sobolev-weighted l2 equivalence for every bounded Lipschitz domain. That's not a standard textbook fact; for Omega=(0,1), the naive periodic coefficients of f(x)=x don't satisfy it. The paper gives no proof and no theorem number, and all experiments are on cubes, so the general-domain claim is untested. Since every rate argument (Props 3.1, 3.2, 5.1) funnels through the m^{-2s} bias decay derived from that lemma, this is a load-bearing gap. Second, the O(n log n) complexity claim assumes CG converges in a number of iterations that is not analyzed. Each matvec is O(m^d log m), but the condition number of the system depends on lambda and n; for lambda = n^{-2s/(2s+d)} it can grow polynomially, so the total cost might be more than O(n log n). The paper says 'total complexity is O(n^{2d/(2s+d)} log n)=o(n)' without justifying the iteration count. That needs a bound or a remark.\n\nWhat is good: the unified use of NUFFT, Toeplitz structure, and CG for the truncated Fourier problem is new in this combination, and the GPU additive model is a first. The ridge-regression proofs are standard and the choices of lambda and m do yield the claimed minimax rates for the finite-basis problem. The empirical scaling to 1e10 points in under a minute is convincing and the comparison against PyGAM is fair. The 'exact' in the abstract is a bit strong—it's exact for the truncated basis, not the full kernel—but they do clarify in Section 2.\n\nMinor: I didn't find the constant-factor slip the reader mentioned; the displayed bound in Prop 3.1 matches the appendix. Some typos (e.g., the additive model's m formula has a confusing negative exponent).\n\nWho this is for: anyone building large-scale kernel or PINN solvers. It deserves a serious referee round, but the authors need to either prove Prop 2.1 (or give a precise theorem number) and address the CG iteration count. If those are fixed, this is a strong paper.","headline":"Useful algorithmic framework, but the theory leans on an unproved Fourier norm-equivalence lemma and an unstated CG iteration assumption.","tokens_in":19605,"tokens_out":5677,"would_cite":true,"duration_ms":62395,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62G08","62G20","65T50"],"pacs":[],"model":"deepseek-v4-flash","headline":"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)}.","keywords":["kernel regression","non-uniform fast Fourier transform (NUFFT)","Sobolev spaces","minimax convergence rates","physics-informed regression","additive models","GPU acceleration","Toeplitz matrices"],"falsifier":"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.","tokens_in":18640,"feed_emoji":"⚡","tokens_out":8114,"duration_ms":79160,"temperature":0.7,"pith_summary":"Kernel ridge regression traditionally costs O(n^3) time and O(n^2) memory, which is why it is rarely used on large data. This paper shows that when the kernel space is represented on a truncated Fourier basis, the entire estimation pipeline—forming the covariance, multiplying by it, and inverting the regularized system—can be done in O(n log n) time and memory with GPU-friendly non-uniform fast Fourier transforms (NUFFT). The result is not an approximation: the Fourier-truncated estimator is fit exactly, and the truncation level is chosen so the statistical error remains minimax-optimal for Sobolev spaces. The authors instantiate the framework for Sobolev regression, physics-informed regression with differential constraints, and additive models, proving convergence rates of O(n^{-2s/(2s+d)}) and O(n^{-2s/(2s+1)}) when known. Empirically the algorithms process up to tens of billions of samples in minutes on one GPU.","feed_headline":"Fourier trick makes kernel regression O(n log n), not O(n^3)","feed_subtitle":"Matching minimax rates with tens of billions of samples processed in minutes on one GPU.","key_machinery":"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.","core_discovery":"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","pith_inferences":["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."],"forward_implications":["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."],"supporting_citations":[{"why":"Supplies the Fourier-series representation and norm equivalence (Proposition 2.1) on which the bias and truncation analysis rests.","marker":"Doumèche et al. (2024)"},{"why":"Supplies the GPU NUFFT implementation used to compute covariance rows and the covariance vector in O(n log n).","marker":"Shih et al. (2021)"},{"why":"Provides the closed-form kernel estimator, the classical O(n^3) baseline, and the ridge-regression lemmas used in the proofs of Propositions 3.1 and 3.2.","marker":"Bach (2024, Chapter 7)"},{"why":"Provides the elementary ridge-regression analysis that the generalization-error decomposition follows.","marker":"Mourtada & Rosasco (2022)"},{"why":"Provides the Toeplitz matrix-vector multiplication result behind the FFT-based fast application of the covariance matrix.","marker":"Golub & Loan (2013)"},{"why":"Generalizes fast recursive block-Toeplitz multiplication to higher dimensions.","marker":"Lee (1986)"},{"why":"Gives the minimax lower bound n^{-2s/(2s+d)} used to identify the rates as optimal.","marker":"Tsybakov (2009)"},{"why":"Establishes existence of minimax Sobolev estimators, the standard against which the Fourier estimator is compared.","marker":"Nemirovski (2000)"},{"why":"Provides the large-data generalized additive model baseline whose cost and memory requirements are compared in Section 5.","marker":"Wood et al. (2014)"}],"fun_headline_variants":["Kernel regression breaks cubic barrier: O(n log n) via Fourier NUFFT","Exact O(n log n) kernel regression on tens of billions of samples","Minimax rates in minutes: kernel regression scales to billions","Sobolev kernel regression at O(n log n) with Fourier-based NUFFT","Kernel methods go billion-scale: Fourier acceleration delivers minimax rates"],"cache_read_input_tokens":2688,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Kernel regression breaks cubic barrier: O(n log n) via Fourier NUFFT","Exact O(n log n) kernel regression on tens of billions of samples","Minimax rates in minutes: kernel regression scales to billions","Sobolev kernel regression at O(n log n) with Fourier-based NUFFT","Kernel methods go billion-scale: Fourier acceleration delivers minimax rates"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000688,"raw_usage":{"total_tokens":2945,"prompt_tokens":727,"completion_tokens":2218,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":471,"completion_tokens_details":{"reasoning_tokens":2117}},"tokens_in":471,"tokens_out":2218,"duration_ms":15109,"temperature":1.0,"reasoning_tokens":2117,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-05T11:44:02.964981+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[],"review_version":1}