Probing the Latent Hierarchical Structure of Data via Diffusion Models

## Abstract

Abstract High-dimensional data must be highly structured to be learnable. Although the compositional and hierarchical nature of data is often put forward to explain learnability, quantitative measurements establishing these properties are scarce. Likewise, accessing the latent variables underlying such a data structure remains a challenge.
In this work, we show that forward-backward experiments in diffusion-based models, where data is noised and then denoised to generate new samples, are a promising tool to probe the latent structure of data.
We predict in simple hierarchical models that, in this process, changes in data occur by correlated chunks, with a length scale that diverges at a noise level where a phase transition is known to take place. Remarkably, we confirm this prediction in both text and image datasets using state-of-the-art diffusion models.
Our results show how latent variable changes manifest in the data and establish how to measure these effects in real data using diffusion models.

## 1 Introduction

Generative artificial intelligence (AI) systems have demonstrated remarkable capabilities in synthesizing data across various modalities, including images (Betker et al., ; Rombach et al., ) and text (Brown, ; Ouyang et al., ; Touvron et al., ). The underlying reasons behind these achievements remain poorly understood. Indeed, natural data are often high-dimensional and thus generically intractable due to the curse of dimensionality (Luxburg & Bousquet, ; Bach, ).
Hence, to be learnable, the distribution of the data must be highly structured. Characterizing this structure is a fundamental challenge central to any theory of learning.

Hierarchical compositionality (Patel et al., ; Mossel, ; Poggio et al., ; Malach & Shalev-Shwartz, ; Schmidt-Hieber, ; Cagnetta et al., ) is a candidate property put forward to rationalize the success of deep architectures. In this view, data can be decomposed into features organized hierarchically. It is well-established that the grammatical structure of most languages is approximately context-free and hierarchical (Chomsky, ; Jäger & Rogers, ),
although it is unclear how well this structure can capture language semantics (Goldberg, ; ).
Likewise, pattern theory (Grenander, ) posits that images have a hierarchical structure. In both cases, obtaining quantitative evidence characterizing this hierarchy and building tools to determine the associated latent variables remain a challenge.

Generative denoising diffusion probabilistic models (DDPMs) offer a new handle to tackle this challenge, particularly through forward-backward experiments, where a controlled level of noise is added to a starting image and then removed to generate a new one. (Ho et al., ; Sclocchi et al., ; Behjoo & Chertkov, ). For small amounts of noise, low-level features of the image change (Ho et al., ; Sclocchi et al., ). Passed a transition point, the class is likely to change (Ambrogioni, ; Sclocchi et al., ; Biroli et al., ; Li & Chen, ), but remarkably
some of the low-level features of the original image are still retained, as predicted in simple hierarchical models of data structure (Sclocchi et al., ).
However, the empirical tests in these works were limited to images and did not explore other data modalities. Moreover, the geometrical structure of the changes occurring in such a process is not known.

In this work, we derive the length scale associated with changes occurring in the forward-backward protocol in a synthetic generative model of hierarchical data, and we show experimentally that our predictions hold in both language and image datasets.
The synthetic model we consider belongs to the class of probabilistic context-free grammars (Rozenberg & Salomaa, ); it is thus defined on a tree graph, and its forward-backward diffusion experiments can be done exactly with a message-passing algorithm.
Specifically, our contributions are as follows:

- •
In the generative model of hierarchically structured data,
using a mean-field description of the forward-backward diffusion process, we show theoretically that changes in the tokens are correlated over a length scale that diverges at the class transition. This phenomenology is a signature of the hierarchy in the data structure, indicating changes in deep latent variables.
- •
We validate our theoretical predictions performing numerical experiments on our synthetic data with a diffusion process used in practice for discrete data, showing the same phenomenology predicted by our theory. To do so, we measure the dynamical susceptibility, an observable used to study the dynamics in physical systems.
- •
We perform forward-backward experiments with state-of-the-art masked diffusion language models (MDLM) (Sahoo et al., ) on ${\tt WikiText}$. We show the presence of a peaking correlation length in the token changes at a finite inversion time, consistently with our theoretical model.
- •
We perform the same experiments with vision Denoising Diffusion Probabilistic Models (DDPM) (Nichol & Dhariwal, ) on ${\tt ImageNet}$. We tokenize the resulting images using the patch embeddings of a contrastively pre-trained vision encoder (Radford et al., ) and show that the correlations of token changes display a qualitative agreement with our analysis.

Overall, our results show how changes in latent variables affect visible data, and directly support the idea that a hierarchical latent structure is central to both language and vision modalities. Moreover, our work puts forward the forward-backward protocol as a tool to probe the latent hierarchical structure of real data.

## 2 Background

### 2.1 Diffusion models

Denoising diffusion models are generative models designed to sample from a distribution by reversing a step-by-step noise addition process (Sohl-Dickstein et al., ; Ho et al., ; Song & Ermon, ; Song et al., ).
Let $t$ indicate the time step in a sequence $[0,\dots,T]$, $q(\cdot)$ the data distribution we wish to sample from and ${\bm{x}}_{0}\sim q({\bm{x}}_{0})$ a sample drawn from this distribution.
Diffusion models consist of: a forward process generating a sequence of increasingly noised data {${\bm{x}}_{t}\}_{1\leq t\leq T}$,
$q({\bm{x}}_{1},\dots,{\bm{x}}_{T}|{\bm{x}}_{0})=\prod_{t=1}^{T}q({\bm{x}}_{t}|
{\bm{x}}_{t-1})$, where at the final time $T$, ${\bm{x}}_{T}$ corresponds to pure noise; a backward process, which reverts the forward one by gradually removing noise. This process is typically obtained by learning the backward transition kernels $p({\bm{x}}_{t-1}|{\bm{x}}_{t})$ using a neural network. Sampling from $q(\cdot)$ is achieved by sampling noise ${\bm{x}}_{T}\sim q({\bm{x}}_{T})$ and then applying the learned backward process to obtain a new sample ${\bm{x}}_{0}\sim q({\bm{x}}_{0})$. Different diffusion models correspond to different choices of the forward process, depending on the data space under consideration (see Yang et al. () for a review).

##### Discrete data

For discrete data, like text, ${\bm{x}}_{0}$ consists of a sequence of tokens $x_{0,i}$, $i\in[d]$, each corresponding to a symbol belonging to a vocabulary ${\mathbb{V}}$.
In this case, we consider masked diffusion with an absorbing state by introducing an additional $[{\sf MASK}]$ symbol (Austin et al., ).
At time step $t$, each non-masked token either stays unchanged or transitions to $[{\sf MASK}]$ with some probability $\beta_{t}$.
Using a one-hot-encoding representation of these $|{\mathbb{V}}|+1$ states, the forward transition matrix reads

$$ $\displaystyle q(x_{t,i}|x_{t-1,i})=(1-\beta_{t})\mathbb{I}+\beta_{t}\mathbf{1} \mathbf{e}_{M}^{\top},$ (1) $$

where $\mathbb{I}$ is the identity matrix, $\mathbf{1}$ a column vector of all ones and $\mathbf{e}_{M}$ the one-hot-encoding vector of the $[{\sf MASK}]$ symbol.
At the final time $T$, $x_{T,i}=[{\sf MASK}]$ for every $i\in[d]$.
In the following, we consider the noise schedule $\beta_{t}=(T-t+1)^{-1}$ such that $q(x_{t,i}=[{\sf MASK}]\,|\,{\bm{x}}_{0})=t/T$ (Austin et al., ).

##### Continuous data

For continuous data, like images, corresponding to ${\bm{x}}_{0}\in\mathbb{R}^{d}$, we consider the time-discretized Gaussian diffusion introduced in (Ho et al., ). In this case, the forward transition matrix reads

$$ $\displaystyle q({\bm{x}}_{t}|{\bm{x}}_{t-1})=\mathcal{N}({\bm{x}}_{t};\sqrt{1- \beta_{t}}{\bm{x}}_{t-1},\beta_{t}\mathbb{I}),$ (2) $$

where $\mathcal{N}$ represents the Gaussian probability distribution and the sequence {$\beta_{t}\}_{1\leq t\leq T}$ is the variance schedule. At the final time $T$, ${\bm{x}}_{T}\sim\mathcal{N}(0,\mathbb{I})$.

##### Forward-backward experiments

Forward-backward experiments involve inverting the diffusion process at some intermediate time $t\leq T$. Starting from ${\bm{x}}_{0}$, the forward process up to time $t$ produces a noisy sample ${\bm{x}}_{t}\sim q({\bm{x}}_{t}|{\bm{x}}_{0})$. The backward process, obtained from a trained diffusion model, produces a new sample $\hat{{\bm{x}}}_{0}(t)\sim p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{t})$.
Thus, we refer to $t$ as the *inversion time* of the forward-backward process.

### 2.2 The Random Hierarchy Model (RHM)

The RHM is a generative model of hierarchically structured data introduced by (Cagnetta et al., ). It belongs to the class of context-free grammars in the field of language theory (Rozenberg & Salomaa, ), and assumes that production rules are random. In its simplest version, it consists of:

- •
A regular tree graph of depth $L$ and branching factor $s$. Each node of the tree corresponds to a discrete random variable.
- •
A discrete vocabulary ${\mathbb{V}}^{(\ell)}$ of size $v$ for each level $\ell=0,1,\dots,L$ of the tree. We call $\ell=0$ the level of the leaves and $\ell=L$ the root level.
- •
A set of production rules defining how each symbol belonging to ${\mathbb{V}}^{(\ell)}$ can be represented at the lower level with the symbols of $({\mathbb{V}}^{(\ell-1)})^{\otimes s}$. For each element of ${\mathbb{V}}^{(\ell)}$, there are $m$ equivalent lower-level representations, which are all distinct and chosen randomly.

We use the notation ${\textnormal{h}}^{(\ell)}_{i}$ to indicate the variable at level $\ell$ and position $i\in[s^{L-\ell}]$.
The leaf nodes ${\textnormal{h}}^{(0)}_{1},\dots,{\textnormal{h}}^{(0)}_{s^{L}}$ correspond to the visible tokens, while the upper-level nodes represent latent variables.
See Appendix A for an example.
We define the tree distance $\tilde{\ell}$ between two visible tokens as the number of edges between them and their lowest common ancestor. Their corresponding real space distance $r$ is $r=s^{\tilde{\ell}}-1$.
Because of the hierarchical structure generating the data, the visible tokens have non-trivial spatial correlations, which depend on the tree distance between them (Cagnetta & Wyart, ).

#### 2.2.1 Bayes-optimal denoising of the RHM using Belief Propagation

Knowing the production rules and the tree structure of the RHM, the probabilities of the latent variables, conditioned on some observation, can be reconstructed exactly (Sclocchi et al., ) using the Belief Propagation (BP) algorithm (Mezard & Montanari, ).
Specifically, if an RHM datum ${\bm{x}}_{0}$ is corrupted by some noise, e.g. via masking a fraction of tokens, resulting in a noisy observation ${\bm{x}}_{t}$, then BP can be used to:

- •
compute the marginal probabilities of any latent or visible variable, conditioned on the noisy observation ${\bm{x}}_{t}$: $p({\textnormal{h}}^{(\ell)}_{i}|{\bm{x}}_{t})$;
- •
sample directly from the posterior $p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{t})$.

If the noisy observation ${\bm{x}}_{t}$ is produced by a forward diffusion process, then sampling from $p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{t})$ is equivalent to integrating exactly the backward diffusion process starting from ${\bm{x}}_{t}$.

Notice that, in this case, since we know the rules that generated the data, BP gives access to $p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{t})$ directly. This is a different situation with respect to real data, like images and text, where $p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{t})$ is unknown.
In the case of the RHM, therefore, sampling from $p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{t})$ with BP corresponds to having access to a neural network achieving perfect generalization in the backward diffusion process.

#### 2.2.2 Diffusion processes in the RHM

For the RHM data, we consider two different diffusion processes.

- •
$\epsilon$-process This is a simplified diffusion process where the noise acts uniformly on the BP priors at the visible tokens (see Section A.1.2 for details). The noising process is controlled by a noise-to-signal ratio $\epsilon\in[0,1]$, which plays the role of time in the standard diffusion processes, such that $\epsilon=0$ at $t=0$ and $\epsilon=1$ at $t=T$. Starting from an RHM datum ${\bm{x}}_{0}$, we indicate with ${\bm{x}}_{\epsilon}$ the noisy observation at the leaf priors. Therefore, BP computes the marginals $p({\textnormal{h}}^{(\ell)}_{i}|{\bm{x}}_{\epsilon})$ and samples from $p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{\epsilon})$.
We study theoretically this process in section 3.1 through a mean-field approach.
- •
Masking diffusion with an absorbing state This is the diffusion process described in section 2.1 for discrete data, which is commonly used in practice.
We study it numerically with BP in section 3.2.

##### Phase transition in the class reconstruction of the RHM

Sclocchi et al. () showed that there exists a regime of the RHM parameters where the probability of reconstructing the class in the $\epsilon$ diffusion process, that is $p({\textnormal{h}}_{1}^{(L)}|{\bm{x}}_{\epsilon})$, undergoes a sharp phase transition at a critical noise level $\epsilon^{*}$ in the limit of large $L$.
Therefore, sampling $\hat{{\bm{x}}}_{0}(\epsilon)\sim p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{\epsilon})$, for $\epsilon<\epsilon^{*}$, $\hat{{\bm{x}}}_{0}(\epsilon)$ and ${\bm{x}}_{0}$ share the same latent ${\textnormal{h}}_{1}^{(L)}$ (i.e. they belong to the same class), while, for $\epsilon>\epsilon^{*}$, the probability that $\hat{{\bm{x}}}_{0}(\epsilon)$ and ${\bm{x}}_{0}$ share the same class corresponds to the random chance $1/v$.

In Figure 9, we show numerically that also in masking diffusion the probability of reconstructing the class $p({\textnormal{h}}_{1}^{(L)}|{\bm{x}}_{t})$ undergoes a phase transition at a specific inversion time $t^{*}$.

## 3 Correlated Blocks of Tokens Change Near the Transition

In this section, we characterize the statistics of how the input tokens of the RHM change in the forward-backward experiments.
Let $x_{0,i}$ denote the $i$-th input token, $i\in[s^{L}]$, and $\hat{x}_{0,i}(t)$ the same token after undergoing a forward-backward experiment with inversion time $t$. We seek to compute the correlations between changes in the tokens as a function of the inversion time $t$. For each token position $i$, we introduce a quantity characterizing the dynamics, a spin variable $\sigma_{i}(t)$. It takes the value $-1$ if after the forward-backward process the symbol of $x_{0,i}$ and $\hat{x}_{0,i}(t)$ remains the same, and $+1$ if the token changed symbol. Using this definition, we can compute the dynamical correlation function between the changes of tokens at positions $i$ and $j$, i.e.,

$$ $\mathcal{C}_{{\bm{x}}_{0},ij}(t)=\langle\sigma_{i}(t)\sigma_{j}(t)\rangle- \langle\sigma_{i}(t)\rangle\langle\sigma_{j}(t)\rangle,$ (3) $$

where $\langle\cdot\rangle$ denotes averaging over different stochastic trajectories. By further averaging over the initial point ${\bm{x}}_{0}$, we define the average dynamical correlation function as $\mathcal{C}_{ij}(t)=\overline{\mathcal{C}_{{\bm{x}}_{0},ij}(t)}$.
Given the correlations, we compute the dynamical susceptibility $\chi(t)$. This quantity is used to study the dynamics in physical systems (Donati et al., ; Toninelli et al., ) and is defined as

$$ $\chi(t)=\frac{\sum_{i=1}^{s^{L}}\sum_{j=1}^{s^{L}}\mathcal{C}_{ij}(t)}{\sum_{i =1}^{s^{L}}\mathcal{C}_{ii}(t)},$ (4) $$

where we normalized by the sum of auto-correlations. Intuitively, the susceptibility measures the volume of the blocks of tokens that change together.

In the case of the $\epsilon$-process, where $\hat{{\bm{x}}}_{0}(\epsilon)$ is sampled from $p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{\epsilon})$, the same definitions hold for the quantities $\mathcal{C}_{ij}(\epsilon)$ and $\chi(\epsilon)$.

### 3.1 Mean-field theory of the ϵ italic-ϵ \epsilon italic_ϵ -process of the RHM

The average correlation function $\mathcal{C}_{ij}(\epsilon)$ can be computed for the $\epsilon$-process of the RHM through a mean-field approximation. This mean-field approach consists of computing the average BP messages at each layer $\ell$, where the average is performed over the possible realizations of the RHM rules.
Let’s consider two leaf nodes ${\textnormal{h}}_{i}^{(0)}$ and ${\textnormal{h}}_{j}^{(0)}$ connected to the common ancestor ${\textnormal{h}}_{k}^{(\ell)}$ at layer $\ell$ through the nodes ${\textnormal{h}}_{l}^{(\ell-1)}$ and ${\textnormal{h}}_{m}^{(\ell-1)}$ (see Figure 1 for an illustration).
Their associated spin variables are therefore $\sigma_{i}^{(0)}$, $\sigma_{j}^{(0)}$, $\sigma_{k}^{(\ell)}$, $\sigma_{l}^{(\ell-1)}$ and $\sigma_{m}^{(\ell-1)}$, where we omit the $\epsilon$ dependence to lighten the notation. Given the tree structure, the joint probability distribution $P(\sigma_{i}^{(0)},\,\sigma_{j}^{(0)})$ can be written as

$$ $\displaystyle\begin{aligned} P(\sigma_{i}^{(0)},\,\sigma_{j}^{(0)})=\sum_{ \sigma_{l}^{(\ell-1)},\sigma_{m}^{(\ell-1)}}P(\sigma_{i}^{(0)}|\sigma_{l}^{( \ell-1)})\ P(\sigma_{j}^{(0)}|\sigma_{m}^{(\ell-1)})\ P(\sigma_{l}^{(\ell-1)}, \sigma_{m}^{(\ell-1)}).\end{aligned}$ (5) $$

Figure: Figure 1: Example of leaf nodes ${\textnormal{h}}_{i}^{(0)}$, ${\textnormal{h}}_{j}^{(0)}$ connected to the common ancestor ${\textnormal{h}}_{k}^{(\ell)}$ through ${\textnormal{h}}_{l}^{(\ell-1)}$ and ${\textnormal{h}}_{m}^{(\ell-1)}$.
Refer to caption: x1.png

Each element in the sum of Equation 5 can be written in terms of BP messages, and its average value can be computed by averaging over the realizations of RHM rules.
The average of $P(\sigma_{i}^{(0)}|\sigma_{l}^{(\ell-1)})$ and $P(\sigma_{j}^{(0)}|\sigma_{m}^{(\ell-1)})$ can be written as a $2\times 2$ matrix ${\bm{T}}^{(\ell-1)}$ only depending on the layer $\ell-1$.
Similarly, also the average of the joint probability $P(\sigma_{l}^{(\ell-1)},\,\sigma_{m}^{(\ell-1)})$ can be represented as a $2\times 2$ matrix ${\bm{C}}^{(\ell-1)}$.
In the mean-field approximation, we neglect the fluctuations of these quantities around their means. Therefore, we compute the average joint probability ${\bm{P}}(\sigma_{i}^{(0)},\,\sigma_{j}^{(0)})$ by substituting the elements in the product of Equation 5 with their means.
For spin variables $i$ and $j$ at tree distance $\ell$, we have

$$ ${\bm{P}}(\sigma_{i}^{(0)},\,\sigma_{j}^{(0)})={\bm{T}}^{(\ell-1)}\ {\bm{C}}^{( \ell-1)}\ {{\bm{T}}^{(\ell-1)}}^{\top}.$ (6) $$

All the expressions for the above quantities are reported in Section A.2.
With a similar procedure, we can compute the average marginal probability ${\bm{p}}(\sigma^{(0)})$.
From these quantities, we obtain the average correlation function $\mathcal{C}_{ij}(\epsilon)$ at each noise level $\epsilon$.

#### 3.1.1 Dynamical correlation length

In the mean-field approach, the average belief $p_{\ell}$ in the original value of a latent variable at layer $\ell$ can be computed through an iterative map

$$ $p_{\ell}=F(p_{\ell-1}),$ (7) $$

where the functional form of $F(p)$ is reported in Section A.2 and the initial condition $p_{0}$ depends on the noise level $\epsilon$ as $p_{0}=1-\epsilon+\epsilon/v$.
In the limit of large depth $L\rightarrow\infty$, the probability $p_{L}$ of reconstructing the class is given by the fixed points of $F(p)$. In the regime where $p_{L}$ undergoes a phase transition, $F(p)$ has three fixed points: two attractive ones, corresponding to $p=1/v$ and $p=1$, and a repulsive one, corresponding to the non-trivial solution of $p^{*}=F(p^{*})$ with $p^{*}\in(\frac{1}{v},1)$. $p^{*}$ corresponds to a critical noise level $\epsilon^{*}=\frac{1-p^{*}}{1-1/v}$.

In the vicinity of $\epsilon^{*}$ and the limit $L\rightarrow\infty$, we can estimate the typical distance over which token changes are correlated by computing the number of layers $\tilde{\ell}$ after which the probability of reconstructing the latent variables $p_{\tilde{\ell}}$ approaches one of the two trivial fixed points $p=1$ and $p=1/v$. This corresponds to the number of layers required to escape the repulsive fixed point $p^{*}$.

Given the iterative map of Equation 7, we can linearize it around the fixed point $p^{*}$ and iterate for $\ell$ layers,

$$ $\displaystyle\Delta p_{\ell}=\left(\frac{dF(p)}{dp}\Bigg{|}_{p^{*}}\right)^{ \ell}\Delta p_{0},$ (8) $$

where $\Delta p_{\ell}=p_{\ell}-p^{*}$. We have that $\frac{dF(p)}{dp}\big{|}_{p^{*}}>1$ and we use the shorthand notation $F^{\prime}_{*}=\frac{dF(p)}{dp}\big{|}_{p^{*}}$. We want to compute the depth $\tilde{\ell}$ at which
${F^{\prime}_{*}}^{\tilde{\ell}}\ |\Delta p_{0}|=\mathcal{O}(1)$.
In terms of the corruption noise $\epsilon$, we have
${F^{\prime}_{*}}^{\tilde{\ell}}\ |\Delta\epsilon|=\mathcal{O}(1)$,
where $\Delta\epsilon=\epsilon-\epsilon^{*}$. Hence, $\tilde{\ell}\sim-\log|\epsilon-\epsilon^{*}|/\log F^{\prime}_{*}$. From the depth $\tilde{\ell}$, we can compute the correlation length in input space as

$$ $\xi\simeq s^{\tilde{\ell}}\sim|\epsilon-\epsilon^{*}|^{-\nu}\quad\text{with } \nu=\frac{\log s}{\log F^{\prime}_{*}},$ (9) $$

that diverges at the critical point: $\lim_{\epsilon\to\epsilon^{*}}\xi=+\infty$. This divergence of the correlation length at the class transition indicates that large blocks of tokens change in concert. In fact, these large correlated changes are caused by the modifications of deeper and deeper latent variables near the transition (see Figure 3 for an illustration). At both smaller and larger noise levels, the correlation length decays. This behavior of the dynamical correlation functions implies that the dynamical susceptibility also peaks at the transition, a hallmark of critical behavior.

Figure: Figure 2: Correlation measures on diffusion samples of the Random Hierarchy Model (RHM). (a-I) In the $\epsilon$-process, the average correlation function shows a correlation length that is maximal for $\epsilon^{*}\simeq 0.74$, corresponding to the class phase transition, with a system-spanning power-law behavior. The full lines are experiments run with Belief Propagation, while the dashed lines are the corresponding mean-field theory description (Section 3.1), showing excellent agreement. (a-II) Correspondingly, also the average susceptibility shows a peak at the transition $\epsilon^{*}$. (b) The same behavior is observed for the correlation function (b-I) and the susceptibility (b-II) for masking diffusion. In this case, the phase transition is observed for inversion time $t^{*}\simeq 0.3\ T$, where both the correlation length and the susceptibility peak. Data for RHM parameters $v=32$, $m=8$, $s=2$, $L=9$, averaged over $256$ starting data and $256$ diffusion trajectories per starting datum.
Refer to caption: x2.png

### 3.2 Numerical experiments

To test our theoretical predictions for the $\epsilon$-process, in Figure 2 (a-I), we present the average correlation functions $\mathcal{C}(r,\epsilon)$, corresponding to $\mathcal{C}_{ij}(\epsilon)$ averaged on all pairs $ij$ such that their real space distance is $r$, and normalized by the auto-correlation $\mathcal{C}(0,\epsilon)$.
We observe that the correlation function displays a system-spanning power-law behavior at a critical value $\epsilon^{*}\approx 0.74$, while it decays faster with distance when $\epsilon\neq\epsilon^{*}$. This observation implies a correlation length that peaks at the critical value $\epsilon^{*}$. Consistently, also the susceptibility $\chi(\epsilon)$ in Figure 2 (a-II) peaks at this critical value. We compare both the correlation functions and the susceptibility measures with the theoretical predictions obtained by the mean-field theory of the $\epsilon$-process (dashed lines in Figure 2 (a-I) and (a-II)), showing excellent agreement. Moreover, in Figure 8 of Appendix A, we test the prediction for the critical exponent of the correlation length of Equation 9, also showing very good agreement.

In the panels (b-I) and (b-II) of Figure 2, we report the average correlation functions $\mathcal{C}(r,t)$ and susceptibility $\chi(t)$ for masking diffusion at different inversion times $t$. Also in this case, the correlation length and the susceptibility are maximal at a specific critical time $t^{*}\approx 0.3\ T$. From Figure 9, we observe that this critical time $t^{*}$ corresponds to the phase transition in the class reconstruction probability.
Although there is not a simple mapping between the masking probability $t/T$ and the noise level $\epsilon$ in the simplified $\epsilon$-process, the qualitative behaviors in the two settings show a remarkable agreement, validating the relevance of our theoretical predictions for both kind of diffusion process.

Figure: Figure 3: Masking diffusion in the RHM for masking fraction (a) $t/T=0.3$ and (b) $t/T=0.5$. The bottom sequence represents the starting datum ${\bm{x}}_{0}$. The blue (green) symbols are the masked ones in ${\bm{x}}_{t}$ that (do not) change feature in $\hat{{\bm{x}}}_{0}(t)$. The leaves of the tree represent the sampled sequence $\hat{{\bm{x}}}_{0}(t)\sim p(\hat{{\bm{x}}}_{0}|{\bm{x}}_{t})$. In the corresponding tree, the red nodes are those that changed features with respect to the generating tree of ${\bm{x}}_{0}$. We observe that larger blocks of changed tokens reflect changes in deeper latent variables.
Refer to caption: x6.png

### 3.3 Spatial correlations in data are not sufficient to get a susceptibility peak

In the RHM, the latent hierarchical structure induces spatial correlations both between the input tokens and in their changes in the forward-backward diffusion.
Therefore, it is natural to ask whether any model of data displaying spatial correlations, even without a latent hierarchical structure, exhibits the same phenomenology of the RHM in the forward-backward diffusion.

In Appendix B, we show that this is not the case. In particular, we consider a Gaussian random field model with a covariance having a power-law decaying spectrum. This induces spatial correlations in the field that decay algebraically with the distance.
We show that performing forward-backward diffusion at different inversion times $t$ results in a variation field having a correlation length that increases monotonically with $t$ and is maximal at the final time $t=T$. This behavior contrasts sharply with the hierarchical data studied here, where the growing length scale occurs in correspondence with a phase transition at a finite inversion time $t^{*}$.

In fact, the mechanisms behind the dynamical correlations are different. For Gaussian random fields, the noise of the diffusion acts as a low-pass filter, which defines a characteristic scale below which a mode is reconstructed in the backward process. For hierarchically structured data, instead, the spatial correlations in the changes are associated with the changes of latent variables at different levels of the hierarchy. Therefore, a diverging correlation length is present only when the reconstruction probability of the root node (i.e., the class) undergoes a phase transition.

## 4 Experiments on Natural Language and Image Data

This section extends our analysis to real-world scenarios, demonstrating that language and vision diffusion models exhibit the same phenomenology as observed in the RHM.

Figure: Figure 4: Forward-backward experiments with language diffusion models. (a) Forward-backward examples for different masking fractions. The words in blue (green) are those that were masked and changed (did not change), while the words in red changed following the backward process. (b) Normalized correlations as a function of index distance $r=|i-j|$ for different fractions of masked tokens. (c) Susceptibility $\chi(t)$ as a function of masking fraction. The results are averaged over $N_{S}=300$ textual samples, each consisting of $N_{T}=128$ tokens, with $N_{R}=50$ noise realizations for each masking fraction. The susceptibility is given by integrating over the domain $r\in[0,10]$.
Refer to caption: x8.png

##### Language diffusion models

We consider Masked Diffusion Language Models (MDLM) (Sahoo et al., ) utilizing the ${\tt GPT2}$ tokenizer. We perform forward-backward experiments starting from samples from the ${\tt WikiText2}$ dataset. In Figure 4 (a), we illustrate how an initial paragraph changes with different inversion times $t$. At small $t$, only a few isolated words are modified. At intermediate $t$, we clearly observe clusters of words changing in a correlated manner. At large $t$, only a small fraction of the initial sentence remains unchanged (see Appendix C for a presentation of the same data in their larger context). In Figure 4 (b-c), we quantify these observations by measuring the average correlation functions and susceptibility(^1^11To avoid finite size effects due to imposing a fixed masking fraction, we integrate the average correlation function up to the maximal correlation length $r\sim\mathcal{O}(10)$.). Strikingly, in line with the phenomenology obtained for the RHM, we find a growing correlation length as $t$ increases, reaching a maximum of $7\div 8$ tokens at a critical inversion time $t^{*}\approx 0.6\,T$, followed by a subsequent decline. Additionally, the susceptibility peaks at $t^{*}$, establishing the existence of a phase transition for the language modality.

##### Vision diffusion models

Figure: Figure 5: Examples of images generated at different inversion times $t$ with forward-backward diffusion. The first column represents the starting images ${\bm{x}}_{0}$, while the other columns represent the generated ones $\hat{{\bm{x}}}_{0}(t)$. The grid indicates the tokens represented inside the CLIP vision encoder. For inversion time $t>0$, the red patches indicate the token embeddings that have a variation magnitude larger than a fixed threshold. These patches of variation appear in domains of growing size around the class transition, observed for $t^{*}\approx 0.6\div 0.7T$ (Figure 6).
Refer to caption: extracted/5935334/Figs/n02009912_ILSVRC2012_val_00005314_t0.png

We extend our analysis to computer vision by considering Improved Denoising Diffusion Probabilistic Models (Nichol & Dhariwal, ), trained on the ${\tt ImageNet}$ dataset. To compute the correlation between changes in the image tokens, we follow recent trends in multimodal LLMs (Liu et al., ; Dai et al., ). Specifically, we divide each image into $7\times 7$ patches and use the last-layer embeddings for each patch from a CLIP ViT-B32 (Radford et al., ) to tokenize the image. Let ${\bm{x}}_{0,i}$ denote the embedding of the $i$-th patch, where $i=(k,l)$ with $k,l\in[7]$. After the forward-backward process, the variation of each patch embedding is given by $\Delta{\bm{x}}_{i}(t)={\bm{x}}_{0,i}-\hat{{\bm{x}}}_{0,i}(t)$.
We then compute the average correlations between the norms of these variations:

$$ $\mathcal{C}_{ij}(t)=\overline{\langle\|\Delta{\bm{x}}_{i}(t)\|\,\|\Delta{\bm{x }}_{j}(t)\|\rangle-\langle\|\Delta{\bm{x}}_{i}(t)\|\rangle\langle\|\Delta{\bm{ x}}_{j}(t)\|\rangle}.$ (10) $$

The susceptibility is subsequently obtained as $\chi(t)=\sum_{ij}\mathcal{C}_{ij}(t)/\sum_{ii}\mathcal{C}_{ii}(t)$.
In Figure 5, we present some examples of starting images and generated ones at different inversion times $t$, together with the grid representing their tokenization. We observe that, for increasing $t$, new semantic elements appear in the generated images, corresponding to blocks of tokens changing in concert.
In Figure 6, we present the average correlation functions and the susceptibility for vision DDPMs, starting from samples of the ${\tt ImageNet}$ validation set (Deng et al., ). At a critical inversion time $t^{*}\approx 0.6\div 0.7\,T$, we observe a peak in susceptibility, signaling the class phase transition in these models. This finding highlights the compositional semantic structure of image data, similar to the phase transitions observed in language diffusion models and the RHM.

Figure: Figure 6: Correlation measures on the variation of CLIP embeddings of images generated with forward-backward diffusion. (a) The average correlation function displays a system spanning power-law behavior for $t^{*}\approx 0.6\div 0.7\ T$, corresponding to the class phase transition (Sclocchi et al., ). (b) In correspondence with the phase transition, the average susceptibility displays a peak. Data obtained with $344$ starting images and $128$ diffusion trajectories per starting image.
Refer to caption: x11.png

## 5 Related Work

##### Phase transitions in diffusion models

Several works have studied phenomena related to phase transitions in diffusion models.
Biroli et al. (); Ambrogioni () show the presence of different dynamical regimes in the diffusion process separated by a ‘speciation’ cross-over where a bimodal distribution merges into a mono-modal one.
Li & Chen () provide bounds for critical time windows appearing in the diffusion of mixtures of strongly log-concave densities.
These works do not consider hierarchical data, and thus do not present growing dynamical susceptibility or length scale at the transition.

##### Hierarchical models of images and text

Generative models have been used to describe the structure of data in several context, including in linguistics and signal processing.
For languages, hierarchically-structured formal grammars are often used as a model of their syntactic structure (Rozenberg & Salomaa, ).
Likewise, pattern theory formalizes the semantic decomposition of visual scenes through a hierarchy of features (Stoyan, ; Jin & Geman, ; Siskind et al., ; Li et al., ).
More recently, images have been described trough a hierarchical decomposition in multi-scale wavelet coefficients (Marchand et al., ; Kadkhodaie et al., ), although the underlying structure, in this case, is not tree-like.

##### Hierarchical models in machine learning theory

Several studies (Mossel, ; Shalev-Shwartz et al., ; Malach & Shalev-Shwartz, ; ) have shown that the correlations between input and task are crucial for determining learnability in hierarchical data. Furthermore, the internal representations of trained deep networks reflect the data’s latent structure (Cagnetta et al., ; Allen-Zhu & Li, ).
Sclocchi et al. () showed that diffusion on hierarchical data can be solved using Belief Propagation. Mei () showed that U-Nets can efficiently approximate the Belief Propagation algorithm on hierarchical data, while Garnier-Brun et al. () provided evidence that transformers can implement the same algorithm.

## 6 Conclusion

In this work, we showed that when data exhibit a hierarchical structure, the changes induced by forward-backward experiments in diffusion models reveal a growing correlation length and susceptibility near a phase transition. At this critical point, changes in the data become highly correlated, reflecting changes in deep latent variables.
In particular, we focused on understanding how modifications in the latent variables manifest in the data, in contrast with common approaches which attempt to reconstruct the latent representations from visible data.

Our predictions for a hierarchical model were confirmed through experiments across different natural data modalities, showing a remarkable level of universality.
This supports the hypothesis that hierarchical and compositional structures are fundamental, universal properties underlying natural data as diverse as images and text.

Future work may include interpreting the large, correlated changes in text in terms of grammatical structure and context variables, possibly sharpening these concepts through the data-driven method introduced in this study.

#### Acknowledgments

We thank Umberto Maria Tomasini for the discussions at the early stages of this project.
We thank Francesco Cagnetta, Daniel Korchinski, and Jack Parley for providing feedback on the manuscript. NL is supported by the EPFL AI4science program.
This work was supported by a grant from the Simons Foundation (#454953 Matthieu Wyart).

## Appendix A The Random Hierarchy Model

##### Data generation in the RHM

A datum of the RHM is generated with the following procedure:

- •
a symbol at the root ${\textnormal{h}}^{(L)}_{1}$ is chosen uniformly at random from ${\mathbb{V}}^{(L)}$;
- •
one of the $m$ production rules associated with the symbol sampled at the root is chosen uniformly at random. This choice defines the values, belonging to ${\mathbb{V}}^{(L-1)}$, of the $s$ children nodes ${\textnormal{h}}^{(L-1)}_{1},\dots,{\textnormal{h}}^{(L-1)}_{s}$;
- •
the procedure is iterated sampling the symbols at layer $\ell-1$ according to the symbols at layer $\ell$ and their respective production rules;
- •
as a result, a string of $s^{L}$ symbols belonging to ${\mathbb{V}}^{(0)}$ is generated at the leaf nodes. This string is the generated datum ${\bm{x}}$.

Figure: Figure 7: Example of RHM with $L=2$, $s=2$, $v=3$, $m=3$. Left: example of production rules with vocabularies ${\mathbb{V}}^{(2)}=\{a,b,c\}$, ${\mathbb{V}}^{(1)}=\{d,e,f\}$, ${\mathbb{V}}^{(0)}=\{g,h,i\}$. Right: one possible datum generated by the production rules, with the hierarchical levels indicated on the right.
Refer to caption: x13.png

### A.1 Denoising the RHM with Belief Propagation

In the RHM, knowing the production rules of its tree structure, the Bayes-optimal denoising of its data can be done exactly using the Belief Propagation (BP) algorithm (Sclocchi et al., ).

In the factor graph of the RHM tree, all latent and visible variables ${\textnormal{h}}^{(\ell)}_{i}$ represent variable nodes, while the RHM rules connecting them are factor nodes.
Each variable node is associated to two BP messages, one coming from below, $\nu_{\uparrow}({\textnormal{h}}^{(\ell)}_{i})$, and one coming from above, $\nu_{\downarrow}({\textnormal{h}}^{(\ell)}_{i})$. The starting point of the BP algorithm is the definition of the messages at the boundaries of the tree, which are the upward messages at the leaves $\nu_{\uparrow}({\textnormal{h}}^{(0)}_{i})$ and the downward message at the root node $\nu_{\downarrow}({\textnormal{h}}^{(L)}_{1})$. Since we consider class-unconditional diffusion processes, we consider the latter as being uniform over the values of ${\mathbb{V}}^{(L)}$. The initialization at the leaves, instead, corresponds to the prior belief on the values of the single visible tokens, which is given by the noisy observation. In the case of diffusion processes, the noisy observation ${\bm{x}}_{t}$ gives prior beliefs $\nu_{\uparrow}({\textnormal{h}}^{(0)}_{i})=p(\hat{x}_{0,i}|x_{t,i})$, which can be computed for the single token by Bayes’ rule and depends on the specific diffusion process under consideration.

#### A.1.1 BP iteration

The initialization of BP is given by the leaf messages $\nu_{\uparrow}({\textnormal{h}}^{(0)}_{i})$, $i\in[s^{L}]$.
For each $s$-patch at level $\ell$, e.g., $\{{\textnormal{h}}^{(\ell)}_{i}\}_{i=1,\dots,s}$, having a common parent node at layer $\ell+1$, e.g. ${\textnormal{h}}^{(\ell+1)}_{1}$, the upward message in the upper level is computed as:

$$ $\displaystyle\tilde{\nu}_{\uparrow}\left({\textnormal{h}}^{(\ell+1)}_{1}=y \right)=\sum_{a_{1},\dots,a_{s}\in{{\mathbb{V}}^{(\ell)}}^{\otimes s}}\psi^{( \ell+1)}\left(y,a_{1},\dots,a_{s}\right)\prod_{i=1}^{s}\ \nu_{\uparrow}\left({ \textnormal{h}}^{(\ell)}_{i}=a_{i}\right),$ (11) $$

$$ $\displaystyle\nu_{\uparrow}({\textnormal{h}}^{(\ell+1)}_{1}=y)=\frac{\tilde{ \nu}_{\uparrow}({\textnormal{h}}^{(\ell+1)}_{1}=y)}{\sum_{y^{\prime}\in{ \mathbb{V}}^{(\ell+1)}}\tilde{\nu}_{\uparrow}({\textnormal{h}}^{(\ell+1)}_{1}= y^{\prime})},$ (12) $$

where the factor ${\psi^{(\ell+1)}}\left(y,a_{1},\dots,a_{s}\right)$ reads

$$ $\displaystyle\begin{aligned} {\psi^{(\ell+1)}}(y,a_{1},...,a_{s})=\begin{cases }1,\quad\text{if }y\rightarrow(a_{1},...,a_{s})\text{ is a rule at layer }( \ell+1)\rightarrow\ell\\ 0,\quad\text{otherwise}.\end{cases}\end{aligned}$ (13) $$

This upward process is iterated from the leaf nodes at $\ell=0$ until the root node at $\ell=L$.
Afterward, BP computes the downward messages.
The initialization at the root node is given by a uniform prior over the symbols of ${\mathbb{V}}^{(L)}$, i.e.

$$ $\nu_{\downarrow}\left({\textnormal{h}}^{(L)}_{1}=a\right)=\frac{1}{v},\quad \forall a\in{\mathbb{V}}^{(L)}.$ (14) $$

For the same $s$-patch at layer $\ell$ and parent node at layer $\ell+1$ as before, the downward message for ${\textnormal{h}}^{(\ell)}_{1}$ is given by

$$ $\tilde{\nu}_{\downarrow}({\textnormal{h}}^{(\ell)}_{1}=a_{1})=\sum_{\begin{ subarray}{c}a_{2},...,a_{s}\in\mathcal{{\mathbb{V}}^{(\ell)}}^{\otimes(s-1)}\\ y\in{\mathbb{V}}^{(\ell+1)}\end{subarray}}{\psi^{(\ell+1)}}(y,a_{1},...,a_{s}) \ \nu_{\downarrow}({\textnormal{h}}^{(\ell+1)}_{1}=y)\prod_{i=2}^{s}\nu_{ \uparrow}({\textnormal{h}}^{(\ell)}_{i}=a_{i}),$ (15) $$

$$ $\nu_{\downarrow}({\textnormal{h}}^{(\ell)}_{1}=a)=\frac{\tilde{\nu}_{ \downarrow}({\textnormal{h}}^{(\ell)}_{1}=a)}{\sum_{a^{\prime}\in{\mathbb{V}}^ {(\ell)}}\tilde{\nu}_{\downarrow}({\textnormal{h}}^{(\ell)}_{1}=a^{\prime})},$ (16) $$

with the same factor node of Equation 13.

At the end of the upward-downward iteration, each variable node ${\textnormal{h}}^{(\ell)}_{i}$ is associated with two BP messages for each symbol of the vocabulary ${\mathbb{V}}^{(\ell)}$: $\nu_{\uparrow}({\textnormal{h}}^{(\ell)}_{i})$ and $\nu_{\downarrow}({\textnormal{h}}^{(\ell)}_{i})$.
Their product gives the marginal probability of the value of the node:

$$ $p({\textnormal{h}}^{(\ell)}_{i}=a)\propto\nu_{\uparrow}({\textnormal{h}}^{( \ell)}_{i}=a)\ \nu_{\downarrow}({\textnormal{h}}^{(\ell)}_{i}=a),\quad a\in{ \mathbb{V}}^{(\ell)}.$ (17) $$

These marginal probabilities are conditioned on the BP messages at the leaf nodes, which can come from a noisy observation of an RHM datum, as is the case for denoising diffusion.

Similarly, sampling from the posterior probabilities given by BP is done by starting sampling from the marginal probability at the root and then iteratively updating the marginal probabilities every time a new node is sampled (Mezard & Montanari, ).

#### A.1.2 Priors at the leaves

##### Masking diffusion

Let’s consider a datum ${\bm{x}}_{0}$ of the RHM undergoing masking diffusion. At any time $t$, the tokens of ${\bm{x}}_{t}$ can have value

$$ $\displaystyle\begin{aligned} x_{t,i}&=x_{0,i},\qquad\ \ \text{if token $i$ has not yet been masked};\\ x_{t,i}&=[{\sf MASK}],\quad\text{if token $i$ has already been masked}.\end{aligned}$ (18) $$

Therefore, given the noisy observation ${\bm{x}}_{t}$, the prior belief $\nu_{\uparrow}({\textnormal{h}}_{i}^{(0)}=a)$ on the value of the token $i$ being equal to $a$ is given by $p(x_{0,i}|x_{t,i})$, that is:

$$ $\displaystyle\begin{aligned} \nu_{\uparrow}\left({\textnormal{h}}_{i}^{(0)}=a \right)&=\delta_{a,\overline{a}},\quad\text{if }x_{t,i}=\overline{a}\in{ \mathbb{V}}^{(0)};\\ \nu_{\uparrow}\left({\textnormal{h}}_{i}^{(0)}=a\right)&=1/v,\quad\forall a\in {\mathbb{V}}^{(0)}\ \text{if }x_{t,i}=[{\sf MASK}].\end{aligned}$ (19) $$

##### ϵ italic-ϵ \epsilon italic_ϵ -process

In this process, instead of running a forward diffusion process, we act directly on the leaf priors.
We introduce a noise-to-signal ratio $\epsilon\in[0,1]$, which controls the noise level instead of the diffusion time $t$.
Starting from a datum ${\bm{x}}_{0}$, whose $i$-th token has value $\overline{a}\in{\mathbb{V}}^{(0)}$, the prior beliefs at the leaf node $i$ taking values in ${\mathbb{V}}^{(0)}$ are defined as

$$ $\displaystyle\begin{aligned} \begin{cases}\nu_{\uparrow}\left({\textnormal{h}} ^{(0)}_{i}=\overline{a}\right)&=1-\epsilon+\epsilon/v,\quad\quad\text{ for }x_ {0,i}=\overline{a};\\ \nu_{\uparrow}\left({\textnormal{h}}^{(0)}_{i}=a\right)&=\epsilon/v,\quad\quad \quad\quad\qquad\forall a\in{\mathbb{V}}^{(0)}\setminus{\overline{a}}.\end{ cases}\end{aligned}$ (20) $$

The role of $\epsilon$ is decreasing the prior belief on the starting value of a token.
This process can be interpreted as an averaged forward diffusion process, where the average is made over different forward trajectories. In the example of masking diffusion, calling $1-\alpha_{t}$ the probability of a token being masked at time $t$, the average prior beliefs at the leaves read

$$ $\displaystyle\begin{aligned} \begin{cases}\left\langle\nu_{\uparrow}({ \textnormal{h}}^{(0)}_{i}=\overline{a})\right\rangle&=\alpha_{t}+\frac{1- \alpha_{t}}{v},\quad\quad\text{ where }x_{0,i}=\overline{a};\\ \left\langle\nu_{\uparrow}({\textnormal{h}}^{(0)}_{i}=a)\right\rangle&=\frac{1 -\alpha_{t}}{v},\quad\quad\quad\quad\qquad\forall a\in{\mathbb{V}}^{(0)} \setminus{\overline{a}},\end{cases}\end{aligned}$ (21) $$

which have the same functional form as Equation 20 by identifying $\epsilon=1-\alpha_{t}$. Both $\epsilon$ and $1-\alpha_{t}$ vary between 0 0 and $1$ and play the role of noise-to-signal ratio in their respective processes. However, the fluctuations of the upward beliefs around their mean in the masking diffusion change the statistics of the BP messages propagating upwards and make the mapping $\epsilon=1-\alpha_{t}$ inaccurate. For example, in the experimental data of Figure 2, the phase transition in the $\epsilon$-process is located at $\epsilon^{*}\simeq 0.74$, while it is found at $1-\alpha_{t^{*}}=t^{*}/T\simeq 0.3$ for masking diffusion.

### A.2 Mean-field theory of the ϵ italic-ϵ \epsilon italic_ϵ -process

##### Computation of the marginal probabilities

Starting from Equation 20 and the BP iterative Equations 11, 12, 15 and 16, Sclocchi et al. () computed the average BP messages at each layer $\ell$, where the average is performed over the possible choices of the RHM rules.
The result consists in the average messages associated with reconstructing the starting value $\overline{a}\in{\mathbb{V}}^{(\ell)}$ of a latent node ${\textnormal{h}}^{(\ell)}_{i}$,

$$ $\displaystyle\left\langle\nu_{\uparrow}\left({\textnormal{h}}^{(\ell)}_{i}= \overline{a}\right)\right\rangle_{\psi}=p_{\ell},\quad\left\langle\nu_{ \downarrow}\left({\textnormal{h}}^{(\ell)}_{i}=\overline{a}\right)\right \rangle_{\psi}=q_{\ell},$ (22) $$

where the average $\langle\dots\rangle_{\psi}$ is performed over the factor nodes $\psi$ representing the randomly chosen rules of the RHM.
The values of $p_{\ell}$ and $q_{\ell}$ can be computed layer-by-layer through the following iterative maps:

$$ $p_{\ell+1}=F(p_{\ell}),\quad q_{\ell-1}=G(q_{\ell},p_{\ell-1}),$ (23) $$

where

$$ $\displaystyle F(p)$ $\displaystyle=\frac{p^{s}+f\frac{m-1}{mv-1}(1-p^{s})}{p^{s}+f(1-p^{s})},$ (24) $\displaystyle G(q,p)$ $\displaystyle=\frac{q\ p^{s-1}+f\frac{m-q}{mv-1}(1-p^{s-1})}{q\ p^{s-1}+f\frac {m-q}{mv-1}(1-p^{s-1})+(v-1)f\frac{m-q}{mv-1}},$ (25) $$

and $f=\frac{mv-1}{v^{s}-1}$.
The initial conditions are given by

$$ $p_{0}=1-\epsilon+\epsilon/v,$ (26) $$

$$ $q_{L}=1/v.$ (27) $$

Notice that the expectation values $p_{\ell}$ and $q_{\ell}$ only depend on the layer $\ell$ and not on the specific position of the node $i$ inside the layer.
Once $p_{\ell}$ and $q_{\ell}$ have been computed for every layer $\ell=0,\dots,L$, the average marginal probability of reconstructing the original value $\overline{a}\in{\mathbb{V}}^{(\ell)}$ of the variable ${\textnormal{h}}^{(\ell)}$ is given by

$$ $P({\textnormal{h}}^{(\ell)}=\overline{a})=\frac{p_{\ell}q_{\ell}}{p_{\ell}q_{ \ell}+\frac{(1-p_{\ell})(1-q_{\ell})}{v-1}}.$ (28) $$

This marginal probability is conditioned on the initialization of the leaf nodes (Equation 20) and only depends on the layer $\ell$, not on the position of the node inside the layer.
Given the initialization of $q_{L}$, the probability of reconstructing the root node $P({\textnormal{h}}^{(L)}=\overline{a})$, that is the class of the datum, is given by

$$ $P({\textnormal{h}}^{(L)}=\overline{a})=p_{L}.$ (29) $$

Therefore, in the limit of large depth $L\rightarrow\infty$, the value of $p_{L}$ is given by one of the fixed of the iterative map $F(p)$.
When $F^{\prime}(1)>1$, $F(p)$ has two fixed points: $p=1$, which is repulsive, and $p=1/v$, which is attractive. This implies that, in this regime, for any noise level $\epsilon>0$ at the leaf nodes, it is impossible to reconstruct the value of the class better than random chance.
When $F^{\prime}(1)<1$, instead, a third non-trivial fixed point $p^{*}=F(p^{*})$ appears, which is repulsive, while both $p=1$ and $p=1/v$ are now attractive. This implies the presence of a phase transition at a specific noise level $\epsilon^{*}=\frac{1-p^{*}}{1-1/v}$. For $\epsilon<\epsilon^{*}$, the class is reconstructed, for $\epsilon>\epsilon^{*}$ it is not.

##### Computation of the correlation functions

Similar to the marginal probabilities, the average correlation function can also be computed through an annealed average over the RHM rules. Let’s consider two leaf nodes ${\textnormal{h}}_{i}^{(0)}$ and ${\textnormal{h}}_{j}^{(0)}$ connected to the common ancestor ${\textnormal{h}}_{1}^{(\tilde{\ell})}$ at layer $\tilde{\ell}$ through the nodes ${\textnormal{h}}_{1}^{(\tilde{\ell}-1)}$ and ${\textnormal{h}}_{2}^{(\tilde{\ell}-1)}$. Given the tree structure, their joint probability distribution can be written as

$$ $\displaystyle\begin{aligned} &P({\textnormal{h}}_{i}^{(0)},{\textnormal{h}}_{j }^{(0)})=\\ &\sum_{{\textnormal{h}}_{l}^{(\tilde{\ell}-1)},{\textnormal{h}}_{m}^{(\tilde{ \ell}-1)}}P\left({\textnormal{h}}_{i}^{(0)}|{\textnormal{h}}_{l}^{(\tilde{\ell }-1)}\right)P\left({\textnormal{h}}_{j}^{(0)}|{\textnormal{h}}_{m}^{(\tilde{ \ell}-1)}\right)\sum_{{\textnormal{h}}_{k}^{(\tilde{\ell})}}P\left({ \textnormal{h}}_{l}^{(\tilde{\ell}-1)},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)} |{\textnormal{h}}_{k}^{(\tilde{\ell})}\right)P\left({\textnormal{h}}_{k}^{( \tilde{\ell})}\right)\end{aligned}$ (30) $$

In the mean-field approach, the average joint probability only depends on the tree-distance $\tilde{\ell}$ between $i$ and $j$ and not their precise location. Moreover, we are only interested in the probability that both the starting values of ${\textnormal{h}}_{i}^{(0)}$, ${\textnormal{h}}_{j}^{(0)}$ are reconstructed, and the probability of only one of the two is reconstructed. In the following, we use an overline $\overline{\dots}$ to indicate the starting value of a variable to be reconstructed.
We need to compute

$$ $\displaystyle\left\langle P({\textnormal{h}}_{i}^{(0)}=\overline{a}_{i},{ \textnormal{h}}_{j}^{(0)}=\overline{a}_{j})\right\rangle_{\psi},$ (31) $\displaystyle\left\langle P({\textnormal{h}}_{i}^{(0)}=\overline{a}_{i},{ \textnormal{h}}_{j}^{(0)}\neq\overline{a}_{j})\right\rangle_{\psi}=\left \langle P({\textnormal{h}}_{i}^{(0)}\neq\overline{a}_{i},{\textnormal{h}}_{j}^ {(0)}=\overline{a}_{j})\right\rangle_{\psi},$ (32) $$

where the average $\langle\dots\rangle_{\psi}$ is done over the possible choices of RHM rules.
Using the same strategy for the computation of the marginal probabilities, we compute the average of each term in Equation 30 by substituting the BP messages with their averages.
For this purpose, we first define the average marginal conditioned on the downward messages at layer $\hat{\ell}$, $P({\textnormal{h}}^{(\ell)}=\overline{a}^{\ell}|q_{\hat{\ell}}=c)$, with $\ell<\hat{\ell}$. This is computed with Equation 28 by iterating the equations 23 between layers 0 0 and $\hat{\ell}$ and using the initial conditions of Equation 26 and $q_{\hat{\ell}}=c$. Therefore, the marginals of Equation 28 correspond to $P({\textnormal{h}}^{(\ell)}=\overline{a}^{\ell}|q_{L}=1/v)$.
For the marginals in Equation 30 we have:

$$ $\langle P\left({\textnormal{h}}_{k}^{(\tilde{\ell})}=\overline{a}^{(\tilde{ \ell})}_{k}\right)\rangle_{\psi}=P({\textnormal{h}}^{(\tilde{\ell})}=\overline {a}^{(\tilde{\ell})}|q_{L}=1/v),$ (33) $$

that is the average marginal computed in Equation 28;

$$ $\left\langle P\left({\textnormal{h}}_{i}^{(0)}=\overline{a}_{i}|{\textnormal{h }}_{l}^{(\tilde{\ell}-1)}=\overline{a}^{(\tilde{\ell}-1)}_{l}\right)\right \rangle_{\psi}=P({\textnormal{h}}^{(0)}=\overline{a}|q_{\tilde{\ell}-1}=1),$ (34) $$

$$ $\left\langle P\left({\textnormal{h}}_{i}^{(0)}=\overline{a}_{i}|{\textnormal{h }}_{l}^{(\tilde{\ell}-1)}\neq\overline{a}^{(\tilde{\ell}-1)}_{l}\right)\right \rangle_{\psi}=P({\textnormal{h}}^{(0)}=\overline{a}|q_{\tilde{\ell}-1}=0).$ (35) $$

The probability terms of the type $P({\textnormal{h}}^{(0)}\neq\overline{a}|\dots)$ are given by $1-P({\textnormal{h}}^{(0)}\neq\overline{a}|\dots)$.
Since these averages only depend on the layer level $\tilde{\ell}$, they are the same for $\left\langle P\left({\textnormal{h}}_{j}^{(0)}|{\textnormal{h}}_{m}^{(\tilde{
\ell}-1)}\right)\right\rangle_{\psi}$.
The last term to compute is the joint $P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)},{\textnormal{h}}_{m}^{(\tilde{
\ell}-1)}|{\textnormal{h}}_{k}^{(\tilde{\ell})}\right)$ which can be expressed in terms of BP messages:

$$ $\displaystyle\begin{aligned} &P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}=a _{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)}=a_{m}|{\textnormal{h}}_{k}^{( \tilde{\ell})}=y\right)\propto\\ &\sum_{\begin{subarray}{c}a_{m+1},...,a_{s}\in\mathcal{{\mathbb{V}}^{(\ell)}}^ {\otimes(s-2)}\end{subarray}}{\psi^{(\ell)}}(y,a_{l},a_{m},\dots,a_{s})\ \nu_{ \uparrow}({\textnormal{h}}^{(\ell)}_{l}=a_{l})\ \nu_{\uparrow}({\textnormal{h} }^{(\ell)}_{m}=a_{m})\prod_{i\neq l,m}^{s}\nu_{\uparrow}({\textnormal{h}}^{( \ell)}_{i}=a_{i}).\end{aligned}$ (36) $$

Computing the averages over the rules, we have:

$$ $\displaystyle\left\langle P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}= \overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)}=\overline{a}_{m}|{ \textnormal{h}}_{k}^{(\tilde{\ell})}=\overline{y}\right)\right\rangle_{\psi}=p _{\ell-1}^{2}\ /\ Z_{\overline{y}}^{(\tilde{\ell}-1)},$ (37) $\displaystyle\left\langle P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}= \overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)}\neq\overline{a}_{m}|{ \textnormal{h}}_{k}^{(\tilde{\ell})}=\overline{y}\right)\right\rangle_{\psi}=f \ p_{\ell-1}(1-p_{\ell-1})\frac{m-1}{mv-1}\ /\ Z_{\overline{y}}^{(\tilde{\ell} -1)},$ (38) $\displaystyle\left\langle P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}\neq \overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)}\neq\overline{a}_{m}|{ \textnormal{h}}_{k}^{(\tilde{\ell})}=\overline{y}\right)\right\rangle_{\psi}=f \ (1-p_{\ell-1})^{2}\frac{m-1}{mv-1}\ /\ Z_{\overline{y}}^{(\tilde{\ell}-1)},$ (39) $\displaystyle Z_{\overline{y}}^{(\tilde{\ell}-1)}=p_{\ell-1}^{2}+f\frac{m-1}{ mv-1}(1-p_{\ell-1}^{2})$ (40) $$

$$ $\displaystyle\left\langle P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}= \overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)}=\overline{a}_{m}|{ \textnormal{h}}_{k}^{(\tilde{\ell})}\neq\overline{y}\right)\right\rangle_{\psi }=0,$ (41) $\displaystyle\left\langle P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}= \overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)}\neq\overline{a}_{m}|{ \textnormal{h}}_{k}^{(\tilde{\ell})}\neq\overline{y}\right)\right\rangle_{\psi }=f\ p_{\ell-1}(1-p_{\ell-1})\frac{m}{mv-1}\ /\ Z_{y}^{(\tilde{\ell}-1)},$ (42) $\displaystyle\left\langle P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}\neq \overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)}\neq\overline{a}_{m}|{ \textnormal{h}}_{k}^{(\tilde{\ell})}\neq\overline{y}\right)\right\rangle_{\psi }=f\ (1-p_{\ell-1})^{2}\frac{m}{mv-1}\ /\ Z_{y}^{(\tilde{\ell}-1)},$ (43) $\displaystyle Z_{y}^{(\tilde{\ell}-1)}=f\frac{m}{mv-1}(1-p_{\ell-1}^{2})$ (44) $$

We can combine these terms with the marginals Equation 33 to obtain $\left\langle P\left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)},{\textnormal{h}}_{
m}^{(\tilde{\ell}-1)}\right)\right\rangle_{\psi}$.
We can write this probabilities in a $2\times 2$ matrix ${\bm{C}}^{(\tilde{\ell}-1)}$ such that:

$$ $\displaystyle C_{11}^{(\tilde{\ell}-1)}=\left\langle P\left({\textnormal{h}}_{ l}^{(\tilde{\ell}-1)}=\overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}-1)} =\overline{a}_{m}\right)\right\rangle_{\psi},$ (45) $$

$$ $\displaystyle C_{12}^{(\tilde{\ell}-1)}=C_{21}^{(\tilde{\ell}-1)}=\left\langle P \left({\textnormal{h}}_{l}^{(\tilde{\ell}-1)}=\overline{a}_{l},{\textnormal{h} }_{m}^{(\tilde{\ell}-1)}\neq\overline{a}_{m}\right)\right\rangle_{\psi},$ (46) $$

$$ $\displaystyle C_{22}^{(\tilde{\ell}-1)}=\left\langle P\left({\textnormal{h}}_{ l}^{(\tilde{\ell}-1)}\neq\overline{a}_{l},{\textnormal{h}}_{m}^{(\tilde{\ell}- 1)}\neq\overline{a}_{m}\right)\right\rangle_{\psi}.$ (47) $$

Similarly, also the conditional marginals of Equations 34 and 35 can be written as a $2\times 2$ matrix ${\bm{T}}^{(\tilde{\ell}-1)}$:

$$ $T^{(\tilde{\ell}-1)}_{11}=\left\langle P\left({\textnormal{h}}_{i}^{(0)}= \overline{a}_{i}|{\textnormal{h}}_{l}^{(\tilde{\ell}-1)}=\overline{a}^{(\tilde {\ell}-1)}_{l}\right)\right\rangle_{\psi},$ (48) $$

$$ $T^{(\tilde{\ell}-1)}_{12}=\left\langle P\left({\textnormal{h}}_{i}^{(0)}= \overline{a}_{i}|{\textnormal{h}}_{l}^{(\tilde{\ell}-1)}\neq\overline{a}^{( \tilde{\ell}-1)}_{l}\right)\right\rangle_{\psi},$ (49) $$

$$ $T^{(\tilde{\ell}-1)}_{21}=\left\langle P\left({\textnormal{h}}_{i}^{(0)}\neq \overline{a}_{i}|{\textnormal{h}}_{l}^{(\tilde{\ell}-1)}=\overline{a}^{(\tilde {\ell}-1)}_{l}\right)\right\rangle_{\psi},$ (50) $$

$$ $T^{(\tilde{\ell}-1)}_{22}=\left\langle P\left({\textnormal{h}}_{i}^{(0)}\neq \overline{a}_{i}|{\textnormal{h}}_{l}^{(\tilde{\ell}-1)}\neq\overline{a}^{( \tilde{\ell}-1)}_{l}\right)\right\rangle_{\psi}.$ (51) $$

Collecting the values of $\langle P({\textnormal{h}}_{i}^{(0)},{\textnormal{h}}_{j}^{(0)})\rangle_{\psi}$ into a $2\times 2$ matrix ${\bm{P}}({\textnormal{h}}_{i}^{(0)},{\textnormal{h}}_{j}^{(0)})$, we finally obtain

$$ ${\bm{P}}({\textnormal{h}}_{i}^{(0)},{\textnormal{h}}_{j}^{(0)})={\bm{T}}^{( \tilde{\ell}-1)}\ {\bm{C}}^{(\tilde{\ell}-1)}\ {{\bm{T}}^{(\tilde{\ell}-1)}}^{ \top}.$ (52) $$

In the language of the spin variables introduced in Section 3, the probability of reconstructing a variable ${\textnormal{h}}_{i}^{(0)}=\overline{a}_{i}$ is the probability that $\sigma^{0}_{i}=+1$, while ${\textnormal{h}}_{i}^{(0)}\neq\overline{a}_{i}$ corresponds to $\sigma^{0}_{i}=-1$.

Figure: Figure 8: $\epsilon$-process in the RHM ($v=32$, $m=8$, $s=2$, $L=9$): correlation function with respect to the token distance $r$, for noise levels $\epsilon$ close to the transition $\epsilon^{*}\simeq 0.74$. (Inset) The correlation function displays system-spanning power-law decay at the transition $\epsilon^{*}\simeq 0.74$, while it decays faster for noise values $\epsilon\neq\epsilon^{*}$. The length scale at which it departs from the critical behaviour defines the correlation length $\xi$. (Main) Rescaling the distance $r$ with $\xi$ given by Equation 9 and $C(r,\epsilon)$ with $\xi^{a}$, $a=1$, the different correlation functions collapse on a single curve. This implies that the power-law scaling $\xi\sim|\Delta\epsilon|^{-\nu}$ of Equation 9 describes well the peaking of the correlation length around the class transition. For this choice of RHM parameters, $\nu\simeq 1.78$. The exponent $a=1$ is obtained by fitting the critical decay $C(r,\epsilon^{*})\sim r^{-a}$ from the data.
Refer to caption: x15.png

Figure: Figure 9: Masking diffusion in the RHM: probability of reconstructing a (latent) node as a function of the inversion time $t$. This is proportional to the masking fraction $t/T$. The probability is averaged over the nodes at a given layer. The probability of reconstructing a leaf node (layer 0 0) decreases smoothly with the inversion time, while the probability of reconstructing the root node (layer $10$), that is the datum class, undergoes a sharper decay from $1$ to $1/v$ at a critical time $t^{*}\simeq 0.2\div 0.3\ T$. This sharp decay is expected to become a step-like transition in the limit of infinite depth $L\rightarrow\infty$. Data for RHM parameters $v=32$, $m=8$, $s=2$, $L=10$, averaged over $10$ diffusion trajectories per $10$ starting data ${\bm{x}}_{0}$.
Refer to caption: x16.png

## Appendix B Gaussian Random Field Model

Consider $u\in[-1,1]^{d}$. Let ${\bm{x}}(u)$ denote a centered Gaussian random field defined over this domain with translational-invariant isotropic covariance function $K(u,u^{\prime})$. Specifically, the field satisfies $\mathbb{E}[{\bm{x}}(u)]=0$ and $\mathbb{E}[{\bm{x}}(u){\bm{x}}(u^{\prime})]=K(u,u^{\prime})=c(\|u-u^{\prime}\|)$, where $c$ is a function depending only on the Euclidean distance $\|u-u^{\prime}\|$.

Assume that the Fourier coefficients $C(k)$ of $c(\|u-u^{\prime}\|)$ satisfy, for large $\|k\|$, $C(k)=\gamma\|k\|^{-a}+o(\|k\|^{-a}),\;\|k\|\to\infty$, with $0{<}a{<}d$. This implies that the Fourier coefficients $\mathbf{X}(k)$ are independent Gaussian random variables, $\mathbf{X}(k)\sim\mathcal{N}(0,\sigma_{k}^{2})$ with $\sigma_{k}^{2}\asymp\|k\|^{-a}$.

### B.1 Forward-backward experiments in Fourier space

Given the independence of the Fourier coefficients $\mathbf{X}(k)$, we apply the diffusion dynamics to each Fourier coefficient independently. The noising process is given by:

$$ $\mathbf{X}(k)_{t}=\sqrt{1-\beta_{t}}\mathbf{X}(k)_{t-1}+\sqrt{\beta_{t}}\eta, \quad\eta\sim\mathcal{N}(0,1),$ (53) $$

for $t=1,2,\dots,T$, where $\beta_{t}\in(0,1)$ are the diffusion coefficients and $\eta$ are independent standard Gaussian variables.

By unrolling the recursion, the forward dynamics can be expressed as

$$ $\mathbf{X}(k)_{t}=\sqrt{\overline{\alpha_{t}}}\mathbf{X}(k)_{0}+\sqrt{1- \overline{\alpha_{t}}}\eta,\quad\eta\sim\mathcal{N}(0,1),$ (54) $$

where $\overline{\alpha_{t}}=\prod_{t^{\prime}=1}^{t}(1-\beta_{t^{\prime}})$.

We then reverse the process at time $t$, following the backward dynamics:

$$ $\mathbf{X}(k)_{t-1}=\frac{1}{\sqrt{\overline{\alpha_{t}}}}\left(\mathbf{X}(k)_ {t}+\beta_{t}\nabla_{\mathbf{X}(k)}\log q(\mathbf{X}(k)_{t})\right)+\sqrt{ \beta_{t}}z,\quad z\sim\mathcal{N}(0,1),$ (55) $$

where $q(\mathbf{X}(k)_{t})$ is marginal probability density of $\mathbf{X}(k)_{t}$ in the forward process and $\log q(\mathbf{X}(k)_{t})$ is the corresponding score function.

Given the forward process, $q(\mathbf{X}(k)_{t})$ is Gaussian and the score function can be computed explicitly:

$$ $\nabla_{\mathbf{X}(k)}\log q(\mathbf{X}(k)_{t})=-\frac{\mathbf{X}(k)_{t}}{ \overline{\alpha_{t}}\sigma_{k}^{2}+1-\overline{\alpha_{t}}}.$ (56) $$

### B.2 Mode retrieval

Our goal is to determine which Fourier coefficients are retrieved after the reverse process. Specifically, we want to compute the modes $k$ for which the distance between the coefficient obtained at the end of the backward process $\widehat{\mathbf{X}}(k,t)_{0}\sim p(\cdot|\mathbf{X}(k)_{t})$ with the starting coefficient $\mathbf{X}(k)_{0}$ is small:

$$ $|\widehat{\mathbf{X}}(k,t)_{0}-\mathbf{X}(k)_{0}|\ll 1.$ (57) $$

Thus, we consider the signal-to-noise ratio (SNR) for each mode $k$

$$ ${\rm SNR}(\kappa,t)=\frac{\kappa^{-a}}{\overline{\alpha_{t}}^{-1}-1},$ (58) $$

where $\kappa=\|k\|$.

Define the critical wavevector magnitude $\kappa^{*}$ where ${\rm SNR}(\kappa^{*},t)=1$:

$$ $\kappa^{*}=\left(\overline{\alpha_{t}}^{-1}-1\right)^{-1/a}$ (59) $$

Modes with $\kappa<\kappa^{*}$ (low-frequency modes) have ${\rm SNR}>1$ and can be retrieved, while modes with $\kappa>\kappa^{*}$ (high-frequency modes) have ${\rm SNR}>1$ are dominated by the noise in the forward dynamics and cannot be reconstructed.

Figure: Figure 10: Gaussian random field model. (a) Relative modal errors as a function of wave-vector magnitude $|k|$. For $|k|>\kappa^{*}$, errors remain large, indicating unsuccessful retrieval of the Fourier coefficients, while for $|k|<\kappa^{*}$, the errors decrease, signifying successful recovery. (b) Spatial correlation function $\mathcal{C}(r,t)$, showing a power law decay at short distances and a cutoff at long distances. The correlation length increases with inversion time $t$. (c) The susceptibility $\chi(t)$ increases monotonically and reaches its maximum at the inversion time $t=T$.
Refer to caption: x17.png

### B.3 Correlation analysis

We seek to compute the correlation of the changes after reverting the process at time $t$. Let ${\bm{x}}(u,t)$ denote the field obtained after reverting the diffusion process at time $t$, at position $u$. In particular, ${\bm{x}}(\cdot,0)$ denotes the starting random field. Define the difference field ${\bm{z}}(u,t)={\bm{x}}(u,t)-{\bm{x}}(u,0)$. Since the two fields are Gaussian, also ${\bm{z}}(\cdot,t)$ is Gaussian.

We are interested in the following spatial correlation function:

$$ $\mathcal{C}(r,t)=\mathbb{E}[{\bm{z}}(u,t)^{2}{\bm{z}}(0,t)^{2}],$ (60) $$

where $r=\|u\|$.
Using Wick’s theorem, we have

$$ $\mathcal{C}(r,t)=\mathbb{E}[{\bm{z}}(u,t){\bm{z}}(u,t)]\,\mathbb{E}[{\bm{z}}(0 ,t){\bm{z}}(0,t)]+2\mathbb{E}[{\bm{z}}(u,t){\bm{z}}(0,t)]^{2}.$ (61) $$

The first term is a constant independent of $r$, while the second term captures the spatial dependence.

To compute $\mathbb{E}[{\bm{z}}(u,t){\bm{z}}(0,t)]$, we express ${\bm{z}}(u,t)$ in terms of its Fourier coefficients $\textbf{Z}(k,t)$. For modes with $\kappa<\kappa^{*}(t)$, we can assume $\textbf{Z}(k,t)\approx 0$. For modes with $\kappa>\kappa^{*}(t)$, $\widehat{\mathbf{X}}(k,t)_{0}$ is approximately independent of $\mathbf{X}(k)_{0}$. Thus, $\textbf{Z}(k,t)$ for $\kappa>\kappa^{*}(t)$ is a Gaussian random variable with zero mean and variance $2\sigma_{k}^{2}$.

Thus, the covariance of ${\bm{z}}$ is

$$ $\mathbb{E}[{\bm{z}}(u,t){\bm{z}}(0,t)]\approx\int_{\|k\|>\kappa^{*}(t)}e^{ik^{ \top}u}\,2\sigma_{k}^{2}d^{d}k.$ (62) $$

Substituting $\sigma_{k}^{2}\asymp\|k\|^{-a}$, we have:

$$ $\mathbb{E}[{\bm{z}}(u,t){\bm{z}}(0,t)]\approx\int_{\|k\|>\kappa^{*}(t)}e^{ik^{ \top}u}\,2\|k\|^{-a}d^{d}k.$ (63) $$

To evaluate the integral, we consider the asymptotic behavior for different regimes of $r$. At short distances $r\ll 1/\kappa^{*}$, the integral over $k$ is dominated by large $\kappa$ and behaves as $\mathbb{E}[{\bm{z}}(u,t){\bm{z}}(0,t)]\approx C_{1}\,r^{a-d}$, where $C_{1}$ is a constant. At long distances $r\gg 1/\kappa^{*}(t)$, the lower limit $\kappa^{*}(t)$ introduces an effective cutoff and the covariance decays faster than any power law.

Therefore, the correlation function $\mathcal{C}(r,t)$ exhibits algebraic decay with exponent $2(a-d)$ for $r\ll 1/\kappa^{*}$ and faster than any power law for $r\gg 1/\kappa^{*}$.

### B.4 Discussion

For the Gaussian random field model, the correlation length $\xi\sim 1/\kappa^{*}(t)$ is a monotonically increasing function of the inversion time $t$, or noise-to-signal ratio (NSR). As a result, the susceptibility $\chi(t)$ – calculated by integrating the correlation function over space – also increases monotonically and reaches its maximum at the inversion time $t=T$, where the ${\rm NSR}=\infty$. This behavior contrasts sharply with the hierarchical data studied here, where a phase transition occurs at a finite time/${\rm NSR}$. As discussed in the main text, this divergence arises due to the geometry of correlations induced by the hierarchical tree structure, which is absent in the Gaussian random field model.

### B.5 Numerical experiments

In Figure 10 (a), we plot the relative modal errors $\mathcal{E}_{k}=\sigma_{k}^{-1}|\widehat{\mathbf{X}}(k,t)_{0}-\mathbf{X}(k)_{0}|$. For $\|k\|>\kappa^{*}$, the errors remain $\mathcal{O}(1)$, indicating that the coefficients are not retrieved, as predicted by our analysis. Conversely, for $\|k\|>\kappa^{*}$, the errors decay, indicating successful recovery of the coefficients. In panel (b), we present the correlations $\mathcal{C}(r,t)$, which exhibit a power law decay followed by a cutoff. Notably, the correlation length increases monotonically with the inversion time $t$. Finally, in panel (c), we plot the susceptibility $\chi(t)$, which reaches its maximum at $t=T$.

## Appendix C Language Diffusion

### C.1 Setup

Here, we briefly describe the particular realization of discrete diffusion used in the MDLM setting, which is detailed in (Sahoo et al., ).

MDLMs are a form of discrete diffusion model tailored for language generation. Unlike autoregressive (AR) models, MDLMs generate text by gradually unmasking tokens, allowing for non-sequential generation. This process is governed by a forward masking and reverse unmasking process, parameterized using a Rao-Blackwellized objective to improve performance.

##### Forward Process:

The forward process is defined by progressively noising a clean input sequence $x$ using a categorical distribution:

$$ $\displaystyle q(z_{t}|x)=\text{Cat}(z_{t};\alpha_{t}x+(1-\alpha_{t})m),$ (64) $$

where $z_{t}$ is the latent variable at time $t$, representing the noisy version of the input sequence, $x$ is the original, clean sequence of tokens, $\text{Cat}(\cdot;\cdot)$ is a categorical distribution over the possible states, $\alpha_{t}$ is the noise schedule function, strictly decreasing from $1$ to 0 0 as $t$ increases, and $m$ is a one-hot vector representing the special masked token. At each time step, a fraction of the data transitions into the masked state.

##### Reverse Process and Rao-Blackwellization:

The reverse diffusion process reconstructs the original data from noisy observations. It is parameterized using a neural network approximation $x_{\theta}(z_{t},t)$, which predicts clean tokens from noisy inputs:

$$ $\displaystyle p_{\theta}(z_{s}|z_{t})=\begin{cases}\text{Cat}(z_{s};z_{t}),& \text{if }z_{t}\neq m,\\ \text{Cat}\left(z_{s};\frac{(1-\alpha_{s})m+(\alpha_{s}-\alpha_{t})x_{\theta}( z_{t},t)}{1-\alpha_{t}}\right),&\text{if }z_{t}=m.\end{cases}$ (65) $$

where $z_{s}$ is the latent variable at a prior time step $s$ (with $s<t$), $x_{\theta}(z_{t},t)$ is a neural network approximation of $x$ given the noisy observation $z_{t}$ at time $t$, and $p_{\theta}(\cdot|\cdot)$ is the model distribution approximating the true reverse process.

The training objective is a negative evidence lower bound (NELBO), expressed as:

$$ $\displaystyle L_{\text{diffusion}}=\sum_{i=1}^{T}\mathbb{E}_{q}\left[\frac{ \alpha_{t(i)}-\alpha_{s(i)}}{1-\alpha_{t(i)}}\log\langle x_{\theta}(z_{t(i)}), x\rangle\right],$ (66) $$

where $T$ is the number of diffusion steps, $\alpha_{t(i)}$, $\alpha_{s(i)}$ is the noise schedules evaluated at time steps $t(i)$ and $s(i)$, respectively, $\mathbb{E}_{q}$ is the expectation over the forward process defined by $q$, and $\langle x_{\theta}(z_{t(i)}),x\rangle$ is the dot product between the neural network output $x_{\theta}(z_{t(i)})$ and the original input $x$.

##### Continuous-Time Likelihood Bounds:

To achieve a tighter approximation to the ELBO, the discrete objective is extended to continuous time as:

$$ $\displaystyle L_{\infty\text{NELBO}}=\mathbb{E}_{q}\int_{0}^{1}\frac{\alpha^{ \prime}_{t}}{1-\alpha_{t}}\log\langle x_{\theta}(z_{t},t),x\rangle\,dt.$ (67) $$

where $\alpha^{\prime}_{t}$ is the time derivative of the noise schedule $\alpha_{t}$.
The integral evaluates the objective over continuous time, providing a tighter bound on the likelihood.
This formulation is invariant to the specific functional form of the noise schedule $\alpha_{t}$, highlighting the robustness of the MDLM approach.

##### Connection to Masked Language Models:

MDLMs leverage a masked diffusion approach where the training objective is a weighted average of classical masked language modeling (MLM) losses:

$$ $\displaystyle L_{\infty\text{NELBO}}=\mathbb{E}_{q}\int_{0}^{1}\frac{\alpha^{ \prime}_{t}}{1-\alpha_{t}}\sum_{\ell}\log\langle x_{\theta}^{\ell}(z_{t}),x^{ \ell}\rangle\,dt,$ (68) $$

where $x^{\ell}$: The $\ell$-th token in the original sequence, $x_{\theta}^{\ell}(z_{t})$: The neural network’s prediction for the $\ell$-th token given the noisy sequence $z_{t}$.
The summation runs over all tokens in the sequence, effectively establishing a connection between MDLMs and BERT-style encoders while equipping them with generative capabilities.

We employ the MDLM proposed in (Sahoo et al., ) to conduct the forward-backward experiments described in Section 4, by first drawing random texts of a fixed token length from the ${\tt WikiText2}$ database, masking a fixed fraction of the tokens $t$, and then performing the backward diffusion process by using the masked sequence as the initial point for the MDLM model.

### C.2 Examples of Text Samples for the Forward-Backward Experiments

Below, we provide examples of texts generated by the forward-backward process using MDLM seeded from ${\tt WikiText2}$ examples for different masking fractions. Selected samples were shown in the main text in Figure 4 (a). We dub the text results after the forward-backward process as U-turn samples. As can be seen by the color coding, correlated blocks of words change together along the denoising process, as described in Section 3, and the semantic meaning of the paragraphs themselves change along the phase transition. In blue we denote masked tokens which have returned to their initial value after the backward process, while red are masked tokens which returned to different values. Green indicates the changes in the final texts. It can be seen that for small masking fractions such as $t/T=0.1$, most of the tokens do not change after masking, while the amount of changed tokens far exceeds the unchanged ones near the phase transition at $t/T=0.5$, hinting at the long range correlations emerging.

## Appendix D Image Diffusion

For image diffusion, we use the publicly available models from Improved Denoising Diffusion Probabilistic Models (Nichol & Dhariwal, ), trained on the ${\tt ImageNet}$ dataset at resolution $256\times 256$. We use the class-unconditional model to ensure a class phase transition at an intermediate diffusion time.
To tokenize the images in a semantically meaningful manner, we use the last-layer embeddings from a CLIP ViT-B32 (Radford et al., ) encoder.
This procedure crops the images to the size $224\times 224$, which get tokenized in $7\times 7$ patches, each of dimension $32\times 32$. The embeddings at the last layer of the CLIP encoder have dimension $768$.

Figure: Figure 11: Examples of images generated at different inversion times $t$. The grid indicates the tokens represented inside the CLIP vision encoder. For inversion time $t>0$, the red patches indicate the token embeddings that have a variation magnitude larger than a fixed threshold. These patches of variation appear in domains of growing size.
Refer to caption: extracted/5935334/Figs/leo_t0.png

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 5
presentation: 2
questions: See weaknesses.
soundness: 3
strengths: (1) The paper introduces novel approaches for analyzing the structure of inputs using pretrained diffusion and language models.

(2) The authors offer a thorough analysis and derivation, with experimental results closely aligning with theoretical expectations.

(3) Multiple schematic diagrams and data visualizations are included, providing valuable insights into the methods.
summary: The paper examines the hierarchical correlation structures among input tokens using a dynamic correlation function and dynamical susceptibility within a forward-backward experimental framework. These variables reveal how two input tokens respond to perturbations when attempting to recover data from noisy inputs. Analyzing diffusion and language models, the study demonstrates an anticipated correlation aligned with spatial structures.
weaknesses: (1) The paper’s presentation could be improved. While there are numerous figures to aid understanding, the main text is somewhat challenging to follow.

(2) Why is the σ in Equation 3 binary? Wouldn’t a continuous measurement be more appropriate? For instance, a small difference in pixel values might not alter the semantic structure of the images, but it would be captured by binary measurement.

(3) Shouldn’t the spatial correlation structures be content-dependent? For example, if the bird and the laptop in Figure 5 were moved slightly farther from the camera, would this change affect the result shown in Figure 2?

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: Please refer to the weaknesses section.
soundness: 3
strengths: 1. The hierarchical perspective provides novel insights into the diffusion model's mechanism and the application of physics is also refreshing. I feel that the community can benefit from these insights, which may give rise to empirical advancements.
2. The paper is well written and clearly communicates the main ideas.
3. The experiments on natural data (image/text) support the theoretical claims.
summary: This paper aims to understand diffusion models through a hierarchical latent variable model. Through this framework, this paper demonstrates the connection between the noise level and the hierarchical levels, as evidenced by a transition phrase. This paper builds on the tools from physics and illustrate their theoretical model with empirical results on practical models.
weaknesses: 1. It'd be great to see attempts at utilizing the theoretical/empirical observations to advance practical model design. Some discussions along this direction would also be appreciated.
2. The tree model seems overly simplified for real-world data like images and languages. For example, one would imagine two high-level variables could become co-parents for some low-level variables, thus breaking the tree structure. I would appreciate a discussion on this limitation and the applicability of the theoretical framework to more general latent models.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. From the analysis, BP denoising appears to be a one-step method that directly samples $\hat{x_0}$ from the noisy observation $x_t$, differing from typical diffusion denoising that iteratively samples $x_{t-1}$ from $x_t$ throughout the process. Does this discrepancy exist, or are the authors also using a denoising schedule similar to real diffusion models?

2. Can we interpret the maximal correlated length achieved at an intermediate noise level (time step) as the model generating class information or lower-level features? If so, this would contrast with existing observations that diffusion processes follow a coarse-to-fine generation pattern (e.g., https://arxiv.org/abs/2303.02490), where lower-level features are generated at the beginning, not in the middle.

3. Figure 4(a) is somewhat unclear. Combined with (c), it seems the authors are suggesting that the largest correlated changing chunk appears at a masking fraction $t/T \in [0.5, 0.7]$. However, this is not immediately evident from (a) alone.
soundness: 3
strengths: 1. The authors aim to capture hidden hierarchical structures within discrete data using the RHM model, with their RHM+BP framework supporting both discrete and continuous diffusion processes.

2. By applying BP for denoising, the authors rigorously analyze phase transitions in the denoising results and identify the critical noise level needed to induce a change in the data class (or low-level feature).
summary: In this paper, the authors examine the hierarchical structure in high-dimensional data by conducting forward-backward experiments within diffusion-based models. They employ a Random Hierarchy Model (RHM) for the data where the tokens of data are generated from a tree structure of latent, they also use Belief Propagation (BP) as the score function to denoise samples.

The authors focus on the phase transition of the average belief $p_L$ of the RHM's root node by analyzing an iterative mapping (Equation 7) and identifying a critical noise level $\epsilon^*$ at which the transition occurs. Based on that, they also compute the minimum layers $\tilde {l}$ needed for the transition, beyond which $p_l$ would collapse to trivial fixed points $\{1/v,1\}$, indicating either a complete reconstruction or randomization of upper latent variables. At this specific noise level, BP can modify the deepest latent layer $\tilde {l}$, yielding the maximum correlation length (i.e. big "chunks" of data token), which is the distance over which token changes remain correlated.

To characterize this effect, the authors introduce **dynamical susceptibility** which exhibits a first increase then decrease curve as expected. They further demonstrate that the dynamical susceptibility curve has the same trend for forward-backward experiments with diffusion models and synthetic RHM experiments.
weaknesses: 1. The paper is somewhat disorganized and hard to follow, as definitions, derivations, and experimental results are heavily interwoven. To improve clarity, consider using theorems or structured definitions to better organize the content (e.g. by moving some derivations, such as Equations 8 and 9 to appendix and summarizing them as a main theorem).

2. In practice, people use real data + score-based denoising; however, the authors use RHM data + BP denoising instead. This discrepancy is insufficiently justified, making the claim that real-world data shares the same hierarchy as RHM unconvincing. While the authors show a similar phase transition phenomenon between the RHM case and real-world diffusion case, they do not rigorously establish a connection between them. Verification by testing real-world diffusion on RHM data may strengthen this claim.

3. The results are somewhat vague and lack practical insights, as it appears that neither the RHM setup nor the forward-backward experiment has direct practical applications. Although the authors mention interpreting the "chunks" that emerge during the forward-backward experiments, they do not provide further discussion or related work on that.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: L57: Does such a divergence definitely indicate a hierarchical structure or are there other ways/reasons divergence could occur? i.e. is this divergence a “proof” of (or very strong evidence for) hierarchy?

L210: Clarification: so, the epsilon-process is itself a mean-field approximation of the discrete diffusion process, but then you use another mean-field on top of that to compute the correlation?

L356: Can you elaborate on why a susceptibility peak is a ‘smoking gun’ for hierarchy? Just because one nonhierarchical example doesn’t have a susceptibility peak doesn’t mean there might not be others that do?

L511: How (or does) this relate to the diffusion-as-spectral-autoregression point of view? Also there is a typo, 'trough'.

General: 

Does the susceptibility divergence tell us anything about how many levels of hierarchy are likely present? Or just that there is at least one level?

The RHM is discrete, and discrete vs continuous diffusion are rather different; can you justify why RHM should be a good model for continuous data/diffusion as well?

Do the MDLM and ImageNet expts actually confirm that a phase transition occurs? Or do we just observe the susceptibility peak and infer a phase transition by analogy to RHM? In particular, it seems that for ImageNet it might actually be possible to run a classifier to determine whether the class changed.
soundness: 4
strengths: The paper is clearly written and the RHM model is nicely studied. The text and image experiments are well-designed (although I still have some questions, below, about how well the conclusions from RHM transfer to real data). I appreciate the application of ideas from Physics to ML problems.
summary: This paper suggests that forward-backward diffusion experiments can be helpful in uncovering hierachical structure in data. They first study a synthetic Random Hierarchical Model and show that a peak of the dynamical susceptibility (related to correlations between blocks of tokens) occurs at a noise level where a phase transition is known to occur in the RHM (i.e. the latent class at the root changes). They then show peaks in the susceptibility in text and image experiments.
weaknesses: Please see Questions.

Meta-review: This paper investigates the hierarchical structure of data using forward-backward experiments with diffusion models. The authors propose that changes in data occur in correlated chunks, with a characteristic correlation length that diverges at a critical noise level associated with a phase transition. These predictions are supported by experiments on synthetic hierarchical data (using a Random Hierarchy Model, RHM) and real-world datasets. The results demonstrate that susceptibility peaks observed at specific noise levels correspond to transitions in the latent structure.

Decision: Accept (Poster)