marginparsep has been altered.
topmargin has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

A Phase Transition in Diffusion Models Reveals the Hierarchical Nature of Data

Antonio Sclocchi 1  Alessandro Favero 1  Matthieu Wyart 1

††footnotetext:  1Institute of Physics, École polytechnique fédérale de
Lausanne (EPFL). Correspondence to<EMAIL_ADDRESS>

Preprint.

###### Abstract

Understanding the structure of real data is paramount in advancing modern
deep-learning methodologies. Natural data such as images are believed to be
composed of features organised in a hierarchical and combinatorial manner,
which neural networks capture during learning. Recent advancements show that
diffusion models can generate high-quality images, hinting at their ability to
capture this underlying structure. We study this phenomenon in a hierarchical
generative model of data. We find that the backward diffusion process acting
after a time $t$ is governed by a phase transition at some threshold time,
where the probability of reconstructing high-level features, like the class of
an image, suddenly drops. Instead, the reconstruction of low-level features,
such as specific details of an image, evolves smoothly across the whole
diffusion process. This result implies that at times beyond the transition,
the class has changed but the generated sample may still be composed of low-
level elements of the initial image. We validate these theoretical insights
through numerical experiments on class-unconditional ImageNet diffusion
models. Our analysis characterises the relationship between time and scale in
diffusion models and puts forward generative models as powerful tools to model
combinatorial data properties.

##  1 Introduction

Figure 1: Images generated by a denoising diffusion probabilistic model
starting from the top-left image and inverting the dynamics at different times
$t$. $T$ corresponds to the time scale when the forward diffusion process
converges to an isotropic Gaussian distribution. At small $t$, the class of
the generated image remains unchanged, with only alterations of low-level
features, such as the eyes of the leopard. After a characteristic time $t$,
the class undergoes a phase transition and changes. However, some low-level
attributes of the original image are retained to compose the new image. For
instance, the wolf is composed with eyes, nose and ears similar to those of
the leopard, and the butterfly inherits its colours and black spots.

Understanding which data are learnable by algorithms is key to machine
learning. Techniques such as supervised, unsupervised, or self-supervised
learning are most often used with high-dimensional data. However, in large
dimensions, for generic data or tasks, learning should require a number of
training examples that is exponential in the dimension Luxburg & Bousquet
(2004); Bach (2017), which is never achievable in practice. The success of
these methods with limited training set sizes, including for images or texts,
implies that such data are highly structured. In particular, the puzzling
phenomenon central to this work is the ability of generative models to compose
a whole new datum by assembling features learned from examples. This feat is
remarkable in foundational models generating text, or in diffusion models
generating images Sohl-Dickstein et al. (2015); Ho et al. (2020); Song & Ermon
(2019); Song et al. (2020) – such as DALL$\cdot$E Betker et al. (2023) and
StableDiffusion Rombach et al. (2022) – investigated here.

Many theoretical studies in machine learning assume that the data distribution
is Gaussian, which is ill-suited to study compositionality. Conversely, we
consider hierarchical generative models of data, inspired by models of formal
grammar and statistical physics, for which the ability of diffusion-based
models to compose can be described exactly.
Diffusion models add noise to images as time increases and learn the reverse
denoising process to generate new samples. If some finite amount of noise is
added, and the process is then reversed, we predict that (i) for small noise,
only low-level features of the image change, as already noticed empirically Ho
et al. (2020), (ii) there exists a threshold noise at which the probability of
remaining in the same class suddenly drops to near-random chance and (iii)
surprisingly, beyond that point, low-level features may persist and compose
the element of a new class. Remarkably, these results appear already evident
in examples such as Figure 1 and can be systematically quantified by
considering the change of internal representations in state-of-the-art CNNs,
whereby a remarkable correspondence to our predictions is found.

Overall, our analysis reveals that diffusion models act at different
hierarchical levels of the data at different time scales within the diffusion
process. It also puts forward hierarchical generative models as a valuable
theoretical tool to address several unanswered questions in machine learning.

###  1.1 Our contributions

On the theoretical side, we consider models belonging to the class of
$L$-level context-free grammars where the structure of the latent variables is
tree-like. We show that the denoising dynamics of diffusion models can be
solved exactly for these models of data using belief propagation (BP). In
particular, the analytical treatment of BP describes the evolution of the
probability of latent variables for different levels of corruption noise and
reveals a phase transition for a threshold noise value when $L$ is large. We
show numerically that the dynamics of latent variables is reflected in the
hidden representation of deep networks previously trained on a supervised
classification task on these data.

On the empirical side, we test these insights by performing a systematic study
of the denoising diffusion dynamics on ImageNet. Following Ho et al. (2020),
we invert the noising process at some time $t$, leading to novel noiseless
images. We then analyse how the representation of a state-of-the-art
convolutional architecture changes between the initial and newly generated
images as a function of both time $t$ and depth of the representation. Our
results confirm the presence of a sharp transition in the class at a given
time or noise level.

###  1.2 Related work

##### Theory of diffusion models

Most of the theoretical work on diffusion models considers simple models of
data. Under mild assumptions on the data distribution, diffusion models
exhibit a sample complexity that scales exponentially with the data dimension
Block et al. (2020); Oko et al. (2023). This curse of dimensionality can be
mitigated through stronger distributional assumptions, such as considering
data lying within a low-dimensional latent subspace De Bortoli (2022); Chen et
al. (2023); Yuan et al. (2023), Gaussian mixture models Biroli & Mézard
(2023); Shah et al. (2023); Cui et al. (2023), graphical models Mei & Wu
(2023), or data distributions that can be factorised across scales Kadkhodaie
et al. (2023). For multimodal distributions such as Gaussian mixtures, the
backward dynamics exhibits a cross-over time when it concentrates toward one
of the modes Biroli & Mézard (2023); Ambrogioni (2023); Raya & Ambrogioni
(2024). Closer to our work, Okawa et al. (2023) consider synthetic
compositional data to empirically show how diffusion models learn to
generalise by composing different concepts. In contrast, we study data that
are not only compositional but also hierarchically structured and make
quantitative predictions on how diffusion models compose features at different
scales.

##### Hierachical models of data structure

Generative models of data have a long history in linguistics Rozenberg &
Salomaa (1997) and have more recently received attention in the context of
machine learning theory. In supervised learning, deep networks can represent
hierarchical tasks more efficiently than shallow networks Poggio et al. (2017)
and can efficiently learn them from an information theory viewpoint Schmidt-
Hieber (2020). For hierarchical models of data, correlations between the input
data and the task are critical for learning Mossel (2016); Shalev-Shwartz et
al. (2017); Malach & Shalev-Shwartz (2018; 2020) and the representations
learnt by neural networks with gradient descent reflect the hidden latent
variables of such models both in CNNs Cagnetta et al. (2023) and transformers
Allen-Zhu & Li (2023). In this work, we use these hierarchical generative
models of data to study the denoising dynamics of diffusion models
theoretically.

##  2 Diffusion models and feature hierarchies

${t=0}$${t=0.5\,T}$${t=0.6\,T}$${t=0.7\,T}$${t=T}$

Figure 2: Left panel. Examples of images generated by reverting the diffusion
process at different times $t$. Starting from the left images $x_{0}$ at time
$t=0$, we generate samples $\hat{x}_{0}(t)\sim p_{\theta}(\hat{x}_{0}|x_{t})$
by first running the diffusion process up to time $t$ and then reverting it,
as described in Section 2.2. At time $t=T$, $x_{T}$ corresponds to isotropic
Gaussian noise and the generated image $\hat{x}_{0}(T)$ is uncorrelated from
$x_{0}$. At intermediate times, instead, a sudden change of the image class is
observed, while some lower-level features are retained. Right panel. Cosine
similarity between the post-activations of the hidden layers of a ConvNeXt
Base Liu et al. (2022) for the initial images $x_{0}$ and the synthesised ones
$\hat{x}_{0}(t)$. Around $t\approx T/2$, the similarity between logits
exhibits a sharp drop, indicating the change in class, while the hidden
representations of the first layers change more smoothly. This indicates that
certain low-level features from the original images are retained for composing
the sampled images also after the class transition. To compute the cosine
similarity, all activations are standardised, i.e., centred around the mean
and scaled by the standard deviation computed on the 50000 images of the
ImageNet-1k validation set. At each time, the values of the cosine similarity
correspond to the maximum of their empirical distribution over $10000$ images
($10$ per class of ImageNet-1k).

This section introduces denoising diffusion probabilistic models and
demonstrates how class-unconditional ImageNet diffusion models operate on
image features across different hierarchical levels at different time scales.
111The code for running the experiments on ImageNet is available at
github.com/pcsl-epfl/forward-backward-diffusion.

###  2.1 Background on denoising diffusion models

Denoising diffusion probabilistic models (DDPMs) Ho et al. (2020) are
generative models designed to sample from a distribution by reversing a step-
by-step noise addition process. In particular, let $q(\cdot)$ represent the
data distribution, and let $x_{0}$ be a sample drawn from this distribution,
i.e., $x_{0}\sim q(x_{0})$. First, DDPMs consist of a forward process which is
a Markov chain generating a sequence of noised data {$x_{t}\\}_{1\leq t\leq
T}$ by introducing isotropic Gaussian noise at each time step $t$ with a
variance schedule {$\beta_{t}\\}_{1\leq t\leq T}$ as follows:

$\displaystyle q(x_{1},\dots,x_{T}|x_{0})$
$\displaystyle=\prod_{t=1}^{T}q(x_{t}|x_{t-1})$
$\displaystyle=\prod_{t=1}^{T}\mathcal{N}(x_{t};\sqrt{1-\beta_{t}}x_{t-1},\beta_{t}\mathbb{I}).$
(1)

Thus, at each time step $t$, we have

$x_{t}=\sqrt{\overline{\alpha_{t}}}x_{0}+\sqrt{1-\overline{\alpha_{t}}}\eta$
(2)

with $\overline{\alpha_{t}}=\prod_{t^{\prime}=1}^{t}(1-\beta_{t^{\prime}})$
and $\eta\sim\mathcal{N}(0,\mathbb{I})$. By selecting the noise schedule such
that $\overline{\alpha}_{t}\to 0$ as $t\to T$, the distribution of $x_{T}$
becomes an isotropic Gaussian distribution. Subsequently, DDPMs reverse this
process by gradually removing noise in a backward process. In this process,
the models learn Gaussian transition kernels $q(x_{t-1}|x_{t})$ by
parametrising their mean and variance using a neural network with parameters
$\theta$ as follows:

$p_{\theta}(x_{t-1}|x_{t})=\mathcal{N}(x_{t-1};\mu_{\theta}(x_{t},t),\Sigma_{\theta}(x_{t},t)).$
(3)

After training, the learnt $p_{\theta}$ can be used to generate novel examples
by initiating the process with $x_{T}\sim\mathcal{N}(0,\mathbb{I})$ and
running it in reverse to obtain a sample from $q$. We refer the reader to Ho
et al. (2020); Nichol & Dhariwal (2021); Dhariwal & Nichol (2021) for more
details regarding the formulation of DDPMs and the technical aspects of the
reverse transition kernels parameterisation with neural networks.

###  2.2 Forward-backward experiments

Previous studies on DDPMs noted that inverting the diffusion process at
different times $t$ starting from an image $x_{0}$ results in samples
$\hat{x}_{0}(t)\sim p_{\theta}(\hat{x}_{0}|x_{t})$ with distinct
characteristics depending on the choice of $t$. Specifically, when
conditioning on the noisy samples $x_{t}$’s obtained by diffusing images from
the CelebA dataset, one finds that for small values of $t$, only fine details
change Ho et al. (2020). We conduct a similar experiment using a class-
unconditional DDPM introduced by Dhariwal & Nichol (2021) on the ImageNet
dataset with 256x256 resolution.

In the left panel of Figure 2, we present some images resulting from this
experiment. For each row, the initial image $x_{0}$ is followed by images
generated by initiating the diffusion process from $x_{0}$, running the
forward dynamics until time $t$, with ${0<t\leq T=1000}$, and ultimately
running the backward dynamics to produce a sample image $\hat{x}_{0}(t)$. Our
observations from these synthetic images are as follows:

  1. (i)

Similarly to the findings in Ho et al. (2020), at small inversion times $t$,
only local features change. Furthermore, the class of the sampled images
remains consistent with that of the corresponding starting images, i.e.,
${{\rm class}(\hat{x}_{0}(t))=\rm class}(x_{0})$ with high probability.

  2. (ii)

There exists a characteristic time scale $t^{*}$ at which the class of the
sampled images undergoes a sudden transition.

  3. (iii)

Even after the class transitions, some low-level features composing the images
persist and are reincorporated into the newly generated image. For instance,
looking at the left panel of Figure 2, in the second row, the jaguar is
composed with the paws and the ears of the dog in the starting picture, or in
the third row, the sofa’s armrests inherit the shape of the car headlights.

Our theory, presented in Section 4 and 5, predicts how features at different
hierarchical levels vary at different time scales of the diffusion dynamics in
accordance with observations (i), (ii), and (iii).

###  2.3 ImageNet hidden representations

To quantify the qualitative observations mentioned earlier, we design an
experiment using the empirically known fact that deep learning models learn
hierarchical representations of the data, with complexity increasing as the
architecture’s depth grows. This phenomenon holds true in both real Olah et
al. (2020); LeCun et al. (2015); Zeiler & Fergus (2014) and synthetic
scenarios Cagnetta et al. (2023); Allen-Zhu & Li (2020). Therefore, we use
these internal representations as a proxy for the compositional structure of
the data. We investigate how the hidden representations of a deep ConvNeXt
Base model Liu et al. (2022), achieving 96.9% top-5 accuracy on ImageNet,
change as a function of the inversion time $t$ and depth $\ell$ of the
representation. In the right panel of Figure 2, we illustrate the value of the
cosine similarity between the post-activations of every hidden layer of the
ConvNeXt for the initial and generated images. We observe that:

  1. (i)

The representations of early layers of the network, corresponding to low-level
and localised features of the images, are the first to change at short
diffusion times and evolve smoothly.

  2. (ii)

At a specific time and noise scale, the similarity between logits experiences
a sharp drop, indicating a transition in the class.

  3. (iii)

Around the class transition, there is an inversion of the similarity curves.
Indeed, the hidden representations in the first layers for the new and
generated images now display the largest alignment. This indicates that low-
level features from the original images can be reused in composing the sampled
images, as qualitatively observed in Figure 2.

To study the robustness of our results with respect to the architecture
choice, in Appendix D, we report the same measurements using ResNet
architectures with varying width and depth He et al. (2016). We find the same
qualitative behaviour as the ConvNeXt in Figure 2.

We now present our theory, which predicts these observations. In Section 3, we
introduce a generative model of data which mimics the hierarchical and
compositional structure of images while being analytically tractable. In
Section 4, we solve exactly the denoising process for this model with
different levels of noising. Finally, in Section 5, we prove the existence of
a phase transition at a critical noise value.

##  3 Hierarchical generative model of data

Figure 3: Sketch of the generative hierarchical model we study in the paper.
In this example, depth $L=3$ and branching factor $s=2$. Different values of
the input and latent variables are represented with different colours.

To model the combinatorial and hierarchical nature of images, we consider
hierarchical generative models Mossel (2016); Shalev-Shwartz et al. (2017);
Malach & Shalev-Shwartz (2018; 2020); DeGiuli (2019); Allen-Zhu & Li (2020);
Cagnetta et al. (2023). Specifically, we study an $L$-level context-free
grammar. It consists of a set of class labels
$\mathcal{C}\equiv\left\\{1,\dots,v\right\\}$ and an alphabet
$\mathcal{A}\,{\equiv}\,\left\\{a_{1},\dots,a_{v}\right\\}$ of $v$ features.
Once the class label $\gamma$ is picked uniformly at random from
$\mathcal{C}$, the data is generated iteratively from a set of production
rules with branching factor $s$ at each layer $\ell$ (see Figure 3, for an
illustration):

$\displaystyle\gamma\mapsto\mu^{(L-1)}_{1},\dots,\mu^{(L-1)}_{s}\text{ for
}\gamma\in\mathcal{C}\text{ and }\mu^{(L-1)}_{i}\in\mathcal{A},$
$\displaystyle\mu^{(\ell)}\mapsto\mu^{(\ell-1)}_{1},\dots,\mu^{(\ell-1)}_{s}\text{
for }\mu^{(\ell)}\in\mathcal{A},\mu^{(\ell-1)}_{i}\in\mathcal{A},$
$\displaystyle\ell\in\\{L-1,\dots,1\\}.$

Since the total size of the data increases by a factor $s$ at each level, the
input data are made of $d\equiv s^{L}$ input features $\mu^{(0)}$. We adopt a
one-hot encoding of these features, ultimately leading to a data vector
$X\in\mathbb{R}^{dv}$. Note that for $\ell\geq 1$, the node variables
correspond to latent variables, and there is no need to specify any choice of
encoding.

For each level $\ell$, we consider that there are $m$ distinct production
rules originating from the same higher-level feature $\mu^{(\ell)}$, i.e.,
there are $m$ equivalent lower-level representations of $\mu^{(\ell)}$. In
addition, we assume that two distinct classes or latent variables cannot lead
to the same low-level representation. This condition ensures, for example,
that two distinct classes never lead to the same data.

We consider the case of the Random Hierarchy Model (RHM) Cagnetta et al.
(2023), for which the $m$ production rules of any latent variable or class are
sampled uniformly at random among the $v^{s}$ possible ones without
replacement. In this case, the total number of possible data produced per
class is $m\cdot m^{s}\cdots m^{s^{L-1}}=m^{\frac{d-1}{s-1}}$, where the
dimension $d$ is defined as $d=s^{L}$. In the following, we use the notation
$X^{(\ell)}_{i}$ to indicate the variable at layer $\ell$ and position
$i\in\\{1,\dots,s^{L-\ell}\\}$.

In the context of unsupervised learning, a key parameter for this model is
$f=m/v^{s-1}$. When $f=1$, all strings of latent variables of size $s$ can be
produced at any level of the hierarchy. This implies that all possible $v^{d}$
input strings are generated, and the data distribution has little structure.
When $f<1$, however, only a small fraction $\sim f^{(d-1)/(s-1)}$ of all
possible strings is generated by the production rules. This implies that
spatial correlations between different input positions appear, reflecting the
hierarchy generating the data.

##  4 Exact denoising of the RHM with message passing

$X_{1}^{(2)}$$\psi^{(2)}$$\nu_{\uparrow}(X_{1}^{(2)})$$X_{1}^{(1)}$$\nu_{\uparrow}(X_{1}^{(1)})$$X_{2}^{(1)}$$\nu_{\uparrow}(X_{2}^{(1)})$$\psi^{(1)}$$X_{1}^{(0)}$$\nu_{\uparrow}(X_{1}^{(0)})$$X_{2}^{(0)}$$\psi^{(1)}$$X_{3}^{(0)}$$X_{4}^{(0)}$$\nu_{\uparrow}(X_{4}^{(0)})$

(Up)

$X_{1}^{(2)}$$\psi^{(2)}$$\nu_{\downarrow}(X_{1}^{(2)})$$X_{1}^{(1)}$$\nu_{\downarrow}(X_{1}^{(1)})$$X_{2}^{(1)}$$\nu_{\uparrow}(X_{2}^{(1)})$$\psi^{(1)}$$X_{1}^{(0)}$$\nu_{\downarrow}(X_{1}^{(0)})$$X_{2}^{(0)}$$\psi^{(1)}$$X_{3}^{(0)}$$X_{4}^{(0)}$$\nu_{\uparrow}(X_{4}^{(0)})$

(Down)

Figure 4: Illustration of the flow of messages in the Belief Propagation
algorithm for the case $s=2$, $L=2$ of the Random Hierarchy Model. The factor
nodes (squares) represent the rules that connect the variables at different
levels of the hierarchy. The downward process is represented only for the
leftmost branch.

In this section, we consider the denoising diffusion process for the RHM.
Given a noisy observation $X^{(0)}=x(t)$ of the input variables at time $t$,
we compute $p(x(0)|x(t))$ exactly, obtaining full control of the statistics of
the ba diffusion process from time $t$ to time $0$. In particular, given the
tree structure of the model, we can compute the marginal probability of the
values of all latent variables conditioned on $x(t)$ by using a message-
passing algorithm. Therefore, we obtain the probability that a latent variable
at level $\ell$ has changed when performing the forward-backward diffusion
process for a duration $t$, a central quantity to interpret Figure 2.

For computing the marginal distributions, we use Belief Propagation (BP)
Mossel (2001); Mezard & Montanari (2009), which gives exact results for a tree
graph such as the Random Hierarchy Model. In this case, the leaves of the tree
correspond to the input variables at the bottom layer, and the root
corresponds to the class variable at the top of the hierarchy. Each rule
connecting variables at different levels corresponds to a factor node, as
shown in Figure 4.

The forward process adds noise to the variables in the input nodes. Each of
these nodes sends its belief on its value at $t=0$ to its parent latent node.
These beliefs, or messages, represent probabilistic estimates of the state of
the sender node. Each latent node receives messages from all its children,
updates its belief about its state, and sends its upward message to its parent
node. This process is repeated iteratively until the root of the tree.
Subsequently, starting from the root, each node sends a downward message to
its children. Finally, the product of the upward and downward beliefs received
at a given node represents the marginal probabilities of its state conditioned
on the noisy observation. Hence, we can use these conditional marginals to
compute the mean values of the variables at all levels of the hierarchy. We
assume that the production rules of the model are known by the inference
algorithm, which corresponds to the optimal denoising process.

The input variables $X^{(0)}$, in their one-hot-encoding representation,
undergo the forward diffusion process of Equation 2, which can be defined in
continuous time and constant $\beta_{t}$ by redefining
$\overline{\alpha_{t}}=e^{-2t}$ and taking the limit $T\rightarrow\infty$ Song
et al. (2020).

The denoising is made in two steps: the initialisation of the messages at the
leaves and the BP iteration.

##### Initialisation of the upward messages

In its one-hot-encoding representation, $X^{(0)}_{i}$ is a $v$-dimensional
vector: taking the symbol $a_{\gamma}\in\\{a_{1},\dots,a_{v}\\}=\mathcal{A}$
corresponds to $X^{(0)}_{i}=e_{\gamma}$, with $e_{\gamma}$ a canonical basis
vector. Its continuous diffusion process takes place in $\mathbb{R}^{v}$:
given the value $X^{(0)}_{i}=x_{i}(t)$, we can compute the probability of its
starting value $p(x_{i}(0)|x_{i}(t))$ using Bayes formula. As derived in
Appendix A, we obtain

$\displaystyle
p(x_{i}(0)=e_{\gamma}|x_{i}(t))=\frac{1}{Z}\,e^{x_{i,\gamma}(t)/\Delta_{t}},$
(4)

with $\Delta_{t}=(1-\overline{\alpha_{t}})/\sqrt{\overline{\alpha_{t}}}$ and
$Z=\sum_{\mu=1}^{v}e^{x_{i,\mu}(t)/\Delta_{t}}$. This computation is performed
independently for each input variable $i$, and therefore does not take into
account the spatial correlations given by the generative model. The
probabilities of Equation 4 are used to initialise the BP upward messages
${\nu_{\uparrow}^{(0)}}=p(x_{i}(0)|x_{i}(t))$ at the input variables.

##### BP iteration

Let ${{\psi}^{(\ell)}}$ be any factor node connecting an $s$-tuple of low-
level variables at layer $\ell-1$, $\\{X_{i}^{(\ell-1)}\\}_{i\in[s]}$, to a
high-level variable $X_{1}^{(\ell)}$ at layer $\ell$. Without loss of
generality, to lighten the notation, we rename the variables as
$Y=X_{1}^{(\ell)}$, taking values $y\in\mathcal{A}$, and
$X_{i}=X_{i}^{(\ell-1)}$, each taking values $x_{i}\in\mathcal{A}$. For each
possible association $y\rightarrow x_{1},\dots,x_{s}$, the factor node
${{\psi}^{(\ell)}}(y,x_{1},...,x_{s})$ takes values

${{\psi}^{(\ell)}}(y,x_{1},...,x_{s})=\begin{cases}1,\;\text{if
}y\rightarrow(x_{1},...,x_{s})\text{ is rule at layer }\ell\\\
0,\;\text{otherwise}.\end{cases}$

The BP upward and downward iterations for the (unnormalised) upward and
downward messages respectively read

${\tilde{\nu}_{\uparrow}^{(\ell+1)}}(y)=\sum_{x_{1},...,x_{s}\in\mathcal{A}^{\otimes
s}}{{\psi}^{(\ell+1)}}(y,x_{1},...,x_{s})\prod_{i=1}^{s}{\nu_{\uparrow}^{(\ell)}}(x_{i}),$
$\displaystyle{\tilde{\nu}_{\downarrow}^{(\ell)}}(x_{1})$
$\displaystyle=\sum_{\begin{subarray}{c}x_{2},...,x_{s}\in\mathcal{A}^{\otimes(s-1)}\\\
y\in\mathcal{A}\end{subarray}}{{\psi}^{(\ell+1)}}(y,x_{1},...,x_{s})$
$\displaystyle\phantom{=}\times{\nu_{\downarrow}^{(\ell+1)}}(y)\prod_{i=2}^{s}{\nu_{\uparrow}^{(\ell)}}(x_{i}),$
(5)

where
$\nu_{\rho}^{(\ell)}(x)=\frac{\tilde{\nu}_{\rho}^{(\ell)}(x)}{\sum_{x^{\prime}}\tilde{\nu}_{\rho}^{(\ell)}(x^{\prime})}$,
$\rho\in\\{\uparrow,\downarrow\\}$. The downward iteration, reported for
$x_{1}$, can be trivially extended to the other variables $x_{i}$ by permuting
the position indices. The values of ${\nu_{\uparrow}^{(0)}}(x_{i})$ and
${\nu_{\downarrow}^{(L)}}(y)$ are set by the initial conditions. In
particular, we initialise ${\nu_{\uparrow}^{(0)}}(x_{i})$ as described in the
previous paragraph and ${\nu_{\downarrow}^{(L)}}(y)=1/v$, which corresponds to
a uniform prior over the possible classes $\mathcal{C}$.222This assumption
corresponds to unconditioned diffusion, where the DDPM is not biased towards
any specific class.

Figure 5: Probability that the latent has not changed in the denoising
process, corresponding to the largest marginal probability computed by BP,
averaged for each layer, for varying inversion times of the diffusion process
$t$. Data for the RHM with $v=32$, $m=8$, $s=2$, $L=10$. Each level of the
tree, indicated in the legend, is represented with a different colour. We
observe the same behaviour of the curves for ImageNet data in Figure 2: the
probability of the correct class has a sharp transition at characteristic time
scale, while the probabilities corresponding to latent variables in the lower
levels change smoothly.

Exact solution We solve the BP upward and backward iterations numerically. In
Figure 5, we show the probability corresponding to the correct symbol for each
node of the tree. Remarkably, we note that (i) the probability for the correct
class at layer $L$ displays a transition at a characteristic time which
becomes sharper for increasing $L$, and (ii) the messages for the correct
input variables and the correct latent variables at low levels of the tree
change smoothly. In particular, the curves for messages at layer $L$ and
layers $\ell<L$ invert their order at the transition, as in our observations
on DDPMs and ImageNet data in Figure 2. This transition is our key finding,
which we explain below.

##  5 Mean-field theory of denoising diffusion

In this section, we make a simplifying assumption for the initial noise acting
on the input and adopt a mean-field approximation to justify the existence of
a phase transition. Remarkably, this approximation turns out to be of
excellent quality for describing the diffusion dynamics. Specifically,
consider a reference configuration at the leaves variables
$X_{i}^{(0)}={\overline{x}}_{i}$ that we would like to reconstruct, given a
noisy observation of it. We assume that for each leaf variable, the noise is
uniformly spread among the other symbols.333This is a mild approximation, as
documented in Appendix C In other words, our belief in the correct sequence is
corrupted by $\epsilon\in[0,1]$:

$\displaystyle\begin{cases}X_{i}^{(0)}={\overline{x}}_{i}\quad\text{with
belief }1-\epsilon,\\\ X_{i}^{(0)}\text{ uniform over alphabet with belief
}\epsilon.\end{cases}$ (6)

Hence, the initialisation condition of the upward BP messages at a leaf node
$X_{i}^{(0)}$ becomes

$\displaystyle\begin{cases}{\nu_{\uparrow}^{(0)}}\left({\overline{x}}_{i}\right)\qquad=1-\epsilon+\epsilon/v,\\\
{\nu_{\uparrow}^{(0)}}\left(x_{i}\neq{\overline{x}}_{i}\right)=\epsilon/v,\end{cases}$
(7)

where $v$ is the alphabet cardinality.

Given these initial conditions and since the production rules are known, if
$\epsilon=0$ – i.e., in the noiseless case – BP can reconstruct all the values
of the latent variables exactly. Conversely, if $\epsilon=1$ – i.e., when the
input is completely corrupted and the belief on the leaves variables is
uniform – the reconstruction is impossible. In general, for a value of
$\epsilon$, one is interested in computing the probability of recovering the
latent structure of the tree at each layer $\ell$ and, as $L\to\infty$, to
decide whether the probability of recovering the correct class of the input
remains larger than $1/v$.

$\displaystyle{0.25}$$\displaystyle{0.50}$$\displaystyle{0.75}$$\displaystyle{1.00}$$\displaystyle{1.25}$$\displaystyle
s\
f$$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle\epsilon$InferenceNo
Inferenceexperimentstheory

Figure 6: Phase diagram for inferring the class node using the upward
iteration of BP. When $sf<1$, BP can infer the class if
$\epsilon<\epsilon^{*}(sf)$. This transition is very well predicted by our
theory. The inference region in the figure corresponds to the phase wherein
the probability of the correct class is larger than the initialisation belief
in the correct values of the leaves, that is $1-\epsilon+\frac{\epsilon}{v}$.
Experimental data are for a single realisation of the RHM with $v=32$, $s=2$,
$L=10$.

##### Upward process

We begin by studying the upward process from the leaves. Consider a true input
tuple ${\overline{x}}_{1},\dots,{\overline{x}}_{s}$ which is associated with
the higher-level feature ${\overline{y}}$. Given the randomness of the
production rules, the messages are random variables depending on the specific
realisation of the rules. We adopt a mean-field or annealed approximation that
neglects the fluctuations coming from the random choice of rules.
Specifically, we approximate the upward message by the average upward message
exiting the corresponding factor node
$\langle{\nu_{\uparrow}^{(1)}}(y)\rangle_{\psi}$ over the possible
realisations of ${\psi}$. In Appendix B, we show that
$\langle{\nu_{\uparrow}^{(1)}}(y)\rangle_{\psi}$ can take only two values: one
for $y={\overline{y}}$ and one for $y\neq{\overline{y}}$, as expected by
symmetry considerations. Therefore, mean messages have the same structure as
Equation 7 and we can define a new $\epsilon^{\prime}$. Introducing the
probability of reconstructions $p=1-\epsilon+\epsilon/v$ and
$p^{\prime}=1-\epsilon^{\prime}+\epsilon^{\prime}/v$, we have

$p^{\prime}=\frac{p^{s}+f\frac{m-1}{mv-1}\left(1-p^{s}\right)}{p^{s}+f\left(1-p^{s}\right)}=F(p).$
(8)

Iterating this procedure across all the levels of the tree, we can compute the
probability of recovering the correct class of the input. In particular, for
large $L$, we are interested in studying the fixed points $p^{*}=F(p^{*})$ of
the iteration map in Equation 8. As derived in Section B.1.1, when $sf>1$,
this map has a repulsive fixed point $p^{*}=1$, which corresponds to
$\epsilon=0$, and an attractive fixed point $p^{*}=\nicefrac{{1}}{{v}}$,
corresponding to $\epsilon=1$. Thus, in this regime, inferring the class from
the noisy observation of the input is impossible. In contrast, when $sf<1$,
$p^{*}=1$ and $p^{*}=\nicefrac{{1}}{{v}}$ are both attractive fixed points,
and a new repulsive fixed point $\nicefrac{{1}}{{v}}<p^{*}<1$ separating the
other two emerges. Therefore, in this second regime, there is a phase
transition between a phase in which the class can be recovered and a phase in
which it cannot. These theoretical predictions are numerically confirmed in
the phase diagram in Figure 6.

Physically, $sf<1$ corresponds to a regime in which errors at lower levels of
the tree do not propagate: they can be corrected using information coming from
neighbouring nodes, thanks to the fact that only a small fraction of the
strings is consistent with the production rules of the generative model.
Conversely, when $sf>1$, even small corruptions propagate through the entire
tree up to the root node and BP cannot infer the class correctly.

##### Downward process

The same calculation can be repeated for the downward process, with the
additional difficulty that the downward iteration mixes upward and downward
messages. We refer the reader to Appendix B for the theoretical treatment.

##### Probabilities of reconstruction

Combining the mean upward and downward messages, we obtain a theoretical
prediction for the probabilities of reconstructing the correct values of the
variables at each layer. We compare our theoretical predictions with numerical
experiments in Figure 7. In these experiments, BP equations are solved exactly
for a given RHM starting with the initialisation of Equation 7. Our theory
perfectly captures the probability of reconstruction for the input nodes and
the class. Moreover, Figure 14 in Appendix B shows that our theory predicts
the probabilities of reconstruction of latent nodes at all layers.

##### Experiment on CNN’s activations

Similarly to our experiment on the ConvNeXt in Section 2, we investigate how
the hidden representation of a model trained to classify the RHM changes when
its input is denoised starting from a corruption noise $\epsilon$. We consider
an instantiation of the RHM with $L=7$, $s=2$, $v=16$, and $m=4$. First, we
train a convolutional neural network with $L=7$ layers, matching the tree
structure of the model, with $n=300k$ training examples up to interpolation.
The resulting architecture has 99.2% test accuracy. To sample new data from
noisy observations of held-out data, we start by sampling the root using the
marginal probability computed with BP. Then, we update the beliefs and the
marginals conditioning on the sampled class, and sample one latent variable at
layer $L-1$. We iterate this procedure node-by-node, descending the tree until
we obtain a sampled configuration at the bottom layer Mezard & Montanari
(2009). For each corrupting noise $\epsilon$ and each layer of the CNN, we
compute the cosine similarity between post-activations for the initial and
generated configurations. Figure 8 shows the obtained curves. Remarkably, we
observe the same qualitative behaviour as in Figure 7, ultimately explaining
the empirical observation of Figure 2.

Figure 7: Probability that the latent has not changed in the denoising
process, corresponding to the largest marginal probability computed by BP, for
varying $\epsilon$. Data for the RHM with $v=32$, $m=8$, $s=2$, $L=10$. Each
level of the tree, indicated in the legend, is represented with a different
colour. The black dashed lines are our mean-field theoretical predictions,
which show excellent agreement with the experiments. In particular, the
inversion between the curves for the top and bottom levels at the phase
transition can be observed. Figure 8: Cosine similarity between the post-
activations following every layer of a deep CNN trained on the RHM ($v=16$,
$m=4$, $s=2$, $L=7$) for the starting and sampled data. Each layer of the
architecture, indicated in the legend, is represented with a different colour.
The curves showcase the same inversion predicted by our theory (cf. Figure 7).

##  6 Conclusions

We have argued that reversing time in denoising diffusion models opens a
window on the compositional nature of data. For synthetic hierarchical
generative data models, low-level features can already change at small times,
but the class remains most often the same. At larger times, a phase transition
is found where the probability of remaining in the same class suddenly drops
to random chance. Yet, low-level features identical to those of the initial
sample can persist and compose the new sample. Strikingly, the same results
are found for ImageNet. These observations may thus be more general, and
testing them in text datasets, using diffusion language models He et al.
(2022), is an interesting venue for future work.

The interplay between the hierarchy in feature space and in time revealed here
may help understand the puzzling success of diffusion models, including the
number of data needed to train such methods, or why they can generalise and
not simply memorise the empirical distribution on which they were trained
Somepalli et al. (2022); Carlini et al. (2023); Yoon et al. (2023). More
generally, our results put forward hierarchical generative models as tools to
understand open questions for other methods, ranging from the emergence of new
skills by the composition of more elementary ones in self-supervised learning
to that of transferable representations in supervised or contrastive learning.

## Impact statement

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none of
which we feel must be specifically highlighted here.

## Acknowledgements

We thank Francesco Cagnetta and Umberto Maria Tomasini for the helpful
discussions. This work was supported by a grant from the Simons Foundation
(#454953 Matthieu Wyart).

## References

  * Allen-Zhu & Li (2020) Allen-Zhu, Z. and Li, Y.  Backward feature correction: How deep learning performs deep learning.  _arXiv preprint arXiv:2001.04413_ , 2020.
  * Allen-Zhu & Li (2023) Allen-Zhu, Z. and Li, Y.  Physics of language models: Part 1, context-free grammar.  _arXiv preprint arXiv:2305.13673_ , 2023.
  * Ambrogioni (2023) Ambrogioni, L.  The statistical thermodynamics of generative diffusion models.  _arXiv preprint arXiv:2310.17467_ , 2023.
  * Bach (2017) Bach, F.  Breaking the curse of dimensionality with convex neural networks.  _The Journal of Machine Learning Research_ , 18(1):629–681, 2017.
  * Betker et al. (2023) Betker, J., Goh, G., Jing, L., Brooks, T., Wang, J., Li, L., Ouyang, L., Zhuang, J., Lee, J., Guo, Y., et al.  Improving image generation with better captions.  _Computer Science. https://cdn. openai. com/papers/dall-e-3. pdf_ , 2(3), 2023.
  * Biroli & Mézard (2023) Biroli, G. and Mézard, M.  Generative diffusion in very large dimensions.  _arXiv preprint arXiv:2306.03518_ , 2023.
  * Block et al. (2020) Block, A., Mroueh, Y., and Rakhlin, A.  Generative modeling with denoising auto-encoders and langevin sampling.  _arXiv preprint arXiv:2002.00107_ , 2020.
  * Cagnetta et al. (2023) Cagnetta, F., Petrini, L., Tomasini, U. M., Favero, A., and Wyart, M.  How deep neural networks learn compositional data: The random hierarchy model.  _arXiv preprint arXiv:2307.02129_ , 2023.
  * Carlini et al. (2023) Carlini, N., Hayes, J., Nasr, M., Jagielski, M., Sehwag, V., Tramer, F., Balle, B., Ippolito, D., and Wallace, E.  Extracting training data from diffusion models.  In _32nd USENIX Security Symposium (USENIX Security 23)_ , pp. 5253–5270, 2023.
  * Chen et al. (2023) Chen, M., Huang, K., Zhao, T., and Wang, M.  Score approximation, estimation and distribution recovery of diffusion models on low-dimensional data.  _arXiv preprint arXiv:2302.07194_ , 2023.
  * Cui et al. (2023) Cui, H., Krzakala, F., Vanden-Eijnden, E., and Zdeborová, L.  Analysis of learning a flow-based generative model from limited sample complexity.  _arXiv preprint arXiv:2310.03575_ , 2023.
  * De Bortoli (2022) De Bortoli, V.  Convergence of denoising diffusion models under the manifold hypothesis.  _arXiv preprint arXiv:2208.05314_ , 2022.
  * DeGiuli (2019) DeGiuli, E.  Random language model.  _Physical Review Letters_ , 122(12):128301, 2019.
  * Dhariwal & Nichol (2021) Dhariwal, P. and Nichol, A.  Diffusion models beat gans on image synthesis.  _Advances in neural information processing systems_ , 34:8780–8794, 2021.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep Residual Learning for Image Recognition.  _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 770–778, June 2016.  doi: 10.1109/CVPR.2016.90.
  * He et al. (2022) He, Z., Sun, T., Wang, K., Huang, X., and Qiu, X.  Diffusionbert: Improving generative masked language models with diffusion models.  _arXiv preprint arXiv:2211.15029_ , 2022.
  * Ho et al. (2020) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  _Advances in neural information processing systems_ , 33:6840–6851, 2020.
  * Kadkhodaie et al. (2023) Kadkhodaie, Z., Guth, F., Mallat, S., and Simoncelli, E. P.  Learning multi-scale local conditional probability models of images.  _arXiv preprint arXiv:2303.02984_ , 2023.
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., and Hinton, G.  Deep learning.  _Nature_ , 521(7553):436, 2015.
  * Liu et al. (2022) Liu, Z., Mao, H., Wu, C.-Y., Feichtenhofer, C., Darrell, T., and Xie, S.  A convnet for the 2020s.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 11976–11986, 2022.
  * Luxburg & Bousquet (2004) Luxburg, U. v. and Bousquet, O.  Distance-based classification with lipschitz functions.  _The Journal of Machine Learning Research_ , 5(Jun):669–695, 2004.
  * Malach & Shalev-Shwartz (2018) Malach, E. and Shalev-Shwartz, S.  A provably correct algorithm for deep learning that actually works.  _arXiv preprint arXiv:1803.09522_ , 2018.
  * Malach & Shalev-Shwartz (2020) Malach, E. and Shalev-Shwartz, S.  The implications of local correlation on learning some deep functions.  _Advances in Neural Information Processing Systems_ , 33:1322–1332, 2020.
  * Mei & Wu (2023) Mei, S. and Wu, Y.  Deep networks as denoising algorithms: Sample-efficient learning of diffusion models in high-dimensional graphical models.  _arXiv preprint arXiv:2309.11420_ , 2023.
  * Mezard & Montanari (2009) Mezard, M. and Montanari, A.  _Information, physics, and computation_.  Oxford University Press, 2009.
  * Mossel (2001) Mossel, E.  Reconstruction on trees: beating the second eigenvalue.  _The Annals of Applied Probability_ , 11(1):285–300, 2001.
  * Mossel (2016) Mossel, E.  Deep learning and hierarchal generative models.  _arXiv preprint arXiv:1612.09057_ , 2016.
  * Nichol & Dhariwal (2021) Nichol, A. Q. and Dhariwal, P.  Improved denoising diffusion probabilistic models.  In _International Conference on Machine Learning_ , pp. 8162–8171. PMLR, 2021.
  * Okawa et al. (2023) Okawa, M., Lubana, E. S., Dick, R. P., and Tanaka, H.  Compositional abilities emerge multiplicatively: Exploring diffusion models on a synthetic task.  _arXiv preprint arXiv:2310.09336_ , 2023.
  * Oko et al. (2023) Oko, K., Akiyama, S., and Suzuki, T.  Diffusion models are minimax optimal distribution estimators.  _arXiv preprint arXiv:2303.01861_ , 2023.
  * Olah et al. (2020) Olah, C., Cammarata, N., Schubert, L., Goh, G., Petrov, M., and Carter, S.  Zoom in: An introduction to circuits.  _Distill_ , 5(3):e00024–001, 2020.
  * Poggio et al. (2017) Poggio, T., Mhaskar, H., Rosasco, L., Miranda, B., and Liao, Q.  Why and when can deep-but not shallow-networks avoid the curse of dimensionality: a review.  _International Journal of Automation and Computing_ , 14(5):503–519, 2017.
  * Raya & Ambrogioni (2024) Raya, G. and Ambrogioni, L.  Spontaneous symmetry breaking in generative diffusion models.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Rombach et al. (2022) Rombach, R., Blattmann, A., Lorenz, D., Esser, P., and Ommer, B.  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 10684–10695, 2022.
  * Rozenberg & Salomaa (1997) Rozenberg, G. and Salomaa, A.  _Handbook of Formal Languages_.  Springer, January 1997.  doi: 10.1007/978-3-642-59126-6.
  * Schmidt-Hieber (2020) Schmidt-Hieber, J.  Nonparametric regression using deep neural networks with relu activation function.  _The Annals of Statistics_ , 48(4):1875–1897, 2020.
  * Shah et al. (2023) Shah, K., Chen, S., and Klivans, A.  Learning mixtures of gaussians using the ddpm objective.  _arXiv preprint arXiv:2307.01178_ , 2023.
  * Shalev-Shwartz et al. (2017) Shalev-Shwartz, S., Shamir, O., and Shammah, S.  Failures of gradient-based deep learning.  In _International Conference on Machine Learning_ , pp. 3067–3075. PMLR, 2017.
  * Sohl-Dickstein et al. (2015) Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _International conference on machine learning_ , pp. 2256–2265. PMLR, 2015.
  * Somepalli et al. (2022) Somepalli, G., Singla, V., Goldblum, M., Geiping, J., and Goldstein, T.  Diffusion art or digital forgery.  _Investigating Data Replication in Diffusion Models_ , 2022.
  * Song & Ermon (2019) Song, Y. and Ermon, S.  Generative modeling by estimating gradients of the data distribution.  _Advances in neural information processing systems_ , 32, 2019.
  * Song et al. (2020) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020.
  * TorchVision maintainers and contributors (2016) TorchVision maintainers and contributors.  Torchvision: Pytorch’s computer vision library.  https://github.com/pytorch/vision, 2016.
  * Yoon et al. (2023) Yoon, T., Choi, J. Y., Kwon, S., and Ryu, E. K.  Diffusion probabilistic models generalize when they fail to memorize.  In _ICML 2023 Workshop on Structured Probabilistic Inference & Generative Modeling_, 2023.
  * Yuan et al. (2023) Yuan, H., Huang, K., Ni, C., Chen, M., and Wang, M.  Reward-directed conditional diffusion: Provable distribution estimation and reward improvement.  _arXiv preprint arXiv:2307.07055_ , 2023.
  * Zeiler & Fergus (2014) Zeiler, M. D. and Fergus, R.  Visualizing and understanding convolutional networks.  In _Computer Vision – ECCV 2014_ , Lecture Notes in Computer Science, pp. 818–833, 2014.

##  Appendix A Belief Propagation initialization for the denoising of the RHM

As discussed in Section 4, we define the diffusion process for the input
variable $X_{i}^{(0)}$ in the space $\mathbb{R}^{v}$. In particular, its value
$x(t)$ at time $t$ is

$\displaystyle
x(t)=\sqrt{\overline{\alpha_{t}}}x(0)+\sqrt{1-\overline{\alpha_{t}}}\eta,$ (9)

with $\eta\sim\mathcal{N}(0,\mathbb{I}_{v})$ and $x(0)$ its starting value at
time $t$, which is a one-hot-encoding vector of the form $x(0)=e_{\mu}$. Given
the value $x(t)$, the conditional probabilities for the values of $x(0)$ are
given by Bayes rule

$\displaystyle
p\left(x(0)=e_{\mu}|x(t)\right)=\frac{p\left(x(t)|x(0)=e_{\mu}\right)\,p\left(x(0)=e_{\mu}\right)}{\sum_{\lambda}p\left(x(t)|x(0)=e_{\lambda}\right)\,p\left(x(0)=e_{\lambda}\right)}.$
(10)

The prior probabilities on $x(0)$ are taken to be uniform over the alphabet,
i.e., $p\left(x(0)=e_{\lambda}\right)=1/v$, $\forall\lambda$, while
$p\left(x(t)|x(0)=e_{\mu}\right)$ is given by the diffusion process of Eq. 9:

$\displaystyle\begin{split}p\left(x(t)|x(0)=e_{\mu}\right)&=C_{t}\exp\left[-\frac{1}{2(1-\overline{\alpha_{t}})}\sum_{\gamma}\left(x_{\gamma}(t)-\sqrt{\overline{\alpha_{t}}}e_{\mu}\right)^{2}\right]=\\\
&=C_{t}\exp\left[-\frac{\left(\|x(t)\|^{2}-\overline{\alpha_{t}}\right)}{2(1-\overline{\alpha_{t}})}\right]\exp\left[\frac{\sqrt{\overline{\alpha_{t}}}}{1-\overline{\alpha_{t}}}x_{\mu}(t)\right],\end{split}$
(11)

where $C_{t}$ is the normalization constant. Putting Eq. 11 into Eq. 10, we
obtain

$\displaystyle
p\left(x(0)=e_{\mu}|x(t)\right)=\frac{1}{Z}e^{\frac{\sqrt{\overline{\alpha_{t}}}}{1-\overline{\alpha_{t}}}x_{\mu}(t)},$
(12)

with
$Z=\sum_{\lambda=1}^{v}e^{\frac{\sqrt{\overline{\alpha_{t}}}}{1-\overline{\alpha_{t}}}x_{\lambda}(t)}$.

##  Appendix B Belief Propagation equations

Given a factor tree-graph, the Belief Propagation (BP) equations compute
iteratively the messages going from the variable nodes to the factor nodes and
vice-versa, starting from the initialization conditions at the leaves and root
of the tree-graph Mezard & Montanari (2009). For the generative model defined
in Sec. 3, the leaves correspond to the variables at the bottom layer while
the root is the class variable at the top of the hierarchy. Each rule,
connecting variables at different layers, corresponds to a factor node. The BP
messages that flow from the variable nodes to the factor nodes, therefore,
correspond to upward messages, while those going from factor nodes to
variables correspond to downward messages (Fig. 9).

To each variable node $X^{(\ell)}_{i}$ at level $\ell$, we associate the
upward messages ${\nu_{\uparrow}^{(\ell)}}$ and downward messages
${\nu_{\downarrow}^{(l)}}$, one for each possible value of the alphabet it can
take. To simplify the notation, here we consider how messages propagate from
one level to the other and we call $Y$ the variable corresponding to the
higher level and $X_{i=1,...,s}$ the lower level ones connected to it. The
factor node connecting them is such that, for each possible association
$y\rightarrow x_{1},\dots,x_{s}$, it takes values

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

The BP upward and downward iterations are defined as follows.

  * •

Upward iteration:

${\tilde{\nu}_{\uparrow}^{(\ell+1)}}(y)=\sum_{x_{1},...,x_{s}\in\mathcal{A}^{\otimes
s}}{{\psi}^{(\ell+1)}}(y,x_{1},...,x_{s})\prod_{i=1}^{s}{\nu_{\uparrow}^{(\ell)}}(x_{i})\
,$ (13)
${\nu_{\uparrow}^{(\ell)}}(y)=\frac{{\tilde{\nu}_{\uparrow}^{(\ell)}}(y)}{\sum_{y^{\prime}}{\tilde{\nu}_{\uparrow}^{(\ell)}}(y^{\prime})}.$
(14)

  * •

Downward iteration:

${\tilde{\nu}_{\downarrow}^{(\ell)}}(x_{1})=\sum_{\begin{subarray}{c}x_{2},...,x_{s}\in\mathcal{A}^{\otimes(s-1)}\\\
y\in\mathcal{A}\end{subarray}}{{\psi}^{(\ell+1)}}(y,x_{1},...,x_{s})\
{\nu_{\downarrow}^{(\ell+1)}}(y)\prod_{i=2}^{s}{\nu_{\uparrow}^{(\ell)}}(x_{i})\
$ (15)
${\nu_{\downarrow}^{(l)}}(x)=\frac{{\tilde{\nu}_{\downarrow}^{(\ell)}}(x)}{\sum_{x^{\prime}}{\tilde{\nu}_{\downarrow}^{(\ell)}}(x^{\prime})}.$
(16)

${\nu_{\uparrow}^{(\ell)}}(y)$ and ${\nu_{\downarrow}^{(\ell)}}(x)$ are
fluctuating quantities that depend on the position of the node.

$Y$$\psi$$\nu^{\uparrow}(Y)$$X_{1}$$\nu^{\uparrow}(X_{1})$$X_{2}$$X_{s}$$\nu^{\uparrow}(X_{s})$

(a)

$Y$$\psi$$\nu^{\downarrow}(Y)$$X_{1}$$\nu^{\downarrow}(X_{1})$$X_{2}$$X_{s}$$\nu^{\uparrow}(X_{s})$

(b)

Figure 9: Factor tree-graph connecting the higher-level feature $Y$ to the
lower-level features $X_{i=1,...,s}$ according to the rules ${\psi}$. The
upward messages $\nu^{\uparrow}(y)$ are computed from the upward messages
$\nu^{\uparrow}(x_{i})$ coming from the nodes $X_{i}$, connected to $Y$
through the rule ${\psi}$ (panel (a)). The downward messages
$\nu^{\downarrow}(x_{1})$, instead, are computed from both the downward
messages $\nu^{\downarrow}(y)$ coming from $Y$ and the upward messages
$\nu^{\uparrow}(x_{i})$ coming from the nodes $X_{i=2,...s}$, connected to
$X_{1}$ through the rule ${\psi}$ (panel (b)).

###  B.1 $\epsilon$-process

In this process, we consider a reference configuration at the leaves variables
$X_{i}^{(0)}={\overline{x}}_{i}$ that we would like to reconstruct, given a
noisy observation of it. As a result of this noise addition, our belief in the
correct sequence is corrupted by $\epsilon\in[0,1]$:

$\displaystyle\begin{cases}X_{i}^{(0)}={\overline{x}}_{i}\quad\text{with
belief }1-\epsilon\\\ X_{i}^{(0)}\text{ uniform over alphabet with belief
}\epsilon.\\\ \end{cases}$ (17)

Therefore, the initialization condition of the upward BP messages at a leaf
node $X_{i}^{(0)}$ is

$\displaystyle\begin{cases}\nu_{\uparrow}\left({\overline{x}}_{i}\right)\qquad=1-\epsilon+\epsilon/v,\\\
\nu_{\uparrow}\left(x_{i}\neq{\overline{x}}_{i}\right)=\epsilon/v,\\\
\end{cases}$ (18)

where $v$ is the corresponding alphabet size.
The initialization condition at the root node $X^{(L)}$, that corresponds to
the messages ${\nu_{\downarrow}^{(L)}}$ for that node, is uniform over the
alphabet $\mathcal{A}$, so that the algorithm has no bias on any specific
class.

####  B.1.1 Upward iteration

We consider the upward iteration when going from the bottom layer to the one
above it. Let $X_{1},\dots,X_{s}$ denote a tuple at the bottom level which is
associated with the reference values ${\overline{x}}_{i}$. This tuple is
connected to the higher level variable $Y$ via a set of rules ${\psi}$ (Fig.
9). According to ${\psi}$, the association from the high-level feature to the
reference low-level sequence ${\overline{x}}_{1},\dots,{\overline{x}}_{s}$ is
given by

${\overline{y}}\rightarrow{\overline{x}}_{1},\dots,{\overline{x}}_{s}.$

We call $\Delta_{\mathbf{w},\mathbf{z}}$ the Hamming distance between two
sequences $\mathbf{w}=[w_{1},\dots,w_{s}]$, $\mathbf{z}=[z_{1},\dots,z_{s}]$
of length $s$.
From Eq. 18, at the bottom layer, the belief in a sequence
${\mathbf{x}}=[x_{1},\dots,x_{s}]$ with
$\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}=k\in\\{0,...,s\\}$ from
${\overline{{\mathbf{x}}}}=[{\overline{x}}_{1},\dots,{\overline{x}}_{s}]$ is

$\mathcal{B}(k)=\left(\frac{\epsilon}{v}\right)^{k}\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s-k}$
(19)

The non-normalized upward messages for the variable $Y$ are given by:

$\tilde{\nu}_{\uparrow}(y)=\sum_{x_{1},...,x_{s}}{\psi}(y,x_{1},...,x_{s})\prod_{i=1}^{s}\nu_{\uparrow}(x_{i})=\sum_{{\mathbf{x}}\in\mathcal{S}}{\psi}(y,{\mathbf{x}})\mathcal{B}\left(\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}\right),$
(20)

where we are using the short-hand notation
${\psi}(y,x_{1},...,x_{s})={\psi}(y,{\mathbf{x}})$ and we have restricted the
sum over the set $\mathcal{S}$ of sequences ${\mathbf{x}}$ that appear in the
possible rules $y\rightarrow x_{1},\dots,x_{s}$. In fact, if
${\mathbf{x}}\notin\mathcal{S}$, then ${\psi}(y,{\mathbf{x}})=0$.

For ${\mathbf{x}}\in\mathcal{S}$, the factor ${\psi}(y,{\mathbf{x}})$ is such
that:

  * •

if $\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}=0$:

$\displaystyle\begin{cases}&{\psi}({\overline{y}},{\overline{x}}_{1},...,{\overline{x}}_{s})=1,\\\
&{\psi}(y,{\overline{x}}_{1},...,{\overline{x}}_{s})=0,\quad
y\neq{\overline{y}}.\\\ \end{cases}$ (21)

  * •

if $\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}>0$:

$\displaystyle\begin{cases}&{\psi}(\tilde{y},x_{1},...,x_{s})=1,\quad\text{for
some $\tilde{y}$ independent of ${\overline{y}}$}\\\
&{\psi}(y,x_{1},...,x_{s})=0,\quad y\neq\tilde{y}.\\\ \end{cases}$ (22)

We can decompose Eq. 20 as

$\tilde{\nu}_{\uparrow}(y)=\delta_{y,{\overline{y}}}\mathcal{B}(0)+\sum_{k=1}^{s}\mathcal{B}\left(k\right)\left[\sum_{\begin{subarray}{c}{\mathbf{x}}\in\mathcal{S}\\\
\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}=k\end{subarray}}\psi(y,{\mathbf{x}})\right].$
(23)

##### Annealed average

$\psi$ is a random quantity and we want to compute the average message
$\langle\tilde{\nu}_{\uparrow}(y)\rangle_{\psi}$ over the possible
realizations of ${\psi}$. We can decompose the selection of the rules in two
steps: sampling the set of $mv-1$ sequences
$\\{{\mathbf{x}},{\mathbf{x}}\neq{\overline{{\mathbf{x}}}}\\}$ and then
associating the $v$ higher-level features $y$ to them. Therefore, for a
generic quantity $A$, we indicate the average over the rules realization
$\langle A\rangle_{\psi}$ as $\langle\langle
A\rangle_{\\{y\\}\leftarrow\\{{\mathbf{x}}\\}}\rangle_{\mathcal{S}}$, where
$\langle\dots\rangle_{\mathcal{S}}$ is the average over the sequence sampling
and $\langle\dots\rangle_{\\{y\\}\leftarrow\\{{\mathbf{x}}\\}}$ is the average
over the $y\leftarrow{\mathbf{x}}$ associations:

$\langle\tilde{\nu}_{\uparrow}(y)\rangle_{\psi}=\delta_{y,{\overline{y}}}\mathcal{B}(0)+\sum_{k=1}^{s}\mathcal{B}\left(k\right)\langle\langle\sum_{\begin{subarray}{c}{\mathbf{x}}\in\mathcal{S}\\\
\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}=k\end{subarray}}\psi(y,{\mathbf{x}})\rangle_{\\{y\\}\leftarrow\\{{\mathbf{x}}\\}}\rangle_{\mathcal{S}}$
(24)

Since for each sequence ${\mathbf{x}}\neq{\overline{{\mathbf{x}}}}$ the
association $y\leftarrow{\mathbf{x}}$ is done randomly, independently of
$\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}$, then from Eq. 22, we have
$\langle\psi(y,{\mathbf{x}})\rangle_{\\{y\\}\leftarrow\\{{\mathbf{x}}\\}}\simeq
1/v$. More precisely, since we have associated the reference sequence
${\overline{{\mathbf{x}}}}$ to ${\overline{y}}$:

$\overline{\psi_{y}}=\langle\psi(y,{\mathbf{x}})\rangle_{\\{y\\}\leftarrow\\{{\mathbf{x}}\\}}=\frac{m-1}{mv-1}\delta_{y,{\overline{y}}}+\frac{m}{mv-1}\left(1-\delta_{y,{\overline{y}}}\right).$
(25)

Therefore:

$\displaystyle\langle\tilde{\nu}_{\uparrow}(y)\rangle_{\psi}$
$\displaystyle=\delta_{y,{\overline{y}}}\mathcal{B}(0)+\overline{\psi_{y}}\sum_{k=1}^{s}\mathcal{B}\left(k\right)\langle\sum_{\begin{subarray}{c}{\mathbf{x}}\in\mathcal{S}\\\
\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}=k\end{subarray}}1\rangle_{\mathcal{S}}=$
(26)
$\displaystyle=\delta_{y,{\overline{y}}}\mathcal{B}(0)+\overline{\psi_{y}}\sum_{k=1}^{s}\mathcal{B}\left(k\right)\langle
n_{k}\rangle_{\mathcal{S}},$

where $n_{k}$ is the number of sequences ${\mathbf{x}}\in\mathcal{S}$ having
Hamming distance $\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}=k$ from
${\overline{{\mathbf{x}}}}$. Since the sequences are sampled randomly, the
numbers $n_{1},...,n_{s}$ are distributed according to a multivariate hyper-
geometric distribution,

$P\left(n_{1},...,n_{s}\right)=\frac{\prod_{k=1}^{s}\binom{\binom{s}{k}(v-1)^{k}}{n_{k}}}{\binom{v^{s}-1}{mv-1}},$
(27)

which gives the averages

$\langle
n_{k}\rangle_{\mathcal{S}}=\frac{mv-1}{v^{s}-1}\binom{s}{k}(v-1)^{k}=f\binom{s}{k}(v-1)^{k},$
(28)

with

$f=\frac{mv-1}{v^{s}-1}.$ (29)

Therefore:

$\displaystyle\langle\tilde{\nu}_{\uparrow}(y)\rangle_{\psi}=\delta_{y,{\overline{y}}}\mathcal{B}(0)+f\
\overline{\psi_{y}}\sum_{k=1}^{s}\mathcal{B}\left(k\right)\binom{s}{k}(v-1)^{k}.$
(30)

From the beliefs Eq. 19, we see that

$\sum_{k=1}^{s}\mathcal{B}\left(k\right)\binom{s}{k}(v-1)^{k}=\sum_{k=1}^{s}\binom{s}{k}(v-1)^{k}\left(\frac{\epsilon}{v}\right)^{k}\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s-k}=\left[1-\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s}\right]=1-\mathcal{B}(0),$
(31)

which gives

$\begin{aligned}
\langle\tilde{\nu}_{\uparrow}(y)\rangle_{\psi}=\delta_{y,{\overline{y}}}\mathcal{B}(0)+f\
\overline{\psi_{y}}\left[1-\mathcal{B}(0)\right]\end{aligned}.$ (32)

The normalization constant is:

$\langle
Z_{\uparrow}\rangle_{\psi}=\sum_{y}\langle\tilde{\nu}_{\uparrow}(y)\rangle_{\psi}=\mathcal{B}(0)+f\left[1-\mathcal{B}(0)\right].$
(33)

Finally, we obtain the average belief for $Y$

$\langle\nu_{\uparrow}(y)\rangle_{\psi}=\frac{\langle\tilde{\nu}_{\uparrow}(y)\rangle_{\psi}}{\langle
Z_{\uparrow}\rangle_{\psi}}=\frac{\delta_{y,{\overline{y}}}\mathcal{B}(0)+f\
\overline{\psi_{y}}\left[1-\mathcal{B}(0)\right]}{\mathcal{B}(0)+f\left[1-\mathcal{B}(0)\right]}$
(34)

We have that:

  * •

for $y={\overline{y}}$

$\langle\nu_{\uparrow}({\overline{y}})\rangle_{\psi}=\frac{\mathcal{B}(0)+f\frac{m-1}{mv-1}\left[1-\mathcal{B}(0)\right]}{\mathcal{B}(0)+f\left[1-\mathcal{B}(0)\right]},$
(35)

  * •

for $y\neq{\overline{y}}$

$\langle\nu_{\uparrow}(y)\rangle_{\psi}=f\frac{m}{mv-1}\frac{1-\mathcal{B}(0)}{\mathcal{B}(0)+f\left[1-\mathcal{B}(0)\right]}.$
(36)

##### Iterating over layers

The average messages in Eqs. 35, 36 are of two kinds: one for the reference
feature ${\overline{y}}$ and another for the others $y\neq{\overline{y}}$, and
they both depend on the previous beliefs through
$\mathcal{B}(0)=\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s}$. Therefore,
the average messages at the higher level have the same structure as those at
the lower level Eq. 18. We can then define a new $\epsilon^{\prime}$:

$1-\epsilon^{\prime}+\frac{\epsilon^{\prime}}{v}=\frac{\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s}+f\frac{m-1}{mv-1}\left[1-\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s}\right]}{\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s}+f\left[1-\left(1-\epsilon+\frac{\epsilon}{v}\right)^{s}\right]}$
(37)

or, equivalently,

$p^{\prime}=\frac{p^{s}+f\frac{m-1}{mv-1}\left(1-p^{s}\right)}{p^{s}+f\left(1-p^{s}\right)}=F(p)$
(38)

with $p^{\prime}=1-\epsilon^{\prime}+\epsilon^{\prime}/v$ and
$p=1-\epsilon+\epsilon/v$. The derivative of $F(p)$ with respect to $p$ is
given by11footnotemark: 1

$F^{\prime}(p)=\frac{m(v-1)}{mv-1}\
\frac{fsp^{s-1}}{\left[p^{s}+f(1-p^{s})\right]^{2}}=fs\frac{p^{s-1}}{\left[p^{s}+f(1-p^{s})\right]^{2}}+\mathcal{O}\left(\frac{1}{v}\right)$
(39)

We can extend the tree in Fig. 9 iteratively to higher levels of the
hierarchy, where the variables $Y$ take the place of the variables $X_{i}$ and
so on.

The iteration Eq. 38 has fixed points $p=1$ (corresponding to $\epsilon=0$) or
$p=1/v$ (corresponding to $\epsilon=1$). An additional repulsive fixed point
at finite $p$ appears if

$F^{\prime}(1)<1,$ (40)

that is

$\frac{m(v-1)}{mv-1}\ fs<1\quad\Rightarrow\quad fs<1+\frac{1-1/m}{v-1}$ (41)
$\boxed{\Rightarrow\quad fs<1+\mathcal{O}\left(\frac{1}{v}\right)}\,.$ (42)

Iteration map

$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle
p$$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle
F(p)$$\displaystyle s\ f=1.5$$\displaystyle s\ f=1.0$$\displaystyle s\
f=0.5$$\displaystyle p$

Class inference phase diagram

$\displaystyle{0.25}$$\displaystyle{0.50}$$\displaystyle{0.75}$$\displaystyle{1.00}$$\displaystyle{1.25}$$\displaystyle
s\
f$$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle\epsilon$InferenceNo
Inferenceexperimentstheory

Figure 10:  Left panel: Iteration map of Eq. 38. For $sf>1$, there are only
two fixed points $p=F(p)$ corresponding to $p_{1}=1$ (repulsive) and
$p_{2}=1/v$ (attractive). For $sf<1$, there is another (repulsive) fixed point
at finite $p^{*}$, $1/v<p^{*}<1$, separating $p_{1}$ and $p_{2}$ (both
attractive). Right panel: Phase diagram for inferring the class node using the
upwards iteration of BP. When $sf<1$, BP can infer the class of the datum if
$\epsilon<\epsilon^{*}(sf)$. This transition is well predicted by
$p^{*}=1-\epsilon^{*}+\epsilon^{*}/v$, with $p^{*}=F(p^{*})$ from Eq. 38.
Experimental data for $v=32$, $s=2$, $L=10$.

####  B.1.2 Downward iteration

We consider the downward process when we try reconstructing the reference
association
${\overline{y}}\rightarrow{\overline{x}}_{1},\dots,{\overline{x}}_{s}$ from
higher-level variable $Y$ to the corresponding lower-level tuple
$X_{1},\dots,X_{s}$, via the set of rules ${\psi}$. We consider the downward
message received by the variable $X_{1}$ (Fig. 9):

$\displaystyle\begin{split}\tilde{\nu}_{\downarrow}(x_{1})&=\sum_{\begin{subarray}{c}x_{2},...,x_{s}\in\mathcal{A}^{\otimes(s-1)}\\\
y\in\mathcal{A}\end{subarray}}{\psi}(y,x_{1},...,x_{s})\
\nu_{\downarrow}(y)\prod_{i=2}^{s}\nu_{\uparrow}(x_{i})=\\\
&=\delta_{x_{1},{\overline{x}}_{1}}\nu_{\downarrow}({\overline{y}})\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})+\sum_{y}\nu_{\downarrow}(y)\sum_{\begin{subarray}{c}x_{2},\dots,x_{s}\\\
{\mathbf{x}}\in\mathcal{S},\Delta_{{\mathbf{x}},{\overline{{\mathbf{x}}}}}>0\end{subarray}}\psi(y,{\mathbf{x}})\prod_{i=2}^{s}\nu_{\uparrow}(x_{i})\end{split}$
(43)

##### Annealed average

To study the iteration of Eq. 43 analytically, we compute the average message
$\langle\tilde{\nu}_{\downarrow}(x_{1})\rangle_{\psi}$ over the realizations
of the random rules ${\psi}$ as done in Sec. B.1.2 for the upward iteration.

We call $n_{x_{1}}$ the number of sequences, having $X_{1}=x_{1}$, that have
been sampled by the choice of the rules. The numbers $n_{x_{1}}$ are
distributed according to a multivariate hyper-geometric distribution,

$P\left(\\{n_{x_{1}}\\}_{x_{1}\in\mathcal{A}}\right)=\frac{\binom{v^{s-1}-1}{n_{{\overline{x}}_{1}}}\prod_{\tilde{x}_{1}\in\mathcal{A}\setminus{\overline{x}}_{1}}\binom{v^{s-1}}{n_{\tilde{x}_{1}}}}{\binom{v^{s}-1}{mv-1}},$
(44)

which gives averages

$\langle
n_{{\overline{x}}_{1}}\rangle=\frac{mv-1}{v^{s}-1}\left(v^{s-1}-1\right)=f\left(v^{s-1}-1\right),$
(45) $\langle
n_{\tilde{x}_{1}\neq{\overline{x}}_{1}}\rangle=\frac{mv-1}{v^{s}-1}v^{s-1}=fv^{s-1}.$
(46)

Averaging the downward messages over the choices of rules ${\psi}$, we obtain:

  * •

for $x_{1}\neq{\overline{x}}_{1}$:

$\displaystyle\begin{split}\langle\tilde{\nu}_{\downarrow}(x_{1})\rangle_{\psi}&=\sum_{y}\nu_{\downarrow}(y)\langle\sum_{\begin{subarray}{c}x_{2},\dots,x_{s}\\\
{\mathbf{x}}\in\mathcal{S}\end{subarray}}\langle\psi(y,{\mathbf{x}})\rangle_{\\{y\\}\leftarrow\\{{\mathbf{x}}\\}}\
\prod_{i=2}^{s}\nu_{\uparrow}(x_{i})\rangle_{\mathcal{S}}=\\\
&=\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}\langle\delta_{{\mathbf{x}}\in\mathcal{S}}\rangle_{\mathcal{S}}\prod_{i=2}^{s}\sum_{x_{i}}\nu_{\uparrow}(x_{i})=\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}f=\\\
&=\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}f,\end{split}$ (47)

where
$\langle\delta_{{\mathbf{x}}\in\mathcal{S}}\rangle_{\mathcal{S}}=\frac{\langle
n_{x_{1}\neq{\overline{x}}_{1}}\rangle}{v^{s-1}}$;

  * •

for $x_{1}={\overline{x}}_{1}$:

$\displaystyle\begin{split}\langle\tilde{\nu}_{\downarrow}({\overline{x}}_{1})\rangle_{\psi}&=\nu_{\downarrow}({\overline{y}})\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})+\sum_{y}\nu_{\downarrow}(y)\langle\sum_{\begin{subarray}{c}x_{2},\dots,x_{s}\\\
{\mathbf{x}}\in\mathcal{S},{\mathbf{x}}\neq{\overline{{\mathbf{x}}}}\end{subarray}}\langle\psi(y,{\mathbf{x}})\rangle_{\\{y\\}\leftarrow\\{{\mathbf{x}}\\}}\
\prod_{i=2}^{s}\nu_{\uparrow}(x_{i})\rangle_{\mathcal{S}}=\\\
&=\nu_{\downarrow}({\overline{y}})\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})+\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}\langle\delta_{{\mathbf{x}}\in\mathcal{S}}\rangle_{\mathcal{S}}\sum_{\begin{subarray}{c}x_{2},\dots,x_{s}\\\
{\mathbf{x}}\in\mathcal{S},{\mathbf{x}}\neq{\overline{{\mathbf{x}}}}\end{subarray}}\prod_{i=2}^{s}\nu_{\uparrow}(x_{i})=\\\
&=\nu_{\downarrow}({\overline{y}})\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})+\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}f\left[1-\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})\right]\end{split}$
(48)

where
$\langle\delta_{{\mathbf{x}}\in\mathcal{S}}\rangle_{\mathcal{S}}=\frac{\langle
n_{{\overline{x}}_{1}}\rangle}{v^{s-1}-1}$.

The normalization factor is

$\langle
Z_{\downarrow}\rangle_{\psi}=\sum_{x_{1}}\langle\tilde{\nu}_{\downarrow}(x_{1})\rangle_{\psi}=\nu_{\downarrow}({\overline{y}})\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})+f\
\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}\left[1-\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})\right]+(v-1)f\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}$
(49)

which gives the normalized average messages:

  * •

for $x_{1}={\overline{x}}_{1}$

$\langle\nu_{\downarrow}({\overline{x}}_{1})\rangle_{\psi}=\frac{\nu_{\downarrow}({\overline{y}})\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})+f\
\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}\left[1-\prod_{i=2}^{s}\nu_{\uparrow}({\overline{x}}_{i})\right]}{\langle
Z_{\downarrow}\rangle_{\psi}};$ (50)

  * •

for $x_{1}\neq{\overline{x}}_{1}$

$\langle\nu_{\downarrow}(x_{1})\rangle_{\psi}=f\frac{\frac{m-\nu_{\downarrow}({\overline{y}})}{mv-1}}{\langle
Z_{\downarrow}\rangle_{\psi}}.$ (51)

##### Iterating over layers

As for the upward process, the average messages for the downward process are
of two kinds, one for the correct value ${\overline{x}}_{1}$ and one for the
other values $x_{1}\neq{\overline{x}}_{1}$. To obtain a mean-field description
of the BP process, we combine the average downward messages with the average
upward ones, by substituting
$\nu_{\uparrow}({\overline{x}}_{i})\rightarrow\langle\nu_{\uparrow}({\overline{x}}_{i})\rangle$
in Eq. 50. We use the notation

$\displaystyle\langle{\nu_{\uparrow}^{(\ell)}}({\overline{x}}_{i})\rangle$
$\displaystyle=p_{\uparrow}^{(\ell)},$ (52)
$\displaystyle\langle{\nu_{\downarrow}^{(\ell)}}({\overline{x}}_{i})\rangle$
$\displaystyle=p_{\downarrow}^{(\ell)},$ (53)

where the upwards and downwards beliefs $p_{\uparrow}^{(\ell)}$,
$p_{\downarrow}^{(\ell)}$ in the correct value for the latent variable
$X^{(\ell)}_{i}$ depend only on the layer $\ell$ and not on the specific
position $i$ inside the layer. Putting together Eqs. 38 and 50, we obtain

$\displaystyle\begin{split}&p_{\uparrow}^{(\ell+1)}=F_{\uparrow}\left(p_{\uparrow}^{(\ell)}\right),\\\
&p_{\downarrow}^{(\ell)}=F_{\downarrow}\left(p_{\downarrow}^{(\ell+1)},p_{\uparrow}^{(\ell)}\right),\end{split}$
(54)

with

$\displaystyle
F_{\uparrow}(p)=\frac{p^{s}+f\frac{m-1}{mv-1}\left(1-p^{s}\right)}{p^{s}+f\left(1-p^{s}\right)},$
(55) $\displaystyle F_{\downarrow}(q,p)=\frac{q\
p^{s-1}+f\frac{m-q}{mv-1}\left(1-p^{s-1}\right)}{q\
p^{s-1}+f\frac{m-q}{mv-1}\left(1-p^{s-1}\right)+(v-1)f\frac{m-q}{mv-1}},$ (56)

and the initialization condition

$\displaystyle p_{\uparrow}^{(0)}=1-\epsilon+\epsilon/v,$ (57) $\displaystyle
p_{\downarrow}^{(L)}=1/v.$ (58)

From $p_{\uparrow}^{(\ell)}$ and $p_{\downarrow}^{(\ell)}$, at layer $\ell$,
the average marginal probability of the correct value $p^{(\ell)}$ is given by

$p^{(\ell)}=\frac{p_{\uparrow}^{(\ell)}p_{\downarrow}^{(\ell)}}{p_{\uparrow}^{(\ell)}p_{\downarrow}^{(\ell)}+\frac{(1-p_{\uparrow}^{(\ell)})(1-p_{\downarrow}^{(\ell)})}{v-1}}.$
(59)

####  B.1.3 Validity of the mean-field theory

Due to the randomness of the production rules, the messages
$\nu_{\uparrow}(x)$, $\nu_{\downarrow}(x)$ are random variables which depend
on the specific realization of the rules. Although their fluctuations are not
captured by the averages computed in Eq. 54, we observe that
$p_{\uparrow}^{(\ell)}$ and $p_{\downarrow}^{(\ell)}$ capture well the average
behaviour of the messages at a given layer. In Fig. 11, the values of
${\nu_{\uparrow}^{(\ell)}}(X_{i}^{(\ell)})$ are reported for the bottom $5$
levels of a Random Hierarchical Model with $L=10$, $s=2$, $v=32$, $m=8$, and
noise level $\epsilon=0.5$. At each layer $\ell$, the index $i$ of the nodes
goes from $1$ to $s^{L-\ell}$ and for each of them there are $v=32$ messages
${\nu_{\uparrow}^{(\ell)}}$, one for each entry of the alphabet. At layer
$\ell=0$, the messages ${\nu_{\uparrow}^{(0)}}$ are initialized according to
Eq. 18. After one iteration, at layer $\ell=1$ we observe that the largest
messages at each node, those corresponding to the most probable features
$x_{i}$, are spread around some mean value that is well captured by the
theoretical prediction of Eq. 54. We observe that also for the upper layers,
the average behaviour of the largest messages is well captured by the theory.
The comparison between the theory and the BP algorithm for every $\epsilon$ is
reported in the Figs. 12 and 13. The upward iteration is reported in Fig. 12
and shows an excellent agreement with the prediction of Eq. 54. In particular,
going from the input layer $\ell=0$ to the class variable $\ell=L$ ($L=10$ in
Fig. 12), the messages for the most probable features show a sharper
transition at a threshold value $\epsilon^{*}$, which corresponds to the phase
transition of the theoretical iteration map in Eq. 38 when
$L\rightarrow\infty$. The downward iteration in Fig. 13 shows that the theory
also captures the trend in $\epsilon$ of the downward messages. However, for
small values of $\epsilon$, we observe that the messages have large
fluctuations around their mean value. The reason for this behaviour is that,
in the Random Hierarchical Model, there is a number $m$ of synonyms
$(x_{1},\dots,x_{s})$ that codes for the same higher level feature $y$.
Therefore, having perfect information on $y$ and on $x_{2},\dots,x_{s}$ is not
enough to perfectly reconstruct the value of $x_{1}$, thereby resulting in
large fluctuations of the messages at small noise level $\epsilon$. This is
different from the upward process where, having perfect information on
$(x_{1},\dots,x_{s})$ allows the perfect reconstruction of $y$. As a result,
the messages in the upward process are more concentrated around their mean
than the downward messages.

Figure 11: Upward BP messages for layers $0$ to $4$ for $\epsilon=0.5$,
$v=32$, $s=2$, $L=10$, $sf=0.5$. Each node has $v$ messages, one per possible
feature. At the input layer (layer $0$), messages have value $1-\epsilon=0.5$
or $\epsilon/v=0.5/32$. Going upward, the values of the messages fluctuate,
but they stay separated into two distinct groups: large messages (i.e., the
most probable feature for each node) and small ones. The annealed mean-field
computation (represented with a black line) captures well the mean value of
the large messages (red dashed line). Figure 12: Largest upward BP messages,
averaged for each layer, for varying $\epsilon$. Data for the Random
Hierarchical Model with $v=32$, $s=2$, $L=10$, $sf=0.5$. Each layer, indicated
in the legend, is represented with a different colour and the black dashed
line is the theoretical prediction from Eq. 54, which shows excellent
agreement with the experiments. The top left panel represents all the layers
together for comparison. Starting from the initialization
$\nu_{\uparrow}=1-\epsilon+\epsilon/v$ at layer $0$, we observe that the
largest upward messages increase as we go to higher levels in the hierarchy
only if $\epsilon$ is smaller than some threshold value. For $\epsilon$ larger
than this threshold, instead, the messages become smaller, and it is not
possible to reconstruct the highest levels in the hierarchy better than random
chance. Figure 13: Largest downward BP messages, averaged for each layer, for
varying $\epsilon$. Data for the Random Hierarchical Model with the same
parameters as Fig. 12. Each layer, indicated in the legend, is represented
with a different colour, while the theoretical prediction from Eq. 54 is
represented with the black dashed line. We observe that the messages in the
downward process have large fluctuations, as represented by their standard
deviations, especially for small $\epsilon$. Still, the theory correctly
captures the trend and becomes more accurate for increasing $\epsilon$. The
top left panel represents all the layers together for comparison. Starting
from the initialization $\nu_{\downarrow}=1/v$ at the top layer ($10$), we
observe that the largest downward messages increase as we go to lower levels
in the hierarchy only if $\epsilon$ is smaller than some threshold value.
Figure 14: Largest marginal probabilities computed by BP, averaged for each
layer, for varying $\epsilon$. Data for the Random Hierarchical Model with the
same parameters as Fig. 12. Each layer, indicated in the legend, is
represented with a different colour and the black dashed line is the
theoretical prediction from Eq. 59, which shows excellent agreement with the
experiments. The top left panel represents all the layers together for
comparison, where the inversion between the top and bottom layers can be
observed (same curves as Fig. 7 in the main text).

##  Appendix C Mapping from time diffusion to $\epsilon$ noise

In the diffusion process for the Random Hierarchy Model defined in Section 4,
the beliefs ${\nu_{\uparrow}^{(0)}}$ at the input variables vary
stochastically in time, according to Equation 4. Instead, in the simplified
model of noise considered in Section 5, at a given noise level $\epsilon$,
these beliefs are fixed to two possible values (cf. Equation 7). To study
whether the $\epsilon$-process is an effective approximation of the time
diffusion process, we define an effective $\epsilon(t)$ depending on the
reverse time of diffusion. At each input node $X^{(0)}_{i}$, we consider the
upward messages ${\nu_{\uparrow}^{(0)}}(x)$ associated to the values $x$ that
are different from the value of $X^{(0)}_{i}$ at time $t=0$. Denoting them as
$\nu_{t}$, we define

$\frac{\epsilon(t)}{v}=\langle\nu_{t}\rangle,$ (60)

where the average $\langle\nu_{t}\rangle$ is performed over all the leaves
variables $i$ and the realizations of the dynamics. $\epsilon(t)$ increases
exponentially in time, according to the noise schedule used in the diffusion
process, as shown in the left panel of Figure 15. The probability of the
correct reconstruction of a given node in the diffusion process is reported as
a function of $\epsilon(t)$ in the right panel of Figure 15. We observe that
the curves for different layers have similar behaviour to those of the
$\epsilon$-process presented in Figure 7 of the main text, confirming that the
latter is an effective approximation of denoising diffusion.

Figure 15: (Left panel) Mapping between the $\epsilon$ values and the
diffusion process. From the average values of the beliefs at the leaves
variables during the diffusion process at time $t$, we compute an effective
$\epsilon(t)$ as $\epsilon(t)/v=\langle{\nu_{\uparrow}^{(0)}}(x)\rangle$, for
the values $x$ different from the starting one, averaging over the
realizations of the diffusion process. (Right panel) Probability of
reconstructing the initial values for the nodes at a given layer during the
diffusion process in time, using the effective $\epsilon$ computed in the left
panel. We observe that the shape of the curves with respect to the effective
noise $\epsilon(t)$ is qualitatively similar to that of the simplified
$\epsilon$-process reported in Fig. 7, supporting that it represents a good
approximation for studying the diffusion process in time.

##  Appendix D Hidden activations for different architectures

We perform the experiments described in Section 2.3 using the internal
representations of different deep convolutional architectures trained for
image classification on ImageNet-1k. We consider the ResNet architecture He et
al. (2016) with varying width and depth: a ResNet 50 achieving 95.4% top-5
accuracy, a Wide ResNet 50 having 95.8% top-5 accuracy, and a ResNet 152
having 96.0% top-5 accuracy TorchVision maintainers and contributors (2016).
The results of the experiments performed with the hidden representations of
these architectures are reported in Figure 16 and show the same qualitative
behaviour as the one observed for the ConvNeXt architecture in Figure 2: the
cosine similarity exhibits a sharp transition for the logits, while it decays
smoothly for the hidden representations at early layers.

(a)

(b)

(c)

Figure 16: Cosine similarity between the post-activations of the convolutional
blocks of different ResNet architectures for the initial images $x_{0}$ and
the synthesized ones $\hat{x}_{0}(t)$. Specifically, panels (a), (b), and (c)
correspond to ResNet50, Wide ResNet50, and ResNet152, respectively He et al.
(2016). As in Figure 2 for the ConvNeXt architecture, the similarity between
logits exhibits a sharp drop around $t\approx T/2$, indicating the change in
class, while the hidden representations of the early layers change more
smoothly. To compute the cosine similarity, all activations are standardised,
i.e., centred around the mean and scaled by the standard deviation computed on
the 50000 images of the ImageNet-1k validation set. At each time, the cosine
similarity values correspond to the maximum of their empirical distribution
over $10000$ images ($10$ per class of ImageNet-1k).