Smoothness-Based Rates for Conditional Density Estimation

Why Smoothness Determines the Rate

A conditional density estimator must learn how the outcome distribution changes with the covariates and must also resolve variation across the outcome space. These are distinct smoothing tasks. The covariate dimension determines how many observations contribute locally; the outcome dimension determines how sharply each observation is smoothed as a density. Smoothness in each set of coordinates determines how much bias those choices introduce.

We study rates in integrated squared error. This norm averages error over both covariate and outcome values, so it is not the same criterion as a pointwise error at one fixed pair. The distinction matters: integrating squared kernel fluctuations over the outcome domain still leaves an inverse power of the outcome bandwidth.

Setup and Estimator

Let $X\in[0,1]^d$ and $Y\in[0,1]^q$, with the domains treated as tori to avoid boundary effects. We observe $n$ independent copies $(X_i,Y_i)$. Write $p_X(x)$ for the covariate density and $f(y\mid x)$ for the conditional density of $Y$ given $X=x$. Assume $0<c_p\le p_X(x)\le C_p<\infty$ and $0\le f(y\mid x)\le B<\infty$.

Let $K_x$ and $K_y$ be bounded, compactly supported, nonnegative kernels on the covariate and outcome spaces, respectively, each integrating to one. For bandwidths $h_x,h_y\in(0,1]$, define

\widehat D_{h_x}(x) = \frac{1}{n h_x^d}\sum_{i=1}^n K_x\!\left(\frac{x-X_i}{h_x}\right), \qquad \widehat N_{h_x,h_y}(x,y) = \frac{1}{n h_x^d h_y^q}\sum_{i=1}^n K_x\!\left(\frac{x-X_i}{h_x}\right) K_y\!\left(\frac{y-Y_i}{h_y}\right).

The raw conditional density estimate is $\widehat N_{h_x,h_y}(x,y)/\widehat D_{h_x}(x)$ when the denominator is sufficiently large. To make the risk argument valid even when it is not, define the estimate to be zero wherever $\widehat D_{h_x}(x)<c_p/2$, and then clip the result to $[0,B]$. Denote this final estimator by $\widetilde f(y\mid x)$. Clipping is pointwise projection onto a range known to contain the true density; it need not preserve exact integration to one.

Definition: The integrated squared-error risk considered here is $$ \mathcal R(\widetilde f) = E\!\left[ \int_{[0,1]^d}\int_{[0,1]^q} \{\widetilde f(y\mid x)-f(y\mid x)\}^2\,dy\,dx \right]. $$

Separate Smoothness: The Main Bound

Suppose $f(y\mid x)$ and $p_X(x)$ have separate Hölder smoothness orders $0<s_x,s_y\le2$ in $x$ and $y$, respectively, with bounded Hölder constants. Use the stated nonnegative kernels with zero first moments in each coordinate whose smoothness order exceeds 1; for orders at most 1, no moment cancellation beyond integral one is needed. For $0<s\le2$, these conditions cancel the nonconstant Taylor terms of degree less than $s$, and the usual Hölder Taylor remainder bounds smoothing bias by the corresponding bandwidth to that smoothness power. Assume $n h_x^d\to\infty$.

Theorem: Under the boundedness, positivity, smoothness, and kernel conditions just stated, there is a constant $C$, independent of $n,h_x,h_y$, such that $$ \mathcal R(\widetilde f) \le C\left\{ h_x^{2s_x}+h_y^{2s_y} +\frac{1}{n h_x^d h_y^q} +\exp(-c n h_x^d) \right\}, $$ for some $c>0$. The exponential term controls the small-denominator event; in particular, this bound does not require an unproved claim that clipping alone makes that event negligible.

Here is the argument. Let $D(x)=E[\widehat D_{h_x}(x)]$, $N(x,y)=E[\widehat N_{h_x,h_y}(x,y)]$, and $r_h(x,y)=N(x,y)/D(x)$. Periodization and the lower bound on $p_X$ imply $D(x)\ge c_p$. The population ratio is a locally weighted average of $f$: its numerator averages $p_X(u)f(v\mid u)$, while its denominator averages $p_X(u)$. Taylor expansion in the $x$ and $y$ coordinates, with the kernel moment cancellations, gives

\int\!\!\int \{r_h(x,y)-f(y\mid x)\}^2\,dy\,dx \le C\{h_x^{2s_x}+h_y^{2s_y}\}.

To bound the random fluctuation, define the good event at $x$ by $G_x=\{\widehat D_{h_x}(x)\ge c_p/2\}$. On $G_x$, subtraction of the two ratios gives the identity

\frac{\widehat N}{\widehat D}-\frac{N}{D} = \frac{\widehat N-N}{\widehat D} + \frac{N(D-\widehat D)}{\widehat D D}.

On this event both denominators are bounded below. Squaring the identity and using $(a+b)^2\le 2a^2+2b^2$ reduces the integrated fluctuation to the integrated variances of $\widehat N$ and $\widehat D$. Independence of the observations gives

\int\!\!\int \operatorname{Var}\{\widehat N(x,y)\}\,dy\,dx \le \frac{C}{n h_x^d h_y^q}, \qquad \int \operatorname{Var}\{\widehat D(x)\}\,dx \le \frac{C}{n h_x^d}.

For the first inequality, the variance of an average is at most $n^{-1}$ times the second moment of one summand. Integrating the squared covariate kernel contributes a factor of order $h_x^d$, and integrating the squared outcome kernel contributes a factor of order $h_y^q$; the estimator's squared normalization contributes $h_x^{-2d}h_y^{-2q}$. Thus the resulting order is $1/(n h_x^d h_y^q)$. The denominator term is multiplied by the bounded population ratio $r_h$, and is of order $1/(n h_x^d)$, which is no larger than $1/(n h_x^d h_y^q)$ because $h_y\le1$. These calculations bound the good-event contribution.

It remains to handle the event omitted from that calculation. For each fixed $x$, the summands defining $\widehat D_{h_x}(x)$ are bounded, their mean is at least $c_p$, and their variance is of order at most $h_x^{-d}$ after the displayed normalization. A Bernstein inequality therefore gives $P(G_x^c)\le C\exp(-c n h_x^d)$. Crucially, the clipped estimate and the true density both lie in $[0,B]$. Hence their squared difference is at most $B^2$, regardless of $h_y$ and regardless of the raw ratio's size. Integrating over the unit-volume $x,y$ domains bounds the exceptional-event risk by $B^2 C\exp(-c n h_x^d)$. Finally, clipping cannot increase squared error when $f\in[0,B]$, because projection onto that interval moves any value no farther from $f$. Combining the bias, good-event, and exceptional-event bounds proves the theorem. $\square$

Choosing Bandwidths and Reading the Rate

The three main terms in the theorem express a direct trade-off: larger bandwidths increase smoothing bias, while smaller bandwidths reduce the number of observations and increase variance. Ignoring the exponentially small term, balance the two bias terms and the variance term by writing $h_x=n^{-a}$, $h_y=n^{-b}$. If all three terms have order $n^{-\rho}$, then

2s_xa=\rho,\qquad 2s_yb=\rho,\qquad 1-da-qb=\rho.

Substituting $a=\rho/(2s_x)$ and $b=\rho/(2s_y)$ into the last equation yields

\rho= \left(1+\frac{d}{2s_x}+\frac{q}{2s_y}\right)^{-1}, \qquad h_x\asymp n^{-\rho/(2s_x)}, \qquad h_y\asymp n^{-\rho/(2s_y)}.

These bandwidths make $n h_x^d$ grow as a positive power of $n$, so the exponential remainder is smaller than the displayed polynomial rate. This is an upper-bound rate under the stated assumptions, not a universal claim of minimax optimality for every model class.

Worked Example: Separate Smoothness in Unequal Dimensions

Suppose $d=2$, $q=1$, $s_x=2$, and $s_y=1$. Then

\rho = \left(1+\frac{2}{4}+\frac{1}{2}\right)^{-1} = \frac12.

Thus $h_x\asymp n^{-1/8}$, $h_y\asymp n^{-1/4}$, and the risk bound has polynomial order $n^{-1/2}$. Checking the terms: $h_x^{2s_x}=n^{-4/8}=n^{-1/2}$; $h_y^{2s_y}=n^{-2/4}=n^{-1/2}$; and $1/(n h_x^d h_y^q)=1/(n\cdot n^{-2/8}\cdot n^{-1/4})=n^{-1/2}$. For $n=4096=2^{12}$, these bandwidths are $h_x=2^{-3/2}=1/\sqrt 8$ and $h_y=2^{-3}=1/8$. The three polynomial terms are each $1/64$, confirming the balance numerically.

Joint Smoothness and Dimension

A common special case assumes isotropic smoothness of order $0<s\le2$ jointly across all $d+q$ coordinates. If the joint density $p_X(x)f(y\mid x)$ and the covariate density are sufficiently smooth, and the product kernel has zero first moments when $s>1$, the population conditional smoother has bias of order $h^s$ when $h_x=h_y=h$. The same variance calculation then gives the familiar integrated-risk balance

\mathcal R(\widetilde f) \lesssim h^{2s}+\frac{1}{n h^{d+q}}+\exp(-c n h^d).

Balancing the first two terms yields $h\asymp n^{-1/(2s+d+q)}$ and risk order $n^{-2s/(2s+d+q)}$. The covariate dimension still appears in the denominator event, whose probability decays exponentially when $n h^d$ grows. This exponential term is not an additional polynomial dimension penalty in the balanced bound.

Worked Example: Joint Smoothness Rate

Let $d=2$, $q=1$, and suppose the joint smoothness order is $s=2$. Then $d+q=3$, so the balanced bandwidth and risk orders are

h\asymp n^{-1/(4+3)}=n^{-1/7}, \qquad \mathcal R(\widetilde f)\lesssim n^{-4/7}.

Indeed, $h^{2s}=h^4=n^{-4/7}$, and $1/(n h^{d+q})=1/(n h^3)=n^{-1+3/7}=n^{-4/7}$. Also $n h^d=n^{1-2/7}=n^{5/7}$, so the denominator-event term is exponentially small in $n^{5/7}$, not merely assumed negligible.

What the Norm and Dimensions Mean

The integrated squared-error norm averages over the entire outcome domain, but it does not erase the outcome dimension from the variance. A single outcome-kernel estimate has squared $L^2$ size of order $h_y^{-q}$; this is why the integrated variance contains $h_y^{-q}$, as well as the covariate factor $h_x^{-d}$. By contrast, a pointwise variance also has order $1/(n h_x^d h_y^q)$, although pointwise risk and integrated risk are different quantities and need separate statements.

The separate-smoothness rate makes the roles explicit. Greater covariate dimension $d$ increases the cost of localizing observations in $x$; greater outcome dimension $q$ increases the cost of resolving a density in $y$. Greater smoothness in either coordinate set allows a larger bandwidth there, reducing that contribution to the variance. If prediction error is averaged under the covariate distribution rather than Lebesgue measure, the upper and lower bounds on $p_X$ make the two integrated risks equivalent up to constants.

Worked Example: Why the Outcome Bandwidth Appears in Integrated Variance

Take a one-dimensional outcome and a normalized kernel $K_y$. For one observation, the outcome factor is $K_y((y-Y_i)/h_y)/h_y$. Its squared integral over $y$ is

\int \left\{\frac{1}{h_y}K_y\!\left(\frac{y-Y_i}{h_y}\right)\right\}^2dy = \frac{1}{h_y}\int K_y(u)^2\,du,

where the equality follows from the change of variable $u=(y-Y_i)/h_y$. Thus even after integrating over $y$, the squared size is proportional to $h_y^{-1}$. With $q$ outcome coordinates, the factor becomes $h_y^{-q}$. Combining it with the covariate localization factor gives $1/(n h_x^d h_y^q)$. Treating the integrated variance as $1/(n h_x^d)$ would omit the cost of estimating a density over the outcome space.

Common Mistake: A Small Failure Probability Is Not Enough

A tempting argument is to show that the denominator is small with exponentially small probability and then declare the ratio estimator's risk contribution negligible. That conclusion does not follow from the probability bound alone: on the failure event, a ratio may be very large, with a size involving inverse powers of $h_y$. Clipping only the denominator keeps division defined; by itself it does not bound the numerator or the ratio. A valid proof must also control the loss on that event. Here, defining the estimate as zero below the denominator threshold and clipping the output to $[0,B]$ bounds squared loss by $B^2$ uniformly over all allowed $h_y$. The exceptional-event contribution is therefore bounded by a constant times its probability, with no hidden outcome-bandwidth factor.

Key takeaway: Under separate smoothness, integrated squared-error risk has bias terms $h_x^{2s_x}+h_y^{2s_y}$ and variance of order $1/(n h_x^d h_y^q)$. The rate depends on both dimensions, and a small-denominator event must be controlled through an actual loss bound, not just a tail probability.

Check Your Understanding

Use the assumptions and rate calculations in this tutorial to answer the following questions.

For $d=1$, $q=2$, $s_x=1$, and $s_y=2$, find $\rho$ and the balanced bandwidth exponents.
Why does integrating the squared outcome-kernel factor over $y$ leave a factor proportional to $h_y^{-q}$?
Under joint smoothness $s$ with equal bandwidths, what balances the squared bias and integrated variance?
Which assumption allows the exceptional-event squared loss to be bounded independently of $h_y$?
Explain why a bound on $P(G_x^c)$, without a bound on squared loss over $G_x^c$, does not by itself establish an integrated risk rate.