Kernel Regression in Structured Non-IID Settings: Theory and Implications for Denoising Score Learning

Thanks for your time and efforts reviewing our paper. We now address raised questions as follows.

• Q1: The target function $f^*$ lies in the interpolation space $[\mathcal{H}]^s$ with $s\geq 1$ may be restrictive. / Discuss the challenges associated with analyzing the case $s < 1$, and whether these challenges are analogous to those in the i.i.d. setting.

  Thanks for your comment. The smoothness condition on $f_\rho^*$ appears in some existing work (Assumption (B) in [1], Lemma A.12 in [4]). In fact, for the derivation of our general theoretical bound (the first statement in Theorem 4.1, Eq.(2)), we can relax our assumption to $s>0$ and obtain exactly the same result, as we leverage assumptions and properties of interpolation spaces for refinements on the technical level.

  While for the specified bounds under conditional orthogonality (the second statement in Theorem 4.1, Eq.(3) and Theorem 4.4), $s\geq 1$ is required to estimate the relevance parameter $r_T$ (Lemma B.11) for providing a concise bound and clear insights, where the relevance parameter is explicitly related to $\alpha_t$. Technically, the smoothness on $f_\rho^*$ enables continuity to convert the relevance in the function space into the relevance $r_0$ in the data space (line 813-815). We consider this specific case to make our conclusion more clear and understandable. Indeed, without the strong smoothness $s\geq 1$, we can also provide estimation (less concise expression) for $r_T$. We merely need to replace $r_0$ with $r_\rho:=\frac{\mathbb{E}\left[\left(f _ \rho^{(r) * }(g_{ij})-f_\rho^{(r)* }(g_{i'j})\right)^2\right]}{4\mathbb{E}\left[f_\rho^{(r)* }(g_{ij})^2\right]}\in [0,1]$, maintaining all convergence guarantees.

• Q2: The experiments involve a three-layer ReLU MLP. Since there is currently no theoretical approximation framework between MLPs and kernel gradient flow (KGF) / KRR in the non-i.i.d. setting, it is better to use KRR instead.

  We thank the reviewer for the insightful comment regarding the use of MLPs. To supplement，we conducted the same experiments using both NTK and RBF kernels. The results show consistent trends with our MLP findings:

  k	NTK Loss at t = 1.0	NTK Loss at t = 0.1	RBF Loss at t = 1.0	RBF Loss at t = 0.1
  1	8.3379e-05	0.7020	0.00268199	0.6412
  16	8.3086e-05	0.6971	0.00175084	0.6532
  32	8.0374e-05	0.7268	0.00123094	0.8061
  64	8.0684e-05	0.7307	0.00102835	1.3014
  128	7.9432e-05	0.8351	0.00069417	2.6334

  These results reinforce both our theoretical insights and the empirical findings from the MLP experiments presented in the paper: increasing $k$ consistently improves performance at $t = 1.0$, where the target function is simpler, but leads to degraded accuracy at $t = 0.1$, where the function exhibits greater complexity. These empirical findings align well with our theory that when $t$ is large (i.e., noise dominates), increasing $k$ is beneficial to generalization.

• Q3: What is the motivation of setting $r \leq d$ in line 172? I notice that $r$ systems are trained independently, and both the theoretical analysis and the example in Section 4.1 consider the special case $r = 1$.

  Thanks for your comment. Considering $r = 1, \dots, d$ allows our framework to handle real-world vector-output tasks like denoising score learning (where outputs correspond to spatial/gradient dimensions). This generalizes prior scalar-output analyses to broader applications. While on the technical level, the vector output analysis can be reduced to $r=1$ for the reason that noise is independent per dimension. This enables provably equivalent analysis via parallel single-dimensional KRR. Focusing on $r=1$ in analysis maintains clarity in proofs and notation, and preserves all key insights.

• Q4: There might be typos in (2) and possibly other equations: The bias convergence rate should be $n^{-\min (s, 2)\theta}$; and the $\Theta$ notation should be $O$ unless both upper and lower bounds are established.

  Thanks for your comment. The bias convergence rate is indeed $n^{-\min (s, 2)\theta}$, we will modify these typos and the usage of asymptotic notations in the revised manuscript.

• Q5: A more comprehensive discussion of relevant prior work.

  Thanks for your suggestion. We will include the following discussion of related prior work in the revised manuscript.

  • (i) Optimal rate of KRR and kernel gradient flow.

    There are rich literature showing the optimal rate of KRR and kernel gradient flow. [3] demonstrates the minimax optimality of KRR when the regression lies in the RKHS and proves the best upper bound of the generalization error is $n^{-\frac{2\beta}{2\beta+1}}$. [4] further extends the result to the misspecified cases when the regression function does not lie in the RKHS. More recently, [5] studies the generalization error curves of a large class of analytic spectral algorithms, including the kernel gradient method, where the regressor $\hat{f}_\lambda$ is obtained by a filter function. Our paper initials the study of the generalization behavior of KRR in the structured non-i.i.d. setting and defers the study of analytic spectral algorithms to future.

  • (ii) Saturation effect of KRR.

    The saturation effect refers to the phenomenon that the kernel ridge regression (KRR) fails to achieve the information theoretical lower bound when the smoothness of the underground truth function exceeds certain level. To be specific, [2] shows that the information theoretical lower bound of the generalization error is $n^{-\frac{s\beta}{s\beta+1}}$ while [3] shows that the best upper bound of the generalization error is $n^{-\frac{2\beta}{2\beta+1}}$. This gap has been widely observed in practices: no matter how carefully one tunes the KRR, the rate of the generalization error can not be faster than $n^{-\frac{2\beta}{2\beta+1}}$. More recently, [1] rigoriously prove the saturation effect by providing a lower bound matching the existing upper bound $n^{-\frac{2\beta}{2\beta+1}}$. In our paper, let $\tilde{r}+\frac{1-\tilde{r}}{k} \asymp n^{-l}$, then the optimal $\theta$ for generalization error is $\theta_{opt}=\frac{(1+l)\beta}{\beta \mathrm{min}(s,2)+1}$ and the optimal generalization error rate is $n^{-\frac{(1+l)\beta \mathrm{min}(s,2)}{\beta \mathrm{min}(s,2)+1}}$. On the one hand, when $l>0$ (implying that $\tilde{r}$ is extremely small and $k$ is extremely large), the saturation effect still holds but with a faster rate. On the other hand, when $l=0$ (e.g. k=1), our result matches the previous result in the i.i.d. setting.

  • (iii) high-dimensional KRR and KGF.

    While kernel regression with a fixed data dimension $d$ has been extensively studied, there are rising attentions to the case when dealing with large-dimensional data. In the large-dimensional setting where $n \asymp d^\gamma$ with $\gamma>0$, [6] finds that for the square-integrable regression function, KRR and kernel gradient flow are consistent if and only if the regression function is a polynomial with a low degree. Another line of work e.g. [7], considers the overfitting behavior of kernel interpolation. On the technical level, [7] utilizes the argument that high-dimensional random kernel matrices can be approximated in spectral norm by linear kernel matrices plus a scaled identity. Regarding the saturation effect, [8] explores this effects for a large class of spectral algorithms (including the KRR, gradient descent, etc.) in large dimensional settings, by providing minimax lower bound and the exact convergence rates. Besides, [9] explores the Pinsker bound problem for kernel regression models that incorporate large-dimensional inner product kernels defined on the sphere $\mathbb{S}^d$, by addressing the scenario where the sample size $n$ is given by $\alpha d^\gamma (1+o_d(1))$ for some $\alpha,\gamma >0$. Our paper defers the analysis in high-dimensional data scnario to future.

We hope above response can address your concern and we are open to discuss more if any question still hold.

[1] ON THE SATURATION EFFECT OF KERNEL RIDGE REGRESSION.

[2] Optimal rates for the regularized learning algorithms under general source condition.

[3] Optimal rates for the regularized least-squares algorithm.

[4] On the optimality of misspecified kernel ridge regression.

[5] GENERALIZATION ERROR CURVES FOR ANALYTIC SPECTRAL ALGORITHMS UNDER POWER-LAW DECAY.

[6] Linearized Two-Layers Neural Networks in High Dimension.

[7] Just Interpolate: Kernel “Ridgeless” Regression Can Generalize.

[8] On the Saturation Effects of Spectral Algorithms in Large Dimensions.

[9] ON THE PINSKER BOUND OF INNER PRODUCT KERNEL REGRESSION IN LARGE DIMENSIONS.

Thanks for your time and efforts reviewing our paper. We now address raised questions as follows.

• Q1: Related works.

  Thanks for your question. We will add statistical learning of SVM [R1] and time-uniform concentration inequalities [R8-R10] in our revised manuscript and focus on discussing the non-i.i.d. analysis [R2,R3,R6,R7] here (deferring the discussion of learning under mixing assumptions [R4,R5,R11-R13] to Q2).

  • Setting

    A line of work established the consistency of SVMs or kernel methods under some specific assumptions e.g. processes satisfying a law of large numbers [R2], or empirical weak convergence [R3]. While the consistency can be characterized under such strong non-i.i.d.-ness, it remains unclear of the convergence speed. Another line of work focuses on the regression over trajectories generated by a dynamic system, both linear [R6] and non-linear [R7]. Though these approaches are without mixing, the surrogate trajectory assumptions limit its applicability to other settings. Our work aim to capture the asymptotic convergence rate for $k$-gap independent data.

  • Technique

    [R6,R7] utilize single-layer block decomposition and handle each blocks by Mendelson’s small-ball method, eschewing the use of standard mixing-time arguments. Regarding the specific $k$-gap independence, our technique relies on a new-designed block scheme which is free of neither Mendelson’s small-ball method nor standard mixing-time arguments. Our novel technique enables us to eliminate the approximation error between dependent blocks and their independent surrogates, revealing the benefit of data relevance (detailed in Q2).

• Q2: Setup.

  • Q2.1: Reason to consider the block-independence.

    Thanks for your question. We clarify the difference from those related works [R4,R5,R11-R13] as follows, whose techniques cannot resolve our setting.

    • Setting

      While many works have derived concentration high-probability bounds for non-i.i.d. data under the geometrically mixing assumption [R4,R5] or the algebraically mixing assumption [R11-R13], our $k$-gap independent case cannot be covered by these works, as the correlation between data points will remain high as long as they are from the same group with size $k$. Therefore, our $k$-gap independent case requires a fundamantally different non-i.i.d. data analysis.

    • Intuition

      Most of the existing results [R11-R13] perform analysis under the intuition to control the approximation error between independent blocks and dependent blocks by mixing coefficients. That being said, the data dependency is treated as an bad effect in their analysis and will finally appear as an error term. However, in our setting, we do not treat the data dependency a bad effect but aim to discover some benefits for reducing the error, otherwise it would be intractable to prove the vanishing generalization error in our setting for general $k$. To this end, we perform refined study in decomposing blocks to eliminate direct error approximation.

    • Insight

      Our analysis demonstrates the benefit of data relevance, which is in totally contrast to current non-i.i.d. data analysis. We believe some of our techniques and methodologies can be potentially extended to a broader settings (e.g., the size of each group could vary).

  • Q2.2: Proof techniques.

    Thanks for your question. Our proof technique differs fundamentally from mixing data approaches [R11-R13] in both decomposition strategy and error handling. While all methods employ block decomposition, the block scheme is different. Prior work does not perform refined study on the intra-block randomness, as they use single-layer decomposition that requires constructing independent surrogate blocks and controlling approximation errors. The crude study of the intra-block randomness fails to characterize the role of data dependence in our $k$-gap independent case. Besides, the size of block in existing methods fails to capture the underlying independence in our specific $k$-gap independent setting, limiting its applicability.

    In contrast, to overcome the above weakness, we introduce a more refined decomposition, yielding multiple partitions with several independent blocks respectively. This enables us to apply standard symmetrization techniques for each partition individually and resolve cross-partition dependencies through our novel incorporating lemma (Lemma C.20). Only by this technique, can we fully utilize the underlying data independence and intra-block dependence, and derive a preciser Bernstein-type high probability bound consistent with the i.i.d. case (up to some logarithm terms).

• Q3: Assumptions.

  • Q3.1: Kernel function / Convention(l. 178).

    Thanks for your good points. Operators related to $T$ are mappings: $L^2\to L^2$. We will provide clear characterizations of other concepts according to [48] in the revised manuscript.

  • Q3.2: The smoothness condition on $f_\rho^*$.

    Thanks for your comment. The smoothness condition on $f_\rho^*$ appears in some existing work (Lemma A.12 in [31]). In fact, In deriving our general theoretical bound Eq.(2), we can relax $s\geq 1$ to $s>0$ and obtain exactly the same result, as we have technically leveraged assumptions and properties of interpolation spaces for refinements.

    While for the specified bounds under conditional orthogonality (Eq.(3) and Theorem 4.4), $s\geq 1$ is required to estimate the relevance parameter $r_T$ (Lemma B.11) for providing a concise bound and clear insights, where the relevance parameter is explicitly related to $\alpha_t$. Technically, the smoothness on $f_\rho^*$ enables continuity to convert the relevance in the function space into the relevance $r_0$ in the data space (line 813-815). We consider this specific case to make our conclusion more clear and understandable. Indeed, without the strong smoothness $s\geq 1$, we can also provide estimation (less concise expression) for $r_T$. We merely need to replace $r_0$ with $r_\rho:=\frac{\mathbb{E}\left[\left(f_\rho^{(r)* }(g_{ij})-f_\rho^{(r)* }(g_{i'j})\right)^2\right]}{4\mathbb{E}\left[f_\rho^{(r)* }(g_{ij})^2\right]}\in [0,1]$, maintaining all convergence guarantees.

    Regarding the terms $i^{-\frac{1}{2}}$ and $a_i\geq c>0$ in this assumption, we follow the expression in the existing result [4] to provide a direct and clear comparison with the i.i.d. setting.

  • Q3.3: The polynomial eigenvalue decay assumption.

    Thanks for your comment. The polynomial spectrum makes our bound clearer and facilitates direct comparison with established results in the i.i.d. setting [31]. This makes the core theoretical insights more accessible. While the assumption simplifies presentation, our framework is readily extensible to general spectra on the technical level. The key terms (e.g. Lemma C.3-C.5) can be expressed directly in terms of individual eigenvalues ($\lambda_1, \dots$) rather than the decay rate $\beta$. For instance, the norm $\Vert f_\lambda\Vert_{[\mathcal{H}]^\gamma}^2$ is fundamentally given by a series (shown below) that depends on the full eigenvalue sequence: $\Vert f_\lambda\Vert_{[\mathcal{H}]^\gamma}^2\asymp\sum_{i=1}^\infty\left(\frac{\lambda_i^p}{\lambda_i+\lambda}\right)^2i^{-1}.$ Therefore, the polynomial decay is primarily a tool for deriving clearer, more interpretable bounds without fundamentally limiting the scope of our technical approach.

• Q4: Application to denoising score learning / Design neural networks.

  Thanks for your question. Our theoretical bounds are derived on KRR to provide insights and inspirations for practical data usage and training. In particular, our theoretical findings that noise multiplicity $k$ should scale with the noise level align well with the intuition that when the noise dominates in the observation, practitioners should add more noise rather than the raw data. Therefore, we believe this finding holds whatever learner is and apply our theory on practical neural network training. / The main contributions and insights of our theoretical results are to provide guidance for data sampling instead of designing neural network to approximate KRR. With our theoretical insights on KRR and empirical validations on neural network training, we believe our work can inspire practical training broadly.

• Q5：Minor Comments

  Thanks for your comments. We apologize for some typos and misunderstandings. For some inappropriate expressions such as "signal", "advanced technique", "roughly", typos "$f_\rho^* \in \mathcal{H}$" and some confused notations such as the asymptotic notations, $\lambda$ in Theorem 1, $\mu_G$, Example 3.2, $g$, we will modify those in the revised manuscript. We clarify other misunderstandings as follows.

  • Q5.1: KRR is applied independently per dimension.

    Thanks for your comment. On the one hand, applying KRR independently per dimension yields cleaner proofs and more interpretable bounds (with standard i.i.d. results [31]) and insights. On the other hand, in practical denoising score learning example, noise is standardly assumed independent per dimension, where multi-output variant case can be reduced to ours.

  • Q5.2: Second moment terms.

    Thanks for your good points. $\sigma_\epsilon$ and $\sigma^2$ control the decay of the population noise conditionally on data, with the former controlling the exponential decay (sub-Gaussian norm) and the latter controlling the quadratic decay (second moment). In parallel, for technical reasons to analyze dependent data, we extend $\sigma_\epsilon$ and $\sigma$ which are single-data conditioning to $\sigma_{\epsilon_{1,2}}$ and $\sigma_G^2$ which are dependent-data conditioning.

We hope above response can address your concern and we are open to discuss more if any question still hold.

We are glad that our rebuttal has addressed most of your other concerns. Thank you once again for your thoughtful follow-up questions. We would like to take this chance to address raised questions as follows.

• Q1: Are you claiming that Theorem 4.1 holds under $s>0$?

  Thanks for your question. Under the relaxation to $s>0$, the first statement (Eq. (2)) holds true and the second statement (Eq.(3)) holds by replacing $r_0$ with $r_\rho:=\frac{\mathbb{E}\left[\left(f_\rho^{(r)* }(g_{ij})-f_\rho^{(r)* }(g_{i'j})\right)^2\right]}{4\mathbb{E}\left[f_\rho^{(r)* }(g_{ij})^2\right]}\in [0,1]$. Previous assumption $s\geq 1$ is used to estimate the relevance parameter $r_T$ (Lemma B.11) for providing a concise bound and clear insights.

• Q2: Can you elaborate on why your $k$-gap independence assumption is not covered by standard mixing assumptions?

  Thanks for your question. For simplicity, we take $\phi$-mixing as an example. We would like to clarify that $\phi$-mixing is a very general assumption, which asssumes $\lim_{i\to \infty} \phi(i)=0$. Our $k$-independent setting can be covered by such general mixing condition with $\phi(i)=0,\ i\geq k$. However, our analysis derives tighter concentration bound over $k$-gap dependent data compared with previous work, yielding a generalization bound that fully captures the benefit of data relevance.

  In the bias-variance excess risk analysis, the key to bound the variance term is to perform concentration for $V=\frac{1}{n^2k^2 }\left[\sum _ {i=1}^{nk} Z_i \right]^2=\frac{1}{n^2k^2 }\sum_ {i,j} Z_iZ_j$, where $Z_i$ are some correlated random variables. To bound the associated differences under $\phi$-mixing conditions, previous concentration techniques treat the data dependency as a bad effect, yielding the following concentration inequality (Theroem 8 in [R12]) for random variables $W_1,...,W_m$ satisfing mixing condition with $\phi(i)$: $P\left[|\sum_{i=1}^m W_i-E[\sum_{i=1}^m W_i]|\geq\epsilon\right]\leq 2\exp(\frac{-2\epsilon^2}{mc^2[1+\sum_{i=1}^{m}\phi(i)]^2}).$ Applying this result on $V$ and let $\phi(i)$ denote the mixing coefficient of the sequence $\\{Z_i Z_j\\} _ {i=1...,nk; \ j=1,...,nk}$, we have $P(|V-E[V]|\geq \epsilon)\leq 2\exp(\frac{-2\epsilon^2n^4k^4}{n^2k^2c^2[1+\sum_{i=1}^{n^2k^2}\phi(i)]^2})$, i.e., $|V-E[V]|\leq O_P(\frac{1+\sum_{i=1}^{n^2k^2}\phi(i)}{nk})$.

  In contrast, to uncover the precise dependency among data points, we develop a new concentration inequality for the random variables $W_1,\dots,W_m$ under our setting (i.e., $k$-gap dependent data) by (1) capturing the underlying data independence and (2) utilizing the intra-block randomness, $P\left[|\sum_{i=1}^m W_i-E[\sum_{i=1}^m W_i]|\geq\epsilon\right]\leq 2\exp(\frac{-\epsilon^2/2}{mc'v^2}),$where $c'$ is some constant and $v^2=\sup_{K\subset\\{1,...,m\\}}\frac{1}{{\rm Card}K}\mathbb{E}\left[\big(\sum_{i\in K} W_i-E[\sum_{i\in K} W_i]\big)^2\right]$(see Proposition D.3). Furthermore, to better use our concentration inequality and develop a sharper bound on the variance error, we decompose $V$ into $V_1,V_2$ and $V_3$ given the specific data structure and apply our technique separately to bound each concentration error $E_i=P(|V_i-E[V_i]|\geq \epsilon), i=1,2,3$. In particular, our results show that $E_1 \leq 2\exp(\frac{-\epsilon^2}{c''/n^3k^2})$, $E_2\leq 2\exp(\frac{-\epsilon^2 }{c''(k-1)^2/n^3k^2})$, $E_3\leq 2\exp(\frac{-\epsilon }{c''(kr_T+1-r_T)/nk})$, where $c''$ is some constant. Overall, the dominating error $2\exp(\frac{-\epsilon }{c''(kr_T+1-r_T)/nk})$ yields the high probability bound $|V-E[V]|\leq O_P(\frac{kr_T+1-r_T}{nk})$, where $r_T$ characterizes the correlation of data (see Definition 3.2).

  For comparison, note that $\phi(i)$ may be difficult to estimate precisely, we could merely bound $\phi(i)$ as a constant level under general condition (see [1]). In this sense, our high probability bound $O_P(\frac{kr_T+1-r_T}{nk})$ is tighter than the one derived from existing technique, i.e., $O_P(nk)$, which is nearly meaningless. Moreover, we notice that even the previous concentration technique is applied to our decomposition method, the final high probability bound $|V-E[V]| \leq O_P(\frac{(1+\sum_{i=1}^{k} \phi(i))^2}{nk})\leq O_P(\frac{k}{n})$ is yet looser than our current result $O_P(\frac{kr_T+1-r_T}{nk})$, as $r_T\in[0,1]$ (here $\phi(i)$ is the mixing coefficient of the sequence $\\{Z_i\\}_{i=1}^{nk}$).

  In conclusion, our analysis provides tighter concentration bound over $k$-gap dependent data compared with previous work and develops a novel decomposion of the variance error analysis, yielding a generalization bound that can precisely captures the data relevance. We believe our techniques can also be potentially extended to a broader settings and we will discuss the comparison between our technique and prior work in the revised version.

• Q3: Theorem 4.4 should mention that the algorithm used is KRR / The theory holds for KRR / The conclusions for NNs are only heuristic based on this understanding.

  Thanks for your comment. We will modify Theorem 4.4 by mentioning that the algorithm used is KRR and the conclusions for NNs are heuristic based on this understanding in the revised manuscript. In numerical experiment, we have added the experiments using both NTK and RBF kernels. The results show consistent trends with our MLP findings:

  k	NTK Loss at t = 1.0	NTK Loss at t = 0.1	RBF Loss at t = 1.0	RBF Loss at t = 0.1
  1	8.3379e-05	0.7020	0.00268199	0.6412
  16	8.3086e-05	0.6971	0.00175084	0.6532
  32	8.0374e-05	0.7268	0.00123094	0.8061
  64	8.0684e-05	0.7307	0.00102835	1.3014
  128	7.9432e-05	0.8351	0.00069417	2.6334

  These results reinforce both our theoretical insights and the empirical findings from the MLP experiments presented in the paper: increasing $k$ consistently improves performance at $t = 1.0$, where the target function is simpler, but leads to degraded accuracy at $t = 0.1$, where the function exhibits greater complexity. We will revise our theorem 4.4 and include these new experimental results in the revision for more consistent explanation.

We hope above response can address your concern and we are open to discuss more if any question still hold.

[1] On the Computation of Mixing Coefficients Between Discrete-Valued Random Variables.

Thanks for your time and efforts reviewing our paper. We now address raised questions as follows.

• Q1: The paper does not propose a new algorithm but focuses mainly on theoretical bounds, which is not a novel direction by itself given that excess bias-variance bounds for kernel methods have been widely studied.

  Thanks for your comment. While excess risk bounds for kernel methods are well-studied in i.i.d. settings, no prior works have studied the excess risk and even learnability in the structured non-i.i.d. setting which is ubiquitous and important in applications like diffusion models. This paper contributes the analysis to the novel structured non-i.i.d. setting. Our theoretical bound explicitly characterizes the learnability in the structured non-i.i.d. setting and demonstrates how multiple noisy observations ($k$) affect generalization.

  Beyond theoretical bounds, our theory also guides noise-sample pairing strategies when training diffusion models: for general $t$, using large $k$ when $t$ increases.

• Q2: The new Bernstein-type concentration inequality mainly builds on existing techniques, including the Intrinsic Dimension Lemma and standard $\phi(a) = e^a - 1$ tricks from related literature, which limits its originality.

  Thanks for your comments. The Intrinsic Dimension Lemma and standard $\phi(a) = e^a - 1$ tricks are used in the standard i.i.d. Bernstein inequality. However, for $k$-gap independent data which is non-i.i.d. in our setting, exisiting techniques cannot handle data dependence. Therefore, we introduce a refined block decomposition, yielding multiple partitions with several independent blocks respectively. The novel technique enables us to apply standard symmetrization techniques for each partition individually and resolve cross-partition dependencies through our novel incorporating lemma (Lemma C.20). Our technique can also extend to more general cases (e.g. varying $k$ per signal, approximately independent block gaps), serving as a new tool in dependent data analysis for machine learning community.

• Q3: Could the authors clarify what key steps / insights make this extension substantially novel compared to existing kernel generalization theory?

  Thanks for your question. Our extension advances kernel generalization theory through two key innovations. One key step lies in the development of the new Bernstein inequality tailored for KRR analysis in the structured non-i.i.d. setting. The other key step is providing additional concentration analysis for $V_3$ (line 547), where we utilize the truncation technique and the novel Bernstein inequality (Lemma B.9). Both the two steps are essential for excess risk analysis in dependent cases. / To existing kernel generalization theory, we provide novel insights for whether data dependencies benefit or hinder the generalization performance of KRR. Our theoretical bounds explicitly blends $\frac{1}{n}$ and $\frac{1}{nk}$, revealing a critical trade-off between relevance and noise sample size: when the correlation level $\tilde{r}$ is large, i.e., the signal dominates in the observed noisy data, increasing $k$ offers little benefits while increasing $k$ helps generalization when the noise component prevails.

• Q4: Clearer notation or explicit definitions for each variance term.

  Thanks for your good points. $\sigma_\epsilon$ and $\sigma^2$ control the decay of the population noise conditionally on data, with the former controlling the exponential decay (sub-Gaussian norm) and the latter controlling the quadratic decay (second moment). In parallel, for technical reasons to analyze dependent data, we extend $\sigma_\epsilon$ and $\sigma$ which are single-data conditioning to $\sigma_{\epsilon_{1,2}}$ and $\sigma_G^2$ which are dependent-data conditioning. We provide explicit definitions as follows. For $(g,y)\in \mathcal{G}\times \mathcal{Y}$, for any dimension $r=1,\dots, d$, denote $\epsilon^{(r)}:=y^{(r)}-f_\rho^{(r)* }(g)$. Then $\sigma_\epsilon$ is the upper bound of the sub-Gaussian norm of $\epsilon^{(r)}$ conditionally on $g$ and $\sigma_{\epsilon_{1,2}}$ is the upper bound of the sub-Gaussian norm of $\epsilon^{(r)}$ conditionally on $g,g'$ where $g\neq g'\in \mathcal{G}$. To be specific, \Vert\epsilon^{(r)}|g\Vert _ {\psi_2}=\inf\left\\{t>0{:}\ \mathbb{E}\left[\exp\left(\left(\epsilon^{(r)}\right)^2/t^2\right)|g\right]\leq2\right\\}\leq\sigma_\epsilon,\ g,g'\in \mathcal{G} \ \text{almost everywhere}. \left\Vert\epsilon^{(r)}|g,g'\right\Vert_{\psi_2}=\inf\left\\{t>0{:}\ \mathbb{E}\left[\exp\left(\left(\epsilon^{(r)}\right)^2/t^2\right)\big|g,g'\right]\leq2\right\\}\leq\sigma_{\epsilon_{1,2}}, \ g,g'\in \mathcal{G} \ \text{almost everywhere}. Similarly, $\sigma^2$ and $\sigma_G^2$ are the upper bounds of the second moment of $\epsilon^{(r)}$ conditionally on $g$ and $g,g'$ respectively. $\mathbb{E}\left[{\epsilon^{(r)}}^2|g\right]\leq\sigma^2,\ \mathbb{E}\left[{\epsilon^{(r)}}^2|g,g'\right]\leq \sigma_G^2, \ g,g' \in \mathcal{G} \ \text{almost everywhere}.$

• Q5: The paper states “Yes” for open access to data and code, but only synthetic data generation details are described; no actual code link or repository is provided, which limits full reproducibility.

  Thanks for your comment. As we can not provide any code link in the rebuttal stage, we will provide the actual code link or repository in the revised manuscript.

• Q6: The paper’s findings on the optimal noise sampling schedule provide some guidance for diffusion models. Could the authors comment on how these results might inspire practical training or data usage strategies more broadly?

  Thanks for your question. The optimal noise sampling schedule provides inspirations for practical training of denoising score learning. Consider training with a fixed batch size $nk$, for varying $t$ (or equivalently, $\alpha_t$), practitioners can adaptively choose $k$ for optimizing the training efficiency: if signal dominates, then setting $k=1$ is enough; while when noise dominates, practitioners are encouraged to increase $k$ to $\Theta\left((1-\alpha_t^{p/2})/\alpha_t^{p/2}\right)$. This adaptive design for noise multiplicity $k$ may advance the empirical study for denoising score learning.

  Additionally, our insights can also extend beyond denoising score learning to flow-based diffusion model (e.g. flow matching, rectified flow) and other noise schedules (e.g. Variance-Exploding (VE), Variance-Preserving (VP), sub-VP, Linear, EDM), unifying noise sampling principles across diffusion paradigms. Overall, our $k$-adaptation principle provides the first theoretical grounding for noise scheduling decisions - enabling systematic design rather than empirical guesswork.

We hope above response can address your concern and we are open to discuss more if any question still hold.

Thanks for your time and efforts reviewing our paper. We now address raised questions as follows.

• Q1: Reliance on strong spectral assumption.

  Thanks for your comment. We consider the polynomial spectrum to make our bound clearer and facilitates direct comparison with established results in the i.i.d. setting (consistent with prior works [1,2]), which can make the core theoretical insights more accessible. While the assumption simplifies presentation, our framework is readily extensible to general spectra on the technical level. The key terms (e.g. $tr(TT_\lambda^{-1})^p$, $\Vert f_\lambda^{(r)}\Vert_{[\mathcal{H}]^\gamma}^2$ and $\Vert f_\lambda^{(r)}-f_\rho^{*(r)}\Vert_{[\mathcal{H}]^\gamma}^2$ in Lemma C.3-C.5) can be expressed directly in terms of individual eigenvalues ($\lambda_1, \lambda_2, \dots$) rather than the decay rate $\beta$. For instance, the norm $\Vert f_\lambda\Vert_{[\mathcal{H}]^\gamma}^2$ is fundamentally given by a series (shown below) that depends on the full eigenvalue sequence: $\Vert f_\lambda\Vert_{[\mathcal{H}]^\gamma}^2\asymp\sum_{i=1}^\infty\left(\frac{\lambda_i^p}{\lambda_i+\lambda}\right)^2i^{-1}.$ Therefore, the polynomial decay is primarily a tool for deriving clearer, more interpretable bounds without fundamentally limiting the scope of our technical approach.

• Q2: Reliance on the Hölder continuity assumption.

  Thanks for your comment. The primary purpose of the Hölder continuity assumption is to eliminate the need for the often unrealistic sub-Gaussian design assumption in deriving our general bounds [1,2]. Technically, it is essential for establishing a uniform concentration bound via covering number estimates (Lemmas C.10 and C.11). Furthermore, this assumption allows us to derive concise estimates for the relevance parameter $r_T$ in specific scenarios, such as when conditional orthogonality holds. We note that Hölder continuity is naturally satisfied by important kernel classes like the Laplace kernel, Sobolev kernels, and Neural Tangent Kernels [1,2].

• Q3: Please clarify whether the blockwise concentration still holds if k varies per signal or if gaps are only approximately independent.

  Thanks for your question. The blockwise concentration remains valid when the number of noisy realizations $k_i$ varies per signal $x_i$. A simple method to handle is replacing $k$ with $k_\mathrm{max}:=\max_i k_i$. The key insight is that the entire training observations $G$ satisfies $k_{\mathrm max}$-gap independent (see definitions in line 146-148), preserving the concentration properties. If gaps are only approximately independent, e.g., some weakly dependent process assuming specific mixing property, our framework can also be readily extended to accomodate the dependency. The core technique is to incorporate the specific dependence structure (e.g., mixing coefficients) to modify the concentration bounds by adding terms quantifying block dependence. Established methods for handling such dependencies can be adapted from [3,4].

• Q4: No comparison with alternative dependent‑data analyses.

  Thanks for your comment. Current literature on dependent-data analyses rely heavily on mixing conditions, which characterize the correlation between random variables from certain distance in the sequence. Differently, we consider the $k$-gap independent data, where the random variables are categorized multiple groups with size $k$, the data point will be highly correlated/dependent with the data in the same group, but is independent of the data from other groups. Then we would like to explain the difference between our work and these related works in the following aspects:

  • Setting and Assumption

    While many works have derived tight concentration high-probability bounds for non-i.i.d. data under the geometrically regular mixing assumption [5,6] or the algebraically mixing assumption [7], which characterize the decay of correlation between random varaiables as the distance increases, our $k$-gap independent case cannot be covered by these works, as the correlation between data points will remain high as long as they are from the same group with size $k$. Therefore, our $k$-gap independent case requires a fundamantally different non-i.i.d. data analysis.

  • Analysis Intuition

    Most of the existing results [5-7] perform analysis under the intuition to control the approximation error between independent blocks and dependent blocks by mixing coefficients. That being said, the data dependency is treated as an bad effect in their analysis and will finally appear as an error term. However, in our setting, we do not treat the data dependency a bad effect but aim to discover some benefits for reducing the error, otherwise it would be intractable to prove the vanishing generalization error in our setting for general $k$. To this end, we perform refined study in decomposing blocks to eliminate direct error approximation.

  • Insight

    Our analysis demonstrating the benefit of data relevance stands in contrast to a long line of work on learning from dependent data, providing insights for understanding data dependence. We believe some of our techniques and methodologies can be potentially extended to a broader settings (e.g., the size of each group could vary).

  On the technical level, prior work [5-7] use single-layer decomposition that requires constructing independent surrogate blocks and carefully controlling approximation errors. On the contrary, we introduce a more refined decomposition, yielding multiple partitions with several independent blocks respectively, which enables us to apply standard symmetrization techniques for each partition individually and resolve cross-partition dependencies through our novel incorporating lemma (Lemma C.20).

• Q5: Could the relevance parameters be estimated from data to guide k adaptively in practice?

  Thanks for your question. In the case when $\alpha_t \to 0$, i.e., the data $x_t$ becomes more noisy, the optimal $k$ can be selected as $k^*=\Theta\left(\frac{1-\alpha_t^{p/2}}{\alpha_t^{p/2}}\right)$. For general $t$, the suggestion from our theory is to use larger $k$ when $t$ increases. Nevertheless, it still requries a bit hyper-parameter tuning for different practical tasks.

• Q6: How sensitive are the bounds to the Hölder index p; can slower continuity be handled?

  Thanks for your question. Under the specific conditional orthogonality, our general upper bound is $R(\lambda) \leq \tilde{\Theta} _ {\mathbb{P}}\left({n}^{-\operatorname*{min}(s,2)\theta}\right) + \tilde{\sigma}^2 O _ \mathbb{P}^{\rm poly}\left(n^{\alpha_0\theta}\left(\frac{\tilde{r} _ 0^{p/2}\vee r_e }{n}+\frac{(1-\tilde{r} _ 0^{p/2}) \wedge (1-r_e)}{nk}\right)\right)$. Overall, $p$ affects both the tightness of our bounds and the benefit of multiple observations ($k$). To be specific, larger $p$ implies the domination of $\frac{1}{nk}$, showing the underlying benefit to generalization by increasing $k$. In contrast, smaller $p$ (slow continuity) implies the generalization error is up to $O_\mathbb{P}^{\mathrm{poly}}(n^{\alpha_0\theta}/n)$.

• Q7: Only toy simulations on a 1‑D mixture that limit empirical support. Practical impact on diffusion training remains speculative without large‑scale experiments. An ablation on real image diffusion training would strengthen claims.

  Thanks for your comment. To demonstrate the applicability of our method beyond toy examples, we conducted an ablation study on the CIFAR-10 dataset. We trained a diffusion model using a dataset of 1024 samples for 100 epochs with a batch size of 1024, optimized using Adam with a learning rate of 2e-3. The model architecture was a two-layer U-Net, and time conditioning was implemented by expanding the time variable t and concatenating it as an additional input channel to the image.

  We report the diffusion loss on the test set with a size of 1024 at $t = 1.0$ and $t = 0.1$ across different values of $k: number of noisy realizations per data$. Each configuration was evaluated over 100 parallel runs to ensure robustness:

  k	Loss at t = 1.0	Loss at t = 0.1
  1	0.00268199	0.09571435
  16	0.00175084	0.09752702
  32	0.00123094	0.10304873
  64	0.00102835	0.12863263
  128	0.00069417	0.18158743

  As shown, increasing $k$ consistently improves performance at $t=1$, indicating better fitting of the complex score function. However, it also leads to degradation at $t=0.1$, consistent with our empirical findings on the MoG settings in the paper. Both our empirical findings on the MoG settings and CIFAR-10 align well with our theory that when $t$ is large (i.e., noise dominates), increasing $k$ is beneficial to generalization.

• Q8: What computational overhead would varying k incur?

  Thanks for your question. Varying $k$ will not incur any computational overhead, as we always conside a fixed batchsize in experiments. No matter how $k$ varies, the computation is the same as the case $k=1$.

We hope above response can address your concern and we are open to discuss more if any question still hold.

[1] On the asymptotic learning curves of kernel ridge regression under power-law decay.

[2] Kernel interpolation generalizes poorly.

[3] Bernstein-type inequality for a class of dependent random matrices.

[4] Bernstein inequality and moderate deviations under strong mixing conditions.

[5] Fast learning from non-iid observations.

[6] Fast learning from $\alpha$-mixing observations.

[7] Rates of convergence for empirical processes of stationary mixing sequences.