Evaluating the design space of diffusion-based generative models

Authors: yuqing  wang
ye  he
molei  tao
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 15:31:55

Abstract: Most existing theoretical investigations of the accuracy of diffusion models, albeit significant, assume the score function has been approximated to a certain accuracy, and then use this a priori bound to control the error of generation. This article instead provides a first quantitative understanding of the whole generation process, i.e., both training and sampling. More precisely, it conducts a non-asymptotic convergence analysis of denoising score matching under gradient descent. In addition, a refined sampling error analysis for variance exploding models is also provided. The combination of these two results yields a full error analysis, which elucidates (again, but this time theoretically) how to design the training and sampling processes for effective generation. For instance, our theory implies a preference toward noise distribution and loss weighting in training that qualitatively agree with the ones used in [Karras et al., 2022]. It also provides perspectives on the choices of time and variance schedules in sampling: when the score is well trained, the design in [Song et al., 2021] is more preferable, but when it is less trained, the design in [Karras et al., 2022] becomes more preferable.

Main paper: 
 introduction: 
 Diffusion models became a very popular generative modeling approach in various domains, including computer vision [20, 7, 27, 28, 38, 51] , natural language processing [6, 34, 37] , various modeling tasks [15, 41, 55] , and medical, biological, chemical and physical applications [3, 17, 43, 49, 23, 56] (see more surveys in [53, 11, 14] ). Karras et al. [30] provided a unified empirical understanding of the derivations of model parameters, leading to new state-of-the-art performance. Karras et al. [31] further upgraded the model design by revamping the network architectures and replacing the weights of the network with an exponential moving average. As diffusion models gain wider usage, efforts to understand and enhance their generation capability become increasingly meaningful.In fact, a rapidly increasing number of theoretical works already analyzed various aspects of diffusion models [32, 19, 52, 16, 12, 8, 18, 9, 13, 39, 44, 25] . Among them, a majority [32, 19, 52, 16, 12, 8, 18] focus on sampling/inference; more precisely, they assume the score error is within a certain accuracy threshold (i.e. the score function is well trained in some sense), and analyze the discrepancy between the distribution of the generated samples and the true one. Meanwhile, there are a handful of results [44, 25] that aim at understanding different facets of the training process. See more detailed discussions of existing theoretical works in Section 1.1.However, as indicated in Karras et al. [30] , the performance of diffusion models also relies on the interaction between design components in both training and sampling, such as the noise distribution, weighting, time and variance schedules, etc. While focusing individually on either the training or generation process provides valuable insights, a holistic quantification of the actual generation capability can only be obtained when both processes are considered altogether. Therefore, motivated by obtaining deeper theoretical understanding of how to maximize the performance of diffusion models, this paper aims at establishing a full generation error analysis, combining both the optimization and sampling processes, to partially investigate the design space of diffusion models.More precisely, we focus on the variance exploding setting [46] , which is also the foundation of continuous forward dynamics in Karras et al. [30] . Our main contributions are summarized as follows:• For denoising score matching objective, we establish the exponential convergence of its gradient descent training dynamics (Theorem 1). We develop a new method for proving a key lower bound of gradient under the semi-smoothness framework [1, 35, 57, 58 ]. • We extend the sampling error analysis in [8] to the variance exploding case (Theorem 2), under only the finite second moment assumption (Assumption 3) of the data distribution. Our result applies to various variance and time schedules, and implies a sharp almost linear complexity in terms of data dimension under optimal time schedule. • We conduct a full error analysis of diffusion models, combining training and sampling (Theorem 3).• We qualitatively derive the theory for choosing the noise distribution and weighting in the training objective, which coincides with Karras et al. [30] (Section 4.1). More precisely, our theory implies that the optimal rate is obtained when the total weighting exhibits a similar "bell-shaped" pattern used in Karras et al. [30] . • We develop a theory of choosing time and variance schedules based on both training and sampling (Section 4.2). Indeed, when the score error dominates, i.e., the neural network is less trained and not very close to the true score, polynomial schedule [30] ensures smaller error; when sampling error dominates, i.e., the score function is well approximated, exponential schedule [46] is preferred. Conclusions and limitations are in Appendix A.Sampling. There has been significant progress in quantifying the sampling error of the generation process of diffusion models, assuming the score function is already approximated within certain accuracy. Most existing works [e.g., 16, 12, 8] focused on the variance preserving (VP) SDEs, whose discretizations correspond to DDPM. For example, Benton et al. [8] is one of the latest results for the VPSDE-based diffusion models, and it only needs a very mild assumption: the data distribution has finite second moment. The iteration complexity is shown to be almost linear in the data dimension and polynomial in the inverse accuracy, under exponential time schedule. However, a limited amount of works [32, 24, 54] analyzed the variance exploding (VE) SDEs, whose discretizations correspond to Score Matching with Langevin dynamics (SMLD) [45, 46] . To our best knowledge, Yang et al. [54] obtained the best result so far for VE assuming the data distribution has bounded support: the iteration complexity is polynomial in the data dimension and the inverse accuracy, under the uniform time schedule. In contrast, our work only assumed that the data distribution has finite second moment, and by extending the stochastic localization approach in [8] to VESDE, we obtain an iteration complexity that is polynomial in the data dimension and the inverse accuracy, under more general time schedules as well. Note the improved complexity in terms of the inverse accuracy and the data dimension dependencies; in fact, under the exponential time schedule, our complexity is almost linear in the data dimension, which recovers the state-of-the-art result for VPSDE-based diffusion models.Training. To our best knowledge, the only works that quantify the training process of the diffusion models are Shah et al. [44] and Han et al. [25] . Shah et al. [44] employed the DDPM formulation and considered data distributions as mixtures of two spherical Gaussians with various scales of separation, together with K spherical Gaussians with a warm start. Then the score function can be analytically solved, and they modeled it in a teacher-student framework solved by gradient descent. They also provided the sample complexity bound under these specific settings. In contrast, our results work for general data distributions for which the true score is unknown, and training analysis is combined with sampling analysis. Han et al. [25] considered the GD training of a two-layer ReLU neural network with the last layer fixed, and used the neural tangent kernel (NTK) approach to establish a first result on generalization error. They uniformly sampled the time points in the training objective, assumed that the Gram matrix of the kernel is away from 0 (implying a lower bound on the gradient), and lacked a detailed non-asymptotic characterization of the training process. In contrast, we use the deep ReLU network with L layers trained by GD and prove instead of assuming that the gradient is lower bounded by the objective function. Moreover, we obtain a non-asymptotic bound for the optimization error, and our bound is valid for general time and variance schedules, which allows us to obtain a full error analysis.Convergence of neural networks training. The convergence analysis of neural networks under gradient descent has been a longstanding challenge and has been developed into an extensive field.Here we will only focus on results mostly related to the techniques used in this paper. One line of them is approaches directly based on neural tangent kernel (NTK) [22, 21, 5, 47, 36] . However, existing works in this direction focus more on either scalar output, or vector output but with only one layer trained under two-layer networks, which is insufficient for diffusion models. Another line of research also considers overparameterized models in a regime analogous to NTK, though not necessarily explicitly resorting to kernels. Instead, it directly quantifies the lower bound of the gradient [1, 35, 2, 57, 58] and uses a semi-smoothness property to prove exponential convergence.Our results align with the latter line, but we develop a new method for proving the lower bound of the gradient and adopt assumptions that are closer to the setting of diffusion models. See more discussions in Section 3.1.We denote ∥ ⋅ ∥ to be the ℓ 2 norm for both vectors and matrices, and ∥ ⋅ ∥ F to be the Frobenius norm.For the discrete time points, we use t i to denote the time point for forward dynamics and t ← i for backward dynamics. For the order of terms, we follow the theoretical computer science convention to use O(⋅), Θ(⋅), Ω(⋅). We also denote f ≲ g if f ≤ Cg for some universal constant C. 
basics of diffusion-based generative models: 
 In this section, we will introduce the basic forward and backward dynamics of diffusion models, as well as the denoising score matching setting under which a model is trained.Consider a forward diffusion process that pushes an initial distribution P 0 to GaussianEQUATIONwhere dW t is the Brownian motion, X t is a d-dim. random variable, and X t ∼ P t . Under mild assumptions, the process can be reversed and the backward process is defined as followsEQUATIONwhere Y 0 ∼ P T , and p t is the density of P t . Then Y T -t and X t have the same distribution with density p t [4] , which means the dynamics (2) will push (near) Gaussian distribution back to (nearly) the initial distribution P 0 . To apply the backward dynamics for generative modeling, the main challenge lies in approximating the term ∇ log p T -t (Y t ) which is called score function. It is common to use a neural network to approximate this score function and learn it via the forward dynamics (1); then, samples can be generated by simulating the backward dynamics (2).In order to learn the score function, a natural starting point is to consider the following score matching objective [e.g., 29]EQUATIONwhere S(θ; t, X t ) is a θ-parametrized neural network, w(t) is some weighting function, and the subscript means this is the continuous setup. Ideally one would like to minimize this objective function to obtain θ; however, p t in general is unknown, and so is the true score function ∇ x log p t . One of the solutions is denoising score matching proposed by Vincent [48] , where one, instead of directly matching the true score, leverages conditional score for which initial condition is fixed so that p t|0 is analytically known.More precisely, given the linearity of forward dynamics (1), its exact solution is explicitly known: Let µ t = ∫ t 0 f s ds, and σ2 t = 2 ∫ t 0 e 2µs-2µt σ 2 s ds. Then the solution is X t = e -µt X 0 + σt ξ, where ξ ∼ N (0, I). We also have X t |X 0 ∼ N (e -µt X 0 , σ2 t I) and g t (x|y) = (2πσ 2 t ) -d/2 exp(-∥xe -µt y∥ 2 /(2σ 2 t )), which is the density of X t |X 0 . Then the objective can be rewritten asEQUATIONwhereC t = E Xt ∥∇ log p t ∥ 2 -E X0 E Xt|X0 ∥∇ log g t (X t |X 0 )∥ 2 .For completeness, we will provide a detailed derivation of these results in Appendix C and emphasize that it is just a review of existing results in our notation.Throughout this paper, we adopt the variance exploding (VESDE) setting [46] , where f t = 0 and hence µ t = 0, which also aligns with the setup in the classic of EDM [30] . 
error analysis for diffusion-based generative models: 
 In this section, we will quantify both the training and sampling processes, and then integrate them into a more comprehensive generation error analysis.In this section, we consider a practical implementation of denoising score matching objective, represent the score by a deep ReLU network, and establish the exponential convergence of GD training dynamics.Training objective function. Consider a quadrature discretization of the time integral in (4) based on deterministic1 collocation pointsEQUATIONEQUATIONEQUATION)Define β j = w(t j )(t j -t j-1 ) 1 σt j to be the total weighting. Consider the empirical version of L (6) .Denote the initial data to be {x i } n i=1 with x i ∼ P 0 , and the noise to be {ξ ij } N j=1 with ξ ij ∼ N (0, I d ). Then the input data of the neural network is {t j , X ij } n,N i=1,j=1 = {t j , x i + σtj ξ ij } n,N i=1,j=1 and the output data is {ξ ij /σ tj } n,N i=1,j=1 if σtj ≠ 0. Consequently, L(θ) (6) can be approximated by the followingEQUATIONWe will use (7) as the training objective function in our analysis. For simplicity, we also denotef (θ; i, j) = β j ∥σ tj S(θ; t j , x i + σtj ξ ij ) + ξ ij ∥ 2 and then Lem (θ) = 1 2n ∑ n i=1 ∑ N j=1 f (θ; i, j). Note the time dependence can be absorbed into the X dependence. More precisely, because σt is a monotonically increasing function of t, we can replace t j in the inputs by σtj to indicate the time dependence. This is then equivalent to augmenting X ij to be d + 1 dimensional with (x i ) d+1 ∶= 0 and (ξ ij ) d+1 ∶= 1. For simplified presentation, we will slightly abuse notation and still use d as the input dimension rather than d + 1.Architecture. The analysis of diffusion model training is in general very challenging. One obvious factor is the complex score parameterizations used in practice such as U-Net [42] and transformers [40, 34] . In this paper, we simplify the architecture and consider deep feedforward networks. Although it is still far from practical usage, note this simple structure can already provide insights about the design space, as shown in later sections, and is more complicated than existing works [25, 44] related to the training of diffusion models (see Section 1.1). More precisely, we consider the standard deep ReLU network with bias absorbed:EQUATIONwhereEQUATIONL , W L+1 ). We consider the gradient descent (GD) algorithm as followsEQUATIONwhere h > 0 is the learning rate. We fix W 0 and W L+1 throughout the training process and only update W 1 , ⋯, W L , which is a commonly used setting in the convergence analysis of neural networks [1, 10, 25] . Initialization. We employ the same initialization as in Allen-Zhu et al. [1] , which is to set (W(0) ℓ ) ij ∼ N (0, 2 m ) for ℓ = 0, ⋯, L, i, j = 1, ⋯, m, and (W (0) L+1 ) ij ∼ N (0, 1 d ) for i = 1⋯, d, j = 1⋯, m.For this setup, the main challenge in our convergence analysis for denoising score matching lies in the nature of the data. 1) The output data that neural network tries to match is an unbounded Gaussian random vector, and cannot be rescaled as assumed in many theoretical works (for example, Allen-Zhu et al. [1] assumed the output data to be of order o(1)). 2) The input data X ij is the sum of two parts: x i which follows the initial distribution P 0 , and a Gaussian noise σtj ξ ij . Therefore, any assumption on the input data needs to agree with this noisy and unbounded nature, and commonly used assumptions like data separability [1, 35] can no longer be used.To deal with the above issues, we instead make the following assumptions. Assumption 1 (On network hyperparameters and initial data of the forward dynamics). We assume the following holds: 1. Data scaling:∥x i ∥ = Θ(d 1/2 ) for all i. 2. Input dimension: d = Ω(poly(log(nN ))).We remark that the first assumption focuses only on the initial data x i instead of the whole solution of the forward dynamics X ij which incorporates the Gaussian noise. Also, this assumption is indeed not far away from reality; for example, it holds with high (at least 1 -O(exp(-Ω(d))) probability for standard Gaussian random vectors. The requirement for input dimension d is to ensure that d is not too small, or equivalently the number of data points is not exponential in d.We also make the following assumptions on the hyperparameters of the denoising score matching. Assumption 2 (On the design of diffusion models). We assume the following holds:1. Weighting: ∑ N j=1 w(t j )(t j -t j-1 )σ tj = O(N ). 2. Variance: σt0 > 0 and σt N = Θ(1). The first assumption is to guarantee that the weighting function w(t) is properly scaled. This expression w(t j )(t j -t j-1 )σ tj is obtained from proving the upper and lower bounds of the gradient of (7) , and is different from the total weighting β i defined above. In the second assumption, σt0 > 0 ensures the output ξ ij /σ tj is well-defined. The σt N = Θ(1) guarantees that the scales of the noise σtj ξ ij and the initial data x i are of the same order at the end of the forward process, namely, the initial data x i is eventually push-forwarded to near Gaussian with the proper variance. Therefore, Assumption 2 aligns with what has been used in practice (see Section 4 and Karras et al. [30] , Song et al. [46] for examples).The following theorem summarizes our convergence result for the training of the score function. Theorem 1 (Convergence of GD). Define a set of indices to be G (s) = {(i, j)|f (θ (s) ; i, j) ≥ f (θ (s) ; i ′ , j ′ ) for all i ′ , j ′ }. Then given Assumption 1 and 2, for any ϵ train > 0, there exists some M (ϵ train ) = Ω (poly(n, N, d, L, T /t 0 , log(1 ϵ train ))), s.t., when m ≥ M (ϵ train ), h = Θ( nN m minj w(tj )(tj -tj-1)σt j ), and k = O(d 1-a 0 2 n 2 N log( d ϵ train )), with probability at least 1 -O(nN ) exp(-Ω(d 2a0-1 )), we haveLem (θ (k) ) ≤ k-1 ∏ s=0 ⎛ ⎝ 1 -C 5 h w(t j * (s) )(t j * (s) -t j * (s)-1 )σ t j * (s) ⎛ ⎝ md a 0 -1 2 n 3 N 2 ⎞ ⎠ ⎞ ⎠ Lem (θ (0) )where the universal constant C 5 > 0, a 0 ∈ ( 1 2 , 1), and (i * (s), j * (s)) = arg max (i,j)∈G (s) w(t j )(t jt j-1 )σ tj . Moreover, when K = Θ(d1-a 0 2 n 2 N log( d ϵ train )), Lem (θ (K)) ≤ ϵ train .The above theorem implies that for denoising score matching objective Lem (θ), GD has exponential convergence. For example, if we simply take j * = min j w(t j )(t j -t j-1 )σ tj , then Lem (θ (k+1) ) is further upper bounded by (1 -C 6 h w(t j * )(t j * -t j * -1 )σ t j * ( mda 0 -1 2 n 3 N 2 )) k+1 Lem (θ (0)). The rate of convergence can be interpreted in the following way: 1) at the kth iteration, we collect all the indices of the time points into G (k) where f (θ; i, j) has the maximum value; 2) we then choose the maximum of w(t j )(t j -t j-1 )σ tj among all such indices and denote the index to be j * (k), and obtain the decay ratio bound for the next iteration as1 -C 6 h w(t j * (k) )(t j * (k) -t j * (k)-1 )σ t j * (k) ( md a 0 -1 2 n 3 N 2 ). Remark 1 (Can ϵ train be arbitrarily small? Some ramifications of the denoising setting). Let us first see some facts about Lem and L. Under minimal assumption of the existence of score function and in the zero-time-discretization-error limit, the score matching objective can be made zero and therefore the denoising score matching objective is bounded below by -C, which is nonnegative and zero only when the data distribution is extremely special (we thus write -C > 0 from hereon unless confusion arises). That is, min θ L(θ) ≥ min any function S L = -C > 0 according to (4). Since Lem → L as the sample size of the training data set n → ∞, we have Lem ≥ -Cc n > 0 for some constant c n > 0 and c n → 0 as n → ∞.However, Theorem 1 seems to imply Lem (θ (k) ) → 0 as k → ∞ since Lem (θ (K) ) ≤ ϵ train and ϵ train is arbitrary, and it seems to contradict the -C > 0 lower bound. However, there is no contradiction due to the combination of two facts. First, the theorem states that for arbitrary ϵ train > 0, there exists a critical size, such that for overparameterized network beyond this size, GD can render the loss Lem (θ) eventually no greater than ϵ train . If we fix the network size, i.e., with m, L, d given, then K is given, and Theorem 1 says nothing about GD's behavior after K iterations. That is, we do not know whether lim sup k→∞ Lem (θ (k) ) = 0. Second, our optimization setting requires the sample size n to be smaller than the network width m (Assumption 1). Thus, when m is fixed, the sample size n is upper bounded.The above discussion implies, within the validity of our theory, for any fixed network width m, if ϵ train is small, the sample size n cannot be too large, meaning Lem (θ) -L(θ) may not be small. Therefore, we can simultaneously have Lem (θ) close to 0 and L(θ) close to -C > 0.Main technical steps for proving Theorem 1. The proof of Theorem 1 is in Appendix D, where the analysis framework is adapted from Allen-Zhu et al. [1] . Roughly speaking, the key proof in this framework is to establish the lower bound of the gradient. Then by integrating it into the semi-smoothness property of the neural network, we can obtain the exponential rate of convergence of gradient descent. For the lower bound of gradient, we develop a new method to deal with the difficulties in the denoising score matching setting (see the discussions earlier in this section).Our new proof technique adopts a different decoupling of the gradient and leverages a high probability bound based on a high-dimensional geometric idea. See Appendix D.1 for a proof sketch and more details.In this section, we establish a nonasymptotic error bound of the backward process in the variance exploding setting, which is an extension to Benton et al. [8] . For simplified notations, denote the backward time schedule as {t← j } 0≤j≤N such that 0 = t ← 0 < t ← 1 < ⋯ < t ← N = T -δ. Generation algorithm.We consider the exponential integrator scheme for simulating the backward SDE (2) with f t ≡ 02 . The generation algorithm can be piecewisely expressed as a continuous-time SDE: for anyEQUATIONInitialization. Denote q t ∶= Law( Ȳt ) for all t ∈ [0, Tδ]. We choose the Gaussian initialization, q 0 = N (0, σ2 T ). Our convergence result relies on the following assumption. Assumption 3. The distribution P 0 has a finite second moment:E x∼P0 [∥x∥ 2 ] = m 2 2 < ∞.Next we state the main convergence result, whose proof is provided in Appendix E.Theorem 2. Under Assumption 3, for any δ ∈ (0, 1) and T > 1, we haveEQUATIONEQUATION) where γ j ∶= t ← j+1 -t ← j for all j = 0, 1, ⋯, N -1 is the stepsize of the generation algorithm in (10) . Theorem 2 is a VESDE-based diffusion model's analogy of what's proved in Benton et al. [8] for VPSDE-based diffusion model, only requiring the data distributions to have finite second moments, and it achieves the sharp almost linear data dimension dependence under the exponential time schedule. The major differences from [8] are (1) the initialization error in the VESDE case is handled differently (see Lemma 10); (2) Theorem 2 applies to varies choices of time schedules, which enables to investigate the design space of the diffusion model, as we will discuss in Section 4. Worth mentioning is, Yang et al. [54] also obtained polynomial complexity results for VESDE-based diffusion models with uniform stepsize, but under stronger data assumption (assuming compact support). Compared to their result, complexity implied by Theorem 2 has better accuracy and data dimension dependencies. A detailed discussion on complexities is given in Appendix I. 1 . (11) represent the three types of errors: initialization error, discretization error, and score estimation error, respectively. Term E I quantifies the error between the initial density of the sampling algorithm q 0 and the ideal initialization p T , which is the density when the forward process stops at time T . Term E D is the error stemming from the discretization of the backward dynamics. Term E S characterizes the error of the estimated score function and the true score, and is related to the optimization error of Lem . Important to note is, in Theorem 2, population loss is needed instead of the empirical version Lem (7) . Besides this, the weighting γ j σ 2 T -t ← j is not necessarily the same as the total weighting in Lem (7) β j , depending on choices of w(t j ) and time and variance schedules (see more discussion in Section 4). We will later on integrate the optimization error (Theorem 1) into this score error E S to obtain a full error analysis in Section 3.3. Remark 2 (sharpness of dependence in d and m 2 2 ). In one of the simplest cases, when the data distribution is Gaussian, the score function is explicitly known. Hence KL(p δ |q T -δ ) can be explicitly computed as well, which verifies that the dependence of parameters d and m 2 2 is sharp in E I and E D .Terms E I , E D , E S inIn this section, we combine the analyses from the previous two sections to obtain an end-to-end generation error bound.Before providing the main result of this section, let us first clarify some terminologies.Time schedule, variance schedule, and total weighting. The terms time schedule and variance schedule respectively refer to the choice of t ← j and σtj in sampling. Meanwhile, note both the training and sampling processes require the proper choices of time and variance, and these choices are not necessarily the same for both processes. For training, the effect of these two is integrated into the total weighting β j , which is also influenced by an additional weighting parameter w(t j ). In this theoretical paper, when studying the generation error, we aim to apply the optimization result to better understand the effect of optimization on sampling. Therefore, to simplify the analysis and discussions in Section 4, we choose the same time and variance schedules for both training and sampling. The main result is stated in the following. Theorem 3. Under the same conditions as Theorem 1,2, and that K is such that GD reaches ϵ train in at most Kth iterations, we haveKL(p δ |q T -δ ) ≲ E I + E D + max 1≤j≤N σ 2 t N -j w(t N -j ) (ϵ train + ϵ n + ϵ est + ϵ approx )where E I , E D are defined in Theorem 2, ϵ train is defined in Theorem 1,ϵ n = | L(θ (K) )-Lem (θ (K) )+ Lem (θ * ) -L(θ * )|, ϵ est = | L(θ * ) -L(θ F )|, ϵ approx = | L(θ F ) + C|.In these terms, C is defined in (5), θ * = arg min θ,s.t., Lem(θ)=0 L(θ) and θ F = arg inf {θ∶S(θ)∈F } | L(θ) + C| with F = {ReLU network function defined in (8) , with d = Ω(poly(log(nN ))), m = Ω (poly(n, N, d, L, T /t 0 )) }.In this theorem, the discretization error E D and initialization error E I are the same as Theorem 2.For the score error E S , our optimization result is valid for general time schedules and therefore can directly fit into the sampling error analysis, which is in contrast to existing works [25, 44] (see more discussions in Section 1.1). The coefficient max j σ 2 t N -j /w(t N -j ) results from different weightings in E S and Lem , i.e., γ j σ 2 T -t ← j and β j . We will discuss the effect of max j σ 2 t N -j /w(t N -j ) under different time and variance schedules in Section 4.E Y t ← j ∼p T -t ← j [∥S(θ; T -t ← j , Y t ← j ) -∇ log p T -t ← j (Y t ← j )∥ 2 ] in E S (see Theorem 2) is to decompose it into the optimization error ϵ train , statistical error ϵ n , estimation error ϵ est , and approximation error ϵ approx . This gives clear intuition to results, but we also note it may not give a tight bound. In fact, we haveϵ n + ϵ train = | L(θ (K) ) -Lem (θ (K) ) + Lem (θ * ) -L(θ * )| + | Lem (θ (K) ) -Lem (θ * )| ≥ L(θ (K) ) + L(θ * ) ≥ 2 min θ L(θ) ≥ -2 C > 0.ϵ n can still be small if we take n → ∞, but that means ϵ train has to be large, and our generation error bound cannot be made 0. It is unclear yet whether this is due to limitation of our analysis or intrinsic, and will be left for future investigation.Another related note is, in this paper, we focus on ϵ train and the effect of optimization, but the analyses of ϵ n , ϵ est , and ϵ approx are also important and possible [13, 39, 25, 50] . On the other hand, again, whether it is optimal to decompose the full error into these four is unclear.To better see the parameter dependence of the error bound in Theorem 3, the following is an example with simplified results, where we employ the schedules in EDM [30] .Corollary 1 (Full error analysis under EDM [30] designs). Under the same conditions as Theorem 3, we haveKL(p δ |q T -δ ) ≲ m 2 2 T 2 + da 2 T 1 a δ 1 a N + (m 2 2 + d) ( a 2 T 1 a δ 1 a N + a 3 T 2 a δ 2 a N 2 ) + 1 N ⎛ ⎜ ⎝ C 9 + ⎛ ⎝ 1 -C 8 h ⎛ ⎝ md a 0 -1 2 n 3 N 2 ⎞ ⎠ ⎞ ⎠ K ⎞ ⎟ ⎠, where C 8 , C 9 > 0 and a = 7 in [30] .4 Theory-based understanding of the design space and its relation to existing empirical counterpartsThis section theoretically explores preferable choices of parameters in both training and sampling, and shows that they agree with the ones used in EDM [30] and Song et al. [46] in different circumstances. 
choice of total weighting for training: 
 This section develops the optimal total weighting β j for training objective (7) . We qualitatively show in two steps that "bell-shaped" weighting, which is the one used in EDM [30] , will lead to the optimal rate of convergence: Step 1) ∥σ tj S(θ; t j , X ij ) + ξ ij ∥ as a function of j is inversely "bell-shaped";Step 2) f (θ; i, j) = β j ∥σ tj S(θ; t j , X ij ) + ξ ij ∥ should be close to each other for any i, j.4.1.1 Inversely "bell-shaped" loss ∥σ tj S(θ; t j , X ij ) + ξ ij ∥ as a function of time index j Proposition 1. Under the same assumptions as Theorem 1, for any θ and i = 1, ⋯, n, we have1. ∀ϵ 1 > 0, ∃ δ 1 > 0, s.t., when 0 ≤ σtj < δ 1 , ∥σ tj S(θ; t j , x i + σtj ξ ij ) + ξ ij ∥ > ∥ξ ij ∥ -ϵ 1 . 2. ∀ϵ 2 > 0, ∃ M > 0, s.t., when σtj > M , ∥σ tj S(θ; t j , x i + σtj ξ ij ) + ξ ij ∥ ≥ M 2 (∥S(θ; t j , ξ ij )∥ -ϵ 2 ).The above proposition can be interpreted in the following way. Given any network S, when σtj is very small, 1 implies that ∥σ tj S(θ; t j , x i + σtj ξ) + ξ ij ∥ is away from 0 by approximately ∥ξ ij ∥ which is of order √ d with high probability, i.e., it cannot be small. When σtj is large, 2 shows that as it becomes larger and larger, i.e., as M increases, ∥σ tj S(θ; t j , x i + σtj ξ) + ξ ij ∥ will also increase. Therefore, the function ∥σ tj S(θ; t j , X ij ) + ξ ij ∥ has most likely an inversely "bell-shaped" curve in terms of j dependence.Ensuring comparable values of f (θ; i, j) for optimal rate of convergence Corollary 2. Under the same conditions as Theorem 1, for some large K ′ > 0, if |f (θ (k+K ′ ) ; i, j)f (θ (k+K ′ ) ; l, s)| ≤ ϵ holds for all k > 0 and all (i, j), (l, s), with some small universal constant ϵ > 0, then we have, for some constantC 7 > 0, Lem (θ (k+K ′ ) ) ≤ ⎛ ⎝ 1 -C 7 h max j=1,⋯,N w(t j )(t j -t j-1 )σ tj ⎛ ⎝ md a 0 -1 2 n 3 N 2 ⎞ ⎠ ⎞ ⎠ k Lem (θ (K ′ )). The above corollary shows that if f (θ (k) ; i, j)'s are almost the same for any i, j, then the decay ratio of the next iteration is minimized. More precisely, the index set G (k) defined in Theorem 1 is roughly the whole set {1, ⋯, N }, and therefore w(t j * )(t j * -t j * -1 )σ t j * can be taken as the maximum value over all j, which consequently leads to the optimal rate.Combining the above two aspects, the optimal rate of convergence leads to the choice of total weighting β j such that f (θ; i, j) = β j ∥σ tj S(θ; t j , X tj ) + ξ ij ∥ is close to each other; as a result, the total weighting should be chosen as a "bell-shaped" curve as a function of j according to the shape of the curve for ∥σ tj S(θ; t j , X tj ) + ξ ij ∥.Before comparing the preferable weighting predicted by our theory and the intuition-and-empiricsbased one in EDM [30] , let us first recall that the EDM training objective3 can be written asEQUATIONwhere Z 1 is a normalization constant, and we denoteβ EDM (σ) = e - (log σ-Pmean ) 2 2P 2 std σ2 +σ 2 data σσ 2 datato be the total weighting of EDM. Note the dependence on σ and time j can be freely switched due to their 1-to-1 correspondence.Figure 2 plots the total weighting of EDM β EDM as a function of σ. As is shown in the picture, this is a "bell-shaped" curve 4 , which coincides with our choice of total weighting in the above theory. When σ is very small or very large, according to Proposition 1, the lower bound of ∥σ tj S(θ; t j , X tj ) + ξ ij ∥ cannot vanish and therefore needs the smallest weighting over all σ. When σ takes the middle value, the scale of the output data ξ ij /σ j is roughly the same as the input data X ij and therefore makes it easier for the neural network to fit the data, which admits larger weighting.This section will discuss the choice of time and variance schedules based on the three errors E S , E D , E I in the error analysis of Section 3.3. Two situations will be considered based on how well the score function is approximated in training: when the network is less trained, E S dominates and polynomial schedule [30] is preferable; when the score function is well approximated, E D + E I dominates and exponential schedule [46] is better.As is shown in Theorem 3, the main impact of different time and variance schedules on score error E S appears in the term max j σ 2 t N -j /w(t N -j ), when the score function is approximated to a certain accuracy. It remains to compute w(t) under various choices of schedules.To ensure fair comparisons between different time and variance schedules, we maintain a fixed total weighting in the training objective. Additionally, to facilitate comparisons with practical usage, we adopt the total weighting in EDM, i.e., β j = C 3 β EDM (σ tj ), for some universal constant C 3 > 0. The reason for using the EDM total weighting is that according to Section 4.1, our total weighting β j should be "bell-shaped" as a function of j, which agrees qualitatively with the one used in EDM.Polynomial schedule [30] vs exponential schedule [46] . We fix ϵ n , ϵ train and apply the two schedules (Table 1 ) separately to the above total weighting β (hence w). Then, compute max j σ 2 t N -j /w(t N -j ) which is a factor in score error E S (Thm.3) in Table 2 . The Exp.'s result1 2 (σ max -σmax ( σ2 min σ2 max ) 1/N ) is larger 5 than the Poly.'s result (σ max -(σ 1/ρ max - σ1/ρ max -σ 1/ρ min N ) ρ) for large N , meaning the poly. time schedule in EDM is better than the exp. schedule in [46] . Note these two terms are both of order 1/N as N → ∞ and therefore the difference lies in their prefactors. Poly. [30] t (σ1/ρ max -(σ 1/ρ max - σ1/ρ min ) N -j N ) ρ Exp. [46] √ t σ2 max ( σ2 min σ2 max ) N -j N Table 2: Comparisons between different schedules.E S (score error) dominatesE D + E I (sampling error) dominates max j σ 2 tj /w(t j ) Choice N Choice Poly. [30] C 4 (σ max -(σ 1/ρ max - σ1/ρ max -σ 1/ρ min N ) ρ ) Ω ( m 2 2 ∨d d ρ 2 ( σmax σmin ) 1/ρ σ2 max ) Exp. [46] C 4 ⋅ 1 2 (σ max -σmax ( σ2 min σ2 max ) 1/N ) Ω ( m 2 2 ∨d d ln( σmax σmin ) 2 σ2 max )In this section, we compare the two different schedules in Table 1 by studying the iteration complexity of the sampling algorithm, i.e., number of time points N , when E D + E I dominates.General rules of comparison. We consider the case when the discretization and initialization errors are bounded by the same quantity ϵ, i.e., E I + E D ≲ ε. Then according to Theorem 2 and Theorem 3, we compute the iteration complexity for achieving this error using the two schedules in Table 1 .To make the comparison more straightforward, we adopt T = t N = Θ(poly(ε -1 )) and therefore σmax = Θ(ε -1/2 ). More details are provided in Appendix I.1.Polynomial schedule [30] vs exponential schedule [46] . As is shown in the last column of Table 2 , the iteration complexity under exponential schedule [46] has the poly-logarithmic dependence on the ratio between maximal and minimal variance (σ max /σ min ) 6 , which is better than the complexity under polynomial schedule [30] , which is polynomially dependent on σmax /σ min . Both complexities are derived from Theorem 2 by choosing different parameters. Remark 3 (The existence of optimal ρ in the polynomial schedule [30] ). For fixed σmax and σmin , the optimal ρ that minimizes the iteration complexity is ρ = 1 2 ln ( σmax σmin ). In [30] , it was empirically observed that with fixed iteration complexity, there is an optimal value of ρ that minimizes the FID. Our result indicates that, for fixed σmax and σmin , hence the desired accuracy in KL divergence being fixed, there is an optimal value of ρ that minimizes the iteration complexity to reach the fixed accuracy. Even though we consider a different metric/divergence instead of FID, our result still provides a quantitative support to the existence of optimal ρ observed in [30] .Conclusions. In this paper, we provide a first full error analysis incorporating both optimization and sampling processes. For the training process, we provide a first result under a deep neural network and prove the exponential convergence into a neighborhood of minima. At the same time, we extend the current analysis to the variance exploding case for sampling. Moreover, based on the full error analysis, we establish a quantitative understanding of the error bound under the two schedules. Consequently, we conclude with a qualitative illustration of the "bell-shaped" weighting and the choices of schedules under well-trained and less-trained cases.Limitations. The network architecture we used in the model is a deep ReLU network. Although being so far the most complicated architecture for theoretical results, it is still far from what is used in practice like U-Nets and transformers. Moreover, regarding the full error analysis, we only focus on the optimization and sampling error and do not dissect the generalization error. When bridging the theoretical results with practical designs of diffusion models, our results are mostly qualitative and we only compare two existing schedules under two extreme cases, when the network is well-trained and less-trained. Thus, theoretical implications on practical designs remain to be explored. We will leave these perspectives for future exploration. The weight in the ith layer of the network S (8) θ (k) The kth iteration of the weights θ through GD ( 9)W (k) ℓThe kth iteration of the weights W ℓ through GD (9) m Width of the network (8)Depth of the network (8) (i * (s), j * (s)) Index of the largest loss and w(t j )(t j -t j-1 )σ tj at the sth iteration (Theorem In this section, we will derive the denoising score matching objective, i.e. show the equivalence of ( 3) and ( 4). For simplicity, we denote S θ to be the neural network we use S(θ; t, X t ).ConsiderEQUATION)where p t is the density of X t .Since p t (x) = ∫ p 0 (y)q t (x|y)dy, where q t (⋅) is the density of X t |X 0 , then we haveE Xt ⟨S θ , ∇ log p t ⟩ = ∫ S ⊺ θ ∇ log p t ⋅ p t dx t = ∫ S ⊺ θ ∇p t dx t = ∬ S ⊺ θ ∇q t (x|y)p 0 (y) dxdy = ∬ S ⊺ θ ∇ log q t (x|y)p 0 (y)q t (x|y) dxdy = E X0∼P0 E Xt|X0∼Qt ⟨S θ , ∇ log q t (x t |x 0 )⟩ . Then (13) = E X0∼P0 E Xt|X0∼Qt [∥S θ ∥ 2 -2 ⟨S θ , ∇ log q t (x t |x 0 )⟩] + E Xt ∥∇ log p t ∥ 2 = E X0 E Xt|X0 ∥S θ -∇ log q t ∥ 2 (∆)+C,whereC = E Xt ∥∇ log p t ∥ 2 -E X0 E Xt|X0 ∥∇ log q t (x t |x 0 )∥ 2 .Moreover, X t |X 0 ∼ N (e -µt X 0 , σ2 t I), and its density function isq t (x|y) = (2πσ 2 t ) -d/2 exp (- ∥x -e -µt y∥ 2 2σ 2 t) .Then(∆) = E X0 E Xt|X0 ∥S θ -∇ log q t ∥ 2 = E X0 E Xt|X0 ∥S θ -∇ x (- ∥X t -e -µt X 0 ∥ 2 2σ 2 t )∥ 2 = E X0 E Xt|X0 ∥S θ + X t -e -µt X 0 σ2 t ∥ 2 = E X0 E ϵt ∥S θ + ϵ t σ2 t ∥ 2 . Let ξ = ϵt σt ∼ N (0, I). Then (∆) = E X0 E ξ σt ⋅ 1 σ2 t ∥σ t S θ + ξ∥ 2 = 1 σt E X0 E ξ ∥σ t S θ + ξ∥ 2In this section, we will prove Theorem 1.Before introducing the concrete proof, we first redefine the deep fully connected feedforward networkr ij,0 = W 0 X ij , q ij,0 = σ(r ij,0 ), r ij,ℓ = W ℓ q ij,ℓ-1 , q ij,ℓ = σ(r ij,ℓ ), for ℓ = 1, ⋯, L S(θ; t j , X ij ) = W L+1 q ij,Lwhere W 0 ∈ R m×d , W L+1 ∈ R d×m and W ℓ ∈ R m×m ; σ is the ReLU activation. We also denote q ij,-1 to be X ij .We also follow the notation in Allen-Zhu et al. [1] and denote D i,ℓ ∈ R m×m to be a diagonal matrix and(D i,ℓ ) kk = 1 (W ℓ q ij,ℓ-1 ) k >0 for k = 1, ⋯, m. Then q ij,ℓ = D ij,ℓ W ℓ q ij,ℓ-1For the objective (7) , the gradient w.r.t. to the kth row of W ℓ for ℓ = 1, ⋯, L is the following∇ (W ℓ ) k Lem (θ) = 1 n n ∑ i=1 N ∑ j=1 w(t j )(t j -t j-1 ) [(W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 R ij,ℓ+1 ) ⊺ (σ tj W L+1 q ij,L + ξ ij )] k q ij,ℓ-1 1 (W ℓ q ij,ℓ-1 ) k >0Throughout the proof, we use both Lem (θ) and Lem (W ) to represent the same value of the loss function, where W = (W 1 , ⋯, W L ), and we let a = b = 1 2 . Next, we will prove Theorem 1.Proof of Theorem 1. First by Lemma 4,Lem (W (0) ) = O(d 2a ∑ j w(t j )(t j -t j-1 )/σ tj ) Also, ∥∇ Lem (θ)∥ ≤ √ L max ℓ ∥∇ W ℓ Lem (θ)∥. Then we have ∥W (k) -W (0) ∥ ≤ k-1 ∑ i=0 h∥∇ Lem (W (i) )∥ ≤ O( √ md 2a-1 N L max j w(t j )(t j -t j-1 )σ tj )hk max i √ Lem (W (i) ) ≤ O( √ md 2a-1 N L max j w(t j )(t j -t j-1 )σ tj d a )hk √ ∑ j w(t j )(t j -t j-1 )/σ tj ∶= ω Let h = Θ( nN m minj w(tj )(tj -tj-1)σt j ) and k = O(d 1-a 0 2 n 2 N log( d ϵ train )), where ϵ train > 0 is some small constant. Then ω = O(log( d ϵ train ) d 1-a 0 with probability at least 1 -O(nN ) exp(-Ω(d 2a0-1 )), Lem (W (k+1) ) ≤ Lem (W (k) ) -h∥∇ Lem (W (k) )∥ 2 + h √ Lem √ ∑ j w(t j )(t j -t j-1 )σ tj O(ω 1/3 L 2 √ m log md a/2 )∥∇ Lem (W (k) )∥ + h 2 √ Lem √ ∑ j w(t j )(t j -t j-1 )σ tj O(L 2 √ md a )∥∇ Lem (W (k) )∥ 2 ≤ ⎛ ⎝ 1 -hw(t j * )(t j * -t j * -1 )σ t j * ⋅ Ω ⎛ ⎝ md a 0 -1 2 n 3 N 2 ⎞ ⎠ ⎞ ⎠ Lem (W (k) ) + hC m 5/6 d 7/12-a0/6 N 1/6 n 2/3 (log m) 1/6 √ L √ ∑ j w(t j )(t j -t j-1 )σ tj min j w(t j )(t j -t j-1 )σ tj max j w(t j )(t j -t j-1 )σ tj ∑ k w(t j )(t j -t j-1 )/σ tj Lem (W (k) ) ≤ ⎛ ⎝ 1 -hw(t j * )(t j * -t j * -1 )σ t j * ⋅ Ω ⎛ ⎝ md a 0 -1 2 n 3 N 2 ⎞ ⎠ ⎞ ⎠ Lem (W (k) )where C > 0 is some constant, a 0 ∈ (1/2, 1); the second inequality follows from Lemma 7 with∥∇ W L Lem (θ (k) )∥ 2 = Ω ⎛ ⎝ md a 0 -1 2 n 3 N 2 w(t j * )(t j * -t j * -1 )σ t j * ⎞ ⎠ Lem (θ (k) ),which is obtained inductively; the last inequality follows from m =Ω ⎛ ⎝ d 13/2-2a0/3 n 14/3 N 11 L 3 (log m) ( maxj w(tj )(tj -tj-1)σt j ∑ k w(tj )(tj -tj-1)/σt j minj w(tj )(tj -tj-1)σt j √ ∑ j w(tj )(tj -tj-1)σt j ) 6 ⎞ ⎠In this section, we will show the main part of the convergence analysis, which is the following lower bound of the gradient. Lemma 1 (Lower bound). With probability 1 -O(nN ) exp(-Ω(d 2a0-1 )), we have∥∇ Lem (θ (0) )∥ 2 ≥ C 6 ⎛ ⎝ md a 0 -1 2 n 3 N 2 w(t j * )(t j * -t j * -1 )σ t j * ⎞ ⎠ Lem (θ (0) )where(i * , j * ) = arg max ∥ √ w(tj )(tj -tj-1) σt j (σ tj W L+1 q ij,L + ξ ij )∥, 1 2 < a 0 < 1, and C 6 > 0 is some universal constant.Below is the proof sketch of Lemma 1.Proof sketch. We first decompose the gradient of the kth row ofW L ∇ (W L ) k Lem (θ) = 1 n w(t j * )(t j * -t j * -1 )(W L+1 ) k ⊺ (σ t j * W L+1 q i * j * ,L + ξ i * j * )q i * j * ,L-1 1 (W L q i * j * ,L-1 ) k >0 ∇1 + 1 n ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )(W L+1 ) k ⊺ (σ tj W L+1 q ij,L + ξ ij ) q ij,L-1 1 (W L q ij,L-1 ) k >0 ∇2where (i * , j * ) indicates the sample index with the largest loss value.Then we first fix (q ij,L-1 ) s = 1, and prove that the index set of both (q i * j * ,L ) s > 0 and ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )σ tj 1 (W L q ij,L-1 ) k >0 (q ij,L ) s > 0 is order m with high probability.Next, we conditioned on the index set we've found, then we can decouple each element of ∇ (W L ) k Lem with high probability. We prove that with high probability∠(W L+1 σt j * q i * j * ,L + ξ i * j * , W L+1 ∑ (i,j)≠(i * ,j * ) α ij q ij,L + ∑ (i,j)≠(i * ,j * ) ᾱij ξ ij ) ≤ π -cd a 0 -1 2 ,for some constant c > 0 and 1 2 < a 0 < 1. Based on this, we show that with probability at least 1 -O(nN ) exp(-Ω(d)),P ((∇ 1 ) s > 0, (∇ 2 ) s > 0) ≥ cd a 0 -1 2 ,for some c > 0. Then we prove that with probability at least 1 -exp(-Ω(mda 0 -1 2 )) |{k ∶ (W k L+1 ) ⊺ v ≥ 0, (W k L+1 ) ⊺ (u + ξ) ≥ 0}| = Θ(md a 0 -1 2 ).with high probability, the event (∇ 1 ) s > 0 and (∇ 2 ) s > 0 has probability at least of order d (a0-1)/2 where a 0 ∈ (1/2, 1).Now, we deal with (q ij,L-1 ) s and prove that if the above results hold for (q ij,L-1 ) s = 1, then there exists an index set with cardinality of order m/(nN ) such that (∇ 1 ) s > 0 and (∇ 2 ) s > 0 also hold in this index set.In the end, combining all the steps above yields the lower bound.Here is the complete proof.Proof. The main idea of the proof of lower bound is to decouple the elements in the gradient and incorporate geometric view. We focus on ∇ W L Lem (θ).Step 1: Rewrite ∇ (W L ) k Lem (θ) to be the (i * , j * )th term g 1 plus the rest nN -1 terms g 2 .Let (i * , j * ) = arg max ∥ √ w(tj )(tj -tj-1) σt j (σ tj W L+1 q ij,L + ξ ij )∥. Let g ij,L = w(t j )(t j -t j-1 )(W L+1 ) k ⊺ (σ tj W L+1 q ij,L + ξ ij ) q ij,L-1 . Then ∇ (W L ) k Lem (θ) = 1 n w(t j * )(t j * -t j * -1 )(W L+1 ) k ⊺ (σ t j * W L+1 q i * j * ,L + ξ i * j * ) q i * j * ,L-1 1 (W L q i * j * ,L-1 ) k >0 ∇1 + 1 n ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )(W L+1 ) k ⊺ (σ tj W L+1 q ij,L + ξ ij ) q ij,L-1 1 (W L q ij,L-1 ) k >0 ∇2 Also define ∇ 1,s = 1 n w(t j * )(t j * -t j * -1 )σ t j * (W L+1 ) k ⊺ W L+1 q i * j * ,L (q i * j * ,L-1 ) s 1 (W L q i * j * ,L-1 ) k >0 ∇11,s + 1 n w(t j * )(t j * -t j * -1 ) (W L+1 ) k ⊺ ξ i * j * (q i * j * ,L-1 ) s 1 (W L q i * j * ,L-1 ) k >0 ∇12,s ∇ 2,s = 1 n ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )σ tj (W L+1 ) k ⊺ W L+1 q ij,L (q ij,L-1 ) s 1 (W L q ij,L-1 ) k >0 ∇21,s + 1 n ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )(W L+1 ) k ⊺ ξ ij q ij,L-1 (q ij,L-1 ) s 1 (W L q ij,L-1 ) k >0 ∇22,sOur goal is to show that with high probability, there are at least O( mdEQUATIONnN ) number of rows k such that ∇ 11,s ≥ 0, ∇ 12,s ≥ 0, ∇ 21,s ≥ 0, ∇ 22,s ≥ 0. Then we can lower bound ∥∇ (W L ) k Lem (θ)∥ 2 by ∥∇ 1 ∥ 2 , which can be eventually lower bounded by Lem (θ).Step 2: Consider [∇ (W L ) k Lem (θ)] s . For (g 2 ) s , first take (q ij,L-1 ) s = 1 for all (i, j) ≠ (i * , j * ). Then we only need to consider∇ ′ 2,s = 1 n ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )(W L+1 ) k ⊺ (σ tj W L+1 q ij,L + ξ ij ) 1 (W L q ij,L-1 ) k >0which is independent of s. For ∇ 1 , since q i * ,j * ,L-1 ≥ 0 which does not affect the sign of this term, we can also first take (q i * ,j * ,L-1 ) s = 1 for all s.Step 3: We focus on ∇ 11 and ∇ 21 and we would like to pick the non-zero elements in this two terms. More precisely, letN 1 = {s | (q i * j * ,L ) s > 0, s = 1, ⋯, m} , N 2 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ s | ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )σ tj 1 (W L q ij,L-1 ) k >0 (q ij,L ) s > 0, s = 1, ⋯, m ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ Let α ij = w(t j )(t j -t j-1 )σ tj 1 (W L q ij,L-1 ) k >0 ≥ 0. Then ∑ (i,j)≠(i * ,j * ) w(t j )(t j -t j-1 )σ tj 1 (W L q ij,L-1 ) k >0 (q ij,L ) s = ∑ (i,j)≠(i * ,j * ) α ij (q ij,L ) s = ∑ (i,j)≠(i * ,j * ) α ij σ(W L q ij,L-1 ) s . If ∑ (i,j)≠(i * ,j * ) α ij (W L q ij,L-1 ) s = (W L ) s ∑ (i,j)≠(i * ,j * ) α ij q ij,L-1 > 0,then there must be at least one pair of (i, j) s.t. α ij (W L q ij,L-1 ) s = α ij σ(W L q ij,L-1 ) s > 0, which implies ∑ (i,j)≠(i * ,j * ) α ij (q ij,L ) s > 0. Therefore, it suffices to considerN 1 = {s | (q i * j * ,L ) s = (W L ) s q i * j * ,L-1 > 0, s = 1, ⋯, m} , N ′ 2 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ s | (W L ) s ∑ (i,j)≠(i * ,j * ) α ij q ij,L-1 > 0 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ .Since (q ij,L-1 ) s ≥ 0, we haveEQUATIONBy Lemma 2 and Proposition 2, we haveP ⎛ ⎝ (W L ) s q i * j * ,L-1 > 0, (W L ) s ∑ (i,j)≠(i * ,j * ) α ij q ij,L-1 > 0 ⎞ ⎠ = P ⎛ ⎝ (W L ) s ∥(W L ) s ∥ q i * j * ,L-1 > 0, (W L ) s ∥(W L ) s ∥ ∑ (i,j)≠(i * ,j * ) α ij q ij,L-1 > 0 ⎞ ⎠ ≥ 1 4 .Also (W L ) s 's are i.i.d. multivariate Gaussian. By Chernoff bound,P (|N 1 ∩ N ′ 2 | ∈ (δ 1 m 4 , δ 2 m 4 )) ≤ 1 -2e -Ω(m)for some small δ 1 ≤ 1 4 and δ 2 ≤ 4, i.e., |N 1 ∩ N ′ 2 | = Θ(m) with probability at least 1 -2e -Ω(m) .Step 4: Next we condition onN 1 ∩ N ′ 2 and consider (W L+1 ) k ⊺ W L+1 σt j * q i * j * ,L + (W L+1 ) k ⊺ ξ i * j * and (W L+1 ) k ⊺ W L+1 ∑ (i,j)≠(i * ,j * ) α ij q ij,L + (W L+1 ) k ⊺ ∑ (i,j)≠(i * ,j * ) ᾱij ξ ij , where ᾱij = α ij /σ ij .We would like to prove that with high probability∠(W L+1 σt j * q i * j * ,L + ξ i * j * , W L+1 ∑ (i,j)≠(i * ,j * ) α ij q ij,L + ∑ (i,j)≠(i * ,j * ) ᾱij ξ ij ) ≤ π -cd a 0 -1 2 ,for some constant c > 0 and 1 2 < a 0 < 1. First, since ξ ij ∼ N (0, I d ), by Bernstein's inequality, with probability at least 1 -exp(-Ω(d)), wehave ∥ξ ij ∥ 2 = Θ(d). Similarly, since (W L+1 q ij,L ) s ∼ N (0, 2∥q ij,L ∥ 2 m), by Berstein's inequality and Lemma 4, with probability at least 1 -exp(-Ω(d)), we have ∥W L+1 q ij,L ∥ 2 = Θ(d). By union bounds, the above holds for all i, j with probability at least 1 -2nN exp(-Ω(d)).Let v = W L+1 ∑ (i,j)≠(i * ,j * ) αij q ij,L +∑ (i,j)≠(i * ,j * ) ᾱij ξij ∥W L+1 ∑ (i,j)≠(i * ,j * ) αij q ij,L +∑ (i,j)≠(i * ,j * ) ᾱij ξij ∥ and u = W L+1 σt j * q i * j * ,L .For notational simplicity, we use ξ to denote ξ i * j * . Fix v, u and consider the probability of eventA A = {v ⊺ (u + ξ) ≤ - √ 1 -c 0 d a0-1 ∥u + ξ∥}for some c 0 > 0 and 1 2 < a 0 < 1. Then consider the following event that has larger probability than AEQUATIONEQUATIONSince (v ⊺ u) 2 ≤ ∥u∥ 2 where the equality holds when v = u ∥u∥ , we haveLHS ≤ c 0 d a0-1 ∥u∥ 2 + 2(v ⊺ u v -(1 -c 0 d a0-1 )u) ⊺ ξ + (v ⊺ ξ) 2Also, since ∥u∥ 2 = O(d) with probability at least 1 -2nN exp(-Ω(d)), we haveP (|2(v ⊺ u v -(1 -c 0 d a0-1 )u) ⊺ ξ| ≥ d a0 ) ≤ 2 exp (-c d 2a0 ∥u∥ 2 ) = 2 exp(-Ω(d 2a0-1))for some constant c > 0.Therefore, with probability at least 1 -O(nN ) exp(-Ω(d 2a0-1 ))LHS of (15) ≤ cd a0 + (v ⊺ ξ) 2 for some constant c > 0.ThenP(v ⊺ (u + ξ) ≤ - √ 1 -c 0 d a0-1 ∥u + ξ∥) ≤ P((v ⊺ (u + ξ)) 2 ≥ (1 -c 0 d a0-1 )∥u + ξ∥ 2 ) ≤ P((v ⊺ ξ) 2 ≥ (1 -c ′ d a0-1 )∥ξ∥ 2 ) = P (v ⊺ ξ ∥ξ∥ ≥ √ (1 -c ′ d a0-1 )) + P (-v ⊺ ξ ∥ξ∥ ≥ √ (1 -c ′ d a0-1 )) = P (∠(v, ξ ∥ξ∥ ) ≥ arccos(- √ (1 -c ′ d a0-1 ))) + P (∠(-v, ξ ∥ξ∥ ) ≥ arccos(- √ (1 -c ′ d a0-1 ))) = 2P (∠(v, ξ ∥ξ∥ ) ≥ π -c ′′ d a 0 -1 2 ) = C (d 1-a 0 2 ) d-1 √ dwhere the second equality follows from Lemma 2; the third equality follows from series expansion; the forth equality follows from ( 16); c ′ , c ′′ , C > 0 are some constants.Thus with probability at least 1 -C (d 1-a 0 2 ) d-1 √ d , v ⊺ (u + ξ) ≥ - √ 1 -c 0 d a0-1 ∥u + ξ∥ i.e.∠(v, u + ξ) ≤ π -cd a 0 -1 2for some c > 0.Then by Lemma 2, with probability at least 1 -O(nN ) exp(-Ω(d)),P ((W k L+1 ) ⊺ v ≥ 0, (W k L+1 ) ⊺ (u + ξ) ≥ 0) ≥ cd a 0 -1 2 ,for some c > 0.Since (W k L+1 ) are iid Guassian vectors for k = 1, ⋯, m, by Chernoff bound on Bernoulli variable1 {(W k L+1 ) ⊺ v≥0,(W k L+1 ) ⊺ (u+ξ)≥0}, we have, with probability at least 1 -exp(-Ω(mda 0 -1 2 )) |{k ∶ (W k L+1 ) ⊺ v ≥ 0, (W k L+1 ) ⊺ (u + ξ) ≥ 0}| = Θ(md a 0 -1 2 ).Step 5: Combining the above 4 steps, we would like to obtain the lower bound of the gradient.For each k, consider (q ij,L-1 ) s for (i, j) ≠ (i * , j * ) and denote q s = {(q ij,L-1 ) s } (i,j)≠(i * ,j * ) = {σ((W L-1 ) s q ij,L-2 )} (i,j)≠(i * ,j * ) . Let qs = {(W L-1 ) s q ij,L-2 } (i,j)≠(i * ,j * ) and qs ∼ N (0, QQ ⊺ ), where each row of Q is q ⊺ ij,L-2 for (i, j) ≠ (i * , j * ). Thus, q s is qs projected to the nonnegative orthant. Let q K ij,ℓ = (q ij,ℓ ) k∈K . Then by Bernstein's inequality, we can also obtain that ∥q K ij,ℓ ∥ 2 = Θ(d) with probability at least 1 -nN exp(-Ω(d)).Let 1 = (1, 1, ⋯, 1) ∈ R nN -1 . Therefore, if ⟨β k , 1⟩ ≥ 0 for some β k ∈ R nN -1 , then at least half of the nonnegative orthant is contained in {v ∈ R nN -1 ∶ ⟨β k , v⟩ ≥ 0}, i.Combine all of the above and apply the Claim 9.5 in Allen-Zhu et al. [1] , we obtain, with probability at least 1 -O(nN ) exp(-Ω(d 2a0-1 )),∥∇ W L Lem (θ (0) )∥ 2 F ≥ 1 n 2 C 6 w(t j * ) 2 (t j * -t j * -1 ) 2 1 d ∥σ t j * W L+1 q i * j * ,L + ξ i * j * ∥ 2 ∥q K i * j * ,L-1 ∥ 2 1 nN md a 0 -1 2 ≥ C 6 md a 0 -1 2 w(t j * )(t j * -t j * -1 )σ t j * 1 n 3 N 2 Lem (θ (0)),where C 6 > 0 is some universal constant, 1 2 < a 0 < 1, and the second inequality follows from the definition of i * , j * .Proposition 2. Consider w ∼ N (0, I), where w ∈ R n . Then ∥w∥ and w ∥w∥ are independent random variables and w ∥w∥ ∼ Unif(S n-1 ). Lemma 2. Let w ∼ Unif (S n-1 ), whereS n-1 = {x ∈ R n |∥x∥ = 1}. Then for two vectors v 1 , v 2 ∈ R n , P(w ⊺ v 1 ≥ 0, w ⊺ v 2 ≥ 0) = π -∠(v 1 , v 2 ) 2π .Proof. Since w ∼ U nif (S n-1 ), we only need to consider the area of the event. It is obvious that the set {w ∈ S n-1 |w ⊺ v i } is a semi-hypersphere. Therefore, we only need to consider the intersection of two semi-hypersphere, i.e.,P(w ⊺ v 1 ≥ 0, w ⊺ v 2 ≥ 0) = area of {w ∈ S n-1 |w ⊺ v 1 ≥ 0} ∩ area of {w ∈ S n-1 |w ⊺ v 2 ≥ 0} area of the hypersphere = π -∠(v 1 , v 2 ) 2π .Next we follow the notations and definitions in Lee and Kim [33] . Consider the unit hypersphere inR d , S d-1 = {x ∈ R d | ∥x∥ = 1}. The area of S d-1 is A d (1) = 2π d/2 Γ(d/2) . Lemma 3. Fix ξ 1 ∈ S d-1 and let ξ 2 ∼ Unif (S d-1 ), where S d-1 = {x ∈ R d |∥x∥ = 1}.Then with probability at least 1 -exp(-Ω(d)), we have ∠(ξ 1 , ξ 2 ) ≤ 3π 4 .Proof. For any fixed ξ 1 , all the ξ 2 's that satisfy ∠(ξ 1 , ξ 2 ) ≥ πθ are on a hyperspherical cap. By Lee and Kim [33] , the area of the hypersperical cap isEQUATIONThenEQUATION)Let θ = π 4 < 1. Then with probability at least 1 -exp(-Ω(d)), we have ∠(ξ 1 , ξ 2 ) ≤ 3π 4 .Consider W i = W (0) i in this section.Lemma 4. If ϵ ∈ (0, 1), with probability at least 1 -O(nN )e -Ω(min(ϵ2 d 4b-1 ,ϵd 2b )) , ∥X ij ∥ 2 ∈ [∥e -µt j x i ∥ 2 + σ2 tj d -ϵσ 2 tj d 2b , ∥e -µt j x i ∥ 2 + σ2 tj d + ϵσ 2 tj d 2b] for all i = 1, ⋯, n and j = 0, ⋯, N -1. Moreover, with probability at least 1 -O(L)e -Ω(mϵ 2 /L) over the randomness of W s for s = 0, ⋯, L, we have ∥q ij,ℓ ∥ ∈ [∥X ij ∥(1 -ϵ), ∥X ij ∥(1 + ϵ)] for fixed i, j. Therefore, with probability at least 1 -O(nN L)e -Ω(min(mϵ 2 /L,ϵ 2 d 4b-1 ,ϵd 2b )) , we haveΩ(d b ) = ∥q ij,ℓ ∥ = O(d a ). Proof. Consider 1 σt j X ij = e -µ t j σt j x i +ξ ij . Since ξ ij ∼ N (0, I), ∥ 1 σt j X ij ∥ 2 follows from the noncentral χ 2 distribution and E∥ 1 σt j X ij ∥ 2 = d+∥ e -µ t j σt jx i ∥ 2 (this includes the time variable at the dth dimension).1 σtj X ij ∥ 2 -E∥ 1 σtj X ij ∥ 2 | ≥ t) ≤ 2 exp ( -c min ( t 2 d , t))i.e., P(|∥e-µt j x i + σtj ξ ij ∥ 2 -(σ 2 tj d + ∥e -µt j x i ∥ 2 )| ≥ σ2 tj t) ≤ 2 exp ( -c min ( t 2 d , t))Therefore, with probability at least 1 -O(nN )e -Ω(min(ϵ2 d 4b-1 ,ϵd 2b )) , ∥X ij ∥ 2 ∈ [∥e -µt j x i ∥ 2 + σ2 tj d - ϵσ 2 tj d 2b , ∥e -µt j x i ∥ 2 + σ2 tj d + ϵσ 2 tj d 2b] for all i = 1, ⋯, n and j = 0, ⋯, N -1, where ϵ ∈ (0, 1). The second part of the Lemma follows the similar proof in Lemma 7.1 of Allen-Zhu et al. [1] . The last part follows from union bound and Assumption 1.Lemma 5 (Upper bound). Under the random initialization of W i for i = 0, ⋯, L, with probability at least 1 -O(nN L)e -Ω(min(mϵ 2 /L,ϵ 2 d 4b-1 ,ϵd 2b )) , we have∥∇ W ℓ Lem (θ (0) )∥ 2 = O (md 2a-1 N max j w(t j )(t j -t j-1 )σ tj ) Lem (θ (0)).Proof. For any ℓ = 1, ⋅, L, we have∥∇ W ℓ Lem (θ)∥ 2 F = m ∑ k=1 ∥∇ (W ℓ ) k Lem (θ)∥ 2 = m ∑ k=1 ∥ 1 n n ∑ i=1 N ∑ j=1 w(t j )(t j -t j-1 ) × [(W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 ) ⊺ (σ tj W L+1 q ij,L + ξ ij )] k q ij,ℓ-1 1 (W ℓ q ij,ℓ-1 ) k >0 ∥ 2 ≤ N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 m ∑ k=1 ∥(W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 ) ⊺ k (σ tj W L+1 q ij,L + ξ ij )∥ 2 ⋅ ∥q ij,ℓ-1 ∥ 2 ≤ C 7 d 2a m d N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 ⋅ ∥σ tj W L+1 q ij,L + ξ ij ∥ 2 ≤ C 7 d 2a mN d max j w(t j )(t j -t j-1 )σ tj Lem (θ)where the first inequality follows from Young's inequality; the second inequality follows from Lemma 4 and Lemma 7.4 in Allen-Zhu et al. [1] ; C 7 > 0Consider W per i = W (0) i + W ′ i for i = 1, ⋯, L in this section. We follow the same idea in Allen-Zhu et al. [1] to consider the network value of perturbed weights at each layer. We use the superscript "per" to denotes the perturbed version, i.e., r per ij,0 = W 0 X ij , q per ij,0 = σ(r per ij,0 ), r per ij,ℓ = W per ℓ q per ij,ℓ-1 , q per ij,ℓ = σ(r per ij,ℓ ), for ℓ = 1, ⋯, L S(θ per ; t j , X ij ) = W L+1 q per ij,LWe also similarly define the diagonal matrix D per ij,ℓ for the above network. The following Lemma measures the perturbation of each layer. The lemma differs from Lemma 8.2 in Allen-Zhu et al. [1] by a scale of d a . For sake of completeness, we state it in the following and the proof can be similarly obtained. Lemma 6. Let ω ≤ 1 C7L 9/2 (log m) 3 d a for some large C > 1. With probability at least 1exp(-Ω(d a mω 2/3 L)), for any ∆W s.t. ∥∆W ∥ ≤ ω, we have1. r per ij,ℓ -r ij,ℓ can be decomposed to two part r per ij,ℓ -r ij,ℓ = r ′ ij,ℓ,1 + r ′ ij,ℓ,2 , where ∥r ′ ij,ℓ,1 ∥ = O(ωL 3/2 d a ) and ∥r ′ ij,ℓ,2 ∥ ∞ = O(ωL 5/2 √ log md a m -1/2 ). 2. ∥D per ij,ℓ -D ij,ℓ ∥ 0 = O(mω 2/3 L) and ∥(D per ij,ℓ -D ij,ℓ )r per ij,ℓ ∥ = O(ωL 3/2 d a ).3. ∥r per ij,ℓ -r ij,ℓ ∥ and ∥q per ij,ℓ -q ij,ℓ ∥ are O(ωL 5/2 √ log md a ).Lemma 7 (Upper and lower bounds of gradient after perturbation). Letω = O ⎛ ⎝ L * em L 9 (log m) 2 n 3 N 3 d 1-a 0 2⋅ min j w(t j )(t j -t j-1 )σ tj max{max j w(t j )(t j -t j-1 )σ tj , ∑ j (w(t j )(t j -t j-1 )σ tj ) 2 } ⎞ ⎠ . Consider θ per s.t. ∥θ per -θ∥ ≤ ω, where θ follows from the Gaussian initialization. Then with probability at least 1 -O(nN )e -Ω(d 2a 0 -1 ) ,∥∇ W ℓ Lem (θ per )∥ 2 = O (md 2a-1 N max j w(t j )(t j -t j-1 )σ tj ) Lem (θ per ), ∥∇ W L Lem (θ per )∥ 2 = Ω ⎛ ⎝ md a 0 -1 2 n 3 N 2 w(t j * )(t j * -t j * -1 )σ t j * ⎞ ⎠ min{ Lem (θ), Lem (θ per )}, for ℓ = 1, ⋯, L.Proof. Consider the following termsEQUATIONEQUATIONEQUATIONThen∥∇ W ℓ Lem (θ) -∇ per W ℓ Lem (θ)∥ 2 F = m ∑ k=1 ∥∇ (W ℓ ) k Lem (θ) -∇ per (W ℓ ) k Lem (θ)∥ 2 ≤ N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 m ∑ k=1 ∥(W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 ) ⊺ k (σ tj W L+1 (q ij,L -q per ij,L ))∥ 2 ⋅ ∥q ij,ℓ-1 ∥ 2 ≤ C 8 d 2a m d N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 ⋅ ∥σ tj W L+1 (q ij,L -q per ij,L )∥ 2 ≤ C ′ 8 d 2a m d N N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 (ωL 5/2 √ log md a ) 2 ≤ C8 ⎛ ⎝ md a 0 -1 2 n 3 N 2 w(t j * )(t j * -t j * -1 )σ t j * ⎞ ⎠ L * emwhere the first two inequalities follow the same as the proof of Lemma 5; the third inequality follows from Lemma 6; the last inequality follows from the definition of ω.Also, we have∥∇ per (W ℓ ) k Lem (θ) -∇ (W ℓ ) k Lem (θ per )∥ ≤ ∥ 1 n n ∑ i=1 N ∑ j=1 w(t j )(t j -t j-1 ) × [(W L+1 D per ij,L W per L ⋯D per ij,ℓ W per ℓ+1 -W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 ) ⊺ (σ tj W L+1 q per ij,L + ξ ij )] k × q per ij,ℓ-1 1 (W per ℓ q per ij,ℓ-1 ) k >0 ∥ + ∥ 1 n n ∑ i=1 N ∑ j=1 w(t j )(t j -t j-1 )[(W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 ) ⊺ (σ tj W L+1 q per ij,L + ξ ij )] k × (q ij,ℓ-1 1 (W ℓ q ij,ℓ-1 ) k >0 -q per ij,ℓ-1 1 (W per ℓ q per ij,ℓ-1 ) k >0 ) ∥ Then ∥∇ per W ℓ Lem (θ) -∇ W ℓ Lem (θ per )∥ 2 F ≤ 2 N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 ⋅ ∥q per ij,ℓ-1 ∥ 2 × m ∑ k=1 ∥(W L+1 D per ij,L W per L ⋯D per ij,ℓ W per ℓ+1 -W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 ) ⊺ k (σ tj W L+1 q per ij,L + ξ ij )∥ 2 + 2 N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 ⋅ ∥q ij,ℓ-1 1 (W ℓ q ij,ℓ-1 ) k >0 -q per ij,ℓ-1 1 (W per ℓ q per ij,ℓ-1 ) k >0 ∥ 2 × m ∑ k=1 ∥(W L+1 D ij,L W L ⋯D ij,ℓ W ℓ+1 ) ⊺ k (σ tj W L+1 q per ij,L + ξ ij )∥ 2 ≤ C 9 (ω 2 L 5 log md 2a )(L 3/2 log m m 1/2 ) m d N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 ⋅ ∥σ tj W L+1 q per ij,L + ξ ij ∥ 2 + C ′ 9 (ω 2 L 5 log md 2a ) m d N n n ∑ i=1 N ∑ j=1 w(t j ) 2 (t j -t j-1 ) 2 ⋅ ∥σ tj W L+1 q per ij,L + ξ ij ∥ 2 ≤ C9 ⎛ ⎝ md a 0 -1 2 n 3 N 2 w(t j * )(t j * -t j * -1 )σ t j * ⎞ ⎠ Lem (θ per )where the first inequality follows from Young's inequality and the above decomposition; the second inequality follows from Lemma 7.4, 8.7 in Allen-Zhu et al. [1] (with s = O(d a mω 2/3 L)) and Lemma 6; the last inequality follows from the definition of ω.For upper bound, we only need to consider ∇ per (W ℓ ) k Lem (θ) and ∇ (W ℓ ) k Lem (θ per ). By similar argument as Lemma 5, with probability at least 1 -O(nN L)e -Ω(min(mϵ 2 /L,ϵ 2 d 4b-1 ,ϵd 2b )) , we have∥∇ per (W ℓ ) k Lem (θ)∥ 2 = O (md 2a-1 N max j w(t j )(t j -t j-1 )σ tj ) Lem (θ per ).Then∥∇ W ℓ Lem (θ per )∥ 2 ≤ 2∥∇ per W ℓ Lem (θ)∥ 2 F + 2∥∇ per W ℓ Lem (θ) -∇ W ℓ Lem (θ per )∥ 2 F = O (md 2a-1 N max j w(t j )(t j -t j-1 )σ tj ) Lem (θ per ).Also,∥∇ W Lem (θ per )∥ 2 ≥ ∥∇ W L Lem (θ per )∥ 2 ≥ 1 3 ∥∇ W L Lem (θ)∥ 2 -∥∇ W ℓ Lem (θ) -∇ per W ℓ Lem (θ)∥ 2 F -∥∇ per W ℓ Lem (θ) -∇ W ℓ Lem (θ per )∥ 2 F = Ω ⎛ ⎝ md a 0 -1 2 n 3 N 2 w(t j * )(t j * -t j * -1 )σ t j * ⎞ ⎠ min{ Lem (θ), Lem (θ per )}.Note when interpolation is not achievable, this lower bound is always away from 0, which means the current technique can only evaluate the lower bound outside a neighbourhood of the minimizer. More advanced method is needed and we leave it for future investigation.Lemma 8 (semi-smoothness). Let ω = Ω. With probability at least 1-e -Ω(log m) over the randomness of θ (0) , we have for all θ s.t. ∥θ -θ (0) ∥ ≤ ω, and all θ per s.t. ∥θ per -θ∥ ≤ ω,Lem (θ per ) ≤ Lem (θ) + ⟨∇ Lem (θ), θ per -θ⟩ + √ Lem (θ) √ ∑ j w(t j )(t j -t j-1 )σ tj O(ω 1/3 L 2 √ m log md a/2 )∥θ per -θ∥ + √ Lem (θ) √ ∑ j w(t j )(t j -t j-1 )σ tj O(L 2 √ md a )∥θ per -θ∥ 2Proof. By definition,Lem (θ per ) -Lem (θ) -⟨∇ Lem (θ), θ per -θ⟩ = 1 n n ∑ i=1 N ∑ j=1 w(t j )(t j -t j-1 )(σ tj W L+1 q ij,L + ξ ij ) ⊺ W L+1 × (q per ij,L -q ij,L - L ∑ ℓ=1 D ij,L W ij,L ⋯W ij,ℓ+1 D ij,ℓ (W per ij,ℓ -W ij,ℓ )q ij,ℓ ) + 1 2σ tj w(t j )(t j -t j-1 )∥σ tj W L+1 (q per ij,L -q ij,L )∥ 2 .Similar to the proof of Theorem 4 in Allen-Zhu et al. [1] , we obtain the desired bound by using Cauchy-Schwartz inequality. Note, in our case, due to the order of input data, we choose s = O(d a mω 2/3 L) in Allen-Zhu et al. [1] and therefore the bound is slightly different from theirs.In this section, we prove Theorem 2. The proof includes two main steps: 1. decomposing KL(p δ |q T -δ ) into the initialization error, the score estimation errors and the discretization errors; 2. estimating the initialization error and the discretization error based on our assumptions. In the following context, we introduce the proof of these two steps separately.Proof of Theorem 2.Step 1: The error decomposition follows from the ideas in [12] of studying VPSDE-based diffusion models. According to the chain rule of KL divergence, we haveKL(p δ |q T -δ ) ≤ KL(p T |q 0 ) + E y∼p T [KL(p δ|T (⋅|y)|q T -δ|0 (⋅|y))],Apply the chain rule again for at across the time schedule (T -t ← j ) 0≤j≤N -1 , the second term can be written asE y∼p T [KL(p δ|T (⋅|y)|q T -δ|0 (⋅|y))] ≤ N -1 ∑ j=0 E yj ∼p T -t ← j [KL(p T -t ← j+1 |T -t ← j (⋅|y j )|q t ← j+1 |t ← j (⋅|y j ))] ≤ 1 2 N -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t E[∥s(θ; T -t ← j , Y t ← j ) -∇ log p T -t (Y t )∥ 2 ]dt ≤ N -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t E[∥s(θ; T -t ← j , Y t ← j ) -∇ log p T -t ← j (Y t ← j )∥ 2 ]dt + N -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t E[∥∇ log p T -t ← j (Y t ← j ) -∇ log p T -t (Y t )∥ 2 ]dt,where the second inequality follows from Lemma 9. Therefore, the error decomposition writes asEQUATIONwhere the three terms in (20) quantify the initialization error, the score estimation error and the discretization error, respectively.Step 2: In this step, we estimate the three error terms in Step 1. First, recall that p T = p * N (0, σ2 T I d ) and q 0 = N (0, σ2 T I d ), hence the initialization error KL(p T |q 0 ) can be estimated as follows,EQUATIONwhere the inequality follows from Lemma 10. Hence we recover the term E I in (11) .Next, since σ t is non-decreasing in t, the score estimation error can be estimated asEQUATIONHence, we recover the term E S in (11) .Last, we estimated the discretization error term. Our approach is motivated by analyses of VPSDEs in [8, 26] . We defines a process L t ∶= ∇ log p T -t (Y t ). Then we can relate discretization error to quantities depending on L t , and therefore bound the discretization error via properties of {L t } 0≤t≤T . According to Lemma 12, we haveN -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t E[∥∇ log p T -t ← j (Y t ← j ) -∇ log p T -t (Y t )∥ 2 ]dt ≤ 2d N -1 ∑ j=0 ∫ t ← j+1 t ← j ∫ t t ← j σ 2 T -t σ 2 T -u σ-4 T -u dudt N1 + N -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t dtσ -4 T -t ← j E[tr(Σ T -t ← j (X T -t ← j ))] N2 - N -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t σ-4 T -t E[tr(Σ T -t (X T -t ))]dt N3 . Since t ↦ σ t is non-decreasing and t ↦ E[tr(Σ T -t (X T -t ))] is non-increasing, we haveN 1 = 2d N -1 ∑ j=0 ∫ T -t ← j T -t ← j+1 ∫ T -u T -t ← j+1 σ 2 t dtσ 2 u σ-4 u du ≤ 2d N -1 ∑ j=0 γ j ∫ T -t ← j T -t ← j+1 σ 4 t σ-4 t dt, N 2 = N -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t dtσ -4 T -t ← j E[tr(Σ T -t ← j (X T -t ← j ))], N 3 ≤ - N -1 ∑ j=0 ∫ t ← j+1 t ← j σ 2 T -t σ-4 T -t dtE[tr(Σ T -t ← j+1 (X T -t ← j+1 ))].Therefore, we obtainEQUATIONThe above bound depends on E[tr(Σ t (X t ))], hence we estimate E[tr(Σ t (X t ))] for different values of t.First, we haveE[tr(Σ t (X t ))] = E[E[∥X 0 ∥ 2 |X t ] -∥E[X 0 |X t ]∥ 2 ] ≤ E[∥X 0 ∥ 2 ] = m 2 2 . Meanwhile, E[tr(Σ t (X t ))] = E[tr(Cov(X 0 -X t |X t ))] = E[E[∥X 0 -X t ∥ 2 |X t ]] -∥E[X 0 -X t |X t ]∥ 2 ≤ E[∥X 0 -X t ∥ 2 ] = σ2 t d Therefore, E[tr(Σ t (X t ))] ≤ min(m 2 2 , σ2 t d) ≲ (1 -e -σ 2 t )(m 2 2 + d).Plug this estimation into (23) and we getEQUATIONT -t ← j )( ∫ t ← j+1 t ← j σ 2 T -t σ-4 T -t ← j dt -∫ t ← j t ← j-1 σ 2 T -t σ-4 T -t dt) ≲ d N -1 ∑ j=0 γ j ∫ T -t ← j T -t ← j+1 σ 4 t σ-4 t dt + m 2 2 ∫ t ← 1 0 σ 2 T -t dt σ4 T + +(m 2 2 + d) N -1 ∑ k=1 (1 -e -σ 2 T -t ← j ) σ4 T -t ← j -σ2 T -t ← j+1 σ2 T -t ← j-1 σ2 T -t ← j-1 σ4 T -t ← j ,where the last inequality follows from the definition of σt and integration by parts. The proof of Theorem 2 is completed.Lemma 9. Let {Y t } 0≤t≤T be the solution to (2) with f t ≡ 0 and p ← t+s|s (⋅|y) be the conditional distribution of Y s+t given {Y s = y}. Let { Ȳt } 0≤t≤T be the solution to (11) q t+s|s (⋅|y) be the conditional distribution of Ȳs+t given { Ȳs = y}. Then for any fixed t ∈ (0, γ j ], we haveE y∼p ← t ← j KL(p ← t ← j +t|t ← j (⋅|y)|q t ← j +t|t ← j (⋅|y)) ≤ 1 2 σ 2 T -t E[∥s(θ; T -t ← j , Y t ← j ) -∇ log p T -t (Y t )∥ 2 ]Proof of Lemma 9. According to [12, Lemma 6] , we haveKL(p ← t ← j +t|t ← j (⋅|y)|q t ← j +t|t ← j (⋅|y)) ≤ -2σ 2 T -t ∫ p ← t ← j +t|t ← j (x|y)∥∇ log p ← t ← j +t|t ← j (x|y) q t ← j +t|t ← j (x|y) ∥ 2 dx + 2σ 2 T -t E p ← t ← j +t|t ← j (x|y) [⟨∇ log p T -t (x) -s(θ; T -t ← j , Y t ← j ), ∇ log p ← t ← j +t|t ← j (x|y) q t ← j +t|t ← j (x|y) ⟩] ≤ 1 2 σ 2 T -t E p ← t ← j +t|t ← j (x|y) [∥∇ log p T -t (x) -s(θ; T -t ← j , Y t ← j )∥ 2 ],where the last inequality follows from Young's inequality. Therefore, Lemma 9 is proved after taking another expectation.Lemma 10. For any probability distribution p satisfying Assumption 3 and q being a centered multivariate normal distribution with covariance matrix σ 2 I d , we haveKL(p * q|q) ≤ m 2 2 2σ 2 .Proof of Lemma 10.KL(p * q|q) ≤ ∫ KL(q(⋅ -y)|q(⋅))p(dy) = ∫ KL(N (y, σ 2 I d )|N (0, σ 2 I d ))p(dy) = 1 2 ∫ ln(1) -d + tr(I d ) + ∥y∥ 2 σ -2 p(dy) = m 2 2 2σ 2 ,where the inequality follows from convexity of KL(⋅|q) and the second identity follows from KLdivergence between multivariate normal distributions.Lemma 11. Let {X t } 0≤t≤T be the solution to (1) with f t ≡ 0 and p 0|t (⋅|x) be the conditional distribution of X 0 given {X t = x}. DefineEQUATIONLet {Y t } 0≤t≤T be the solution to (2) with f t ≡ 0 and q 0|t (⋅|x) be the conditional distribution ofEQUATIONThen we have for all t ∈ (0, T ),EQUATIONd dt E[Σ t (X t )] = 2σ 2 t σ-4 t E[Σ t (X t ) 2 ].Proof of Lemma 11. We first represent ∇logp t (X t ) and ∇ 2 log p t (X t ) via m t (X t ) and Σ t (X t ). Since {X t } 0≤t≤T solves (1), X t = X 0 + σt ξ with (X 0 , ξ) ∼ p ⊗ N (0, I d ). Therefore, according to Bayes rule, we haveEQUATIONwhere the second identity follows from the fact that p t|0 (⋅|x) = N (x, σ2 t I d ). The last identity follows from the definition of m t (X t ) in Lemma 11. Similarly, according to Bayes rule, we can computeEQUATIONwhere the second identity follows from the fact that p t|0 (⋅|x) = N (x, σ2 t I d ) and the definition of Σ t (X t ) in Lemma 11.According to Bayes rule, we haveEQUATION∥X t -x∥ 2 σ2 t )p(dx)andEQUATIONEQUATIONwhere Z t = ∫ exp(h t (x))p(dx) is a (random) normalization constant. From the above computations, we can see that q 0|t (dx|Y t ) = p 0|T -t (dx|X T -t ) for all t ∈ [0, T ]. Therefore, we havemt (Y t ) = E X∼q 0|t (⋅|Yt) [X] = m T -t (X T -t ), Σt (Y t ) = Cov X∼q 0|t (⋅|Yt) (X) = Σ T -t (X T -t ),where the identities hold in distribution. Therefore, to prove the first statement, it suffices to compute d mt (Y t ). To do so, we first compute dh t (x), d[h(x), h(x)] t , dZ t and d log Z t .EQUATIONAccording to the definition of Y t and (26), we haveEQUATIONApply ( 29) and ( 30) and we getEQUATIONandEQUATIONIf we further define R t (Y t ) ∶=EQUATIONWith ( 29), ( 30), ( 31), ( 32) and ( 33), we haveEQUATIONwhere most terms cancel in the last identity. Therefore, the first statement is proved. Next, we prove the second statement. We haved dt E[Σ T -t (X T -t )] = d dt E[ Σt (Y t )] = d dt E[Σ T -t (X T -t )] = d dt E Yt∼p T -t [E q 0|t (⋅|Yt) [x ⊗ 2] -mt (Y t ) ⊗ 2] = d dt E q0 [x ⊗ 2] - d dt E[ mt (Y t ) ⊗ 2] = -E[-2 mt (Y t )d mt (Y t ) ⊺ + d[ mt (Y t ), mt (Y t ) ⊺ ]] = 2σ -3 T -t σT -t E[ Σt (Y t ) 2 ]dt = -2σ 2 T -t σ-4 T -t E[Σ t (X T -t ) 2 ],where the second last identity follows from (34) and the last identity follows from the definition of σt . Last, we reverse the time and getd dt E[Σ t (X t )] = 2σ 2 t σ-4 t E[Σ t (X t ) 2 ].The proof is completed.Lemma 12. Under the conditions in Lemma 11, let {Y t } 0≤t≤T be the solution to (2) with f t ≡ 0. Define L t ∶= ∇ log p T -t (Y t ), then for any t ∈ [t ← j , t ← j+1 ), we have E[∥L t -L t ← j ∥ 2 ] = 2d ∫ t t ← j σ 2 T -u σ-4 T -u du + σ-4 T -t ← j E[tr(Σ T -t ← j (X T -t ← j ))] -σ-4 T -t E[tr(Σ T -t (X T -t ))]Proof of Lemma 12. First, according to the definition of L t and Y t , it follows from Itô's lemma thatEQUATIONEQUATIONwhere the last step follows from applying the Fokker Planck equation of (1) with f t ≡ 0, i.e., ∂ t p t = σ 2 t ∆p t . Most of the terms are cancelled after applying the Fokker Planck equation. Now, for fixed s > 0 and t > s, define E s,t ∶= E[∥L t -L s ∥ 2 ]. Apply Itô's lemma and (35) , we haveEQUATIONwhere ∥A∥ F denotes the Frobenius norm of any matrix A. According to (27) , we havedE s,t dt = 2σ 2 T -t E[∥∇ 2 log p T -t (Y t )∥ 2 F ] = 2σ 2 T -t E[∥∇ 2 log p T -t (X T -t )∥ 2 F ] = 2σ 2 T -t E[∥-σ -2 T -t I d + σ-4 T -t Σ T -t (X T -t )∥ 2 F ] = 2dσ 2 T -t σ-4 T -t -4σ 2 T -t σ-6 T -t E[tr(Σ T -t (X T -t ))] + 2σ 2 T -t σ-8 T -t E[tr(Σ T -t (X T -t ) 2 )] = 2dσ 2 T -t σ-4 T -t -4σ 2 T -t σ-6 T -t E[tr(Σ T -t (X T -t ))] -σ-4 T -t d dt E[tr(Σ T -t (X T -t ))],where the last identity follows from the proof of Lemma 11. Therefore, for any t ∈ [t ← j , t ← j+1 ), we haveE t ← j ,t = 2d ∫ t t ← j σ 2 T -u σ-4 T -u du -4 ∫ t t ← j σ 2 T -u σ-6 T -u E[tr(Σ T -u (X T -u ))]du -∫ t t ← j σ-4 T -u d du E[tr(Σ T -u (X T -u ))] = 2d ∫ t t ← j σ 2 T -u σ-4 T -u du -4 ∫ t t ← j σ 2 T -u σ-6 T -u E[tr(Σ T -u (X T -u ))]du -σ-4 T -t E[tr(Σ T -t (X T -t ))] + σ-4 T -t ← j E[tr(Σ T -t ← j (X T -t ← j ))] + 4 ∫ t t ← j σ 2 T -u σ-6 T -u E[tr(Σ T -u (X T -u ))] = 2d ∫ t t ← j σ 2 T -u σ-4 T -u du -σ-4 T -t E[tr(Σ T -t (X T -t ))] + σ-4 T -t ← j E[tr(Σ T -t ← j (X T -t ← j ))]The proof is completed.In this section, we consider a special case when the data distribution is a mixture of Gaussians, i.e.,EQUATIONwhere N (x; m, σ 2 I d ) is the density of Gaussian random vector with mean m and covariance σ 2 I d . In this case, the score function ∇ log p t (x) can be explicitly calculated from any t > 0, see Lemma 14. Therefore, the sampling process (11) can be implemented with zero score estimation error via the following piecewise SDE: for any t ∈ [t ← j , t ← j+1 ),EQUATIONThe iterates, ( Ȳt ← j ), are all Gaussians with explicit means and covariance matrices, see Lemma 15. As a consequence, we can quantify the quantity, KL(p δ |q T -δ ) in Theorem 2 explicitly since both distributions are Gaussians.Lemma 13 (KL-divergence error for Gaussian data distribution). Assume the data distribution has density given by (39) . Let ( Ȳt ← j ) N j=0 be defined by (40) with initial condition Ȳ0 ∼ N (0, σ2EQUATIONwhere E σ is a positive constant depending on the variance schedule (σ T -t ← j ) N j=0 , given byE -1 σ = (σ 2 + σ2 δ )σ 2 T (σ 2 + σ2 T ) 2 + N -1 ∑ j=0 (σ 2 + σ2 δ )(σ 2 T -t ← j -σ2 T -t ← j+1 ) (σ 2 + σ2 T -t ← j+1 ) 2.Proof of Lemma 13. p δ = p * N (0, σ2 δ ) = N (m, (σ 2 + σ2 δ )I d ) and q T -δ = Law( Ȳt ← N ) = N (m N , Σ N ) with (m N , Σ N ) given in Lemma 15. Therefore, we haveEQUATION) -d 2 + d 2 ( (σ 2 + σ2 δ )σ 2 T (σ 2 + σ2 T ) 2 + N -1 ∑ j=0 (σ 2 + σ2 δ )(σ 2 T -t ← j -σ2 T -t ← j+1 ) (σ 2 + σ2 T -t ← j+1 ) 2 ) -1 + (σ 2 T + (σ 2 + σ2 T ) 2 N -1 ∑ j=0 (σ 2 T -t ← j -σ2 T -t ← j+1 ) (σ 2 + σ2 T -t ← j+1 ) 2 ) -1 ∥m∥ 2 .Lemma 14 (Explicit score function for mixture of Gaussian target). Assume the data distribution has density given by (39) , then the score function is given byEQUATIONProof. Since the forward process (1) with f t ≡ 0 is the just a process that keeps adding noise, the density p t along the process is a convolution between data density and a Gaussian density with mean zero and covariance σ2 t I d :p t (x) = p * N (⋅ ; 0, σ2 t I d )(x) = N (x; m, (σ 2 + σ2 t )I d ).Therefore, we have∇p t (x) = (2π(σ 2 + σ2 t )) -d/2 exp ( - ∥x -m∥ 2 2(σ 2 + σ2 t ))( -x -m σ 2 + σ2 t ) = - x -m σ 2 + σ2 t N (x; m, (σ 2 + σ2 t )I d ).(42) follows directly from the above computations.Lemma 15 (Gaussian iterates along the trajectory). Assume the data distribution has density given by (39) . Let ( Ȳt ← j ) be defined by (40) with initial condition Ȳ0 ∼ N (0, σ2 T I d ). Then for all 0 ≤ j ≤ N , Ȳt ← j ∼ N (m j , Σ j ) withEQUATIONEQUATIONProof of Lemma 15. According to Lemma 14, (40) can be written asd Ȳt = -2σ 2 T -t Ȳt ← j -m σ 2 + σ2 T -t ← j dt + √ 2σ 2 T -t d Wt ,which implies that for any j ∈ 0, 1, ⋯, N -1:EQUATIONwhere (U j ) N j=1 are i.i.d. standard Gaussian vectors in R d . Since Ȳt ← 0 is Gaussian, by induction, we prove that Ȳt ← j is Gaussian for all j = 1, ⋯, N . Denote Ȳt ← j ∼ N (m j , Σ j ). According to (45) and the independence between U j+1 and Ȳt ← j , we havem j+1 = (1 - σ2 T -t ← j -σ2 T -t ← j+1 σ 2 + σ2 T -t ← j )m j + σ2 T -t ← j -σ2 T -t ← j+1 σ 2 + σ2 t ← j m, ⇒ m j+1 -m = σ 2 + σ2 T -t ← j+1 σ 2 + σ2 T -t ← j (m j -m), ⇒ m j = σ 2 + σ2 T -t ← j σ 2 + σ2 T -t ← 0 (m 0 -m) + m = σ2 T -t ← 0 + σ2 T -t ← j σ 2 + σ2 T -t ← 0 m.Again, according to (45) and the independence between U j+1 and Ȳt ← j , we get a relation between consecutive covariance matrices:EQUATIONT -t ← j+1 (σ 2 + σ2 T -t ← j+1 ) 2 I d , ⇒ Σ j (σ 2 + σ2 T -t ← j ) 2 = Σ 0 (σ 2 + σ2 T -t ← 0 ) 2 + j-1 ∑ l=0 σ2 T -t ← l -σ2 T -t ← l+1 (σ 2 + σ2 T -t ← l+1 ) 2 I d , ⇒ Σ j = ( (σ 2 + σ2 T -t ← j ) 2 σ2 T (σ 2 + σ2 T -t ← 0 ) 2 + j-1 ∑ l=0 (σ 2 + σ2 T -t ← j ) 2 (σ 2 T -t ← l -σ2 T -t ← l+1 ) (σ 2 + σ2 T -t ← l+1 ) 2 )I d .Proof of Theorem 3. We only need to deal with E S . By applying the same schedules to training objective, we obtainE S = N -1 ∑ j=0 σ 2 t N -j w(t N -j ) ⋅ w(t N -j )(t N -j -t N -j-1 ) 1 σt N -j E X0 E ξ ∥σ t N -j s(θ; t N -j , X t N -j ) + ξ∥ 2 + N -1 ∑ j=0 σ 2 t N -j w(t N -j ) ⋅ w(t N -j )(t N -j -t N -j-1 ) ⋅ C ≤ max j σ 2 t N -j w(t N -j ) ⋅ ( L(W ) + C).Together withL(W (k) ) + C ≤ | L(W (K) ) -Lem (W (K) ) + Lem (θ * ) -L(θ * )| + | Lem (W (K) ) -Lem (θ * )| + | L(θ * ) -L(θ F )| + | L(θ F ) + C|, we have the result.H Proofs for Section 4.1 H.1 Proof of "bell-shaped" curve Proof of Proposition 1. Fix x i , ξ ij , σt N . By definition of the network S(θ; t j , X ij ), it is continuous with respect to X ij .For 1, X ij = x i +σ tj ξ ij , and thus S(θ; t j , X ij ) is also continuous w.r.t. σtj . Also, since σtj ∈ [0, σt N ], there exists M 0 > 0, s.t., S(θ; t j , x i + σtj ξ ij ) ∈ [-M 0 , M 0 ] d . Then for any ϵ 1 > 0, there exists δ = ϵ1√ dM0 > 0, s.t., ∀ 0 ≤ σtj < δ 1 , we have ∥σ tj S(θ; t j , x i + σtj ξ ij ) + ξ ij ∥ ≥ ∥ξ ij ∥ -∥σ tj S(θ; t j , x i + σtj ξ ij )∥ ≥ ∥ξ ij ∥ - √ dM 0 σtj ≥ ∥ξ ij ∥ -ϵ 1 .For 2, by the positive homogeniety of ReLU,S(θ; t j , x i + σtj ξ ij ) = σtj S (θ; t j , x i σtj + ξ ij ) .Consider σtj ≥ M , for some M > 0. ThenEQUATIONFor any y ∈ D(M ), whereD(M ) = {y ∈ R d ∶ y s ∈ [(ξ ij ) s -|(x i ) s |/M, (ξ ij ) s + |(x i ) s |/M ], ∀ s = 1, ⋯, d},∥S(θ; t j , y)∥ ≥ ∥S(θ; t j , ξ ij )∥ -∥S(θ; t j , y) -S(θ; t j , ξ ij )∥.(47) Since S is differentiable a.e., by the fundamental theorem of calculus, S(θ; t j , y) -S(θ; tEQUATIONwhere 46), (47) , and (48), we haveM 1 = max s (x i ) s ⋅ ess sup x∈D(M2) ∥S ′ x (θ; t j , x)∥ < +∞ for some fixed 0 < M 2 < M . Combining (∥σ tj S(θ; t j , x i + σtj ξ ij ) + ξ ij ∥ ≥ M 2 (∥S(θ; t j , ξ ij )∥ - M 1 M ) -∥ξ ij ∥ = M 2 (∥S(θ; t j , ξ ij )∥ - M 1 M - ∥ξ ij ∥ M 2 ) . Then ∀ ϵ 2 > 0, there exists M = max { M1+ √ M 2 1 +4ϵ2∥ξij ∥ 2ϵ2, M 2 } > 0, s.t., when σtj > M , we have∥σ tj S(θ; t j , x i + σtj ξ ij ) + ξ ij ∥ ≥ M 2 (∥S(θ; t j , ξ ij )∥ -ϵ 2 ).Proof of Corollary 2. If |f (θ (k) ; i, j) -f (θ (k) ; l, s)| ≤ ϵ for all i, j, l, s and k > K, then by Lemma 1 and 7, we choose the maximum f (θ (k) ; i, j) for the lower bound, which is of order O(ϵ) away from the other f (θ (k) ; i, j) ′ s. Therefore, we can take j * (k) = arg max j f (θ (k) ; i, j) and absorb the O(ϵ) error in constant factors. Then the result naturally follows.Recall that the training objective of EDM is defined in the followingE σ∼p train E y,n λ(σ)∥D θ (y + n; σ) -y∥ 2 = 1 Z 1 ∫ 1 σ e - (log σ-Pmean ) 2 2P 2 std ⋅ σ2 + σ 2 data σ2 σ 2 data ⋅ σ2 E X0,ξ ∥σs(θ; t, X t ) + ξ∥ 2 dσ. Let β j = C 1 β EDM , i.e., w(t j ) σtj (t j -t j-1 ) = C 1 ⋅ e - (log σt j -Pmean ) 2 2P 2 std ⋅ σ2 tj + σ 2 data σ2 tj σ 2 data ⋅ σtj w(t j ) = C 1 ⋅ σtj t j -t j-1 ⋅ e - (log σt j -Pmean ) 2 2P 2 std ⋅ σ2 tj + σ 2 data σ2 tj σ 2 data ⋅ σtj EDM. Consider σt = t and t j = (σ 1/ρ max -(σ 1/ρ max - σ1/ρ min ) N -j N ) ρ for j = 0, ⋯, N . Then w(t j ) = C 1 ⋅ t j t j -t j-1 ⋅ e - (log t k -Pmean ) 2 2P 2 std ⋅ t j 2 + σ 2 data t j σ 2 data = C 1 ⋅ 1 (σ 1/ρ max -(σ 1/ρ max - σ1/ρ min ) N -j N ) ρ -(σ 1/ρ max -(σ 1/ρ max - σ1/ρ min ) N -j+1 N ) ρ ⋅ e - (log (σ 1/ρ max -(σ 1/ρ max -σ 1/ρ min ) N -j N ) ρ -Pmean ) 2 2P 2 std ⋅ (σ 1/ρ max -(σ 1/ρ max - σ1/ρ min ) N -j N ) 2ρ + σ 2 data σ 2 dataThen the maximum ofσ 2 t j w(tj ) = tj w(tj ) appears at j = N max j σ 2 tj w(t j ) = max j t j w(t j ) = σmax σ 2 data e (Pmean -log σmax) 2 2P 2 std C 1 (σ 2 max + σ 2 data ) ⋅ ⎛ ⎝ σmax - ⎛ ⎝ σ1/ρ max - σ1/ρ max - σ1/ρ min N ⎞ ⎠ ρ ⎞ ⎠Song et al. [46] . Consider σt = √ t and t j = σ2 max ( Under the EDM choice of variance, σt = t for all t ∈ [0, T ], and study the optimal time schedule when E D + E I dominates. First, it follows from Theorem 2 thatσ2 min σ2 max ) N -j N for j = 0, ⋯, N . Then w(t j ) = C 1 ⋅ √ t j t j -t j-1 ⋅ e - (log √ t j -Pmean ) 2 2P 2 std ⋅ t j + σ 2 data √ t j σ 2 data = C 1 ⋅ 1 σ2 max ( σ2 min σ2 max ) N -j N -σ2 max ( σ2 min σ2 max ) N -j+1 N ⋅ e - ⎛ ⎜ ⎝ log σmax( σmin σmax ) N -j N -Pmean ⎞ ⎟ ⎠ 2 2P 2 std ⋅ σ2 max ( σ2 min σ2 max ) N -j N + σ 2 data σ 2 data Then max j σ 2 tj w(t j ) = max j 1 2w(t j ) = σmax σ 2 data e (Pmean -log σmax) 2 2P 2 std C 1 (σ 2 max + σ 2 data ) ⋅ 1 2 (σ max -σmax (E I + E D ≲ m 2and it is easy to see that our theoretical result supports what's empirically observed in EDM that there is an optimal value of ρ that minimizes the FID.Exponential time schedule. we consider γ j = κ(Tt ← j ) with κ = ln(T /δ) N , we haveE I + E D ≲ m 2 2 T 2 + d ln(T /δ) 2 N + (m 2 2 + d)( ln(T /δ) 2 N + ln(T /δ) 3 N 2 ) Therefore, to obtain E I + E D ≲ ε, it suffices to require T = Θ( m2 ε 1 2) and the iteration complexityN = Ω( m 2 2 + d ε ln( m 2 δε 1 2 ) 2 ) When m 2 ≤ O( √ d),the exponential time schedule is asymptotic optimal, hence it is better than the polynomial time schedule when the initilization error and discretization error dominate. Once we let δ = σmin , T = σmax = Θ( m2 ε 1/2 ), the iteration complexity isN = Ω( m 2 2 ∨ d d ln ( σmax σmin ) 2 σ2 max ).Now we adopt the variance schedule in [46] , σt = √ t for all t ∈ [0, T ], it follows from Theorem 2 thatE I + E D ≲ m 2 2 T + d N -1 ∑ j=0 γ 2 j (T -t ← j ) 2 + (m 2 2 + d)( ∑ T -t ← j ≥1 γ 2 j (T -t ← j ) 3 + ∑ T -t ← j <1 γ 2 j (T -t ← j ) 2 )Polynomial time schedule. we consider Tt ← j = (δ 1/a + (N -j)h) a with h = T 1/a -δ 1/a N and a > 1, γ j = a(δ 1/a + (N -jϑ)h) a-1 h for some ϑ ∈ (0, 1). We have γ j /h ∼ a(Tt ← j )a-1 a andE I + E D ≲ m 2 2 T + da 2 T 1 a δ 1 a N + (m 2 2 + d) a 2 T 1 a δ 1 a N Therefore, to obtain E I + E D ≲ ε, it suffices to require T = Θ( m 2 2 ε) and the iteration complexityN = Ω(a 2 ( m 2 2 δε ) 1 a m 2 2 + d ε ) Once we let δ = σ2 min , T = σ2 max = Θ( m 2 2 ε) and a = ρ, the iteration complexity isN = Ω( m 2 2 ∨ d d ρ 2 ( σmax σmin ) 2/ρ σ2 max ).Compared to exponential time schedule with the EDM choice of variance schedule, this iteration complexity is worse up to a factor ( σmax σmin )1/ρ .Exponential time schedule. we consider γ j = κ(Tt ← j ) with κ = ln(T /δ) N , we haveEQUATIONε ) and the iteration complexityEQUATIONε ) and a = ρ, the iteration complexity isN = Ω( m 2 2 ∨ d d ln ( σmax σmin ) 2 σ2 max ).Compared to exponential time schedule with the EDM choice of variance schedule, this iteration complexity has the same dependence on dimension parameters m 2 , d and the minimal/maximal variance σmin , σmax .Optimality of Exponential time schedule. For simplicity, we assume m 2 2 = O(d). Then under both schedules in [30] and [46] , E I s only dependent on T , and are independent of the time schedule. Both E D s satisfy 2 is convex on the domain x ∈ (0, ∞), according the Jensen's inequality, 2 reaches its minimum when τ j are constant-valued for all j, which implies the exponential schedule is optimal to minimize E D , hence optimal to minimize E D + E I .E D ≲ d N -1 ∑ j=0 γ 2 j (T -t ← j ) 2 ≲ ε Let τ j = ln ( T -t ← j T -t ← j+1 ) ∈ (0, ∞). Then γj T -t ← j = 1 -e -τj and ∑ δ<T -t ← j <T τ j = ln(T /δ) is fixed. Since x ↦ (1 -e -x )∑ δ<T -t ← j <T γ 2 j (T -t ← j )• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [NA] .Justification:Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [NA]Justification:Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [NA] Justification:Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper poses no such risks.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.Otherwise it is no longer GD training but stochastic GD.The exponential integrator scheme is degenerate since ft ≡ 0. Time discretization is applied when we evaluate the score approximations {S(θ; t, Ȳt)}.In EDM[30], they use Pmean = -1.2, P std = 1.2, σ data = 0.5, σmin = 0.002, σmax = 80.This horizontal axis is in log-scale and the plot in regular scale is a little bit skewed, not precisely a "bell" shape. However, we remark that the trend of the curve still matches our theory.This holds under parameters used in either Song et al.[46] or Karras et al.[30].The exponential time schedule under the variance schedule in[30] also has the poly-logarithmic dependence on σmax/σmin. Under both variance schedules in[30] and[46], it can be shown that exponential time schedule is optimal. Details are provided in Appendix I.1.nN 5/2 L 1/2 √ m √ maxj w(tj )(tj -tj-1)σt j ∑ k w(tj )(tj -tj-1)/σt j minj w(tj )(tj -tj-1)σt j) and byLemma 8, 
 
 Appendix:  


Decision: This paper presents an analysis of diffusion models, taking into account both the training and sampling processes. The authors derive non-asymptotic convergence bounds for denoising score matching under gradient descent and obtain a refined sampling error analysis for variance exploding models. This unified analysis provides insights into the design choices of diffusion models, including the selection of noise distribution, loss weighting functions, and time and variance schedules.

The paper has perhaps potential to guide practical applications. However, the results presented therein are complex and it will be difficult for readers to understand their practical implications. While the paper aims to connect theoretical results with practical design choices, it somewhat lacks in concrete examples and guidelines. This is reinforced by the lack of empirical results. Additionally the authors rely on strong assumptions which may not capture the practice of diffusion models. 

So while I recommend the paper be accepted, the authors are strongly encouraged to enhance the clarity and accessibility of the presentation. This could include providing more intuitive explanations of the theoretical results and their practical implications, offering informal versions of the theorems, summarizing key features etc. A more concrete discussion of practical implications should be included. and the authors should demonstrate how their theoretical findings can help guiding the selection of noise distributions, weighting functions, and time and variance schedules, providing specific examples and applications. Presenting detailed experiments to validate these guidelines would improve significantly the paper's impact.

Reviews: review: paper_summary: paper_summary: This paper provides a full error analysis (considering both training and sampling) for a class of score-based SDE diffusion models, and using the results to understand why and when certain time and variance schedules are preferred.
strength_weakness: strength_weakness: - The stated theorems are quite dense and are not easy to parse. The paper would benefit from having informal version of these results highlighting the key qualitative features
- Only an SDE based sampler (the exponential integrator scheme) is considered, while comparisons with other samplers (such as the Probability Flow ODE) are missing
questions: questions: - In Remark 1, it was stated that the KL divergence can be explicitly computed when the data distribution is Gaussian, would you provide such result in the appendix?

Minor remarks:
- typo in line 16: variance -> various 
- missing lemma number in line 201
limitations: limitations: Yes
paper_summary: paper_summary: This work studies the training and sampling process and then achieves the end-to-end analysis for diffusion models. For the optimization process, this work uses an over-parameterized NN and gradient descent to prove the convergence rate. For the sampling process, this work provides VESDE results and explains the great performance of SOTA VESDE. After these results, this work explains the “bell shaped” weight function used in application from the theoretical perspective.
strength_weakness: strength_weakness: Weakness 1: The high-dimensional data assumption is strong. It would be better to discuss the application of high-dimensional data.

Weakness 2: The technique novelty of Theorem 1 is needed to be discussed in detail. Section 3.1 claims that this work develops a new method for obtaining a lower bound of gradient. However, it does not discuss this method on the main page. It would be helpful to discuss the technique novelty in detail.

Weakness 3: As shown in Section 3.1, this work mentions that the data separability assumption (corresponds to parameter $\delta$ in [1]) can not be directly used due to the property of diffusion models. This work replaces $\delta$ with $t_0/T$ when choosing $m$, which seems to be related to the $E_I$ term in Theorem 2. The technique challenge is unclear after replacing $\delta$ with $t_0/T$. 

Weakness 4: As shown in [2], their results can be straightforwardly extended to any linear SDE, including VESDE. Hence, it would be better to discuss the technique of Theorem 2.	

[1] Allen-Zhu, Z., Li, Y., & Song, Z. (2019, May). A convergence theory for deep learning via over-parameterization. In International conference on machine learning (pp. 242-252). PMLR.

[2] Benton, J., Bortoli, V. D., Doucet, A., & Deligiannidis, G. (2024). Nearly d-linear convergence bounds for diffusion models via stochastic localization.
questions: questions: Please see the Weakness part.

Question 1. Can you relax the high-dimensional data assumption by using the low-dimensional manifold assumption [3]?

Question 2. The guarantee in [4] is a pure $W_2$ guarantee and can not be directly compared to the $KL$ guarantee in this work. Can you discuss it in detail?

Comment 1: It would be better to add a conclusion and limitation part at the end of this work.

Typo 1: In line 201, the reference Lemma is unclear.

[3] Tang, R., & Yang, Y. (2024, April). Adaptivity of diffusion models to manifold structures. In International Conference on Artificial Intelligence and Statistics (pp. 1648-1656). PMLR.
limitations: limitations: This work does not discuss the limitation and societal impace in a independent paragraph.
paper_summary: paper_summary: This work develops theoretical bounds on the training error and sample complexity of score-based diffusion models, which further imply a complete convergence analysis (by combining the both). In addition, the derived bounds shed light on the efficient hyper-parameters selection, including re-weighting functions, time discretization schemes and variance schedules, supporting former empirical works in a principled manner.
strength_weakness: strength_weakness: 1. The paper is organized in a almost technical way, and some details are not easy to follow. It would be helpful to improve the readability by adding e.g., (key) notations collections and formulation illustrations. Also, there is no conclusion section.
2. The derived bounds in main theorems are too technically complex. Is it possible to simplify them to summarize the core algebraic dependence on hyper-parameters (including the data dimension $d$, sample size $n$, model capacity $(m,L)$, discretization steps $N$)? 
3. There are no experiments to support the theory. Particularly, is it feasible to numerically verify the exponential convergence in Theorem 1? 
4. Are there any differences on the training between over-parameterized neural networks applied to normal supervised learning, and here in the generative modeling (as the score network)? Note that this work selects the *same* feed-forward network (without biases) for *all* time steps (see the equation between Line 173 and 174), ignoring the time information in the architecture of score networks.
questions: questions: 1. Please provide more details of questions raised in the weaknesses section above.
2. The scale of model initializations is near-zero (see  Line 181-182), which seems to be one of the typical configurations in the over-parameterization regime (lazy training?). Is this setting reasonable and representative in practice? It would be better to provide (numerical) examples about this, particularly in this generative modeling case. 
3. Details 1: In Line 201, "... high probability (see Lemma..)." 
4. Details 2: What is the definition of $a$ in Theorem 1 (see $d^{2a-1.5-4c}$)?
5. Details 3: In Proposition 1, there seems to be a confusion on the mixed use of $S$ and $S^*$ (particularly, Point 2). 
6. Details 4: In Line 297-299, it states that $\epsilon_n$ (the distance between the empirical loss and population loss at the GD outputs) can be estimated by central limit theorem. Can the authors provide more details on it (since the learned model parameters are *dependent* of data)? 
7. Details 5: In Table 1, why is the exponential schedule selected as $\bar{\sigma}_t=\sqrt{t}$? In VE-SDE, $\sigma_t$ is an exponential function of $t$ (and hence $\bar{\sigma}_t^2:=2\int_0^t \sigma_s^2 ds$). Note that in this case, Assumption 2 (Point 1) may not hold.
limitations: limitations: Based on the developed theory, it is worthy to further explore new weighting functions, time discretization schemes and variance schedules to better tune diffusion models.
paper_summary: paper_summary: This paper analyzes both the training and sampling errors of diffusion models. The analysis sheds light on two practically relevant design choices, which are the noise level weighting (during training) and discretization (during sampling).
strength_weakness: strength_weakness: All the theoretical results are very hard to even parse, let alone extract intuition or guidance for practice. The authors argue that their results provide justification for the design choices of Karras et al. [31], but it seems to me that they are shoehorned to fit the chosen narrative. In other words, it is not clear to me that the exact same error bounds could not be used to justify any other schedule used in practice. I find Section 4.1 in particular to be very handwavy. I appreciate the attempt to connect theoretical results to practical choices, but I think the connection is too tenuous here.

Some claims made in the text seem unsubstantiated to me:
- The authors claim that the NTK can only be analyzed in very restricted settings (lines 103-104). I don't see why this is true, as the NTK applies to all architectures in the lazy regime.
- I don't see how Theorem 1 implies that GD has an exponential convergence (lines 223-224) as product factors depend on $s$.
- I don't see how Theorem 2 implies that the sampling error is almost linear in $d$ as $\gamma_k$ and $\sigma_t$ may depend on $d$.

Minor:
- line 140 and eq (4): $C$ should be $C_t$ as it depends on $t$
- line 150: aspectS
- line 154: by A deep ReLU
- line 176: introduce $W$ before line 174
- line 364: we maintain a fixed ~the~ total weighting
questions: questions: - I don't understand the sentence in lines 296-299, which mentions empirical and population losses even though the paper does not tackle the generalization error. Also, I don't see how the central limit theorem applies here, as individual terms in the empirical average over the training set are not independent since the network parameters themselves depend on the training set. Could the authors please clarify?
- Do the authors believe that their work can inform practical decisions? If yes, what new predictions can be made using their theoretical results? If no (which does not necessarily imply the results are not significant), I think the paper should be more honest about this.
limitations: limitations: Contrary to what is stated in the checklist, limitations are not discussed in the paper. There is no discussion section.

- A first limitation is that only the approximation and optimization errors are studied, and generalization is not tackled. This means that the score network learns the empirical distribution of the training set, and thus reproduces training data points during sampling. This is only briefly mentioned in passing in the middle of the text but should be stated in the introduction (and potentially discussion).

- Some assumptions are questionable. First, the first and last layer are not trained and left to their random initialization. Second, and more importantly, the chosen asymptotic setting corresponds to extremely wide networks trained on very few data points (probably so that they can memorize their training set). This thus corresponds to a rather unrealistic setting.

All in all, I am divided about this paper. On the one hand, it provides a unified analysis of both training and sampling errors, which opens up understanding how to adapt the discretization schedule to score errors. On the other hand, independently of the mathematical correctness of the results (which were not checked), they seem very artificial (due to the restricted setting of the analysis in section 3, and the very qualitative arguments in section 4), and not clearly presented at all. It is thus not clear to me how they can be useful to the theory and practice communities. I recommend rejection, but I am willing to change my decision if the authors convince me that my assessment is incorrect.

Comments: title: Thank you for the rebuttal
Continued response to the reviewer
Discussion
Response to reviewer comments
Continued response to the reviewer
comment: I am satisfied with the rebuttal and will keep my score.
Thanks for the detailed discussion on the technique novelty and limitation. Since the rebuttal addresses my concerns, I will raise my score to $5$.
Thank you for the clarifications. I have increased my score.
> Regarding the dependence on hyper-parameters of the upper bound, the dependence on model capacity ($m,L$) and sample size $n$ seems still unclear, in the sense that the upper bound is not related to the network depth $L$, and also independent of the network width $m$ and sample size $n$ when the total GD iteration $K$ goes to infinity (see the answer of W2).

Thank you for the expert questions. Here are technical explanations. In the simplified bound (for better seeing hyper-parameter dependences) that was provided in our last reply, namely

$KL(p_\delta|q_{T-\delta}) \lesssim  \frac{\mathrm{m}_2^2}{T^2} +\frac{d a^2 T^{\frac{1}{a}} }{\delta^{\frac{1}{a}}N}+(\mathrm{m}_2^2+d)( \frac{a^2T^{\frac{1}{a}}}{\delta^{\frac{1}{a}} N}+\frac{a^3T^{\frac{2}{a}}}{\delta^{\frac{2}{a}}N^2})$
$+\frac{C}{N}(C_1+(1-C_2 h (\frac{md^{\frac{a_0-1}{2}}}{n^3N^2}))^K),$

the constant $C_1$ actually hides many terms: 
$$ C_1=\epsilon_n+\bar{\mathcal{L}}(\theta_\mathcal{F})+|\bar{\mathcal{L}}(\theta_\mathcal{F})+\bar{C}|, $$
where $\epsilon_{n}=|\bar{\mathcal{L}}(\theta^{(K)})-\bar{\mathcal{L}}\_{\rm em}(\theta^{(K)})|$, $-\bar{C}$ is the true minimum of $\bar{\mathcal{L}}(\theta)$ defined in (5), and $\theta_{\mathcal{F}}=\arg\inf_{\{\theta:S(\theta)\in\mathcal{F}\}} |\bar{\mathcal{L}}(\theta)+\bar{C}|$ with $\mathcal{F}=\big\lbrace$ ReLU network function, with $d=\Omega({\rm{poly}}(\log (nN))), m=\Omega\left(\text{poly}\big(n,N,d,L,T/t_0\big)\right)\big\rbrace$. Then, $\epsilon_n$ is the statistical error, $\bar{\mathcal{L}}(\theta_\mathcal{F})$ is the estimation error, and $|\bar{\mathcal{L}}(\theta_\mathcal{F})+\bar{C}|$ is the approximation error.

Regarding the depth dependence $L$, indeed it is not shown in the rate. However, it was actually reflected in the function approximation error $|\bar{\mathcal{L}}(\theta_\mathcal{F})+\bar{C}|$; for example, if $L$ is too small, this error could be large. We did not analyze the function approximation theory of neural network, which is another profound area. Since that area is extensively studied, one can just leverage the existing results, which will give $L$ dependence.

Regarding the disappearing of $m$ and $n$ dependences as $K\to\infty$, here is how the dependences were hidden: 1) there is a lower bound of the width $m$ in Assumption 1; once $m$ is greater than this bound, i.e. the model is sufficiently overparametrized, it no longer has effect on the optimization error in the infinite $K$ limit. 2) Regarding the sample size $n$, its effect is mainly hidden inside the statistical error $\epsilon_n$. Please kindly see more details in the reply below.

> On the estimation error $\epsilon_n$, the analysis is definitely more complex than a trivial application of the law of large numbers, since the data and parameters are decoupled under this setting. Can authors kindly provide at least (promising) ideas to handle this difficulty?

We sincerely apologize for the confusion caused, which is likely due to our previous typo in the Corollary 1, and we deeply appreciate the opportunity to clarify: we did not intend to establish any estimation of the statistical error $\epsilon_n$, and our updated version actually looks like this: 

>> **Theorem.** Under the same conditions as updated Theorem 1 (global rebuttal) with $k=K$ and Theorem 2,
    we have 
    $$
        KL(p_\delta|q_{T-\delta}) \lesssim E_I+E_D+\max_k \frac{\sigma^2_{t_{N-k}}}{w(t_{N-k})} \bigg(\epsilon_{\rm{train}}+\epsilon_{n}+\bar{\mathcal{L}}(\theta_\mathcal{F})+|\bar{\mathcal{L}}(\theta_\mathcal{F})+\bar{C}|\bigg)
    $$
    where $E_I,E_D$ are defined in Theorem 2, $\epsilon_{\rm train}$ is defined in updated Theorem 1 (global rebuttal), $\epsilon_{n}=|\bar{\mathcal{L}}(\theta^{(K)})-\bar{\mathcal{L}}\_{\rm em}(\theta^{(K)})|$, $\bar{C}$ is defined in (5), and $\theta_{\mathcal{F}}=\arg\inf_{\{\theta:S(\theta)\in\mathcal{F}\}} |\bar{\mathcal{L}}(\theta)+\bar{C}|$ with $\mathcal{F}=\big\lbrace$ReLU network function, with $d=\Omega({\rm{poly}}(\log (nN))), m=\Omega\left(\text{poly}\big(n,N,d,L,T/t_0\big)\right)\big\rbrace$.

We'd deeply appreciate it just in case the reviewer is willing to shed more light on $\epsilon_n$'s estimation. In particular, the reviewer pointed out that the analysis is more complex than a trivial application of the law of large numbers, since the data and parameters are decoupled. May we confirm if "parameters" refer to neural network parameters, or training hyperparameters? If the latter, does the decoupling matter? If the former, we thought after one step of training, the NN parameters become coupled with the data. Or, could "decoupled" be, by any chance, a typo, and "coupled" was intended instead? We thought law of large numbers is actually easier when random variables are independent from each other, and when there are coupling (i.e. correlation) it could be more tricky.


---
In any case, thank you again for your time and consideration!
We sincerely thank the reviewer for discussing with us and helping us improve the quality of our manuscript.

> Although the current manuscript can not be refined, would you please state the plan of notations collections, formulation illustrations, and the conclusion section?

Apologize for not having provided enough plan in the rebuttal. Here it is.

For notation collections, we will add a detailed table in the appendix summarizing all our symbols for readers' easier references.

For the formulation illustrations, we will add more remarks including: 1) more discussions on the special structure of this denoising score matching problem and how it is related and difference from traditional nonlinear regression problems (the former is not interpolating), 2) comparison between our assumptions and results of this model and previous overparameterized neural network models, especially weakened assumption (we do not require data separability as shown in Allen-Zhu et al. [1]), a decomposition of the non-interpolation setting in diffusion models compared to the interpolation ones used in previous models, and improved parameter dependence on the input dimension $d$ in the error bound, 3) better quantitative interpretation of our results under specific schedules with explicit parameter dependences.

The conclusion section will start with the following: We provide a first full error analysis incorporating both optimization and sampling processes. For the training process, we provide a first result under a deep neural network and prove the exponential convergence into a neighborhood of minima, while for sampling, we also extend the current analysis to the variance exploding case. Moreover, based on the full error analysis, we establish quantitative understandings of the error bound under the two schedules. Consequently, we conclude a qualitative illustration of the "bell-shaped" weighting and the choices of schedules under well-trained and less-trained cases.

Then we will discuss limitations and future directions. More precisely, the network architecture we used in the model is a deep ReLU network. Although being so far the most complicated arhitecture for theoretical results, it is still far from what is used in practice like U-Nets and transformers. Regarding the full error analysis, we only focus on the optimization and sampling error, and do not disect the generalization error. When bridging the theoretical results with practical designs of diffusion models, our results are mostly qualitative and we only compare two schedules under well-trained and less-trained cases.

---

**We apologize for exceeding the character limit and would like to kindly direct the reviewer to our next Comment for the continuation of our reply.**
We deeply thank the reviewer for reading our rebuttal and discussing with us.

> The significant improvement of the dependence of  is helpful. However, the proof sketch of the improved theorem is not provided. It would be better to provide an intuitive proof sketch, which is necessary to check the correctness.

The main improvement is on proving the lower bound of the gradient. Therefore, we will provide a proof sketch for it, and the rest of the proof follows the similar framework as Allen-Zhu et al [1].


We first decompose the gradient of the $k$th row of $W_L$ $\nabla\_{(W_L)\_k}\bar{\mathcal{L}}\_{\rm em}(\theta)=\underbrace{\frac{1}{n}w(t_{j^*})(t_{j^*}-t_{j^*-1}){(W_{L+1})^k}^\top(\bar{\sigma}\_{t_{j^*}}W_{L+1} q_{i\^\* j\^\*,L}+\xi_{i\^\* j\^\*})q_{i\^\* j\^\*,L-1} {1}\_{(W_L q_{i\^\* j\^\*,L-1})\_k>0}}\_{\nabla_1}$
    $+\underbrace{\frac{1}{n}\sum_{(i,j)\ne (i^*,j^*)}w(t_j)(t_j-t_{j-1}){(W_{L+1})^k}^\top(\bar{\sigma}\_{t_j}W_{L+1}q_{ij,L}+\xi_{ij})\,q_{ij,L-1} {1}\_{(W_L q_{ij,L-1})\_k>0}}_{\nabla_2}$

where $(i^*,j^*)$ indicates the sample index with the largest loss value.


Then we first fix $(q_{ij,L-1})\_s=1$, and prove that the index set of both $(q\_{i^\*j^\*,L})\_s> 0$ and $\sum\_{(i,j)\ne (i^*,j^*)}w(t_j)(t_j-t_{j-1})\bar{\sigma}\_{t_j}{1}
\_{(W_Lq_{ij,L-1})\_k>0}(q_{ij,L})_s>0$ is order $m$ with high probability.

Next, we conditioned on the index set we've found, then we can decouple each element of $\nabla_{(W_L)_k}\bar{\mathcal{L}}\_{\rm em}$ with high probability. We then prove that with high probability, the event $(\nabla_1)_s> 0$ and $(\nabla_2)_s>0$ has probability at least of order $d^{(a_0-1)/2}$ where $a_0\in(1/2,1)$.

Now, we deal with $(q_{ij,L-1})\_s$ and prove that if the above results hold for $(q_{ij,L-1})_s=1$, then there exists an index set with cardinality of order $m/(nN)$ such that $(\nabla_1)_s> 0$ and $(\nabla_2)_s>0$ also hold in this index set.

In the end, combining all the steps above yields the lower bound.


> Furthermore, I think it would be better to discuss the limitation in the rebuttal phase and add an independent limitation paragraph in the next version.

We thank the reviewer for the suggestion. We will add a limitation paragraph to our paper. Below is a quick summary:

The network architecture we used in the model is a deep ReLU network. Although being so far the most complicated arhitecture for theoretical results, it is still far from what is used in practice like U-Nets and transformers. Regarding the full error analysis, we only focus on the optimization and sampling error, and do not disect the generalization error. When bridging the theoretical results with practical designs of diffusion models, our results are mostly qualitative and we only compare two schedules under well-trained and less-trained cases. Finer analysis techniques are needed for providing a quantitative demonstration of the parameters in diffusion models, as well as for motivating more designs of weightings and schedules. We will leave them for future exploration.
Please participate to the discussion. Thanks
Thanks for your careful response. The significant improvement of the dependence of $d$ is helpful. However, the proof sketch of the improved theorem is not provided. It would be better to provide an intuitive proof sketch, which is necessary to check the correctness. 

Furthermore, I think it would be better to discuss the limitation in the rebuttal phase and add an independent limitation paragraph in the next version.
Thanks for the clarifications by authors. I find most of them adequate, except the following three points: 
 - Although the current manuscript can not be refined, would you please state the plan of notations collections, formulation illustrations, and the conclusion section?
- Regarding the dependence on hyper-parameters of the upper bound, the dependence on model capacity $(m,L)$ and sample size $n$ seems still unclear, in the sense that the upper bound is not related to the network depth $L$, and also independent of the network width $m$ and sample size $n$ when the total GD iteration $K$ goes to infinity (see the answer of W2). 
- On the estimation error $\epsilon_n$, the analysis is definitely more complex than a trivial application of the law of large numbers, since the data and parameters are decoupled under this setting. Can authors kindly provide at least (promising) ideas to handle this difficulty?
We sincerely thank the reviewer for carefully reading our rebuttal and being willing to increase the score.

> This first contribution...

Thank you for the comment and apologize for potential confusions. The main goal of this work is not to revamp the general methodology for analyzing neural network training, but rather to adapt existing tools to obtain a full theoretical guarantee of diffusion model's generation capability. The former is a holy grail, but the latter, we believe, is also rather important. Nevertheless, although the sampling process of diffusion model has been analyzed for ~2 years already and relatively well understood, it is only half of the generation pipeline, and no result existed that combines it with the other half, namely the score training. In fact, denoising score matching training has some significant difference from settings in which existing training analyses can directly apply, and we suspect this is why this problem has not been already solved before. Please allow us to elaborate:

For the denoising score matching problem, there are a lot of special properties and additional issues. Specifically, the input and output data in the denoising score matching problem have specific structures such as: 

1) The input data $X_{ij}=x_i+\bar{\sigma}\_{t_j}\xi_{ij}$ is noisy and therefore the whole model is non-interpolation, which is in contrast to the interpolation assumption in existing techniques.
2) The scales of the input and output data cannot be freely assumed in the theory. More precisely, the variances of the noise $\bar{\sigma}\_{t_j}$ in the input and output data $-\xi_{ij}/\bar{\sigma}\_{t_j}$ changes as $\bar{\sigma}\_{t_j}$ varies. Moreover, to favor the sampling process after training the score, the scale of input data needs to be of order $\sqrt{d}$, and the scale of outputs is at least of order $\sqrt{d}$. This is in sharp contrast to data assumptions in existing literatures. For example, in Allen-Zhu et al. [1], they assume the output data to be $o(1)$ and the input data to have order 1 scale.


Therefore, while still leveraging the framework of existing theoretical analyses, we had to develop additional techniques to deal with the above issues. For 1, we decomposed the model into an interpolation component and a non-interpolation one in order to apply the existing techniques on the interpolation part. For 2, due to the special data structure and scalings, existing techniques, which rely on a commonly used data separability assumption, no longer suffice. Instead, we develop a new geometric technique to overcome this issue, and the only price to pay is that the order of the dimension of the input data $d$ should not be too small, which has now been significantly relaxed in the new version as stated in the global rebuttal.

Hopefully our goal and contributions are better clarified by the above discussions. 

> The second contribution...

We apologize for the confusion caused. Our full error analysis is completely quantitative, and we will now explain why our presentation on preferred schedule(s) is qualitative. This explanation could start with a question: why can't we just substitute the hyperparameter choices of a specific schedule in our error bound, get a concrete number, and then compare that number with other schedules' numbers, or even optimize that number to get better schedules?

That is because our full error bound is a function of not only hyperparameters provided by a schedule, but also many other things (i.e. what the reviewer mentioned as "a lot of wiggle room"). For example, a) there is a dependence on the data distribution, which is reflected in a generalization gap term in the error bound; b) there is dependence on the network structure (for our setting, width $m$ and depth $L$); c) there is dependence on optimization hyperparameters. If all these were known, then a number could indeed be spit out, but then, we feel, it becomes a bit less interesting because we are no longer comparing general schedules per se.

Therefore, we don't compare precise values, and instead adopt some qualitative arguments, so that some generic yet practical guides can still be provided. More precisely, we consider the following two general situations: $E_S\gg E_I+E_D$ and $E_S\ll E_I+E_D$.  Then we substitute the two specific time and variance schedules into the error bound, and compare either the values or the complexities of the time steps under the same error as is shown in Table 2. The schedule with smaller value or complexity will be the preferable choice. The results in Table 2 may appear informal; however, they are indeed rigorously calculated in Appendix E3 and F (see lines 844,846,853-875), and could be formalized into corollaries or propositions. Here we chose not to present the full details in Table 2 merely due to the page limit. We will clarify this part in a revision.
I thank the authors for their very detailed answer. This fully addresses my concerns about the validity of the results. After reading the other reviews and rebuttals, I am willing to increase my score.

The main point that is still not clear for me is what we learn from the results in this paper.

- This first contribution is the quantitative analysis of the learning and sampling errors. While sampling is well understood, learning is not, and the analysis in the paper suffers from the same limitations as the rest of the literature (although relaxing the very-wide assumption in the new theorem 1 above seems to be a major advance?). I suppose the main contribution lies in the novel gradient lower bound, which is relegated to the appendix, but that readers more expert than this reviewer may read and appreciate. 
   * As a side remark, if theoretical progress is slow, then should we keep pushing in the same direction and keep submitting results with the same limitations but slightly more general assumptions every conference cycle? Or rather try to tackle these fundamental limitations with new approaches? I'm not sure that the first option is the most efficient way to drive progress in theoretical research. But this is more of a systematic issue of the field and beside the point here.

- The second contribution is qualitative analysis of noise schedules in the literature. As the authors pointed out, this paper predicts that two different sampling schedules should be used based on the quality of the score approximation, which is an interesting hypothesis. My main qualm here still stands: the authors aim for quantitative statements, but it seems that the error analyses leave a lot of wiggle room in the epsilons and deltas to fit any desired schedule. To be clear, I see this as a negative point: it does not tell us which schedule is the best, and not even how to evaluate a given schedule beyond saying that it should have a few qualitative features. I am all in favor of qualitative analyses, which I are usually underrated in theoretical papers, but I think they should be _simple_. It would greatly strengthen the papers if the authors studied a simplified setting where the error bounds take a simple form, and it can be clearly seen by a non-expert reader what features of a sampling discretization matter. In its current state, the paper does not give _intuition_ about the problem, and readers have to trust the authors that their qualitative results are an accurate description of typical phenomena.
> Q1. I don't understand the sentence in lines 296-299, which mentions empirical and population losses even though the paper does not tackle the generalization error. Also, I don't see how the central limit theorem applies here, as individual terms in the empirical average over the training set are not independent since the network parameters themselves depend on the training set. Could the authors please clarify?

We deeply apologize for the confusions. Yes, this paper focuses on optimization and sampling errors and does not tackle generalization error. The reason we mention the population loss and empirical loss is the following. In training, we need to optimize the empirical loss. However, when estimating the whole generation error, the population loss is necessary (see Theorem 2, the term $E_S$). Therefore, the gap between the population and empirical versions is integrated in the full error analysis (Corollary 1). We will further clarify this point in our next version.

For the $\epsilon_n$, we sincerely apologize; it was our mistake, and we forgot to correct it. The "central limit theorem" part should be the law of large numbers, i.e., this $\epsilon_n$ is the statistical error that converges to 0 as $n\to \infty$. We currently do not have refined statistical error estimation, and this limitation will be clarified in a revision. 

> Q2. Do the authors believe that their work can inform practical decisions? If yes, what new predictions can be made using their theoretical results? If no (which does not necessarily imply the results are not significant), I think the paper should be more honest about this.



Thank you for a very important question. We focused a general theory and did not explore all possible ramifications, but here is one: when the score is very well trained, it is more preferrable to use the schedules in [Song et al. 2020] to sample/generate; otherwise, the design in [Karras et al. 2022] is more preferable.

This important point is already clarified in a revision. However, if the reviewer wants us to state in the limitation section that we presented a general theory but did not explore enough practical implications, we are happy to do so.





> Contrary to what is stated in the checklist, limitations are not discussed in the paper. There is no discussion section.
> * A first limitation is that only the approximation and optimization errors are studied, and generalization is not tackled. This means that the score network learns the empirical distribution of the training set, and thus reproduces training data points during sampling. This is only briefly mentioned in passing in the middle of the text but should be stated in the introduction (and potentially discussion).

We sincerely apologize for not clarifying it in a separate section. In our next version, we will add a discussion section and state clearly that we do not tackle the generalization in the introduction. 

> * Some assumptions are questionable. First, the first and last layer are not trained and left to their random initialization. Second, and more importantly, the chosen asymptotic setting corresponds to extremely wide networks trained on very few data points (probably so that they can memorize their training set). This thus corresponds to a rather unrealistic setting.

We would like to thank the reviewer for this comments. 

First, the assumption that "first and last layer are not trained" is equivalent to applying fixed linear transforms respectively to the input and the output of the network. This will not affect the trainability of the model if the depth $L$ and width $m$ are properly chosen. Also, it is a commonly used assumption in theoretical analysis (see e.g. [2,3]; also lines 179-180 in the paper). 

Second, we agree with the reviewer that "wide networks trained on very few data points" is a little bit "unrealistic". However, this is the limitation of all the current theoretical works about the convergence of neural networks due to technical difficulties. Existing works (e.g. [1]) are making efforts to get closer to practical settings but the progress is very slow. If there are important references we missed, we'd deeply appreciate it; otherwise, we sincerely hope that the reviewer could consider if we could agree about the current stage of theoretical research, and allow us to explore in a sustainable way (a 3:Reject adds one straw instead; is our work really that bad, especially in the context of the existing literature?)



[1] Zou et al. An improved analysis of training over-parameterized deep neural networks. NeurIPS 2019.

[2] Allen-Zhu et al. A convergence theory for deep learning via over-parameterization. ICML 2019

[3] Cai et al. Neural temporal-difference learning converges to global optima. NeurIPS 2019.