Conditional Diffusion Models are Minimax-Optimal and Manifold-Adaptive for Conditional Distribution Estimation

## Abstract

Abstract We consider a class of conditional forward-backward diffusion models for conditional generative modeling, that is, generating new data given a covariate (or control variable). To formally study the theoretical properties of these conditional generative models, we adopt a statistical framework of distribution regression to characterize the large sample properties of the conditional distribution estimators induced by these conditional forward-backward diffusion models. Here, the conditional distribution of data is assumed to smoothly change over the covariate. In particular, our derived convergence rate is minimax-optimal under the total variation metric within the regimes covered by the existing literature. Additionally, we extend our theory by allowing both the data and the covariate variable to potentially admit a low-dimensional manifold structure. In this scenario, we demonstrate that the conditional forward-backward diffusion model can adapt to both manifold structures, meaning that the derived estimation error bound (under the Wasserstein metric) depends only on the intrinsic dimensionalities of the data and the covariate.

## 1 Introduction

Conditional distribution estimation aims to estimate the distribution (or its density if exists) of a response variable $Y$ given some covariate or predictor variable $X$, which is a fundamental problem in statistics with wide applicability in finance, economics [], biology [] and social science, to name just a few.
The conditional distribution provides a full characterization of the dependence structure of the response variable on the predictors, which allows one to gain deeper insights about the data characteristics beyond those from a simple mean regression model, such as capturing uncertainty and addressing multiple-modality.
Conditional density estimation has received significant attention from both statistics and machine learning community with proposed estimators ranging from classical nonparametric estimates such as those based on smoothing techniques [, , , ], Bayesian nonparametric estimates [], to some recent methods that utilize deep neural networks [].

Although there is a rich literature on conditional distribution estimation, many existing methods, such as the classical nonparametric estimators based on kernel smoothing [, , ], suffer from some limitations. One notable drawback of these classical methods is the requirement for the existence of a conditional density function, which is often violated when the response variable $Y$ contains discrete components or itself is a high-dimensional object with low-dimensional structures. As a result, most classical methods can only deal with data with small dimensions, and their performance deteriorates quickly when the dimension increases, thus they suffer from the curse of dimensionality.
In addition, these classical estimators generally do not have the ability to adapt to any potential intrinsic structure, such as the manifold structure of the data, a characteristic of many modern high-dimensional datasets.
For distribution or density estimation in the unconditional setting, estimators based on deep generative models appear to overcome the aforementioned challenges. Constructing a distribution estimator implicitly by specifying its data-generating process naturally allows singular structures in the data.
Additionally, an emerging body of literature on the theoretical understanding of deep generative models, including diffusion-based models [, , , ], demonstrates that such models provide an estimator of the underlying distribution with convergence rates dependent only on the intrinsic dimension of the data.

Motivated by advancements in deep generative modeling, in this work we explore *conditional diffusion models* based on deep neural networks for conditional distribution estimation, accounting for possible low-dimensional manifold structures on either (or both) the covariate $X$ and response $Y$. Unlike other generative model estimation procedures, such as GANs [] and variational auto-encoders [], which explicitly incorporate low-dimensional structures by operating in or maintaining a low-dimensional latent space, diffusion models operate directly in the original ambient data space. Therefore, it is both interesting and important to formally study whether they can still adapt to any low-dimensional structure, if present.
Towards these goals, we consider conditional distribution estimators implicitly defined through a class of conditional forward-backward diffusion models with conditional score matching. We investigate theoretical properties of such estimators through a finite-sample analysis of their statistical error bounds with respect to various metrics and examine their dependence on the intrinsic dimension and certain smoothness characteristics of the data.
The key findings and contributions of our work can be summarized in the following:

- •
Our rates are minimax-optimal under the total variation metric in the classical setting when the conditional distribution admits a smooth density function that also varies smoothly across different covariate values.
- •
Our models encompass unconditional distribution estimation and nonparametric mean regression as special cases. When restricted to the former, our derived estimation error bounds achieve the minimax rate under both the total variation and the Wasserstein metrics. For the latter, our rates recover the classical minimax rate of nonparametric regression under the $L_{2}$ risk.
- •
Our results show that conditional diffusion estimators are adaptive to intrinsic manifold structures when either (or both) the covariate $X$ and response $Y$ are concentrated around some lower-dimensional manifold; thus, our model can handle high-dimensional *distribution regression* with covariates exhibiting low-dimensional structures.

Other related works.
There is vast literature on nonparametric conditional distribution estimation. In addition to smoothing-based methods such as the ones employ kernel smoothing or local polynomial regression [], there are other approached based on mixture model [] , Gaussian processes [, ], and nonparametric Bayes [, ], among others.
There has been a recent line of work that utilizes deep generative approach for conditional sampling such as [] and []. [] utilizes conditional GAN based approach and derived a consistent conditional density estimator but no convergence rates or error bounds are provided. There is also a growing body of theoretical research on diffusion generative models, though most have not considered the conditional setting as we do, such as [, , , , , , , , , , ] and []. Among these, [, ] and [] study the approximation error and generalization ability of diffusion models for estimating (unconditional) distributions when data exhibits low-dimensional structures, which are shown to attain the minimax optimality in the 1-Wasserstein metric for distributions supported on low-dimensional hyperplanes [] and general submanifolds []. Some works, such as  [] and [], explicitly leverage low-dimensional data structures during the generative sampling phase of implementing diffusion models. In scenarios where a covariate is available, some recent works [, ] explore the theoretical properties of conditional diffusion models; however, they do not account for manifold structures when deriving their convergence rates.

## 2 Forward-Backward Diffusion Model and its Conditional Variant

In this section, we will begin by reviewing the forward-backward diffusion model for (unconditional) distribution estimation. After that, we will introduce its adaptations for estimating conditional distributions under the statistical framework of distribution regression.

### 2.1 Forward-backward diffusion model with score matching

Forward-backward diffusion models (see, e.g., [, , , ]) have emerged as a new state-of-the-art class of generative models for estimating and generating samples from an underlying data distribution $\mu^{\ast}$ on a data space $\mathcal{M}_{Y}\subset\mathbb{R}^{D_{Y}}$. In a typical forward-backward diffusion model, two diffusion processes are utilized collaboratively: one process is designed for the estimation of a time-dependent score function describing the direction towards a high data probability region, while the other is for generating samples through a time-inhomogeneous process, based on the estimated score functions. Consequently, this model overcomes the slow convergence issues (e.g., due to the multimodality of $\mu^{\ast}$) commonly observed in models that rely solely on a single diffusion process, such as Langevin diffusion. Throughout the remainder of this paper, the term “diffusion model” specifically refers to the forward-backward diffusion model.

More concretely, the first diffusion process in the diffusion model, often referred to as the forward diffusion, employs a simple diffusion starting from $p_{\rm\scriptsize data}$ that admits a closed-form solution and converges exponentially quickly to its limiting distribution. In this paper, we focus on the commonly used Ornstein–Uhlenbeck (OU) process as the forward process, which gradually injects Gaussian noise into the data and is described as a stochastic differntial equation (SDE)

$$ $\displaystyle{\rm d}\overrightarrow{Y}_{t}=-\delta_{t}\overrightarrow{Y}_{t}\, {\rm d}t+\sqrt{2\delta_{t}}\,{\rm d}B_{t},\ \overrightarrow{Y}_{0}\sim\mu^{*},$ (1) $$

where $\{B_{t}\,:\,t>0\}$ denotes the standard Brownian motion in $\mathbb{R}^{D_{Y}}$ and $\{\delta_{t}:\,t\geq 0\}$ is some (possibly time-dependent) drift coefficient. Note that the OU process admits the closed form solution $\overrightarrow{Y}_{t}=m_{t}\overrightarrow{Y}_{0}+\int_{0}^{t}\frac{m_{t}}{m_
{s}}\,\sqrt{2\delta_{s}}\,{\rm d}B_{s}$; thus, the conditional distribution of $\overrightarrow{Y}_{t}$ given $\overrightarrow{Y}_{0}=y$ is $\mathcal{N}(m_{t}\,y,\,\sigma_{t}^{2}I_{D_{Y}})$, where $m_{t}=\exp\big{(}-\int_{0}^{t}\delta_{s}\,{\rm d}s\big{)}$ and $\sigma_{t}^{2}=1-m_{t}^{2}$.
Therefore, the marginal distribution of $\overrightarrow{Y}_{t}$, denoted as $p_{t}$, converges exponentially quickly to its limiting distribution $p_{\infty}=\mathcal{N}(0,I_{D_{Y}})$ under the Kullback–Leibler divergence.

The second diffusion process in the diffusion model, usually called the backward diffusion, reverses the forward diffusion and can be written as the following SDE,

$$ $\displaystyle{\rm d}\overleftarrow{Y}_{t}=\big{[}\delta_{T-t}\overleftarrow{Y} _{t}\ +$ $\displaystyle\ 2\delta_{T-t}\nabla\log p_{T-t}(\overleftarrow{Y}_{t})\big{]}\, {\rm d}t+\,\sqrt{2\delta_{T-t}}\,{\rm d}B_{t},\ \ \overleftarrow{Y}_{0}\sim p_ {T}.$ (2) $$

Under mild conditions on $\mu^{\ast}$ [, ] (valid for our setting), the distribution of $\overleftarrow{Y}_{t}$ is $p_{T-t}$, so that $\overleftarrow{Y}_{T}\sim p_{0}=\mu^{*}$. Since $p_{T}$ is close to $p_{\infty}=\mathcal{N}(0,I_{D_{Y}})$, one can instead initialize the backward diffusion using the easy-to-sample distribution $p_{\infty}$, i.e. set $\overleftarrow{Y}_{0}\sim\mathcal{N}(0,I_{D_{Y}})$.
The drift term of the backward diffusion depends on the time dependent score function $\nabla\log p_{t}$ defined through the forward diffusion; therefore, the forward and the backward diffusions together constitutes a generative model for sampling from $\mu^{*}$.

Equations (1) and (2) define the forward-backward diffusion model at the population-level. In a standard statistical setting, we utilize independent and identically distributed (i.i.d.) samples $\{Y_{i}\}_{i=1}^{n}$ from $\mu^{*}$ to estimate the time-dependent score function $\nabla\log p_{t}$ in the backward diffusion. The estimation is achieved by the so-called score matching [, ]. Specifically, one first numerically simulates for some sufficiently large time horizon $T$ from a sample-level forward process $\{y_{t}:\,t\in[0,T]\}$, which is SDE (1) initialized at the empirical distribution of the data, that is, $y_{0}\sim\widehat{\mu}_{n}=n^{-1}\sum_{i=1}^{n}\delta_{Y_{i}}$, with $\delta_{y}$ denoting the point mass (Dirac) measure at a point $y$. One then uses a score approximating map $S_{\theta}(y,t)$ over space and time, indexed by a parameter $\theta$, e.g., (deep) neural networks with controlled depth and number of non-zero parameters, to estimate the true underlying score function $\nabla\log p_{t}(y)$, by minimizing the following ($L_{2}$-)score matching risk (over $\theta$):

$$ $\displaystyle\int_{\tau}^{T}\mathbb{E}_{y_{t}\sim p_{t}(\cdot\,|\,y_{0}),\,y_{ 0}\sim\widehat{\mu}_{n}}\big{[}\|S_{\theta}(y_{t},t)-\nabla\log p_{t}(y_{t}\,| \,y_{0})\|^{2}\big{]}\,\lambda(t)\,{\rm d}t,$ $$

where $p_{t}(\cdot\,|\,y)$ denotes the distribution of $\overrightarrow{Y}_{t}$ in forward diffusion (1) initialized at $\overrightarrow{Y}_{0}=y$ for any $y\in\mathcal{M}_{Y}$.
Here, $\lambda(t)$ is a weighting function (over time), and $\tau$ is an early-stopping threshold for preventing the explosion (singularity) of the score function as $t\to 0$ commonly employed in practice [, ].
Equivalently, this score estimation step can be efficiently carried out in practice by simulating a trajectory $\{y_{t}:t\geq 0\}$ from SDE (1) starting from each data point $Y_{i}$, that is, $y_{t}\sim p_{t}(\cdot\,|\,Y_{i})$. One then uses $S_{\theta}(y,t)$ to match the ensemble of all sample score functions $\nabla\log p_{t}(y_{t}\,|\,Y_{i})$ over all $n$ simulated trajectories, by minimizing the following empirical risk function (over $\theta$):

$$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\int_{\tau}^{T}\mathbb{E}_{y_{t}\sim p_{ t}(\cdot\,|\,Y_{i})}\big{[}\|S_{\theta}(y_{t},t)-\nabla\log p_{t}(y_{t}\,|\,Y_ {i})\|^{2}\big{]}\,\lambda(t)\,{\rm d}t.$ (3) $$

We will adopt this statistical formulation of score matching to facilitate our theoretical analysis, leveraging tools from statistical learning theory.

Finally, let $\widehat{S}(x,t)=S_{\widehat{\theta}}(x,t)$ denote the resulting score estimator. The distribution estimator of $\mu^{*}$ based on the forward-backward diffusion model is then $\widehat{p}_{T-\tau}$, where $\widehat{p}_{t}$ represents the distribution of $\overleftarrow{Y}^{\dagger}_{t}$ for $t\in[0,T-\tau]$, and $\overleftarrow{Y}^{\dagger}_{t}$ follows the SDE below with a plugged-in score,

$$ $\displaystyle{\rm d}\overleftarrow{Y}^{\dagger}_{t}=\big{[}\delta_{T-t} \overleftarrow{Y}^{\dagger}_{t}\ +\ 2\delta_{T-t}\widehat{S}(\overleftarrow{Y} ^{\dagger}_{t},T-t)\big{]}\,{\rm d}t+\,\sqrt{2\delta_{T-t}}\,{\rm d}B_{t},\ \ \overleftarrow{Y}^{\dagger}_{0}\sim\mathcal{N}(0,I_{D_{Y}}).$ (4) $$

### 2.2 Conditional forward-backward diffusion model with conditional score matching

A notable characteristic of diffusion models is their flexibility in incorporating a covariate or control variable, denoted as $X\in\mathcal{M}_{X}\subset\mathbb{R}^{D_{X}}$, to guide the generation of new data $Y\in\mathcal{M}_{Y}\subset\mathbb{R}^{D_{Y}}$. This can be equivalently formulated as the statistical problem of generating samples from the conditional distribution $\mu^{*}_{Y|x}$ of $Y$ given $X=x$ for any covariate value $x\in\mathcal{M}_{X}$. To facilitate the borrowing of information across different covariate values, it is commonly assumed that the conditional distribution $\mu^{*}_{Y|x}$ varies smoothly with $x\in\mathcal{M}_{X}$. This assumption underlies a statistical framework often referred to as *distribution (density) regression* in the literature [, , ]. Distribution regression expands the classical (nonparametric) mean regression by estimating not only the conditional expectation $\mathbb{E}[Y|X=x]$ as a smooth function of $x$, but also the entire conditional distribution $\mu^{*}_{Y|x}$ that varies smoothly with $x$. Compared to classical distribution regression methods based on kernel smoothing, which require $\mu^{*}_{Y|x}$ to admit a density function (thus termed density regression in the early literature, see [, ]), conditional diffusion model-based methods are more flexible. They can be more generally applicable to cases where $\mu^{*}_{Y|x}$ is supported on a low-dimensional manifold and is therefore singular; see Section 3.2 for details.

A natural way to convert a diffusion model into a conditional diffusion model for sampling from $\mu^{*}_{Y|x}$ is to replace the (unconditional) score function $\nabla\log p_{t}$ in the backward diffusion with some conditional score function $\nabla\log p_{t}(\cdot\,|\,x)$ satisfying $p_{0}(\cdot\,|\,x)=\mu^{*}_{Y|x}(\cdot)$ and $p_{T}(\cdot\,|\,x)=p_{\infty}(\cdot)$.
Earlier literature considers the so-called classifier guidance method for estimating the conditional score using Bayes’ rule (especially when covariate
$x$ is discrete or categorical; see, e.g., [, ]):

$$ $\displaystyle\nabla\log p_{t}(y_{t}\,|\,x)=\nabla\log p_{t}(y_{t})+\nabla\log c _{t}(x\,|\,y_{t}),$ $$

Here, the first term, $\nabla\log p_{t}(y_{t})$, is the unconditional score function defined in the forward diffusion (1) in the unconditional diffusion model, which can be estimated via score matching. The second term, $\nabla\log c_{t}(x\,|\,y_{t})$, is the likelihood function of an external “classifier” trained to predict
$x$ from $y_{t}$. In other words, the classifier guidance method incorporates information from the covariate value $x$ into the unconditional score function through the gradient of an external classifier to guide the backward diffusion for generating samples from $\mu^{*}_{Y|x}$.

In this work, for the sake of theoretical simplicity and to avoid the need to analyze and quantify the statistical accuracy of an external classifier, we consider a classifier-free method, which at the population level directly applying a *conditional score matching* method [, , , ] based on simulating $\overrightarrow{Y}_{t}$ from the same (marginal) forward diffusion

$$ $\displaystyle{\rm d}\overrightarrow{Y}_{t}=-\delta_{t}\overrightarrow{Y}_{t}\, {\rm d}t+\sqrt{2\delta_{t}}\,{\rm d}B_{t},\ \overrightarrow{Y}_{0}\sim\mu^{*}_ {Y}.$ (5) $$

Here, $\mu^{*}_{Y}$ denotes the marginal distribution of $Y$. Consider a generic conditional score approximating map $S_{\theta}(y,x,t)$ over space, covariate value and time, indexed by a parameter $\theta$. In the conditional score matching step, we minimize the following ($L_{2}$-)conditional score matching risk over $\theta$:

$$ $\displaystyle\int_{\tau}^{T}\mathbb{E}_{y_{t}\sim p_{t}(\cdot\,|\,y_{0}),\,(x, y_{0})\sim\mu^{\ast}_{X,Y}}\big{[}\|S_{\theta}(y_{t},x,t)-\nabla\log p_{t}(y_{ t}\,|\,y_{0})\|^{2}\big{]}\,\lambda(t)\,{\rm d}t,$ (6) $$

where $\mu^{*}_{X,Y}$ denotes the joint distribution of $(X,Y)$, and $p_{t}(\cdot\,|\,y)$ is the distribution of $\overrightarrow{Y}_{t}$ in forward diffusion (5) initialized at $\overrightarrow{Y}_{0}=y$ for any $y\in\mathcal{M}_{Y}$.
Here, the early stopping threshold $\tau$ and the weighting function $\lambda(t)$ are defined as before. It is straightforward to show (see Lemma C.12 in Appendix C.2) that if $S_{\theta}(y,x,t)$ can range over all possible conditional score functions, then the global minimizer of the preceding risk function is precisely the true underlying conditional score function $\nabla\log p_{t}(y_{t}\,|\,x)$. Here, $p_{t}(y_{t}\,|\,x):\,=\mathbb{E}_{y_{0}\sim\mu^{*}_{Y|x}}[p_{t}(y_{t}\,|\,y_{0
})]$ denotes the conditional distribution of $\overrightarrow{Y}_{t}$ given $X=x$ after marginalizing out $\overrightarrow{Y}_{0}$, where $(X,\overrightarrow{Y}_{0})\sim\mu^{*}_{X,Y}$ and $\overrightarrow{Y}_{t}$ follows forward diffusion (5) starting from $\overrightarrow{Y}_{0}$.

The corresponding conditional backward diffusion for sampling from $\mu^{\ast}_{Y|x}$ is then given by

$$ $\displaystyle{\rm d}\overleftarrow{Y}_{t|x}=\big{[}\delta_{T-t}\overleftarrow{ Y}_{t|x}\ +$ $\displaystyle\ 2\delta_{T-t}\nabla\log p_{T-t}(\overleftarrow{Y}_{t|x}\,|\,x) \big{]}\,{\rm d}t+\,\sqrt{2\delta_{T-t}}\,{\rm d}B_{t},\ \ \overleftarrow{Y}_{ 0|x}\sim p_{T}(\cdot\,|\,x).$ (7) $$

Here, we added a subscript $x$ in the notation $\overleftarrow{Y}_{t|x}$ to indicate that, unlike the forward process (5), the backward diffusion process is $x$-dependent. Similar to the (unconditional) diffusion model, choosing a sufficiently large $T$ can guarantee $p_{T}(\cdot\,|\,x)\approx p_{\infty}(\cdot)=N(0,\,I_{D_{Y}})$, which is independent of $x$. We will refer to equations (5) and (7) as the *conditional (forward-backward) diffusion model* for the remainder of the paper.

To estimate the conditional score $\nabla\log p_{t}(\cdot\,|\,x)$ using i.i.d. observations $\{(X_{i},Y_{i})\}_{i=1}^{n}$ sampled from the joint distribution $\mu^{*}_{X,Y}$ of $(X,Y)$ under the statistical framework of distribution regression, one can again minimize the following empirical version of conditional score matching risk (6),

$$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}$ $\displaystyle\,\int_{\tau}^{T}\mathbb{E}_{y_{t}\sim p_{t}(\cdot|Y_{i})}\big{[} \|S_{\theta}(y_{t},X_{i},t)-\nabla\log p_{t}(y_{t}\,|\,Y_{i})\|^{2}\big{]}\, \lambda(t)\,{\rm d}t.$ (8) $$

Finally, let $\widehat{S}(y,x,t)=S_{\widehat{\theta}}(y,x,t)$ denote the corresponding conditional score estimator. For each $x\in\mathcal{M}_{X}$, the conditional distribution estimator of $\mu^{*}_{Y|x}$ based the conditional forward-backward diffusion model is then $\widehat{p}_{T-\tau}(\cdot\,|\,x)$, where $\widehat{p}_{t}(\cdot\,|\,x)$ is the distribution of $\overleftarrow{Y}^{\dagger}_{t|x}$ for $t\in[0,T-\tau]$, and $\overleftarrow{Y}^{\dagger}_{t|x}$ follows the SDE below with a plugged-in conditional score,

$$ $\displaystyle{\rm d}\overleftarrow{Y}^{\dagger}_{t|x}=\big{[}\delta_{T-t} \overleftarrow{Y}^{\dagger}_{t|x}\ +$ $\displaystyle\ 2\delta_{T-t}\widehat{S}(\overleftarrow{Y}^{\dagger}_{t|x},x,T- t)\big{]}\,{\rm d}t+\,\sqrt{2\delta_{T-t}}\,{\rm d}B_{t},\ \ \overleftarrow{Y} ^{\dagger}_{0|x}\sim\mathcal{N}(0,I_{D_{Y}}).$ (9) $$

### 2.3 Neural network class for conditional score function approximation

Definition (neural network class): A class of neural networks $\Phi(H,\,W,\,R,\,B,\,V)$ with height $H$, width vector $W=(W_{1},\,W_{2},\ldots,\,W_{H+1})$, sparsity $R$, norm constraint $B$, and function norm constraint $V$ is defined as
$\Phi(H,W,R,B,V)=\big{\{}f(\cdot)=(A^{(H)}\operatorname{ReLU}(\cdot)+b^{(H)})
\circ\cdots\circ(A^{(2)}\operatorname{ReLU}(\cdot)+b^{(2)})\circ(A^{(1)}x+b^{(
1)}),\ \mbox{such that}\ A^{(i)}\in\mathbb{R}^{W_{i}\times W_{i+1}};\,b^{(i)}
\in\mathbb{R}^{W_{i+1}};\,\sum_{i=1}^{H}(\|A^{(i)}\|_{0}+\|b^{(i)}\|_{0})\leq R
;\,\max_{i}\|A^{(i)}\|_{\infty}\vee\|b^{(i)}\|_{\infty}\leq B;\,\|f\|_{\infty}
\leq V\,\big{\}}$, where $\operatorname{ReLU}(x)=\max\{0,\,x\}$ denotes the rectified linear unit activation function and the $\max$ function is applied elementwise to a vector.

According to [] and [], the smoothing effect of gradually injecting Gaussian noise into the data distribution during the forward diffusion (5) suggests that the optimal size of the neural network for effectively approximating $\nabla\log p_{t}(\cdot,|,x)$ should decrease as $t$ increases. This observation motivates us to consider a neural network class whose size diminishes over time. For technical convenience, we discretize the time and adopt the following piece-wise constant complexity neural network class, as utilized in []:

$$ $\begin{array}[]{r}\mathcal{S}_{NN}=\Big{\{}S(y,x,t)=\sum_{i=1}^{\mathcal{I}}S_ {i}(y,x,t)\cdot\mathbf{1}\left(t_{i-1}\leq t<t_{i}\right)\,\Big{|}\,S_{i}\in \Phi\left(H_{i},W_{i},R_{i},B_{i},V_{i}\right),i\in[\mathcal{I}]\Big{\}},\end{array}$ $$

where $\tau=t_{0}<t_{1}<\cdots<t_{\mathcal{I}}=T$, $\frac{t_{i+1}}{t_{i}}=2$ for any $i\in[\mathcal{I}]$, and $\tau=2^{-\mathcal{I}}T$ for some $\mathcal{I}$ to be determined later.
We have also conducted a simulation study (see Appendix A) to demonstrate the effectiveness of this theoretically guided neural network architecture compared to a standard single ReLU neural network (across both space and time). Recall that $\widehat{S}(y,x,t)$ denotes the conditional score estimator, defined as the minimizer of the conditional score matching risk (8) over the class $\mathcal{S}_{NN}$ with the weight function $\lambda(t)=t$ (although any other weights such as $\lambda(t)\equiv 1$ would also suffice). We define a (truncated) estimator $\widehat{\mu}_{Y|x}$ for $\mu^{*}_{Y|x}$ as the distribution of $\overleftarrow{Y}^{\dagger}_{T-\tau|x}\cdot\mathbb{1}(\|\overleftarrow{Y}^{
\dagger}_{T-\tau|x}\|_{\infty}\leq L)$, where $(L,T)$ are large enough constants so that $\mathcal{M}_{Y}\subset\mathbb{B}_{\mathbb{R}^{D_{Y}}}(0,L/2)$ and $p_{T}(\cdot,|,x)\approx\mathcal{N}(0,I_{D_{Y}})$. Here, we truncate the random variable $\overleftarrow{Y}^{\dagger}_{T-\tau|x}$ to guarantee a bounded support for the induced distribution estimator $\widehat{\mu}_{Y|x}$, which is solely for technical reasons.

## 3 Main Theoretical Results

In this section, we present our main theoretical results characterizing the statistical accuracy of the conditional diffusion model for conditional distribution estimation (or distribution regression) under two scenarios. In the first scenario, we consider the classical density regression setting where the conditional distribution $\mu^{*}_{Y|x}$ admits a density function relative to the Lebesgue measure of the data space $\mathcal{M}_{Y}\subset\mathbb{R}^{D_{Y}}$. We derive the convergence rate of the estimator under both the total variation and the Wasserstein metrics. In particular, our derived convergence rate is minimax-optimal under the total variation metric within the regime covered by the existing literature [] (see Remark 1 for further details), and extends to a broader regime.
In the second scenario, we consider a high-dimensional distribution regression setting where both the response variable $Y\in\mathbb{R}^{D_{Y}}$ and the covariate variable $X\in\mathbb{R}^{D_{X}}$ reside in high-dimensional ambient spaces characterized by large $D_{Y}$ and $D_{X}$. However, the covariate space $\mathcal{M}_{X}$ of $X$ has an intrinsic (or effective) dimension $d_{X}$ that is significantly smaller than $D_{X}$. Furthermore, given any $x\in\mathcal{M}_{X}$, the corresponding data space of $Y$, denoted as $\mathcal{M}_{Y|x}$, can be $x$-dependent and also has a small intrinsic dimension $d_{Y}<D_{Y}$.
We demonstrate that the conditional diffusion model effectively adapts to the underlying manifold structures of both the data and the covariate variable. Specifically, we show that the convergence rate of the estimator depends solely on the intrinsic dimensions $(d_{Y},d_{X})$, rather than the ambient dimensions $(D_{Y},D_{X})$.

In the following, we denote $d_{\rm TV}(\mu,\nu)$ and $W_{1}(\mu,\nu)$ as the respective total variation distance and the 1-Wasserstein distance between two distributions $\mu$ and $\nu$. We denote $\mathcal{M}=\{(x,y)\,:\,x\in\mathcal{M}_{X},y\in\mathcal{M}_{Y|x}\}$ as the joint space of $(X,Y)$ and $\mathcal{M}_{Y}=\bigcup_{x\in\mathcal{M}_{X}}\mathcal{M}_{Y|x}$ as the (marginal) data space. We use the notation $a\vee b$ and $a\wedge b$ to denote the respectively shorthand of $\max\{a,b\}$ and $\min\{a,b\}$. For a sequence $\{a_{n}:\,n\geq 1\}$, we use $\Theta(a_{n})$ to indicate the order of $a_{n}$ up to a multiplicative constant as $n\to\infty$, and $\widetilde{\Theta}(a_{n})$ to indicate the order of $a_{n}$ up to a multiplicative constant and logarithmic terms of $n$. Similarly, we use $\mathcal{O}(a_{n})$ and $\widetilde{\mathcal{O}}(a_{n})$ to indicate at most of order $a_{n}$.

### 3.1 Classical density regression in Euclidean space

In this subsection, we consider the classical density regression setting where both the covariate space $\mathcal{M}_{X}\subset\mathbb{R}^{D_{X}}$ and the data space $\mathcal{M}_{Y}\subset\mathbb{R}^{D_{Y}}$ are compact subsets (with open interiors) of the Euclidean spaces, and the conditional distribution $\mu^{*}_{Y|x}$ admits a density function, denoted as $\mu^{*}(y\,|\,x)$, relative to the Lebesgue measure of $\mathbb{R}^{D_{Y}}$. For simplicity, we assume $\mathcal{M}_{Y}=[-1,1]^{D_{Y}}$ and $\mathcal{M}_{X}=[-1,1]^{D_{X}}$. In order to derive a non-asymptotic bound to the expected total variation distance and Wasserstein distance between the conditional distribution estimator $\widehat{\mu}_{Y|X}$ and the target $\mu^{*}_{Y|X}$, with the expectation taken over $X\sim\mu^{*}_{X}$, we impose certain smoothness condition to the condition density function $\mu^{*}(y\,|\,x)$ relative to $(y,x)$ as in the classical density regression literature [, ]. Specifically, we assume that $\mu^{*}(y,|,x)$, as a function of $(y,x)$, is $C^{\alpha_{Y},\alpha_{X}}$-smooth, where $\alpha_{Y}$ and $\alpha_{X}$ quantify the respective smoothness in the response variable $y$ and the covariate $x$. Note that a function $f(y,x)$ being $C^{\alpha_{Y},\alpha_{X}}$-smooth implies that, around any point $(y_{0},x_{0})$, there exists a local polynomial approximation of $f$, with an approximation error of order $\mathcal{O}(\|y-y_{0}\|^{\alpha_{Y}}+\|x-x_{0}\|^{\alpha_{X}})$; a rigorous definition can be found in Appendix B. Formally, we make the following assumptions.

Assumption A (smoothness and lower boundness of $\mu^{*}_{Y|x}$): For each $x\in\mathcal{M}_{X}$, the conditional distribution $\mu^{*}_{Y|x}$ admits a density function $\mu^{*}(y\,|\,x)$ that is $C^{\alpha_{Y},\alpha_{X}}$-smooth in $(y,x)$. Moreover, there exists a positive constant $c$ so that $\mu^{*}(y\,|\,x)\geq c$ holds for any $x\in\mathcal{M}_{X},y\in\mathcal{M}_{Y}$.

Assumption B (regularity of the drift coefficient): The drift coefficient $\delta_{t}$ is infinitely differentiable and there exist positive constants $c_{1},c_{2}$ so that $c_{1}\leq\delta_{t}\leq c_{2}$ for any $t\geq 0$.

Here the lower bound requirement of $\mu^{*}(y\,|\,x)$ is a commonly made assumption for distribution estimation in the classical density regression literature, and is also imposed in [, ] for analyzing (unconditional) diffusion models.
In practical applications, the drift coefficient $\delta_{t}$ is typically chosen as a positive constant independent of $t$, thus naturally satisfying Assumption B.

###### Theorem 1 (Density regression in Euclidean space) .

Suppose Assumptions A and B are satisfied. Let $\varepsilon_{1}=n^{-1\,\big{/}\,\big{(}2\alpha_{Y}+D_{Y}+\frac{\alpha_{Y}}{
\alpha_{X}}D_{X}\big{)}}$. If we take $\tau=\widetilde{\Theta}\big{(}\varepsilon_{1}^{2(\alpha_{Y}+1)}\big{)}$, $T=\Theta(\log n)$, and neural network sizes satisfying $H_{i}=\Theta(\log^{4}n)$, $\|W_{i}\|_{\infty}=\widetilde{\Theta}\big{(}n^{\frac{D_{X}}{2\alpha_{X}+D_{X}}
}t_{i}^{\frac{-\alpha_{X}D_{Y}}{2\alpha_{X}+D_{X}}}\wedge\varepsilon_{1}^{-D_{
Y}-\frac{\alpha_{Y}D_{X}}{\alpha_{X}}}\big{)}$,
$R_{i}=\widetilde{\Theta}\big{(}n^{\frac{D_{X}}{2\alpha_{X}+D_{X}}}t_{i}^{\frac
{-\alpha_{X}D_{Y}}{2\alpha_{X}+D_{X}}}\wedge\varepsilon_{1}^{-D_{Y}-\frac{
\alpha_{Y}D_{X}}{\alpha_{X}}}\big{)}$, $B_{i}=\exp(\Theta(\log n^{4}))$ and $V_{i}=\Theta\big{(}\sqrt{\frac{\log n}{t_{i}\wedge 1}}\big{)}$ for $i\in[\mathcal{I}]$(^1^11Here we use the notation $[\mathcal{I}]=\{1,2,\cdots,\mathcal{I}\}$.) with $\mathcal{I}=\log_{2}(\frac{T}{\tau})$, then it holds with probability at least $1-n^{-1}$ that

$$ $\displaystyle\,\mathbb{E}_{x\sim\mu^{*}_{X}}\big{[}d_{\rm TV}(\widehat{\mu}_{Y |x},\,\mu^{*}_{Y|x})\big{]}=\widetilde{\mathcal{O}}\bigg{(}n^{-\frac{1}{2+ \frac{D_{X}}{\alpha_{X}}+\frac{D_{Y}}{\alpha_{Y}}}}\bigg{)},$ (10) and $\displaystyle\,\mathbb{E}_{x\sim\mu^{*}_{X}}\big{[}W_{1}(\widehat{\mu}_{Y|x}, \,\mu^{*}_{Y|x})\big{]}=\widetilde{\mathcal{O}}\bigg{(}n^{-\frac{1}{2+\frac{D_ {X}}{\alpha_{X}}}}\vee n^{-\frac{1+\frac{1}{\alpha_{Y}}}{2+\frac{D_{X}}{\alpha _{X}}+\frac{D_{Y}}{\alpha_{Y}}}}\bigg{)}.$ (11) $$

###### Remark 1 .

In the special case when $\alpha_{X}\in[0,1]$, [] shows that a well-designed kernel-based estimator can achieve the same convergence rate (10) in the total variation metric as our conditional diffusion model-based estimator; furthermore, this rate is shown to be minimax-optimal. Therefore, our result implies the minimax-optimality of the conditional diffusion model for density regression in the regime where $\alpha_{X}\in[0,1]$, although our upper bound is also applicable to $\alpha_{X}>1$.

###### Remark 2 .

When specializing to the unconditional case with no covariate (that is, taking $D_{X}=0$ in Theorem 1), our derived estimation error bounds (10) and (11) reduce respectively to the minimax rate of (unconditional) distribution estimation under the total variation metric and the Wasserstein metric [, ]. However, unlike the upper bound proofs in [, ], which rely on generative adversarial network (GAN) type estimators, our proof demonstrates that the diffusion model is also minimax-optimal for distribution estimation. In particular, our results recover those from [] as a special case (by taking $D_{X}=0$).

###### Remark 3 .

The derived $W_{1}$ error bound (11) comprises two terms $n^{-\frac{\alpha_{X}}{2\alpha_{X}+D_{X}}}$ and $n^{-\frac{\alpha_{Y}+1}{2\alpha_{Y}+\frac{D_{X}\alpha_{Y}}{\alpha_{X}}+D_{Y}}}$. The first term resembles the classical minimax rate of nonparametric regression under the $L_{2}$ risk and can be interpreted as mainly capturing the estimation error related to learning the dependence of the response variable $Y$ on the covariate $X$, so that it only depends on the smoothness and intrinsic dimension of $X$. Technically, this term arises from the approximation of the conditional score function for large time $t$, where finer details of the conditional distribution in $Y$ have been smoothed out and only the global dependence on $X$ matters. The second term reflects the estimation error of recovering the entire conditional distribution of $Y$ given $X$, and depends on characteristics related to the response variable $Y$, such as the smoothness $\alpha_{Y}$ of the conditional density function and the dimension of $Y$. Interestingly, the derived rate suggests a phase transition phenomenon: if the dimension of the response variable $D_{Y}$ satisfies $D_{Y}\leq 2+\frac{D_{X}}{\alpha_{X}}$, then the estimation error under the $W_{1}$ metric remains of order $n^{-{\alpha_{X}}/{(2\alpha_{X}+D_{X})}}$ regardless of the smoothness level $\alpha_{Y}$, and the $W_{1}$ estimation error is dominated by the error of capturing the global dependence of the response variable $Y$ on the covariate variable $X$; otherwise, the $W_{1}$ estimation error is influenced by both the smoothness $\alpha_{Y}$ of conditional density on $Y$ and the smoothness $\alpha_{X}$, which captures the finer details of the conditional distribution of $Y$ given $X$.

###### Remark 4 .

A recent related work [] also explores theoretical properties of conditional diffusion model, and show the minimax optimality of diffusion model under the total variation distance. In our work, we allow the conditional distribution of $Y$ given $X$ to have difference smoothness levels $\alpha_{Y}$ and $\alpha_{X}$ on the response $Y$ and covariate $X$; in comparison, [] assumes the two smoothness levels are the same. The varying smoothness levels can allow for the applicability of the results in more general settings. For instance, when specializing to the mean regression case, where the conditional distribution is a Gaussian distribution centered at the evaluation of an $\alpha_{X}$ smooth regression function over $\mathbb{R}^{D_{X}}$, our derived estimation error bound (10) under $d_{\rm TV}$ (that is, taking $\alpha_{Y}\to\infty$ in Theorem 1) can recover the classical minimax rate $n^{-\frac{\alpha_{X}}{2\alpha_{X}+D_{X}}}$ of nonparametric regression under the $L_{2}$ risk.

###### Theorem 1 (Density regression in Euclidean space) .

###### Remark 1 .

###### Remark 2 .

###### Remark 3 .

###### Remark 4 .

### 3.2 High-dimensional distribution regression with low-dimensional manifold structures

In this subsection, we consider the case where both the covariate space $\mathcal{M}_{X}$ and the response space $\mathcal{M}_{Y}$ may have low-dimensional structures in their respective ambient spaces $\mathbb{R}^{D_{X}}$ and $\mathbb{R}^{D_{Y}}$.
For the covariate space $\mathcal{M}_{X}$, the low-dimensional structure is imposed in terms of its upper Minkowski dimension, which is related to the growth of its packing number (see Assumption C). This low-dimensional structure is notably less stringent than a typical manifold assumption, as it does not require any smoothness properties of $\mathcal{M}_{X}$. For the response space $\mathcal{M}_{Y}$, since we allow the conditional distribution $\mu^{*}_{Y|x}$ to have different supports, denoted as $\mathcal{M}_{Y|x}$, for each $x\in\mathcal{M}_{X}$, we can decompose $\mathcal{M}_{Y}$ as $\bigcup_{x\in\mathcal{M}_{X}}\mathcal{M}_{Y|x}$.
In our theory, we require each “section” $\mathcal{M}_{Y|x}$ to be a smooth submanifold in $\mathbb{R}^{D_{Y}}$; additionally, we require $\mathcal{M}_{Y|x}$ to vary smoothly with $x$ (see Assumption D). We list the concrete assumptions as follows.

Assumption C (intrinsic dimension of $\mathcal{M}_{X}$): $\mathcal{M}_{X}$ is compact set in $\mathbb{R}^{D_{X}}$ and there exist constants $(C_{1},C_{2})$ so that for any $\varepsilon>0$, any $\varepsilon_{1}\in(0,\varepsilon)$, and any $x\in\mathcal{M}$, we have(^2^22A set $P\subseteq S$ is a $\varepsilon$-packing of $S$ if for every $x,x^{\prime}\in P$ we have $\|x-x^{\prime}\|>\varepsilon$.) $\mathbb{M}(\mathbb{B}_{\mathcal{M}_{X}}(x,\varepsilon),\|\cdot\|,\varepsilon_{
1}):=\max\{m:\exists\varepsilon_{1}$-packing of $B_{\mathcal{M}_{X}}(x,\varepsilon)$ of size $m\}\leq C_{2}(\frac{\varepsilon_{1}}{\varepsilon})^{-d_{X}}$.

This assumption naturally holds if $\mathcal{M}_{X}$ is a compact subset of a $d_{X}$-dimensional hyperplane, or more generally, a compact $d_{X}$-dimensional submanifold embedded in $\mathbb{R}^{D_{X}}$ with its reach(^3^33The reach of a closed subset $A\subset\mathbb{R}^{D}$ is defined as $\tau_{A}=\inf_{p\in A}{\rm dist}(p,\operatorname{Med}(A))=\inf_{z\in \operatorname{Med}(A)}{\rm dist}(z,A)$, where ${\rm dist}(z,A)=\inf_{p\in A}\|p-z\|$ denotes the distance function to $A$, and $\operatorname{Med}(A)$ is the medial axis of $A$ consisting of the points that have at least two nearest neighbors.) bounded away from zero. Therefore, the constant $d_{X}$ in the assumption can be interpreted as the intrinsic dimension of $\mathcal{M}_{X}$. Next, we introduce our assumption on the conditional distribution $\mu^{*}_{Y|X}$. For easy understanding, we present an informal assumption here and postpone the more rigorous and detailed version to Appendix B in the supplement. Recall that $\mathcal{M}$ denotes the joint space of $(X,Y)$.

Assumption D (smoothness of $\mu^{*}_{Y|X}$, informal version): For any $x\in\mathcal{M}_{X}$, $\mathcal{M}_{Y|x}$ is a $d_{Y}$-dimensional submanifold in $\mathbb{R}^{D_{Y}}$, and $\mu^{*}_{Y|x}$ admits a density with respect to the volume measure of $\mathcal{M}_{Y|x}$, which is uniformly lower bounded away from zero. Moreover, for any $\omega=(x_{0},y_{0})\in\mathcal{M}$ and $x\in\mathbb{B}_{\mathcal{M}_{X}}(x_{0},r_{0})$, there exists an encoder-decoder pair $\big{(}Q^{\omega}_{x}(y),\,G^{\omega}_{x}(z)\big{)}$, such that $Q^{\omega}_{x}(\cdot)$ maps $y\in\mathbb{B}_{\mathcal{M}_{Y|x}}(y_{0},r_{0})$ to a low-dimensional latent variable $z\in\mathbb{R}^{d_{Y}}$, and $G^{\omega}_{x}(\cdot)$ reconstructs the data $y$ through the latent variable $z$. Here, the decoder $G^{\omega}_{x}(z)$ is $C^{\beta_{Y},\beta_{X}}$-smooth in $(z,x)$; and the induced (local) conditional density function $v^{\omega}(z|x)$ of the latent variable $z$, as the pushforward measure through the encoder $Q^{\omega}_{x}(y)$ of the restriction of the measure $\mu^{*}_{Y|x}$ onto $\mathbb{B}_{\mathcal{M}_{Y|x}}(y_{0},r_{0})$, is $C^{\alpha_{Y},\alpha_{X}}$-smooth in $(z,x)$.

Constants $(\beta_{Y},\beta_{X})$ in Assumption D quantify the smoothness of the submanifold $\mathcal{M}_{Y|x}$, which is the image (range) of the decoder $G^{\omega}_{x}(z)$. Specifically, index $\beta_{Y}$ characterizes the smoothness level of the manifold $M_{Y|x}$ supporting the response variable for any fixed $x\in\mathcal{M}_{X}$; and index $\beta_{X}$ characterizes the smoothness level of the section manifold $M_{Y|x}$ in $x$, that is, how similar $\mathcal{M}_{Y|x}$ and $\mathcal{M}_{Y|x^{\prime}}$ are when $x$ is close to $x^{\prime}$ in $\mathcal{M}_{X}$. In contrast, constants $(\alpha_{Y},\alpha_{X})$ in Assumption D quantify the smoothness of the conditional distribution $\mu^{*}_{Y|x}$, or more precisely, the corresponding conditional density function on its supporting manifold $\mathcal{M}_{Y|x}$. Specifically, the index $\alpha_{Y}$ characterizes the smoothness level of the conditional density function in the response variable $Y$ for any fixed $x\in\mathcal{M}_{X}$, while the index $\alpha_{X}$ captures how smoothly the conditional density function changes with $x$. The constant $d_{Y}$ in the assumption can be viewed as the intrinsic dimension of the response variable $Y$.
Similar to [], the requirements on the conditional distribution $\mu^{*}_{Y|x}$ in Assumption D are stated in a local manner since a manifold, as a topological space, is only locally defined. In fact, many common manifolds, such as spheres, do not admit a global parameterization (or encoder-decoder representation).

###### Theorem 2 (Distribution regression on manifolds) .

Suppose Assumptions B, C and D are satisfied with $\beta_{Y}\geq\alpha_{Y}\vee 1+1$ and $\beta_{X}\geq\alpha_{X}+\frac{\alpha_{X}}{\alpha_{Y}}$. Let $\varepsilon_{1}=n^{-1\,\big{/}\,\big{(}2\alpha_{Y}+d_{Y}+\frac{\alpha_{Y}}{
\alpha_{X}}d_{X}\big{)}}$. If we take $\tau=\widetilde{\Theta}(\varepsilon_{1}^{2(\alpha_{Y}+1)})$, $T=\Theta(\log n)$, and neural network sizes satisfying $H_{i}=\Theta(\log^{4}n)$, $\|W_{i}\|_{\infty}=\widetilde{\Theta}\big{(}n^{\frac{d_{X}}{2\alpha_{X}+d_{X}}
}t_{i}^{\frac{-\alpha_{X}d_{Y}}{2\alpha_{X}+d_{X}}}\wedge\varepsilon_{1}^{-d_{
Y}-\frac{\alpha_{Y}d_{X}}{\alpha_{X}}}\big{)}$,
$R_{i}=\widetilde{\Theta}\big{(}n^{\frac{d_{X}}{2\alpha_{X}+d_{X}}}t_{i}^{\frac
{-\alpha_{X}d_{Y}}{2\alpha_{X}+d_{X}}}\wedge\varepsilon_{1}^{-d_{Y}-\frac{
\alpha_{Y}d_{X}}{\alpha_{X}}}\big{)}$, $B_{i}=\exp(\Theta(\log n^{4}))$ and $V_{i}=\Theta\big{(}\sqrt{\frac{\log n}{t_{i}\wedge 1}}\big{)}$ for $i\in[\mathcal{I}]$ with $\mathcal{I}=\log_{2}(\frac{T}{\tau})$, then it holds with probability at least $1-n^{-1}$ that

$$ $\mathbb{E}_{x\sim\mu^{*}_{X}}\big{[}W_{1}(\widehat{\mu}_{Y|x},\,\mu^{*}_{Y|x}) \big{]}=\widetilde{\mathcal{O}}\bigg{(}n^{-\frac{1}{2+\frac{d_{X}}{\alpha_{X}} }}\vee n^{-\frac{1+\frac{1}{\alpha_{Y}}}{2+\frac{d_{X}}{\alpha_{X}}+\frac{d_{Y }}{\alpha_{Y}}}}\bigg{)}.$ $$

###### Remark 5 .

Since $\widehat{\mu}_{Y|x}$ and $\mu^{*}_{Y|x}$ are almost surely mutually singular measures (supporting on different submanifolds), the total variation metric is always $1$ and, therefore, not suitable for quantifying their closeness. Consequently, our error bound is stated only in terms of the Wasserstein metric. In fact, even in the (unconditional) distribution estimation case, [] shows that no estimator can achieve estimation consistency under the total variation metric.

###### Remark 6 .

Theorem 2 demonstrates that the statistical accuracy of the conditional diffusion model depends solely on the intrinsic dimensions $(d_{X},d_{Y})$ rather than the ambient dimensions $(D_{X},D_{Y})$, modulo multiplicative constants and logarithmic terms. This indicates that the conditional diffusion model can adapt to the low-dimensional manifold structures in both the response and the covariate variables.
In particular, when there is no low-dimensional manifold structures (i.e., $D_{X}=d_{X}$ and $D_{Y}=d_{Y}$), the $W_{1}$ error bound in Theorem 2 recovers the $W_{1}$ error bound in Theorem 1 in the classical density regression.

###### Remark 7 .

The same remarks after Theorem 1 in the previous subsection also apply: when specializing to the unconditional case with no covariate (that is, taking $D_{X}=0$ in Theorem 1), our error bound reduces to the minimax rate of (unconditional) distribution estimation under the Wasserstein metric [] with a (sufficiently smooth) manifold structure; when specializing to the mean regression case, our error bound can recover the classical convergence rate $n^{-\frac{\alpha_{X}}{2\alpha_{X}+d_{X}}}$ of nonparametric regression when the covariate $X$ is supported on a $d_{X}$-dimensional submanifold [, ].

###### Remark 8 .

Several works [, ] have also studied the unconditional diffusion model with a low-dimensional structure, where the data lies in a subspace. However, our work addresses a more general setting in which the manifold is unknown and can be highly nonlinear. In particular, for a linear subspace, the (unconditional) diffusion process can be decomposed into the tangent part and orthogonal part, so that the subspace estimation error and the estimation error of the distribution on the subspace can be decoupled and analyzed separately. In comparison, for a nonlinear manifold, such a decomposition does not exist, and the manifold estimation error and distribution estimation error are coupled in a complicated manner. Furthermore, due to the nonlinearity, in our approximation error analysis using neural networks, we have to locally approximate a class of projection operators of the nonlinear manifold that changes cross the manifold, rather than approximating a single global projection operator onto a linear subspace.

###### Theorem 2 (Distribution regression on manifolds) .

###### Remark 5 .

###### Remark 6 .

###### Remark 7 .

###### Remark 8 .

### 3.3 Proof highlights

Since Theorem 2 extends Theorem 1 by incorporating manifold structures, we will only outline the proof for the former in this subsection.
All missing definitions, formal assumptions, and detailed proofs are provided in the appendices of the supplementary material for this paper.

Our strategy for bounding the distribution estimation error mainly follows the pipeline of [, ]. First, we construct a specific neural network within the class $\mathcal{S}_{NN}$ to approximate the true conditional score function $\nabla\log p_{t}(\cdot,|,x)$ with controlled error, which is summarized in the following lemma.

###### Lemma 1 (Score approximation error by neural network class) .

Under the same neural network sizes $\{(H_{i},W_{i},R_{i},B_{i},V_{i})\}_{i=1}^{\mathcal{I}}$ and choices of $\tau$, $T$ as in Theorem 2, there exists neural network $\phi_{i}(w,x,t)\in\Phi(H_{i},W_{i},R_{i},B_{i},V_{i})$ for any $i\in[\mathcal{I}]$ so that

$$ $\displaystyle\mathbb{E}_{\mu^{*}_{X}}\bigg{[}\int_{t_{i-1}}^{t_{i}}\int_{ \mathbb{R}^{D}}\big{\|}\phi_{i}(w,x,t)-\nabla\log p_{t}(\cdot\,|\,x)(w)\big{\| }^{2}p_{t}(\cdot\,|\,x)(w)\,{\rm d}w\,{\rm d}t\bigg{]}$ $\displaystyle=\begin{cases}\widetilde{\mathcal{O}}\bigg{(}n^{-\frac{2\alpha_{Y }}{2\alpha_{Y}+d_{Y}+\frac{\alpha_{Y}}{\alpha_{X}}d_{X}}}+t_{i}^{-1}n^{-\frac{ 2\cdot(\beta_{Y}\wedge\frac{\beta_{X}\alpha_{Y}}{\alpha_{X}})}{2\alpha_{Y}+d_{ Y}+\frac{\alpha_{Y}}{\alpha_{X}}d_{X}}}\bigg{)},&\mbox{if }\tau\leq t_{i}<n^{- \frac{2}{2\alpha_{Y}+d_{Y}+d_{X}\frac{\alpha_{Y}}{\alpha_{X}}}};\\ \widetilde{\mathcal{O}}\Big{(}n^{-\frac{2\alpha_{X}}{2\alpha_{X}+d_{X}}}t_{i}^ {-\frac{\alpha_{X}d_{Y}}{2\alpha_{X}+d_{X}}}\Big{)},&\mbox{if }n^{-\frac{2}{2 \alpha_{Y}+d_{Y}+d_{X}\frac{\alpha_{Y}}{\alpha_{X}}}}\leq t_{i}\leq T.\end{cases}$ $$

The primary technical challenge in proving Lemma 1 arises from the fact that the space $\mathcal{M}_{Y|x}$ is a general (possibly nonlinear) manifold that depends on $x$. To address this, we partition the joint space $\mathcal{M}$ of $(X,Y)$ into small pieces with varying resolution levels in $X$ and $Y$, tailored to the smoothness levels $(\alpha_{Y},\alpha_{X})$, dimensions $(d_{Y},d_{X})$, and times $t_{i}$. Within each pieces, we carefully construct local polynomials to approximate the local charts (i.e., the decoders $G^{\omega}_{x}(z)$ defined in Assumption D) of $\mathcal{M}_{Y|x}$. The actual proof is much more involved and delicate in order to optimally balance between the approximating neural network size and the approximate error; see Appendix C for details.
Based on Lemma 1, we can now utilize the complexity of $\mathcal{S}_{NN}$ to control the generalization error for our conditional score estimator $\widehat{S}$, which minimizes the empirical score matching risk (8). The result is summarized as follows.

###### Lemma 2 (Score matching generalization error) .

It holds with probability at least $1-n^{-1}$ that,

$$ $\displaystyle\mathbb{E}_{\mu^{*}_{X}}\bigg{[}\int_{t_{i-1}}^{t_{i}}\int_{ \mathbb{R}^{D_{Y}}}\big{\|}\widehat{S}(w,x,t)-\nabla\log p_{t}(\cdot\,|\,x)(w) \big{\|}^{2}\,p_{t}(\cdot\,|\,x)(w)\,{\rm d}t\,{\rm d}w\bigg{]}$ $\displaystyle\lesssim\underset{S\in\mathcal{S}_{i}}{\min}\mathbb{E}_{\mu^{*}_{ X}}\bigg{[}\int_{t_{i-1}}^{t_{i}}\int_{\mathbb{R}^{D_{Y}}}\|S(w,x,t)-\nabla \log p_{t}(\cdot\,|\,x)(w)\|^{2}\,p_{t}(\cdot\,|\,x)(w)\,{\rm d}t\,{\rm d}w \bigg{]}$ $\displaystyle\qquad\qquad\qquad+\frac{R_{i}H_{i}\log\big{\{}R_{i}H_{i}\|W_{i} \|_{\infty}(B_{i}\vee 1)\,n\big{\}}\,(\log n)^{2}}{n},\quad\mbox{ for each $i \in[\mathcal{I}]$.}$ $$

[] derived a similar result in the context of (unconditioned) distribution estimation without $x$; in addition, their error bound is not a high probability bound but instead takes another expectation with respect to the randomness of $\widehat{S}$. In contrast, our proof for the conditional distribution estimation requires a *high probability bound*. We utilize more technical tools in empirical process theory, such as the localization and peeling techniques [], to derive such a high probability bound as in Lemma 2. The rest of the analysis is similar to a standard analysis for score-based diffusion models [, , ], where we apply Girsanov’s theorem to relate the distribution estimation error with the obtained $L_{2}$ score estimation error.

###### Lemma 1 (Score approximation error by neural network class) .

###### Lemma 2 (Score matching generalization error) .

## 4 Conclusion

In this study, we investigate the theoretical properties of conditional forward-backward diffusion estimators within the statistical framework of distribution regression. Our results identify the primary sources of error in conditional distribution estimation using conditional diffusion models and include earlier results on unconditional distribution estimation and nonparametric mean regression as special cases. Notably, our findings demonstrate that although (conditional) diffusion models operate directly in the original ambient data space and do not explicitly incorporate low-dimensional structures, the resulting conditional distribution estimators can still adapt to intrinsic manifold structures when either (or both) the covariate
$X$ and response $Y$ are concentrated around a lower-dimensional manifold.
Our analysis also offers practical guidance for designing the neural network approximation family to optimally control different types of errors. This includes recommendations for the architecture of the neural network, as well as how the network’s size (depth, width, sparsity, etc.) should depend on various problem characteristics, such as sample size, smoothness levels, and intrinsic dimensions.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 8
presentation: 3
questions: Please see the weaknesses section.
soundness: 3
strengths: - The contributions are, while possibly not surprising, are meaningful and interesting results extending the theoretical analysis of diffusion models.
- The paper is generally very well written and easy to follow, given the technical nature of the content. In particular the remarks on the theorems make understanding the implications of them more accessible.
summary: Previous work has studied the min-max optimality of score-based diffusion models in estimating the density of an unknown distribution, and in estimating the intrinsic dimension of data lying on a submanifold of a sample space. 

This work extends these results to the conditional modelling setting, where a conditioning variate changes the distribution being modeled. The distribution is assumed to vary smoothly with the conditioning variate.

These results show that very similar rates of convergence of the total variation and Wasserstein distances between the true and  estimated data distributions hold to the unconditional modeling setting, with additional terms depending on the dimensionality of the conditioning variate. In the limit of the conditioning variate having dimension zero, the same rates as previous works on the unconditional setting are obtained.
weaknesses: - The practical implications of the theoretical work presented are not completely clear to me. The authors mention that the results suggest how to design networks, but from my understanding on theorem 2, this suggests simply that larger networks produce better bounds?
- The implications of assumption D are unclear to me. The variable $r_0$ is undefined. Is this condition meant to hold for any $r_0$? If so then this implies a need for the decoding to be smooth globally, with I believe would limit topological structure of the low-dimension data manifold. If this is local condition, the is this not a consequence of the smooth manifold structure of $\mathcal{M}$?

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
questions: 1) Are the covariates X also located in the Besov space support?
2) The marginals M_X and M_Y are assumed to be  [-1,1]  D-X, resp. D_Y dimensional cube. How does the scaling of the data to the [-1, 1] cube impact some other intrinsic features, e.g. latent translations? What type of scaling transformations would you recommend?
3) Can you provide guidelines on how you estimate $\alpha_X$ and $\alpha_Y$ or do you assume these values?
4) How does presented results in Theorem 2 relate to latent diffusion models?

5) p.7 Remark 1: Does this remark imply that the kernel estimator can well approximate the performance of the conditional diffusion based estimator?

6) It would be great if authors can  include experiments for different embedded lower dimensional structures, and demonstrate how the presented theory supports and navigates the experimentally obtained results. Using even simple synthetic geometric objects would be great.

7) Can you provide more discussion on the implications of the different regimes discussed in Remark 3? In particular how would you select the smoothness parameters $\alpha_x$ to investigate the transition phases?

8 )   How does Lemma 2 helps to assess the generalization error?
soundness: 4
strengths: - The paper is mathematically well written. The notation is consistent and precise. Definitions and proofs are provided in the extensive Appendix. These can be useful for deriving and investigating other novel properties and algorithms.
- The problem is relevant and well chosen and the previous work is well mapped. 
- There is an interesting potential applications for measuring the quality of the conditional diffusion estimators and guided  sampling. The paper also  provides a background for the reasoning about the selection of the model architecture.
- The interesting contribution is in analyzing the latent dimension representations of the data and the adaptability of the higher dimensional estimator.
- The results provide insight into the impact of the smoothness of the underlying data manifold on the performance of the estimator.
summary: - The paper presents a derivation of the estimation error bound for the diffusion models with sampling guided by the covariate regressor. The authors closely build upon the methodology from Oko et al 2023 https://proceedings.mlr.press/v202/oko23a/oko23a.pdf who explored the statistical learning theory for diffusion models via score matching estimators and proved nearly minimax optimal estimation rates in total variation and Wasserstein metric of order 1 for fully connected neural-networks (including suggestions on generalisation to U-nets and ConvNets). 

- The presented paper aims to introduce the statistical learning theory for conditional diffusion models using the distribution regression of the covariate structure guiding the sampling. Authors investigate the properties of the estimators using the finite sample analysis of the statistical error bound under various metrics, the relations of latent lower dimensional data structures and the smoothness properties of the data. 

- The main contribution is: 
   1) Given the assumptions on the size of the architecture,  the estimation error bound of the conditional diffusion estimator is minimax optimal under total variation and Wasserstein distance of order 1, authors include the discussion on the generalization to the unconditional case as a special case, 
   2) if there is a low-dimensional latent space structure within the high-dimensional data manifold (say within the feature space or feature and label space), the resulting conditional distribution estimators adapts to the lower-dimensional latent space.
weaknesses: 1) Paper is missing motivational examples and practical experiments to demonstrate how the presented theory can be used on real-life problems. The results are theoretically interesting and there is no issue with the soundness of the arguments. Nevertheless, without a practical demonstration on how to use the derived results, the main contribution of the paper may be harder to access to the wider machine learning community.

2) While it is great to see the methodology related to the proofs as a part of the main paper, it would be great to see particular direct applications in the main body of the paper.

3) The usual approach to estimate the gradients of log p_t(y_t|y) (see line 200) is to use the gradient of unconditional model and add gradient of log-likelihood of the ‘reverse predictor' which detects ‘the covariate values’ x based on the features constructed from whom the unconditional values ‘y_t’. Authors decided to use marginal forward diffusion which results in empirical estimation of the joint density between the covariate X (the guiding covariates) and Y (the unconditional diffusion). The discussion on incorporating computational efficiency of this approach in comparison to the 'reverse predictor' approach is missing.
4) The proven Theorems assume the initial constraints of the sizes of the architecture such as depth, width, norm of the parameters, etc. These constraints are set as functionals of the order of magnitude of the data sample $n$. It is not well discussed what happens if these constraints are not satisfied.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: What is p_data on Line 122?

Authors need to check all the parenthetical and narrative citations throughout the whole paper. Many were mis-used. For example, Line 316 should be parenthetical citation.

Definition of neural network class is already given in Tang and Yang, 2024 (AISTATS). Authors need to reference it or directly call it from the citation.

In remark 1, can authors discuss the case when alpha_X>1? What does this mean to the minimax rate?

In Remark 2, can the authors clarify why D_Y need to satisfy this condition to have the first term dominate the second one? 




Rong Tang and Yun Yang. Adaptivity of diffusion models to manifold structures. 27th International Conference on Artificial Intelligence and Statistics, 2024.
soundness: 3
strengths: 1. Authors introduce SDE in a simple way and made it easy to connect SDE to distribution regression.

2. The paper extends the unconditional distribution estimator to the conditional one and provides theoretical results. Theories are given in an elegant way that they clearly explain the connections to the theories of unconditional distribution estimator, mean regression problems, and nonparametric regression under l2 error. Theoretical results are interesting and inspiring.
summary: This paper uses distribution regression to estimate the conditional distribution induced by conditional forward-backward diffusion models. Authors utilize neural network to estimate the key points and also allow for high-dimensional X and Y with low-dimensional structures. Solid theoretical results are provided.
weaknesses: The paper mentions that it provides practical guidance to design the neural network, including the network size, smoothness level and so on. Though it is a theoretical work, some empirical results would significantly demonstrate the usefulness of the paper. I would recommend authors run some simulations to show how the theoretical results guide the design of neural networks. For example, under what scenarios, the errors are optimally controlled?

The theoretical results are interesting.  Can the authors highlight the difficulties and summarize the new ways of proof if there are? This could be a huge contribution to the statistical learning community. I didn't see the challenges from my end and all seem to be standard.

Meta-review: This paper studies conditional diffusion models for generative modeling, focusing on generating data given a covariate. Using a statistical framework of distribution regression, the authors analyze the large-sample properties of conditional distribution estimators under the assumption that the conditional distribution changes smoothly with the covariate. They derive a minimax-optimal convergence rate under the total variation metric, consistent with existing literature. The analysis is further extended to scenarios where both data and covariates lie on low-dimensional manifolds, showing that the model adapts to these structures. The resulting error bounds under the Wasserstein metric depend only on the intrinsic dimensionalities of the data and covariates. This paper has received unanimous support from the reviewers. Therefore, I recommend acceptance.

Decision: Accept (Poster)