Likelihood Training of Cascaded Diffusion Models via Hierarchical Volume-preserving Maps

We greatly appreciate that the reviewer's close reading of our work, and their thoughtful concerns. Below, we address their reservations on the derivation of Theorem 5.1 and the necessity of the volume preserving map.

Eq. 27 -> Eq. 28. The shift from $\mathbf{z}$ to $\mathbf{x}$ is correct and not a typographical error — this is due to Theorem 2 in Shirdhonkar. We provide a more in-depth treatment of this theorem for the reviewer, see below for details.

On the connection to the EMD metric.

More detailed treatment of Theorem 2 in Shirdhonkar et. al. We have provided two additional statements in Appendix A. Namely, Lemma A.1, which restates Theorem 2 in Shirdhonkar et. al according to our notation, and Lemma A.2, which shows in detail how Lemma A.1 can be applied to Theorem 5.1. To answer the reviewer’s immediate question, observe that for all possible values of $j, s, n$, the term $2^{-j(s + n/2)}$ is upper bounded by 1. Therefore it can be safely ignored in the upper bound, which is all we require for our Theorem 5.1 (Shirdhonkar et. al, 2008 provides both upper and lower bounds on the EMD).

Significance of Theorem 5.1 --- it does not seem like this allows diffusion models to be trained with the EMD measure. Actually, this theorem shows precisely that diffusion models can be trained with EMD. Namely, by training with the loss function Eq. 17, our model minimizes the EMD between the true and predicted scores of the data distribution $\mathcal{X}$, which can be understood as unnormalized histograms in Eq. 19. This is the surprising result of Theorem 5.1 – a model that minimizes Eq. 17 is simultaneously a hierarchical model on $\mathbf{z}^{(s)}$ and an EMD-based model on the original $\mathbf{x}$. We believe that this connection is significant, since the EMD measure is known to have desirable properties for inferring distances between spatially structured data.

On the ELBO and the necessity of volume-preserving maps. Indeed, when the goal is to compute a lower bound on the likelihood (i.e., the ELBO), the reviewer is correct, and volume-preserving maps are not strictly necessary. However, this derivation is different (and we argue, inferior) to our framework (Eq. 15) in two ways:

First, the resulting ELBO is looser than our ELBO (Eq. 15), since it does not utilize the volume preservation property to obtain the equality in Eq. 9. Empirically, this difference is significant: in Appendix C and Tables 4-6 we compare our model to a hierarchical model that does not utilize a volume-preserving $h$ (i.e., a standard cascading hierarchy) derived without Eq. 9. With all other hyperparameters (architecture, noise schedule, and code implementation) held constant, our model improves significantly over the standard hierarchical baseline.

Second, its inability to produce exact likelihoods $p_\theta(x)$ (e.g. Eq. 6) greatly limits its applicability. For example, the derivation does not work with autoregressive models or normalizing flows, while our framework does. Furthermore, a key property of diffusion models is that they can be interpreted both as an ODE (Eq. 6) and a latent variable model (Eqs. 3-5). As such, this derivation would also be useful for diffusion models that use the exact likelihood formula under Eq. 6 (such as Song et. al, 2020b, Song et. al, 2021, Lipman et. al, 2022, or Zheng et. al, 2023). While we derive our hierarchical diffusion model in its latent variable form, we note that all reported likelihoods (i.e., Tables 1-3) can also be computed exactly via Eq. 6.

We greatly appreciate the time and care the reviewer took to study our work, the in-depth discussion of our theoretical contributions, and the encouraging comments. Below are detailed responses to their concerns.

The framework proposed considers all homeomorphisms that satisfy Eq. 8, but only two linear examples of $h$ are empirically evaluated. We believe a major contribution of our work is demonstrating general conditions under which hierarchical models also yield tractable and exact likelihood models, which can greatly inform future research on diffusion models, many of which are hierarchical. We leave further exploration of various $h$, including non-linear cases, to further work. We do recognize that the linear scope of our experiments should be clearly stated early on, and we have added sentences clarifying this point to the abstract, introduction, and Section 4 (introducing the hierarchical volume-preserving map). To quote our added statement in the list of contributions in Section 1: “We show that two well-known multiscale decompositions, the Laplacian pyramid and wavelet transforms, are hierarchical volume preserving maps, where the linearity of the maps reduce the volume-preserving property into more familiar notions of orthogonality and tight frames (Unser et. al, 2011).”

Connection to the EMD Metric.

Transport between scores, rather than images. Indeed, Theorem 5.1 relates our training loss (Eq. 17) to the transport between scores of the images, rather than images themselves. However, this is not so out of the ordinary. The notion of applying Earth Mover’s Distance to features of the image histogram, rather than its pixels directly, is a relatively common strategy in the literature — for example, the seminal work Rubner et. al, 2000 (as cited in our text) finds that the EMD distance on image “signatures” such as texture, color, and joint color/position embeddings work better than the raw pixels themselves.

The Wasserstein-p distance is defined between probability distributions, and the score is not a probability distribution. In Eq. 18, we do in fact consider the setting where $x$ and $y$ are unnormalized and do not sum to one (see the sentence leading up to Eq. 18). This unnormalized case is known as partial matching, and is a setting considered in Rubner et. al, 2000 as well as Shirdhonkar et. al, 2008 which our theory builds off of (see respective references in our work). On the other hand, non-negative entries in the score can indeed be problematic, but this can be circumvented by subtracting the most negative entry from both $x$ and $y$. Namely, let $c$ be such that $c \leq x_{ij}$ and $c \leq y_{ij}$ for all $i, j$. Then $x - c$ and $y - c$ are surely non-negative. Moreover, inspecting Eq. 18, we can see that this transformation does not affect the unnormalized Wasserstein metric, since the $c$ from each unnormalized histogram cancels.

Free energy term. We do not quite follow the reviewer’s argument here. What is the reviewer referring to as the “additional free energy term” included in the $z^{(1)}$ component of the wavelet score? As it stands, the equation provided by the reviewer appears to agree entirely with our original formulation (Eq. 15).

Preserving likelihoods when $z$ and $x$ are not the same dimension. We agree that likelihood computation is much trickier when $z$ and $x$ are not the same dimension. However, we can preserve a pushforward measure when we assume $z \in \mathcal{M}$ and $x \in \mathcal{N}$ are from differentiable manifolds $\mathcal{M}$ and $\mathcal{N}$, and $h$ is a diffeomorphism between the two spaces (we recommend referring to [1, 2]). This is satisfied because $h$ is linear and invertible in the cases we consider. [1] Foundations of differential geometry. S Kobayashi, K Nomizu. 1996. [2] Normalizing flows on riemannian manifolds. M Gemici, D Rezende, S Mohamed. 2015.

Acknowledging the trade-off of a likelihood weighting scheme optimizing FID. Optimizing for FID, we obtain BPD 9.55 and 10.31 respectively for W-PCDM and LP-PCDM, respectively. We have added an acknowledgement of this trade-off.

Why not use $M=1$ for OOD detection? Thank you for this suggestion! We tried using $M=1$, and saw little appreciable change in performance (AUROC of out-of-detection classification of SVHN images was reduced by a hundredth of a point), so we have included this result in our work instead of the $M=2$ case, which we had chosen following the experiments in Nalisnick et. al, 2019. Anomaly detection with $M=1$ is clearly more difficult than $M=2$, so this is a better result. As the reviewer suggests, the superior performance even at $M=1$ may be closely related to the concentration of measure in high dimensions.

Minor suggestions. We thank the reviewer for bringing these errors to our attention, and have incorporated them into the text.

We appreciate the reviewer's thoughtful comments and insightful critiques. We also thank the reviewer for their extensive questions. We acknowledge that our work builds extensively on existing principles in computer vision, statistics, and signal processing, and hope that our point-by-point response has provided some intuition for the principles underlying our work.

Informal and imprecise theoretical presentation. We do take very seriously criticisms of the informality and incorrectness of our text. If the reviewer’s concerns with the theoretical presentation are related to the other points (and minor comments) further in the review, we hope we have addressed them in a satisfactory manner and the reviewer may reconsider this judgment. We are happy to answer any further questions.

Contribution of the paper. The paper uses Laplacian pyramid / Wavelet transforms as is. To summarize our contributions, we have proposed a generalized framework for hierarchical modeling that allows for tractable likelihood computations, which was previously not considered or possible with related works. While two notable maps are the Laplacian pyramid / Wavelet transform, our framework is certainly not limited to these transforms. We are not aware of Laplacian pyramid / wavelet diffusion models that derive tractable likelihoods or even conduct density estimation experiments. Finally, we demonstrate significant improvements to the state-of-the-art with our hierarchical framework, and enable future work with more complex hierarchical volume-preserving maps.

There are few / no comparisons to 2022 and 2023 methods. We did not include more recent results as we did not find significant improvement over existing results in Table 1. As per the reviewer’s request, we have incorporated several more 2022 and 2023 methods (Soft Truncation Kim et. al 2022, INDM Kim et. al 2022, i-DODE Zheng et. al 2023). Crucially, we still obtain significant improvements over these methods. We would be happy to include any other results that the reviewer suggests.

There are no comparisons to diffusion models with super-resolution, wavelet, or laplace transforms. As discussed in our work, there did not previously exist super-resolution, wavelet, or Laplacian pyramid -based diffusion models that preserve likelihood computation. Therefore, we cannot form any comparisons on the experiments in Section 6. That being said, as pointed out by another reviewer, it is possible to fold the latent scales of standard cascading hierarchies (that do not satisfy Definition 1) into a looser variational lower bound. However, this comes at a significant cost to modeling performance, as seen in Table 4-6 (Appendix C).

The experiments are roundabout, and don’t directly substantiate specific contributions made in this paper. Can the reviewer elaborate on this point? If it is due to points in the minor comments below, we hope we have adequately addressed this criticism. The two central claims in this paper are that 1) we enable exact likelihood estimation, which is substantiated by Lemma 4.1 and 2) that hierarchical modeling with wavelets and Laplacian pyramids approximate an EMD loss on the scores, which improves modeling performance and is substantiated by significant improvements in performance (Tables 1-3).

There are no ablations. Our paper claims that hierarchical volume-preserving maps allow for likelihood modeling in cascaded models. Unfortunately, hierarchical maps that are not volume-preserving (such as a standard cascading model) do not allow for exact likelihood computation, therefore direct ablations are not possible. However, as previously mentioned, one can still obtain a loose loser bound on the likelihood. We compare our work against this lower bound in a series of ablations in Appendix C and Tables 4-6, and find that using a standard cascading hierarchy (end of Section 4.1) significantly reduces modeling performance.

Minor comments:

Eq. 3 differs from Ho et. al, 2020. The RHS here is different from Ho et al, 2020 because they are different Markov processes. Eq. 2 in Ho et. al, 2020 describes the forward diffusion process, whereas Eq. 3 here describes the reverse diffusion process. The $q$ in Eq. 5 is also taken over the same reverse diffusion. For more information, see Appendix E.1 of Kingma et. al, 2021 (as cited in our work).

Premise of the paper. We do not understand the reviewer’s concern here. Why must we assume independence of the scales $z^{(1)}, z^{(1)}, \dots, z^{(S)}$ for the invariance to hold? In fact, we consider solely the dependent case, as this is the necessary conditional structure for a hierarchical multiscale model.

Square root in Eq. 8. Indeed, the square root is not strictly necessary here. However, this entire term (including the square root) is mathematically meaningful in normalizing flows and differential geometry as the instantaneous change in volume between differentials (infinitesimal volumes) of the two Riemannian manifolds at $x$ and $h(x)$ (Berger and Gostiaux, 2012), so we have opted to keep the square root.

What does $z^{(1)}, \dots, z^{(S)} = h(x)$ mean? Since $h$ defines a hierarchical map, the range of $h$ is the cartesian product of the $S$ scale spaces, i.e., $h: \mathcal{X} \rightarrow \mathcal{Z}^{(1)} \times \mathcal{Z}^{(2)} \times \dots \times \mathcal{Z}^{(S)}$, where $x \in \mathcal{X}$ and $z^{(s)} \in \mathcal{Z}^{(s)}$. We have swapped sides of the equality and added parentheses around the RHS in the equation $h(\mathbf{x}) = (z^{(1)}, \dots, z^{(S)})$ to emphasize that the functional output of $h$ is a cartesian product, and further clarified the functional definition of $h$ in Section 5.1.

How is $p_\theta$ “trained” in Lemma 4.1? We have removed the word “trained” in this sentence. $p_\theta$ need not be trained in any manner — the theorem holds as long as $p_\theta$ is a valid likelihood function on $\mathbf{z}^{(1)}, \dots, \mathbf{z}^{(S)}$.

What is $p_\theta(x)$? What is the model likelihood? The model likelihood is the joint likelihood of the model and data, rigorously defined as $p_\theta(x)$. $p_\theta$ simply refers to a likelihood function that takes a data point (or its latent representation) as input, and returns its modeled density. We have changed "model likelihood" to simply "likelihood" to reduce confusion.

The statements after Eq. 9 are not true. We state after Eq. 9 that Eq. 9 (i.e., $\log p_\theta(x) = \log p_\theta(h(x))$) demonstrates a form of probabilistic invariance. Formally, letting $g_x(h) = p_\theta(h(x))$, the set of hierarchical volume-preserving maps $\mathcal{H}$ form an equivalence class under the the equivalence relation $g_x(h) = g_x(h’)$ for $h, h’ \in \mathcal{H}$ and $x \in \mathcal{X}$. This equivalence relation can be clearly seen by applying Eq. 9 to $h$ and $h’$, and invoking the transitive property.

$h$ in Eq. 10 is a linear operator. Isn't it a function? Indeed, in Eq. 10, $h$ is a function that is happens to also a linear operator, since we are considering the standard cascading hierarchy as described in Ho et. al, 2020, which can be seen as a bilinear resampling operation. In general, $h$ can be an arbitrary function that satisfies Eq. 9.

The * in (j*2). The asterisk is simply the product, not a convolution. We see how it can be confusing, and have replaced the asterisk with a dot (i.e., $j \cdot 2$) for clarity.

What does Figure 2 show? In this figure, we show the Laplacian Pyramid coefficients, i.e., $z^{(1)}, z^{(2)}, …, z^{(4)}$ from left to right.

What is a tight frame? Informally, a tight frame describes a representation (of a signal) that is norm-preserving. For a more formal treatment, please see Unser et. al, 2011 (as cited in our manuscript).

What is $h$ in the Laplacian Pyramid section in 4.2? What is $h$ in the Wavelet Decomposition section in 4.2? $h$ refers to the Laplacian pyramid or wavelet transform, respectively. Formally, $h(\mathbf{x}) = (z^{(1)}, z^{(2)}, …, z^{(S)})$ where where $(z^{(1)}, z^{(2)}, …, z^{(S)})$ are as defined in the respective subsection (Wavelets or Laplacian Pyramids) in Section 4.2. While the functional form is not explicit to maintain generality w.r.t. $S$, this definition is rigorous.

Wasn’t the goal of this paper to reduce resolution? Resolution vs scale? The goal of the paper is not to reduce resolution. Hierarchical models are generally useful because they enable an explicit multi-scale model of the data. While super-resolution is used to obtain finer scales from coarser scales (e.g. $\mathbf{z}^{(s)}$ from $\mathbf{z}^{(s - 1)}$), all scales are retained in our model, and therefore information is not lost. In fact, in Ho et. al, 2022 (in the paper) the sequence of scales actually contains more parameters (i.e., pixels) than the original image. Moreover, when $h$ is the wavelet transform, the size of the largest scale is 25% the original image.

Relation between Lemma 4.1 and Eq. 15. The reviewer’s derivation is not true, as $x \neq z^{(1)}…z^{(S)}$. $x$ is only equivalent to its scales under the transformation $h$, i.e., $h(x) = z^{(1)}…z^{(S)}$. In general, transformations (like $h$) do not preserve likelihoods. In other words, $p_\theta(\mathbf{x}) \neq p_\theta(h(\mathbf{x}))$ for an arbitrary function $h$. Therefore, to our knowledge Eq. 15 is not true unless Lemma 4.1 holds and $h$ is a volume-preserving map.

What is the q in eq 16? The $q$ in Eq. 16 is the same $q$ as in Eqs. 3 and 5.

What are the $z_{0:T}^{(s)}$ in Eq. 16? $z_{0:T}^{(s)}$ are the variables of the diffusion process corresponding to each scale $s$, where $z_0^{(s)}$ = $z^{(s)}$. This is much like $x_{0:T}$ are variables of standard diffusion process as described in Section 3 where $x_0 = x$, except that the diffusion described by $z_{0:T}^{(s)}$ are additionally conditioned on the prior scales $z^{(<s)}_0 = z^{(<s)} = (z^{(1)},z^{(2)}, \dots, z^{(s-1)})$.

Understanding Eq. 18. This equation defines a distance between two distributions $\mathbf{x}$ and $\mathbf{y}$ based on the cost to transport mass from each pair of locations ($(i,j)$ in $\mathbf{x}$ and $(k, \ell)$ in $\mathbf{y}$) so that $\mathbf{x}$ looks like $\mathbf{y}$. The infimum ensures that this cost is the most reasonable (i.e., efficient) transport plan. Intuitively, we are trying to measure the effort required to shift pixel intensities from the “high mass” areas of $\mathbf{x}$ towards the “high mass” areas of $\mathbf{y}$, until $\mathbf{x}$ resembles $\mathbf{y}$. For a rigorous understanding, we encourage the reader to look into more detailed treatments on optimal transport [1] or Earth Mover’s Distance (Rubner et. al, 2000).

[1] Computational Optimal Transport. https://arxiv.org/abs/1803.00567

Understanding Eq. 19. We agree with the reviewer that Theorem 5.1 is quite surprising! The main theoretical machinery for this result can be found in Theorem 2 of Shirdhonkar and Jacobs, 2008 (as cited in our manuscript). In essence, the multiscale decomposition offered by wavelet and laplacian pyramid transformations allow us to approximate Eq. 18 via differences between the corresponding scales of $\mathbf{x}$ and $\mathbf{y}$. This is significant for two reasons. First, Eq. 19 is much easier to compute than Eq. 18. Second, Eq. 19 is also known as the Earth Mover’s Distance, which is a well-known measure for distances between spatially structured data (e.g. images) that works better than the L2 distance. Therefore, Theorem 5.1 suggests that the performance of our algorithm could be directly attributed to the hierarchical prior.

What is $h$ in Section 6 (experiments)? We discuss $h$ at length in Section 4. In general, $h$ is a hierarchical volume-preserving map — any function that maps data (e.g., images) $\mathbf{x}$ to a sequence of latent variables $\mathbf{z}^{(1)}, \dots, \mathbf{z}^{(S)}$ and satisfies Eq. 8. In our experiments, $h$ refers to either the Laplacian pyramid transform or the wavelet transform (LP-PCDM and W-PCDM, respectively).

Citing simple diffusion. If the reviewer is referring to “simple diffusion: End-to-end diffusion for high resolution images” by Heek et. al, 2023, it appears that this method specifically does not involve any multiscale (i.e., hierarchical, cascading) architecture. (See e.g. Section 1 and Figure 2 caption of their paper.) Our understanding is that they tackle diffusion purely on the pixel space.

Why are there no comparisons to other multiscale diffusion models (Section 6.1, 6.2, 6.3)? There are no other multiscale diffusion models that produce tractable likelihoods — in fact, this is the main contribution of our work. Tractable likelihoods are required for any experiments in Sections 6.1, 6.2, and 6.3. The only known method to obtain a lower bound on the likelihood is via a variational lower bound. We derive and compare against this setting in Appendix C and Tables 4-6, and show that our method performs significantly better. On the other hand, we have added more results for standard diffusion models 2022 and 2023. We note that we still demonstrate significant improvements over these results. We are happy to include any other recent works, or likelihood-based multiscale diffusion models the reviewer suggests.

Comparing OOD benchmark to modern diffusion models. We have included a column in Table 2 against the state-of-the-art diffusion model in likelihood modeling, the variational diffusion model (VDM, Ho et. al, 2021).

What do the rows in Table 2 mean? Why is there no comparison against modern diffusion models? We include detailed information about the experiment in Table 2 in Appendix B.1. We compare against the current state-of-the-art diffusion model, variational diffusion models (VDM, Ho et. al, 2021).

Thank you for your response, we are more than happy to clarify.

On math.

$p(x)$ has two definitions. First, neither Eq. 7 nor Eq. 9 are definitions. Eq. 7 is obtained by marginalization, which is a fundamental property of expectations (https://en.wikipedia.org/wiki/Marginal_distribution). Eq. 9 is obtained via Lemma 4.1, which builds from the continuous change of variables formula (see, e.g., a simpler formulation in [1] or the full manifold-based formulation in Ben-Israel, 1999 in our text), once again derived from fundamental principles in probability theory. Intuitively, the key difference here is the assumptions on $h$. The marginalization formula in Eq. 7 is the "general" case, where $x$ and $z^{(1)}, \dots, z^{(S)}$ are not related by any structured function $h$. On the other hand, Eq. 9 describes our case, where $h$ satisfies Eq. 8. We only obtain Eq. 9 under special assumptions (think: circles have more special properties than squares).

$p_\theta(x)$ admits multiple arguments. Indeed, this is the case. Here, $p_\theta(a)$ just means "the likelihood of the random variable $a$ given the model parameterized by $\theta$." This overloading of $p_\theta$ is common practice. For example, in DDPM (Ho et. al, 2020), $p_\theta$ takes just $x$ (the modeled likelihood), but is also overloaded to take the joint distribution $x_0, \dots, x_T$ (see the first sentence of Section 2) and a single conditional variable $x_{t-1} | x_t$ (Eq. 1). You can see similar overloading for $q$ (throughout Section 2). Ultimately, they are different densities over different variables related by $h$ that happen to coincide when $h$ is a volume-preserving map, which is why Eq. 9 is so useful.

Understanding Eq. 15. Note that we are considering the log probabilities in Eq. 15. The sum of log probabilities is equal to the product of the probabilities. Therefore we equivalently have $p(x) = p(z^{(1)}) p(z^{(2)}|x^{(1)}) \cdots p(z^{(S)} | z^{(1)}, \dots, z^{(S-1)})$. The RHS of this equality becomes the joint likelihood $p(z^{(1)}, \dots, z^{(S)})$ by the probabilistic chain rule. Therefore this is just an application of Eq. 9. Many likelihood models consider log probabilities rather than probabilities for numerical stability.

Volume preservation. We note that the volume preservation property only holds when considering the entire map $h(x) = (z^{(1)}, \dots, z^{(S)})$. It is not clear what precisely this means for just the relationship between $z^{(s)}$ and $z^{(s - 1)}$. Intuitively, converting from $z^s$ and $z^{s-1}$ can be seen as a super-resolution step. Therefore, I do not see evidence that $p(z^{(s)} | z^{(<s)})$ should be equal across $s$, nor that an invertible transformation of the form $h(z^{(s-1)}) = z^{(s)}$ should exist. Super-resolution is an ill-posed task (which we really can only achieve by relying on a very strong prior -- in our case, the prior of natural images), so rigorous relations between different scales may be difficult. However, we are happy that the reviewer has put thought in this direction, and would be happy to discuss further.

Concrete derivation for $p(x)$. Note again that the overloading of $p$ is quite common (see above answer). We do provide a derivation for Eq. 9 in the proof for Lemma 4.1 (Appendix A). Can the reviewer elaborate on why this derivation is not sufficient?

On claims. Our paper does not claim to have an improved loss over other models. Our claim is precisely that we enable a tractable likelihood in hierarchical models, whereas other models do not. This provides many downstream applications related to likelihood modeling, including density estimation, out-of-distribution detection, and neural compression, which we evaluate our models on.

Does this make sense? Please let us know if you have any other concerns.

[1] NICE: Nonlinear Independent Components Estimation. https://arxiv.org/abs/1410.8516

We are grateful for the reviewer's appreciation of the clarity and relevance of our work, and for raising insightful concerns. Below we respond to each comment.

Differences from existing hierarchical diffusion models. Indeed, hierarchical diffusion models exist in the literature, and we thank the reviewer for the references, which we have added to the text. However, as discussed, these models do not support tractable likelihood computation. The main contribution of our work is to show the conditions under which hierarchical models such as these are also feasible likelihood models (Eq. 8), which enable many downstream analyses (anomaly detection, data compression, data enrichment, etc.). To further support the importance of our hierarchical volume-preserving property and the exact likelihood computation it enables (Eqs. 8-9), we provide additional experiments in Appendix C and Tables 4-6 that compare our model against $h$ that does not satisfy Eq. 8.

How is Theorem 5.1 different from Kwon et al., 2022? While the theorems appear similar, Kwon et. al, 2022 differs significantly from our approach. Namely, Kwon et. al, 2022 (as well as several other works - De Bortoli et. al, 2021, Chen et. al, 2021, Lipman et. al, 2022, Shi et. al, 2022, all cited in our work) consider optimal transport in the data space $\mathbf{R}^d$, where $d$ is the image size, e.g., 32x32x3 = 3072 for CIFAR10, whereas we consider optimal transport on the 2D (i.e., discretized $\mathbf{R}^2$) spatial grid of image and score pixels. In particular, Theorem 5.1 relates multiscale score matching to score matching under the Earth Mover’s Distance ($p=1$), which on spatially structured data (e.g. images, audio) has many well-known properties and is regarded as a superior metric to the L2 distance used in original diffusion models (see Rubner et al., 2000). Ultimately, this results in distinctly different theoretical properties from the above-mentioned related work (including Kwon et. al, 2022), and is also the first such connection to EMD on images.

Architecture description in Appendix B. Indeed, there is a separate score model for each scale. We have added a sentence clarifying this point in Appendix B. We have also added the number of parameters for each model (also Appendix B).

Minor edits:

Closing parentheses. Thank you for bringing this to our attention. Indeed, the \left( and \right) latex operations are very helpful, and we have incorporated them here and elsewhere in the text.

$h(x)$ notation. We apologize for the confusion. The range of $h$ is the cartesian product of the $S$ scale spaces, i.e., $h: \mathcal{X} \rightarrow \mathcal{Z}^{(1)} \times \mathcal{Z}^{(2)} \times \dots \times \mathcal{Z}^{(S)}$. We have reversed the equality and added parentheses around the RHS in the equation $h(x) = (z^{(1)}, \dots, z^{(S)})$ to alleviate this confusion, and further clarified the functional definition of $h$ in Sections 3.1 and 5.1.

We thank the reviewer for their continued and thorough efforts in evaluating our work, and for their additional suggestions and questions. Below we provide our responses:

Further concerns.

1. Scaling of the diffusion process. Indeed, any scaling of the coefficients results in a scaling of the respective diffusion processes. However, we maintain that this rescaling would not change the model performance in terms of FID or BPD. In the former, the FID is simply a function of embeddings of the sampled images, so it is agnostic to the sampling process (e.g., agnostic to whether images were sampled autoregressively or all at once, or via scaled or unscaled diffusion processes). In the latter, any alternative scaling of the wavelet transform is no longer volume-preserving, and incurs a penalty in the form of the $\sqrt{\text{det}\left(\left[\frac{\partial}{\partial x} h(x)\right]^T\left[\frac{\partial}{\partial x} h(x)\right]\right)}$ term in the change-of-variables formula for $p_\theta(x)$ (Eq. 20). In terms of the negative log-probability (proportional to BPD), this results in a correction factor that would actually cancel out the change in log probability due to scaling coefficients. We do not consider these non-volume-preserving cases because this correction term can be difficult to compute or derive.

   However, we do understand the reviewer's concern. We interpret the underlying question as whether the wavelet transform (and Laplacian pyramid) results in a "real" improvement in model performance over vanilla diffusion models, or just a "rescaling" of the loss in some trivial manner. We firmly believe that it is the former, and elucidate this in the following experiment. Taking the learned conditional score networks $s_\theta(z^{(s)}|z^{(<s)}, t)$ in a trained W-PCDM model on CIFAR10, we observe that we can reconstruct a conditional "vanilla" score $s_\theta(x, t)$ by applying $h^{-1}$ to the hierarchical conditional scores -- i.e., the inverse of Eqs. 31-32. This holds due to the simple form of the score (i.e., $\frac{x + \epsilon}{\sigma}$, Vincent, 2011) when q(x_k|x_0) is Gaussian (which is our case, see e.g., Ho et. al, 2020). After this transformation, we evaluate the ELBO of this "vanilla" diffusion model on the CIFAR10 test set. We ultimately obtain the same gain in performance as the W-PCDM (i.e. $2.65 \rightarrow 2.35 BPD$), even though the present formulation no longer transforms the pixel scores in any way after they are computed (only to convert from the wavelet scores to the pixel scores). This tells us that the gain modeling performance is unrelated to the scaling of the wavelet coefficients.

2. Architecture of the score model. Please note that we do discuss in-depth the architecture of the score model in Appendix B -- in summary, we use the VDM architecture with no changes except for the input image, which now concatenates the noise-perturbed scale variable $\tilde{z}^{(s)}$ with the conditional information $z^{(<s)}$. We have also added an explicit description of the input-output relationship of our model in Appendix B.

   FLOPs. While our models have more total parameters, our W-PCDM implementation actually has significantly fewer FLOPs than the standard diffusion model, and is therefore both faster and smaller in memory footprint. This is mainly due to the smaller input sizes at each scale. Overall, we approximate that the sum total FLOPs across all scales (versus a vanilla VDM model) is reduced by at least a factor of 2. (For more details, see "Calculating FLOPs of W-PCDM.") As discussed in [1], FLOPs and throughput are as important indicators for model efficiency / complexity as the number of parameters.

   Wall-clock time. On an NVIDIA A6000 GPU, we observe that the reduction in wall-clock time between W-PCDM and a standard VDM follows closely with our approximations. In our CIFAR10 model, there is a 2.3x reduction in wall-clock time during inference versus VDM (26 min to evaluate 30K examples versus 1 hour).

[1] The Efficiency Misnomer. https://arxiv.org/pdf/2110.12894.pdf

Calculating FLOPs of W-PCDM. Note that for an input image of size $C \times H \times W$ (where $C$, $H$, and $W$ are the channels, width, and height respectively) the largest input resolution of any scale in the W-PCDM model is size $3C \times (H / 2) \times (W / 2)$. Each subsequent model then takes an input of size $3C \times (H / 2^s) \times (W / 2^s)$ for $s = 1, \dots S-1$, with the final scale being $C \times (H / 2^S) \times (W / 2^S)$. A U-Net is comprised almost entirely of convolutional layers and attention layers -- plus relatively inexpensive element-wise and group-wise nonlinearities. The FLOPs of a convolutional layer is usually $C_{in} \times C_{out} \times (H_{kernel} \times W_{kernel}) \times (H_{out} \times W_{out})$. Therefore, fixing the network architecture, the FLOPS of a convolutional layer of the W-PCDM model at the first scale level is approximately $1/4$ that of the original model on the original image ($H_{out}$ and $W_{out}$ are halved w.r.t. the original model). Some exceptions are the first and last convolutional layers of the model, and upsampling/downsampling layers. The FLOPs of an attention layer is quadratic in the image size, i.e., $(H \times W) \times (H \times W)$, and thus reduces by 1/16 in the first scale of a wavelet model. Applying the same argument to the remaining scales, one can find that the FLOPs at each scale is approximately $2^{-2s}$ that of the original FLOPs, where $s$ is the scale level. The sum of all FLOPs across scales is thus less than $1/2$ that of the original FLOPs (by applying a crude geometric bound $1/4 \sum_{s=0}^S (1 / 2)^s \leq 1 / 2$).

Clarity of paper.

Mixing up $`\citet`$ and $`\citep`$. Thank you for this catch. We have standardized all citations with $`\citep`$.

Rephrasing page 1 "To regain $p_\theta(x)$ evaluation capabilities", and using $`\paragraph`$ for paragraph titles. Again, thank you kindly for the suggestions -- we have incorporated them.

We greatly appreciate the detailed response of the authors. I now admit the consistency of the optimal transport theorem that connects the EMD loss to the score function with the rebuttal with Reviewer y69M. On the other hand, I am curious about the following points.

1. I have been concerned that the volume-preserving map with respect to the orthogonal wavelets is actually the scaling of the wavelet bases (which is equal to normalization) to preserve volume, and I doubt that the significant gain in terms of NLL (=BPD) is obtained by this aspect of loss design. In a practical point of view, the scaling in each wavelet basis leads to the scaling in the forward (thus reverse process) of the conditional diffusion models with downsampled wavelet coefficients, and in my opinion, the FID gain (and maybe the BPD gain) comes from this rescaling of the diffusion process.

2. As I have concerned, the (conditional) likelihood-based diffusion models that the paper used utilizes multiple neural netowrks, using one network for different scale to work as hierarchical models. The final version of the manuscript should also contain the input-output relationship and the architecture of the score model (including the number of Resblocks, the Attention resolutions, so on....) I confirmed that the number of parameters are considered in Appendix B, but the wall-clock time (or the FLOPs, if possible) is not considered yet. This makes me suspect that the gain of this paper is because of arbitrarily bulked models.

For these unresolved issues, I keep my current review score.

---

Clarity of paper

(Overall) Please do not mix up $`\citet`$ (without parenthesis) and $`\citep`$ (with parenthesis).
(Page 1) To regain $p_\theta(x)$ evaluation capabilities $\to$ For $p_\theta(x)$ to regain evaluation capabilities
(Overall) Please use $`\paragraph`$ instead of bold text for the title of paragraphs.