On the Relation Between Linear Diffusion and Power Iteration

## Abstract

Abstract Recently, diffusion models have gained popularity due to their impressive generative abilities. These models learn the implicit distribution given by the training dataset, and sample new data by transforming random noise through the reverse process, which can be thought of as gradual denoising. In this work, we examine the generation process as a “correlation machine”, where random noise is repeatedly enhanced in correlation with the implicit given distribution.
To this end, we explore the linear case, where the optimal denoiser in the MSE sense is known to be the PCA projection. This enables us to connect the theory of diffusion models to the spiked covariance model, where the dependence of the denoiser on the noise level and the amount of training data can be expressed analytically, in the rank-1 case.
In a series of numerical experiments, we extend this result to general low rank data, and show that low frequencies emerge earlier in the generation process, where the denoising basis vectors are more aligned to the true data with a rate depending on their eigenvalues. This model allows us to show that the linear diffusion model converges in mean to the leading eigenvector of the underlying data, similarly to the prevalent power iteration method.
Finally, we empirically demonstrate the applicability of our findings beyond the linear case, in the Jacobians of a deep, non-linear denoiser, used in general image generation tasks.

## 1 Introduction

Recently, diffusion models have gained much popularity as very successful generative models, showcasing impressive performance in image generation tasks (Dhariwal & Nichol, ; Ho et al., ; Song & Ermon, ; Song et al., ). These models learn the implicit distribution given by the training dataset, and sample new data by transforming random noise inputs through a reverse diffusion process, which can be thought of as gradual denoising. More formally, it has been shown in Kadkhodaie et al. () that learning the underlying distribution is equivalent to optimal denoising at all noise levels.

In order to shed more light onto the mechanism behind the success of diffusion models, in this work we analyze the behavior of denoisers in the context of image generation, where pure noise is gradually processed into a sample from a given (implicit) distribution by gradual denoising. Unlike other works, e.g. Kadkhodaie et al. (), we focus on the denoiser(s) throughout the generation process, and not only on the final generated data.

To this end, we suggest the following simple model to illustrate our point.
Consider the class of linear denoisers, where the optimal denoiser is given by a PCA projection. To simulate the diffusion generation process, we learn a series of projections onto noisy data at different noise levels, and use them to transform pure noise into samples from the underlying distribution. Given this simple model we can inspect the evolution of eigenvectors spanning gradual projections with decreasing noise levels, as well as the distribution of the generated data samples.

We show that the correlation of the noisy basis eigenvectors with their clean version decays as the noise level increases, with a rate determined by the eigenvalues and the size of the training dataset. In other words, we show that low frequencies, corresponding to large eigenvalues, emerge earlier in the reverse process as was empirically observed in (Ho et al., ), and analyze how more training data contributes to generalization (Kadkhodaie et al., ). Analytically, this corresponds to the spiked covariance model (Johnstone, ), in which we bound this decay for the leading eigenvector (corresponding to the largest eigenvalue).

Next, we demonstrate the applicability of our findings to more general, non-linear deep denoisers. Although the network is not linear, its application can be written as a linear operation of the Jacobian calculated on the input image. We empirically show that the aforementioned decay of eigenvector correlations is prevalent also in in the Jacobians of a deep denoiser, in the final stages of image generation, thus showing the relevance of our analysis in a broader context, and not just in a simplified linear case.

## 2 Background and Related Work

Since their introduction by Sohl-Dickstein et al. (), diffusion models have been vastly used in image generation tasks (Dhariwal & Nichol, ; Ho et al., ; Song & Ermon, ; Song et al., ), more general computer vision tasks (Amit et al., ; Baranchuk et al., ; Brempong et al., ; Cai et al., ), and in other domains such as natural language processing (Austin et al., ; Hoogeboom et al., ; Li et al., ; Savinov et al., ; Yu et al., ) and temporal data modeling (Alcaraz & Strodthoff, ; Chen et al., ; Kong et al., ; Rasul et al., ; Tashiro et al., ). On top of their practical success, different flavors of training and sampling have risen based on interesting theoretical reasoning, e.g., considering the statistical properties of the intermediate data (Song et al., ; Sohl-Dickstein et al., ), or by framing the problem in the form of stochastic differential equations (SDEs) (Karras et al., ; Song et al., ; ; Chen et al., ) or score based generative models (Song & Ermon, ; ). In this work, we look at diffusion models in the context of iterative denoising, and focus on the properties of the learned denoiser (Milanfar & Delbracio, ).

Recently Kadkhodaie et al. () showed that the learned denoising functions are equivalent to a shrinkage operation in a basis adapted to the underlying image. In this sense the diffusion denoiser is an adaptive filter (Milanfar, ; Talebi & Milanfar, ; ). While they focus on the analysis of the nonlinear denoiser at the point of the final generated data, we are interested in the evolution (adaptation) of the denoiser throughout the generation process, and its dependence on the noise level. To this end, we suggest a simple linear denoising model, presented in Section 3. In this case, the (optimal) denoiser does depend on the underlying image, and its dependence on the noise level can be traced analytically, as we show hereafter.

Due to their phenomenal empirical success, some attempts have been devoted towards providing theory supporting the sample and iteration complexity of diffusion models. The current body of work can be generally parted to attaining iteration complexity bounds assuming approximately accurate scores (Li et al., ; ; Chen et al., ; Huang et al., ; Benton et al., ), and to assessing the sample complexity to learn the score functions (Chen et al., ; Block et al., ; Biroli & Mézard, ).
Among these works, many assume a low dimensional data distribution (Bortoli, ; Li & Yan, ; Oko et al., ; Chen et al., ; Wang et al., ), which is a reasonable assumption in theoretical works. Yet, it might particularly explain the gap between the current iteration bounds and the much lower complexity apparent in practice (Li & Yan, ).
In our work, we consider linear models and deduce a linear sample complexity bound associated with learning the score function in Sec. 4 and discuss the tradeoffs of the synthesis conversion rate in Sec. 4.1.
The previous works mentioned above mainly develop bounds assuming specific samplers and scaling details, which differ from our setting. In addition, they generally bound the Total Variation distance (under varying assumptions on the target distributions), which is not trivial to translate to the generated covariance matrix that we focus on even in the linear Gaussian case (Devroye et al., ). The difference in our setting enables us to connect the theory of diffusion models to a broad body of work concerning the spiked covariance model (Johnstone, ), and supports the analysis of denoising diffusion as a correlation machine, which is the main purpose of this paper.

In the setting of Statistical Mechanics, Biroli & Mézard () analyse diffusion models in very large dimensions, focusing on the Curie-Weiss model of ferromagnetism. As an introduction to their work, they also discuss a simple linear score model, in the context of the sample complexity of learning the score function. They focus their discussion on the case of Gaussian data, where the eigenvalues of the covariance matrices can be typically characterised. Unlike their work, we consider data that reside in a low dimensional subspace, with no specific distribution, described in Sec. 4.

## 3 Linear Diffusion - Problem Setup

For our analysis, we define the following simple iterative linear generation model. First, define the standard diffusion model.
Let $q$ denote the natural data distribution and let $x_{0}\sim q$ be a sample from the natural data ($x\in\mathbb{R}^{d}$). The forward (diffusion) process is defined (Ho et al., ) by

$$ $q(x_{t}|x_{t-1})=\mathcal{N}(\sqrt{1-\beta_{t}}x_{t-1},\beta_{t}\mathbf{I})$ (1) $$

for some fixed noise schedule $\{\beta_{t}\}_{t=1}^{T}$ and $x_{0}\sim q$. It can be shown that

$$ $q(x_{t}|x_{0})=\mathcal{N}(\sqrt{\bar{\alpha}_{t}}x_{0},(1-\bar{\alpha}_{t}) \mathbf{I}),$ (2) $$

where $\alpha_{t}=1-\beta_{t}$ and $\bar{\alpha}_{t}=\Pi_{s=1}^{t}\alpha_{s}$.
For our simplified model, consider the process (without scaling),

$$ $q(x_{t}|x_{t-1})=\mathcal{N}(x_{t-1},\sigma^{2}_{t}\mathbf{I}).$ (3) $$

This implies that $x_{t}=x_{t-1}+\epsilon_{\sigma_{t}}$, where $\epsilon_{\sigma_{t}}\sim\mathcal{N}(0,\sigma^{2}_{t}\mathbf{I})$ for some fixed noise schedule $\{\sigma_{t}\}_{t=1}^{T}$. We discard the scaling to comply with previous analysis of the spiked covariance model (Nadler, ) (more details in Section 4). This corresponds to the ”Exploding Variance” formulation, used with Langevin dynamics to sample data as a variant of score based diffusion models (Song & Ermon, ; Song et al., ; Song & Ermon, ). We choose to present the ”standard” diffusion models in the setting of denoising diffusion Ho et al. () and not using the score-based approach entirely, as we focus our discussion on the qualities of the denoiser.

The reverse (generation) process is defined using a parameterized distribution model $p_{\theta}$, generally defined by the Markov process

$$ $\displaystyle p_{\theta}(x_{0:T})=p(x_{T})\Pi_{t=1}^{T}p_{\theta}(x_{t-1}|x_{t }),$ (4) $\displaystyle p_{\theta}(x_{t-1}|x_{t})\triangleq\mathcal{N}(\mu_{\theta}(x_{t },t),\Sigma_{\theta}(x_{t},t)),$ (5) $$

where $p(x_{T})=\mathcal{N}(0,\mathbf{I})$. By choices of parametrization and loss manipulations (see (Ho et al., )), one generally learns to estimate the error $\epsilon_{\theta}(x_{t},t)$, where

$$ $\mu_{\theta}(x_{t},t)=\frac{1}{\sqrt{\alpha_{t}}}\Big{(}x_{t}-\frac{\beta_{t}} {\sqrt{1-\bar{\alpha}_{t}}}\epsilon_{\theta}(x_{t},t)\Big{)},$ (6) $$

$\Sigma_{\theta}(x_{t},t)=e^{2}_{t}\mathbf{I}$, and $e_{t}$ is a designed schedule (usually chosen to be equal to $\sigma_{t}$). Thus, the reverse process can be expressed as a denoising chain

$$ $D_{t}(x_{t})=\frac{1}{\sqrt{\alpha_{t}}}\Big{(}x_{t}-\frac{\beta_{t}}{\sqrt{1- \bar{\alpha}_{t}}}\epsilon_{\theta}(x_{t},t)\Big{)}+e_{t}z,$ (7) $$

where $z\sim\mathcal{N}(0,\mathbf{I})$ and $z_{1}=0$. This is a stochastic denoiser, which preserves the Markovian property of the forward process. Later versions suggested similar (non Markovian) deterministic denoisers, e.g., DDIM (Song et al., ), or more general stochastic denoiser chains, for a continuous forward model (InDI (Delbracio & Milanfar, )).

In our case, we restrict the denoisers to be a linear function of $x_{t}$. Thus, the optimal denoiser (in the $\ell_{2}$ sense) is given by the PCA projection onto the target distribution. For the reverse process, we learn a simple PCA denoiser (projection) based on $X_{t-1}$, which is the cleaner version of the training set $X=\{x_{1},...,x_{n}\}$ at time $t-1$. Thus, at each time step we learn

$$ $D^{t}_{PCA}(x_{t})=D^{t}_{PCA}(x_{t};X_{t-1})=P_{t}(x_{t};X_{0}+E_{\bar{\sigma }_{t}}),$ (8) $$

where each column in $E_{\bar{\sigma}_{t}}$ is distributed by $\mathcal{N}(0,\bar{\sigma}^{2}_{t}\mathbf{I})$ and $\bar{\sigma}_{t}$ is a function of $\{\sigma_{s}\}_{s=1}^{t}$. Our simple denoising procedure is based on the sequential application of $P_{t}\in\mathbb{R}^{d\times d}$, which is the projection on perturbed principal components with respect to the clean data distribution $q$. It is a deterministic denoiser given the sampling of training data and noises, which does not depend neither on $x_{t}$ nor on $x_{0}$. Nevertheless, this model is relevant in differentiable environments of more complex settings such as DNN based denoisers, as we show in Section 5.

Empirical Demonstration of a Linear Diffusion Model.
To illustrate the forward and backward processes in the linear case, we perform a numerical simulation using the MNIST dataset, which is simple enough to be estimated via a linear model. We start here with the training and generation procedures, and use the same setting and trained denoiser to demonstrate our findings throughout the paper.

In the following experiment we simulate the process described above using the MNIST dataset (we use the default train / test splits). In the class conditioned case, we learn a PCA denoiser with 30 components for each time step where $x_{t}=x_{t-1}+\epsilon_{\sigma_{t}}$, $\sigma_{t}\propto t$, and $T=65$ iterations. Figure 1 shows a (decimated) example of digit generation from pure noise, where we apply the sequence of denoisers $D^{t}_{PCA}=P_{t}$, which will be more accurately defined in Section 4.

Figure: Figure 1: Digit generation from pure noise (class conditioned). The reverse process runs from left to right, top then bottom.
Refer to caption: extracted/5930746/figures/digit_generation_30.png

In order to understand the reverse process, we now turn to analyze the gradual change of $P_{t}$, that might be expressed by the angle between the clean and noisy components over time.

## 4 Linear Diffusion as Basis Perturbation

We now turn to analyze the linear model presented above and show how the generation process can be seen as a kernel “correlation machine”. Specifically, we are interested in the temporal (i.e., noise level) dependence of the denoiser Equation 8 throughout the generation process.
Recall that at each time step $x_{t}=x_{t-1}+\epsilon_{\sigma_{t}}$, where $\epsilon_{\sigma_{t}}\sim\mathcal{N}(0,\sigma^{2}_{t}\mathbf{I})$ (Equation 3). Since the noise is assumed to be Gaussian, we can write $x_{t}=x_{0}+\epsilon_{\bar{\sigma}_{t}}$, where $\bar{\sigma}_{t}=\sqrt{\sum_{i=0}^{t}\sigma_{i}^{2}}$. Assume that the data distribution is such that its population covariance is given by

$$ $\displaystyle\Sigma_{t}$ $\displaystyle=\mathbb{E}x_{0}x_{0}^{\dagger}+\bar{\sigma}^{2}_{t}\mathbf{I}= \sum_{i=0}^{r-1}\lambda^{2}_{i}u_{i}u_{i}^{\dagger}+\bar{\sigma}^{2}_{t} \mathbf{I}\triangleq\Sigma_{0}+\bar{\sigma}^{2}_{t}\mathbf{I},$ (9) $$

where $r-1<d$, i.e., the data reside in a low dimensional subspace (which is generally true for natural data). This is known as the “spiked model” (Johnstone, ), with a vast body of work covering the distribution and identifiability of the spikes spectrum (e.g., (Nadler, )). Throughout the paper, we use the term ”index” to refer to to the index $i$ in 9, where the eigenvalues $\lambda_{i}$ are ordered largest to smallest.

Given $n$ samples concatenated as columns in the matrix $X_{0}$, at each time step we learn the PCA basis associated with $X_{t}=X_{0}+E_{\bar{\sigma}_{t}}$, by the diagonalization of the sample covariate matrix

$$ $\displaystyle\hat{\Sigma}_{t}$ $\displaystyle=\frac{1}{n}X_{t}X_{t}^{\dagger}=\frac{1}{n}(X_{0}X_{0}^{\dagger} +X_{0}E_{\bar{\sigma}_{t}}^{\dagger}+E_{\bar{\sigma}_{t}}X_{0}^{\dagger}+E_{ \bar{\sigma}_{t}}E_{\bar{\sigma}_{t}}^{\dagger})\triangleq U_{t}S_{t}U_{t}^{ \dagger}.$ (10) $$

Thus, during the reverse process, at each time step we apply the projection

$$ $D^{t}_{PCA}=P_{t}=U_{t}U_{t}^{\dagger}.$ (11) $$

In order to understand the diffusion generation process, we analyze the decay of the product $\langle u_{i}^{t},u_{i}\rangle$ over time. Note, that there are two drivers of change in the perturbation of $u_{i}$ to $u_{i}^{t}$. The first being the added noise, i.e., $\left\lVert\Sigma_{t}-\Sigma_{0}\right\rVert$. This is the key in the diffusion process and our main focus. The second, is in the finite sample approximation $\left\lVert\hat{\Sigma}_{t}-\Sigma_{t}\right\rVert$. This source of error is interesting in the context of sample complexity, as it encompasses the approximation of the denoiser learned from a finite dataset, the equivalent of the sample complexity of learning the score function in (Chen et al., ; Block et al., ; Biroli & Mézard, ).

For the rank-1 case, (Nadler, ) presented a finite sample theorem which holds with high probability for the closeness between the leading eigenvalue and eigenvector of sample and population PCA under a spiked covariance model similar to Equation 9. They bound the angle between the leading empirical eigenvector and its population counterpart with approximately $\mathcal{O}(d)$ sample complexity, and a linear dependence on the noise level:

$$ $\mathbb{E}\sin{\theta_{\text{PCA}}}=\mathbb{E}\sqrt{1-\langle u^{t},u\rangle^{ 2}}\approx\frac{\bar{\sigma}_{t}}{\lambda}\sqrt{\frac{d}{n}},$ (12) $$

where $\bar{\sigma}_{t}$ is assumed to be small and $d\gg 1$. This result shows that the leading eigenvector rotates in a rate proportional to the noise level. Our numerical experiments on the MNIST dataset (detailed in Section 4.1) show that this is a good approximation in practice, also for the rank-r case (Fig. 2).

Figure: Figure 2: The sine of the angle between the clean principal components and their noisy versions, colored by the order of the eigenvalues (the darkest being largest eigenvalue). Low frequencies emerge earlier in the generation process (at higher noise levels). This motivates Assumption 4.1, that extends Equation 12 to higher ranks.
Refer to caption: extracted/5930746/figures/theta_pca_mnist.png

Notice that in Equation 12 the angle is inversely linked to the eigenvalue, inferring a slower change with higher eigenvalues. In the reverse process, we gradually move from pure noise or high noise levels to smaller noise variance. Given the lower slope of the components corresponding to larger eigenvalues, we interpret the result in Fig. 2 as the earlier emergence of low frequencies in the generation process. The first component to be visible in the generated image is the one with the largest eigenvector, as it is the first one that shows a correlation in high noise levels. Throughout the generation process, when the noise level decreases, the next components take presence, by the order of their associated eigenvalue - from the larger to the smaller. Finally, the components with the smallest eigenvalues appear when the noise level is low.

In the linear case, Equation 12 shows that the diffusion model’s sample complexity is determined by the sample complexity of PCA, with a linear dependence on the dimension of the data.
To further enhance our understanding of the relationship between the amount of training data and the generalization of the diffusion model, we repeat the experiment with varying datasets sizes. Figure 3 shows the angle to noise profile for selected principal components, with the indices $0,5,10$ (left to right; index 0 0 corresponds to the largest in a list of ordered eigenvalues). Increasing the amount of training data improves robustness to noise and enables the emergence of higher frequency components at higher noise levels, thereby capturing more nuances in the generated data.

Figure: (a)
Refer to caption: extracted/5930746/figures/sin_pca_0_different_sizes.png

### 4.1 The Generated Distribution

We now turn to discuss the distribution of the generated output, and how it relates to the natural data distribution. First, we analyze a generation process by repetitive denoising without additional noise, and show how it relates to power iteration. Then we discuss a similar process only with the injection of noise, reminiscent of other common sampling methods (e..g (Ho et al., )).
Given our linear model, the generation process is essentially the linear transformation given by the matrix

$$ $\mathcal{P}_{T}=\Pi_{t=0}^{T}P_{t}=P_{0}\cdots P_{t}\cdots P_{T}.$ (13) $$

The generated output can be expressed as

$$ $x_{g}=\mathcal{P}_{T}\xi,$ (14) $$

where $\xi\sim\mathcal{N}(0,\sigma_{T})$.
Other than the visual aesthetic of the generated images, we are interested in their distribution, and how well it represents the natural distribution of train images. Thus, we would like to compare the generated covariance $\mathbb{E}{x_{g}x_{g}^{\dagger}}$ to the natural covariance $\Sigma_{0}$.

In this context, a natural comparison is the power iteration (PI) method, which may be used to estimate the leading eigenvector of a matrix. This can be seen as another iterative form of generating data from random vectors. Unlike our projection, in PI we ”project” a random vector onto the entire matrix, i.e. including the eigenvalues. In this case the denoiser would be $D^{t}_{PI}=\Sigma_{0}\forall t$, where we ignore the normalization and focus on the direction of the final vector, since there is no normalization constraint for generated data in diffusion models.

We now turn to show how the reverse process performed by a repeated denoising as in Equation 13 converges in mean to PI. To this end, we make the following assumptions.

###### Assumption 4.1 .

Assume that Equation 12 holds for all eigenvectors, i.e.,

$$ $\mathbb{E}\sqrt{1-\langle u_{i}^{t},u_{i}\rangle^{2}}\approx\frac{\bar{\sigma} _{t}}{\lambda_{i}}\sqrt{\frac{d}{n}},$ (15) $$

for $i=,\dots,r-1$.

This assumption is the extension of Equation 12 to higher ranks, and is motivated by our simulations (Fig. 2). In addition, we make the following assumption regarding the cross products of components of different indices, at consecutive time steps.

###### Assumption 4.2 .

For each index $i$ there exists a time $\tau_{i}$, where for $t\leq\tau_{i}$ and $j\leq i$,

$$ $\mathbb{E}\langle u_{i}^{t},u_{j}^{t+1}\rangle=0.$ (16) $$

In addition, $\tau_{i}>\tau_{j}$ for $i<j$.

This assumption is supported by our simulations in Fig. 5, and will be further discussed hereafter. Assumptions 4.2, 4.1 are an extension of Nadler () to higher ranks. We leave their explicit derivation to future work, and focus on their implications to linear diffusion.

We are now ready to state our main result.

###### Theorem 4.3 (Convergence to Power Iteration) .

Let $\sigma_{t}=\frac{1}{T}$, $t=0,\dots,T$. Assuming 4.2, 4.1, in the limit $T\to\infty$,

$$ $\mathbb{E}{x_{g}x_{g}^{\dagger}}\propto u_{0}u_{0}^{\dagger}.$ (17) $$

###### Proof.

Let us analyze the product in Equation 13 to show how it relates to the power method. The linear operator representing the reverse process can be written as

$$ $\mathcal{P}_{T}=U_{0}\Pi_{t=0}^{T-1}(U_{t}^{\dagger}U_{t+1})U_{T}^{\dagger}.$ (18) $$

The matrix product $U_{t}^{\dagger}U_{t+1}$ can be analyzed using the extension of Equation 12 to higher ranks. Given 4.1, the expected inner product with the natural data component $u_{i}=u_{i}^{t=0}$ is given by

$$ $\mathbb{E}\langle u_{i}^{t},u_{i}\rangle\approx 1-\frac{\bar{\sigma}^{2}_{t}}{ \lambda_{i}^{2}}\frac{d}{n}.$ (19) $$

Figure: Figure 4: Schematic illustration of the basis perturbation, per index.

The evolution of this product over time is depicted in Figure 4. We are interested in the projection of $u_{t+1}$ onto $u_{t}$, which is the cosine of the angle $\Delta\theta=\theta_{t+1}-\theta_{t}$. This angle is tractable for small noise levels, so we divide our analysis to two parts: $0\leq t\leq\tau$ and $\tau\leq t\leq T$, where the choice of $\tau$ will soon be motivated.

First, we inspect the limit of $t\to 0$ ($0\leq t\leq\tau$). For small angles, we can write

$$ $\displaystyle\Delta\theta=\arccos\Big{(}1-\frac{\bar{\sigma}^{2}_{t+1}}{ \lambda^{2}}\frac{d}{n}\Big{)}-\arccos\Big{(}1-\frac{\bar{\sigma}^{2}_{t}}{ \lambda^{2}}\frac{d}{n}\Big{)}\approx\frac{d}{\lambda^{2}n}(\bar{\sigma}^{2}_{ t+1}-\bar{\sigma}^{2}_{t})=\frac{\sigma_{t+1}^{2}d}{\lambda^{2}n},$ (20) $$

since $\arccos{\theta}\approx\frac{\pi}{2}-\theta$ and $\bar{\sigma}_{t}^{2}=\sum_{\tau=0}^{t}\sigma_{\tau}^{2}$.
The diagonal elements in $U_{t}^{\dagger}U_{t+1}$ are then given by

$$ $\mathbb{E}\langle u_{i}^{t},u^{t+1}_{i}\rangle\approx\cos{\frac{\sigma_{t+1}^{ 2}d}{\lambda_{i}^{2}n}},$ (21) $$

where the off-diagonal elements are negligible, since

$$ $\mathbb{E}\langle u_{i}^{t},u^{t+1}_{j}\rangle\approx\mathbb{E}\langle u_{i}^{ t},u^{t}_{j}\rangle=0,$ (22) $$

which holds for $t\leq\tau_{r-1}$ by Assumption 4.2. Notice, that in small angles, $\langle u_{i}^{t}-u_{i}^{t+1},u\rangle=(\sigma_{t+1}^{2}d)/(\lambda^{2}n)\to 0$, so the vectors $u_{i}^{t}$ are co planar, as depicted in Figure 4.
Thus, the time point basis correlations $U_{t}^{\dagger}U_{t+1}$ form an approximately diagonal matrix with the fraction $c_{i}\triangleq\cos{\frac{\sigma_{t+1}^{2}d}{\lambda_{i}^{2}n}}$ on the diagonal, where $c_{i}>c_{j}$ for $i<j$. We eliminate the dependence of $c_{i}$ on $t$ by choosing the constant schedule $\sigma_{t}=1/T\ \forall t$, to simplify the proof. However, many schedules can be used, as long as $c_{i,t}>c_{j,t}$ remains correct.
Define the partial linear diffusion operator until time $\tau$ by $\mathbb{E}\mathcal{P}_{\tau}=\Pi_{t=0}^{\tau}P_{t}$. Then

$$ $\displaystyle\mathbb{E}\mathcal{P}_{\tau}$ $\displaystyle=U_{0}\begin{pmatrix}c_{0}^{\tau}&&\\ &\ddots&\\ &&c_{r-1}^{\tau}\end{pmatrix}U_{\tau}^{\dagger}=U_{0}c_{0}^{\tau}\begin{ pmatrix}1&&\\ &\big{(}\frac{c_{1}}{c_{0}}\big{)}^{\tau}&\\ &&\ddots\end{pmatrix}U_{\tau}^{\dagger}\underset{\tau\to T}{\to}U_{0}\begin{ pmatrix}c_{0}^{\tau}&&\\ &0\\ &&\ddots\end{pmatrix}U_{\tau}^{\dagger},$ (23) $$

where the diagonal elements decay as $\tau$ grows larger, as $c_{i}>c_{j}$ for $i<j$. Similarly to power iteration, the convergence rate depends on the ratio $c_{1}/c_{0}$. The convergence rate might not be fast enough for the process to converge while the small angles approximation still holds. Thus, we continue with the second phase of our analysis, showing the convergence of the full reverse process.

We now turn to analyse the phase where $\tau\leq t\leq T$. In high noise levels, the correlation with the natural basis is low, and the products $U_{t}^{\dagger}U_{t+1}$ are not exactly diagonal. However, the correlation ”leaks” to a close neighborhood of the original component and the temporal products are still somewhat concentrated around their diagonal. This process happens in accordance with Equation 4.1, where the spreading of the diagonal elements happens for high indices in lower values of $t$ (less noise is needed to spread the correlation). This leads us to Assumption 4.2, claiming that for each index $i$ there exists a time $\tau_{i}$ after which the small angle approximation does not hold; $\tau_{i}>\tau_{j}$ for $i<j$. This is apparent in practice, and depicted in 5 (left image per duo). However, given the decaying diagonal structure of the partial operator $\mathcal{P}_{\tau}$, we will now show that 4.2 is sufficient for the total operator to converge as desired.

Suppose we added one more matrix multiplication to our former analysis, i.e. observe

$$ $\mathbb{E}\mathcal{P}_{\tau}U_{\tau+1}=U_{0}c_{0}^{\tau}\begin{pmatrix}1&&\\ &\big{(}\frac{c_{1}}{c_{0}}\big{)}^{\tau}&\\ &&\ddots\end{pmatrix}U_{\tau}^{\dagger}U_{\tau+1}.$ (24) $$

Assumption 4.2 guarantees $U_{\tau}^{\dagger}U_{\tau+1}$ is diagonal just enough not to spoil the diagonality of the next partial operator $\mathbb{E}\mathcal{P}_{\tau+1}$. Since the elements of the partial product $\mathbb{E}\mathcal{P}_{\tau}$ decay faster with $i$ than any single product $U_{\tau}^{\dagger}U_{\tau+1}$, $\mathbb{E}\mathcal{P}_{\tau+1}$ is also diagonal. Overall, the final product is a diagonal matrix with a spectrum that converges to be concentrated around the first eigenvalue, where we can control the distribution of the generated data by the choice of the diffusion parameters. Figure 5 shows our simulation of the process, supporting both assumption 4.2 and the result stated by this theorem.
∎

Figure: (a)
Refer to caption: extracted/5930746/figures/UU_0.png

Oftentimes, the reverse process includes the injection of noise to the intermediate images (e.g. (Ho et al., )). The overall transformation in this case is given by

$$ $x_{g}=\Sigma_{t=0}^{T}\Pi_{\tau=0}^{t}P_{\tau}\xi_{t}=P_{0}\cdots P_{T}\xi_{T} +\dots+P_{0}\xi_{0}$ (25) $$

for some schedule $\{\xi_{t}\}_{t=0}^{T}$ (for example, $\xi_{T}\sim\mathcal{N}(0,1)$ and $\xi_{t}=\mathcal{N}(0,1/T)$ for $t=0,\dots,T-1$).
In this case, the generated output is a combination of a (purely) noisy image that was repeatedly correlated to converge to $v_{0}$ (as shown above), with generally lower noise levels that are ”lightly” correlated, although to the cleaner projection operators. The generated output can thus be seen as a combination of three conceptual parts, with a different balance of the noise level and the portrayed components.

###### Assumption 4.1 .

###### Assumption 4.2 .

###### Theorem 4.3 (Convergence to Power Iteration) .

###### Proof.

#### The first eigenvector

The first part of the sum in Equation 25 is $P_{0}\cdots P_{T}\xi_{T}$, the estimation of the eigenvector with the largest eigenvalue, as shown above theoretically in Equation 23 and empirically in the right matrix of the bottom right duo in Fig. 5. The ”strongest” noise is repeatedly correlated to be concentrated around the first eigenvector.

#### The entire (clean) spectrum

The last part in Equation 25 is $P_{0}\xi_{0}$, a weak noise level that is spread across all components. This noise is very lightly and not repeatedly correlated, although to a clean version of the natural data basis.

#### In between

The third part consists of all the intermediate products $\Pi_{\tau=0}^{t}P_{\tau}\xi_{t}$. The product operators $\Pi_{\tau=0}^{t}P_{\tau}$ preserve varying parts of the natural spectrum, according to $t$ - as $t$ grows, the total projection tends to retain only the components associated with larger eigenvalues. This can bee seen in Fig. 5. The right matrix in each pair shows the product $\Pi_{\tau=0}^{t}P_{\tau}$ for varying values of $t$. The total projections range from the entire spectrum (top left), to only the leading eigenvalue (bottom right). In between, the products are diagonal matrices where the entries in the indices of the smaller eigenvalues have already diminished, in a similar way to the convergence described in Equation 23.

Thus, we get a combination of a solid estimation of the leading eigenvector, together with a more uniform and week sampling of the components with low eigenvalues in the natural data basis. In between, the intermediate projections are at different levels of convergence to the leading eigenvector, hence tend to be more concentrated on components with large eigenvalues as $t\to T$. The freedom in choice of schedule $\{\xi_{t}\}_{t=0}^{T}$, allows control of the spread of the final distribution on the natural data components.

Figure: (a)
Refer to caption: extracted/5930746/figures/distribution_over_PC_range.png

To inspect this, we can plot the empirical absolute distribution of generated images over the clean PCs, given by

$$ $c_{i}=\frac{1}{n}\sum_{j=1}^{n}\frac{|\langle u_{i},x_{j}\rangle|}{\left\lVert x _{j}\right\rVert_{2}},$ (26) $$

where $v_{i}$ is the clean principal component with index $i$ (defined in Equation 9).
Figure 6 shows the empirical distribution of generated images over the clean principal components. On the left, we plot the distribution without injected noise, i.e. $x_{g}=\mathcal{P}_{T}\xi$ (as in Equation 14), for various values of $T$. As we show above, the distribution tends to be concentrated on the first eigenvector as $T$ increases. The center plot shows the distribution of the process including the injected noise in the intermediate denoising steps. While in the low indices the dominant behavior is similar to the former case, the higher orders do not converge to zero and maintain their presence in the generated distribution. On the right, we picked the best configuration ($T=65$ in this case) to approximate the natural distribution. In addition to the convergence in mean, we included the standard deviation of the natural and generated samples, resulting in a decent fit to the target distribution.

## 5 Empirical Extension to Deep Denoisers

In the linear case described above, the optimal denoiser is given by the PCA projection onto the clean(er) data. These denoisers are computed with the training data, and their principal components do not depend on the input in the reverse process. When the denoiser is nonlinear, and might be implemented using a deep neural network, its input-output mapping can be locally expressed via the network Jacobian, by

$$ $D(x_{t})=\nabla D(x_{t})x_{t}=V_{t}\Lambda_{t}V_{t}^{\dagger}x_{t},$ (27) $$

where $V_{t}\Lambda_{t}V_{t}^{\dagger}$ denotes the eigen decomposition of the Jacobian calculated at $x_{t}$. For simplicity, we assume that the Jacobian is symmetric and non-negative (which is approximately true (Mohan et al., )). Note that in this case, the denoising base depends on the input image (and noise level). While the network is non linear, we can follow the generation path in the sampling process and inspect the basis of the network Jacobians calculated at the intermediate sampled points $x_{t}$. We can then trace $\sin{\theta_{J}}=\sqrt{1-\langle v_{i}^{t},v^{t=0}_{i}\rangle^{2}}$ where $v_{i}^{t}$ is the $i^{th}$ column in $V_{t}$ defined in Equation 27, in a similar way to our simulations of the linear case (Figure 2). This can be calculated per generation path, where $x_{0}$ is the final generated image, and $V_{0}$ is the basis of the Jacobian calculated at this final point.

Figure 7 shows $\sin{\theta_{J}}$ calculated using the Jacobians of a UNet based diffusion model, described in (Ning et al., ). This model was simply chosen as the (^1^11[https://paperswithcode.com/sota/image-generation-on-celeba-64x64](https://paperswithcode.com/sota/image-generation-on-celeba-64x64))state-of-the-art in the task of image generation considering the CelebA dataset at the time of writing this paper. We used the default settings and calculated the Jacobians at the final iterations. We plot the results for the leading 300 Jacobian eigenvectors, where the color is assigned by the index - darker colors for lower indices $i$. We repeated the experiment sampling images from the CelebA dataset (left) and CIFAR 10 (right). Even though the denoising model is far from linear, the decay of the angle between the denoising basis in high noise levels and the natural denoising basis is similar to the decay in the linear case (compare to Figure 2). In this case as well, the correlation of the low indices (and hence low frequencies) withstands higher noise levels, thus appearing first in the generation process. As this is the basis of our analysis comparing the reverse diffusion process to power iteration, this experiment shows that our analysis is relevant in a broader context and not just in the simplified linear case.

Figure: (a)
Refer to caption: extracted/5930746/figures/sin_theta_generation_celeba.png

## 6 Conclusion

In this paper, we discuss a simple diffusion model with a linear denoiser and normalization free sampler, that allows us to cast the diffusion problem as noisy PCA, and make the connection to the spiked covariance model assuming that the natural data distribution reside in a low dimensional subspace. This enables us to show that in the linear case, the generation process acts as a “correlation machine”, where initial random noise is repeatedly correlated to noisy estimations of the natural data basis, to finally embody the true distribution, in a manner similar to the power iteration method. We show that in this process, low frequencies emerge earlier, and more data contributes to a richer representation per the same diffusion configuration. Finally, we demonstrate the relevance of our analysis also in a deep, non-linear diffusion denoiser.

We acknowledge the limitation of admitting a linear model, with its lack of ability to represent the complex data often expected of diffusion models. While our theoretical setting is modest, we empirically demonstrate how our observations deduced from a simple linear model and classic theory (Johnstone, ; Nadler, ) are relevant to more general models and datasets. This enables us to shed light on the internal mechanism powering this technology, and connect it to a rich pool of theory and prevalent methods such as power iteration.

## Acknowledgments

The authors extend their sincere gratitude to Mojtaba Ardakani, Tom Duerig and Shady Abu-Hussein for their valuable feedback. RG and DW were partially supported by the Israeli Innovation Authority, the Israeli council of higher education and the center for AI and Data Science at Tel University. RG and DW are also thankful for KLA and Google for their support.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: (1) what are $x_j$ and $n$ in eq. (26)? 
(2) where $v_j$ right after eq. (26) should be $u_j$.
(3) line 487: what is $J$ in $sin(\theta_J)$? is it $i$?
soundness: 2
strengths: The authors aim to provide an analysis of the diffusion process, an important question for diffusion models. This is achieved using a simple model.

Overall, the writing is easy to follow.
summary: This paper studies the generation process of diffusion models. Focusing on linear models where the diffusion process is similar to performing noisy PCA, the authors show that the generation process of the linear diffusion model is similar to the power iteration method in that it converges to the leading eigenvector of the underlying data distribution.
weaknesses: Providing theoretical studies for complex problems in deep learning and deep generative models is often very challenging. Thus, it is common practice to study a simpler problem that can shed light on the underlying mechanisms of more complex ones. While this work falls into this category, I found that the results for the linear diffusion model diverge significantly from the phenomena observed in real cases.

(1) Theorem 4.3 shows that the diffusion process only converges to the dominant eigenvector (the one corresponding to the largest eigenvalue). Thus, it gives a bad estimation for the true distribution with a dimension larger than 1. While eq. (25) also provides a formulation when noise is injected in the intermediate steps, (i) it does not provide an estimation guarantee for the entire spectrum, (ii) many recent diffusion models such as DDIMDMP-solver do not require the injection of noise in the intermediate steps.

(2) The results in Section 5 (e.g, Figure 7) on real cases only demonstrate the convergence to the entire spectrum, not the dominant eigenvector. This also supports the gap between the analysis in Theorem 4.3 and the phenomenon observed in real cases. 

(3) Initially, I thought Section 5 was intended to support the argument that the diffusion process converges to the "true distribution" as in Theorem 4.3. However, this is not the case, as Figure 7 plots the convergence of the network Jacobian evaluated at $x_t$. Thus, even at $t= 0$, it depends on the learned network and does not represent the true distribution of the data; the convergence probably holds as long as the network is Lipschitz. In other words, the convergence does not indicate the quality of the diffusion process. For example, a network that outputs all zeros would also exhibit fast convergence. I may be misinterpreting the results; if so, I'd appreciate clarification.

---

confidence: 5
contribution: 1
initial_rating: 3
presentation: 2
questions: I'd like to summarize my final assessment of this paper here:

1. First of all, the authors argue a linear diffusion models generate samples that concentrate on the first eigenvector of the covariance matrix. However, the diffusion models considered are not trained with the standard denoising score matching loss, rather, it is trained with a sequential denoising loss that is very different from the standard diffusion model. Besides the training loss, the reverse diffusion process considered in this paper is also not the standard reverse diffusion process. Due to this differences between the model considered in the paper and the standard diffusion model (which is not clearly explained in the paper), the conclusion in this paper can be misleading to the general public. The theorem states the linear diffusion models can only generate images that concentrate on the first eigenvector, thus of limiting diversity. However, as I've shown in my comment, the optimal linear diffusion models trained with standard diffusion denoising score matching loss are equivalent to the score function of noisy Gaussian distributions. Therefore, standard reverse sampling using the optimal linear diffusion models is equivalent to sampling from a Gaussian distribution, which can generate samples with great diversity. I urge the author to carefully take care of this discrepancy between standard diffusion and the model they considered, to avoid any misleading. 

2. Second, in the section 5 of the paper the authors investigate the Jacobian and state the linear diffusion model is similar to the nonlinear diffusion model. However, while there indeed is certain similarity between case and nonlinear case, a significant amount of dissimilarity also exists but is ignored by the authors without any discussion. The implication of this section is also unclear. I expect the authors to discuss the property of the final generated samples from the nonlinear diffusion models with some mathematical formulation. Please notice that for linear model, the Jacobian fully explains the network behavior since Jacobian is the linear weight itself. However, for nonlinear models, the Jacobian $\nabla f(x_t)$only captures how the network's output change w.r.t. to input $x_t$ at a local region. However, to investigate the property of the final samples, we need to study the property of $f(x_t)$ itself. 

3. Third, some assumptions made in this paper is not valid. For example, the author cannot assume $\sigma_t=\frac{1}{T}$ (line 354), which implies $\epsilon_t\sim N(0,\frac{1}{T})$. Under this assumption, more sampling steps result in smaller amount of noise at the end stage of diffusion forward process. However, the noise level at the end of the forward diffusion should be sufficiently large so that the noisy distribution at time step is approximately a Gaussian distribution with zero mean. If this is violated, the reverse sampling process starting from a random Gaussian noise can not recover the clean distribution. Though the authors in the comment argue they use valid schedule in the numerical experiment, I do believe they should derive their theorem from valid schedule as well. Unfortunately since the theorem derived in the paper explicitly depends on $T\rightarrow \infty$, the $\epsilon_t$ will be extremely small, directly contradict the basic assumption of diffusion model, therefore the theorem is less meaningful.

In conclusion, I suggest rejection for the paper and I believe major revision should be made before the next submission.
soundness: 1
strengths: The paper has a clear structure and is easy to read. The topic is also interesting since studying the evolution of the denoisers along diffusion trajectory is important but less well explored in the literature.
summary: This paper explore diffusion reverse sampling process by analyzing a linear diffusion model. The main result is that the reverse sampling process leads to samples that concentrate on the first eigenvector of the covariance matrix of the underlying data.
weaknesses: I have doubt on both the theoretical and practical aspects of this work. The theoretical results do not seem correct (correct me if I am wrong) and the connection to practice is weak. I list my questions as below:

*Major questions:

(i) Can the author explain how do they conclude Theorem 4.3 ? The author claimed that the final projection operator is a diagonal matrix with a spectrum concentrate around the first eigenvalue (line 409-412) and they supported this claim qualitatively in Figure 5. However, I cannot find a close-form expression for this projection matrix as in the limit of $\tau\rightarrow \infty$. From (24) I would assume 
\begin{align}
\mathbb{E}\mathcal{P_{\tau}}=U_0 c_0^\tau\text{diag}[1, (\frac{c_1}{c_0})^\tau,...] U_0^T,
\end{align}
, which do not really hold since equation (20) only holds for small t. If we assume this expression holds, we have 
\begin{align}
\mathbb{E}x_gx_g^T=(U_0c_0^{2\tau}\text{diag}[1, (\frac{c_1}{c_0})^{\tau},...])(\mathbb{E}[\epsilon_T\epsilon_T^T])(\text{diag}[1, (\frac{c_1}{c_0})^{\tau},...]U_0^T).
\end{align}

By the assumption that diffusion stepsize is $\frac{1}{T}$ and there are $T$ steps in total (line 360-361), we have $\epsilon_T\sim\mathcal{N}(0,\frac{1}{T})$. Therefore $\mathbb{E}x_gx_g^T=U_0\text{diag}[c_0^{2\tau}/T, c_1^{2\tau}/T,...]U_0^T$. Then how does this term equal to $u_0u_0^T$ as in equation (17)? Shouldn't all eigenvectors go to zero as $T\rightarrow \infty$ ? 


(ii) In section 5 the author claims that correlation of the low indices (lower frequencies) withstands higher noise levels. However, it is shown in Figure 7 that many of the dark curves remain flat during a wide range of time steps and only drop when t approaches 0, which is very different from linear case (Figure 3) that show a gradual evolution of the sine angles. Therefore I disagree with that the decaying behavior of the nonlinear model  is similar to that in the linear case.

*Minor questions:
(i) The Figure 7 is very difficult to read since all the useful information is concentrated near $t=0$. I suggest the authors to redraw the plot using log scale.

(ii) The paper could benefit from putting more details on the theoretical results in an appendix. Not sure whether authors are allowed to add one.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - Most experiments for linear case are done on MNIST- could the authors clarify to what extent the assumptions made in the spiked covariance model apply to real-world datasets? For example, in cases where the low-dimensional subspace assumption does not hold.
- Moreover, since power iteration is such a major theme of the paper, a small discussion of the power iteration method in the introduction would help enhance readability.
- The term $u_i^t$, has never been defined. Moreover, a consistent notation of keeping time in the superscript and index in the subscript  (or vice-versa) for both vectors and matrices will help avoid confusion.
soundness: 3
strengths: The paper provides a novel perspective on diffusion models for understanding the generation process by drawing parallels with power iteration. The work opens up new research directions for improving the interpretability and performance of generative models by highlighting the eigenvector alignment properties in diffusion models. The paper is well-written with most claims justified theoretically and backed by empirical evidence.
summary: The explores the relationship between the power iteration method and diffusion models with a linear denoiser, that can be modeled as a PCA projection, enabling analysis via the spiked covariance model. The authors demonstrate that, in the reverse process, high eigen-value components emerge earlier, and leading eigenvectors align with the data distribution.
weaknesses: The paper does a good job at conveying the idea and intuition from the simulations, however, the proof of the main result (Theorem 4.3), could be refined more. Particularly, I can not fully understand the proof in the high noise regime i.e., $\tau \leq t \leq T$. A more explicit explanation for the following statement would enhance the rigor: how exactly does Assumption 4.2 guarantee that $U_\tau^\top U_{\tau+1}$ is diagonal just enough not to spoil the diagonality of the next partial operator $\mathbb E \mathcal P_{τ+1}$? Can the authors justify this mathematically?

Meta-review: This paper presents a theoretical study on diffusion models. The main result is the establishment of a connection between power iteration and linear diffusion. It is claimed that this result can provide insight to the theoretical understanding of the generation process of diffusion models. The major weakness pointed out is the linear diffusion setting considered in this work is too restricted and the techniques developed can not be applied to more practical settings. There are also some concerns on the proofs that are not fully addressed.

Decision: Reject