Score Matching for Conditional Density Estimation

Why Match Scores Instead of Densities?

A conditional density model describes how the distribution of an outcome $Y$ changes with predictors $X$. A direct density-based loss may require evaluating the model density at observed outcomes, including its normalizing constant. That constant can be difficult to calculate when the model is specified only up to proportionality. Score matching offers another route: it compares gradients of log densities with respect to the outcome. Since a normalizing constant that depends on $x$ but not on $y$ disappears under this gradient, it need not be evaluated.

The key step is an integration-by-parts identity. It turns a loss involving the unknown true conditional score into an objective involving only the candidate model score and its outcome divergence. The identity is valid only under appropriate smoothness, integrability, common-support, and boundary conditions. In particular, boundary terms cannot be discarded merely because the densities are normalized.

Conditional Scores and the Objective

Definition: Let $p(y\mid x)$ be a conditional density on a region $\Omega_x\subseteq\mathbb{R}^d$, with respect to $d$-dimensional Lebesgue measure. Where $p(y\mid x)>0$ and is differentiable in $y$, its conditional score with respect to the outcome is $$ s_p(y,x)=\nabla_y\log p(y\mid x). $$ For a candidate conditional density $q(y\mid x)$, define $s_q(y,x)=\nabla_y\log q(y\mid x)$. The divergence of this score is $\operatorname{div}_y s_q(y,x)=\sum_{j=1}^d\partial s_{q,j}(y,x)/\partial y_j$.

The gradient here is with respect to $y$, not $x$. Suppose the candidate density is specified as $q(y\mid x)=h(y\mid x)/Z(x)$, where $h$ is positive and $Z(x)=\int_{\Omega_x}h(y\mid x)\,dy$ is finite and positive. Since $Z(x)$ does not vary with $y$,

s_q(y,x)=\nabla_y\log h(y\mid x).

Thus the score can be evaluated without calculating $Z(x)$. This advantage does not remove the need for $h$ to define a valid density: its normalizing integral must exist, and the candidate and true distributions must have the support required by the score-matching identity.

Conditional score-matching identity: For a fixed $x$, let $p$ be the true conditional density and $q$ a candidate density on the same region $\Omega_x$. Define the Fisher divergence $$ D_F(p\|q\mid x)=\frac12\int_{\Omega_x}p(y\mid x)\|s_p(y,x)-s_q(y,x)\|^2\,dy $$ and the score-matching objective $$ J(q\mid x)=\int_{\Omega_x}p(y\mid x)\left\{\frac12\|s_q(y,x)\|^2+\operatorname{div}_y s_q(y,x)\right\}\,dy. $$ If integration by parts is valid and the oriented boundary integral $\int_{\partial\Omega_x}p(y\mid x)s_q(y,x)\cdot n(y)\,dS(y)$ is zero, then $$ D_F(p\|q\mid x)=J(q\mid x)+\frac12\int_{\Omega_x}p(y\mid x)\|s_p(y,x)\|^2\,dy. $$

Here $n(y)$ is the outward unit normal and $dS$ is surface measure. The last term in the identity depends on the true distribution but not on $q$. Consequently, minimizing $J(q\mid x)$ over candidates is equivalent to minimizing the Fisher divergence, provided the stated boundary condition holds. For random predictors, the conditional objective is averaged over the distribution $\mu$ of $X$: $J(q)=\int J(q\mid x)\,d\mu(x)$. The same identity then holds after integration over $x$, when the integrals are finite and the conditions hold for $\mu$-almost every $x$.

Deriving the Identity and Its Conditions

For the boundary-free identity, a sufficient setting is a region $\Omega_x$ with piecewise smooth boundary, a positive density $p(\cdot\mid x)$ that is continuously differentiable in the outcome, and a candidate score $s_q$ that is continuously differentiable there. Require the squared scores, the divergence term, and the boundary integral to be integrable. The divergence theorem or its one-dimensional integration-by-parts version must apply. These hypotheses can be weakened, but they cannot simply be omitted.

To see what happens when the boundary contribution is not zero, define $$ B(x)=\int_{\partial\Omega_x}p(y\mid x)s_q(y,x)\cdot n(y)\,dS(y). $$ Under the same integration-by-parts regularity conditions, but without setting $B(x)$ to zero, the exact relation is

D_F(p\|q\mid x)=J(q\mid x)-B(x)+\frac12\int_{\Omega_x}p(y\mid x)\|s_p(y,x)\|^2\,dy.

This version makes the role of the boundary explicit. The usual objective is equivalent to Fisher divergence up to a candidate-independent constant only when $B(x)=0$, or when a boundary contribution is otherwise known to be candidate-independent over the model class.

Worked Example: Integration by Parts on a Bounded Interval

Let $Y$ have a positive continuously differentiable density $p$ on $(a,b)$, and let $q$ be a positive twice continuously differentiable candidate density there. In one dimension, $s_q(y)=q'(y)/q(y)$. The boundary integral is the oriented endpoint expression $$ B=p(b)s_q(b)-p(a)s_q(a). $$ Expanding the squared-score divergence gives $$ D_F(p\|q)=\frac12\int_a^b p(s_p-s_q)^2\,dy =\frac12\int_a^b ps_q^2\,dy-\int_a^b p's_q\,dy+\frac12\int_a^b ps_p^2\,dy, $$ because $ps_p=p'$. Integration by parts yields $$ \int_a^b p's_q\,dy=[ps_q]_a^b-\int_a^b p s_q'\,dy =B-\int_a^b p s_q'\,dy. $$ Substitution therefore gives $$ D_F(p\|q)=\int_a^b p\left(\frac12s_q^2+s_q'\right)\,dy-B+\frac12\int_a^b ps_p^2\,dy. $$ If $p(a)s_q(a)=p(b)s_q(b)=0$, then $B=0$ and the usual score-matching identity follows. If either endpoint product is nonzero, dropping it changes the objective's relationship to the Fisher divergence.

Result: Under the regularity conditions above, the conditional score-matching identity follows by expanding the squared difference and integrating the cross term by parts: $$ \begin{aligned} D_F(p\|q\mid x) &=\frac12\int_{\Omega_x}p\|s_q\|^2\,dy-\int_{\Omega_x}\nabla_y p\cdot s_q\,dy +\frac12\int_{\Omega_x}p\|s_p\|^2\,dy\\ &=\frac12\int_{\Omega_x}p\|s_q\|^2\,dy-B(x) +\int_{\Omega_x}p\,\operatorname{div}_y s_q\,dy +\frac12\int_{\Omega_x}p\|s_p\|^2\,dy\\ &=J(q\mid x)-B(x)+\frac12\int_{\Omega_x}p\|s_p\|^2\,dy. \end{aligned} $$ The second equality uses $\nabla_y p=p\,s_p$ and the divergence theorem applied to $p s_q$, for which $\operatorname{div}_y(p s_q)=\nabla_y p\cdot s_q+p\,\operatorname{div}_y s_q$. If $B(x)=0$, this proves the boundary-free identity. $\square$

Three Worked Applications

Worked Example: Gaussian Conditional Regression

Suppose $q(y\mid x)=N(\beta x,4)$, so the conditional variance is known to be $4$. Its log density, up to terms not depending on $y$, is $-(y-\beta x)^2/8$, and hence $$ s_q(y,x)=-\frac{y-\beta x}{4},\qquad \frac{\partial s_q}{\partial y}=-\frac14. $$ Take two observations, $(x_1,y_1)=(1,1)$ and $(x_2,y_2)=(2,3)$. The empirical score-matching objective, the average of the model loss at these observations, is $$ \widehat J(\beta)=\frac12\left[\frac12\left\{\left(\frac{\beta-1}{4}\right)^2+\frac14\left(\frac{2\beta-3}{4}\right)^2\right\}-\frac14-\frac14\right] =\frac{(\beta-1)^2+(2\beta-3)^2}{64}-\frac14. $$ The residual sum of squares in the numerator is $$ (\beta-1)^2+(2\beta-3)^2=5\beta^2-14\beta+10, $$ whose derivative is $10\beta-14$. Thus the objective is minimized at $\hat\beta=14/10=1.4$. Substitution gives a residual sum of squares $0.2$, so $\widehat J(1.4)=0.2/64-1/4=-0.246875$. The negative value is not a problem: the divergence-equivalent objective can be shifted by a constant. The minimizer here agrees with least squares because the score loss reduces to a positive constant times squared residuals plus a term independent of $\beta$.

Worked Example: Fitting an Unnormalized Gaussian Model

At a fixed predictor value $x$, consider the unnormalized model $h_\theta(y\mid x)=\exp(-\theta y^2/2)$ on the real line, with $\theta>0$. Its normalizing constant is $Z(x)=\sqrt{2\pi/\theta}$, but the outcome score can be obtained directly: $$ s_q(y,x)=\frac{\partial}{\partial y}\log h_\theta(y\mid x)=-\theta y,\qquad \frac{\partial s_q}{\partial y}=-\theta. $$ For observations $y_1=1$ and $y_2=2$, the average empirical objective is $$ \widehat J(\theta)=\frac12\left[\frac12\theta^2(1^2+2^2)-2\theta\right] =\frac54\theta^2-\theta. $$ Differentiating gives $\widehat J'(\theta)=\frac52\theta-1$, so the minimizer over $\theta>0$ is $\hat\theta=2/5=0.4$. The calculation never evaluates $Z(x)$; the term $-\log Z(x)$ in the log density has zero derivative with respect to $y$. This example illustrates the computational benefit, not a claim that two observations provide a generally reliable estimate.

Worked Example: A Half-Line Boundary Term That Cannot Be Dropped

Let the true density on $(0,\infty)$ be $p(y)=\lambda e^{-\lambda y}$, where $\lambda>0$, and let the candidate be $q(y)=\alpha e^{-\alpha y}$, where $\alpha>0$. Their scores are $s_p=-\lambda$ and $s_q=-\alpha$, with $s_q'=0$. Therefore $$ J(q)=\int_0^\infty \lambda e^{-\lambda y}\frac{\alpha^2}{2}\,dy=\frac{\alpha^2}{2}, \qquad D_F(p\|q)=\frac12(\lambda-\alpha)^2. $$ The outward unit normal at $0$ is $n=-1$; there is no contribution at infinity. Consequently the oriented boundary integral is $$ B=p(0)s_q(0)n=\lambda(-\alpha)(-1)=\lambda\alpha. $$ The endpoint product itself is $p(0)s_q(0)=-\lambda\alpha$, whereas the oriented boundary contribution is $+\lambda\alpha$. Equivalently, the integration-by-parts boundary term is $[p s_q]_0^\infty=0-(-\lambda\alpha)=+\lambda\alpha$. When $\alpha=\lambda$, the endpoint product is $-\lambda^2$ and the oriented boundary contribution is $+\lambda^2$.

For every $\alpha>0$, the true-score constant is $\frac12\int p s_p^2\,dy=\lambda^2/2$. The boundary-aware identity checks exactly: $$ J(q)-B+\frac{\lambda^2}{2} =\frac{\alpha^2}{2}-\lambda\alpha+\frac{\lambda^2}{2} =\frac12(\lambda-\alpha)^2=D_F(p\|q). $$ At the true candidate $\alpha=\lambda$, the Fisher divergence is zero but $J(q)=\lambda^2/2$. Minimizing $J(q)=\alpha^2/2$ over $\alpha>0$ instead drives $\alpha$ toward zero. This failure is caused by the nonzero boundary term, not by an error in the integration-by-parts identity.

From the Population Objective to Data

If $(X_i,Y_i)$, $i=1,\ldots,n$, are observations from a joint distribution with conditional density $p(y\mid x)$, the natural empirical criterion for a candidate $q_\theta$ is $$ \widehat J(\theta)=\frac1n\sum_{i=1}^n\left\{\frac12\|s_{q_\theta}(Y_i,X_i)\|^2+\operatorname{div}_y s_{q_\theta}(Y_i,X_i)\right\}. $$ This is the sample average corresponding to $E[\frac12\|s_{q_\theta}(Y,X)\|^2+\operatorname{div}_y s_{q_\theta}(Y,X)]$. It can be minimized over $\theta$ without evaluating conditional density normalizers, as long as the outcome score and its divergence can be computed. The empirical criterion is an estimation procedure; its finite-sample minimizer is not guaranteed to equal the population minimizer.

For an unbounded region such as $\mathbb{R}^d$, a common sufficient boundary requirement is that the boundary flux over expanding regions vanishes, for example $$ \lim_{R\to\infty}\int_{\partial B_R}p(y\mid x)s_q(y,x)\cdot n(y)\,dS(y)=0, $$ where $B_R$ is the ball of radius $R$, together with integrability of the terms in the identity. Pointwise decay alone is not always enough: the surface area also grows, so the integrated flux must vanish. With conditional supports that vary with $x$, the boundary condition must be checked on each relevant $\partial\Omega_x$, for almost every predictor value under consideration.

Common Pitfalls and Interpretation

Differentiating with respect to the wrong variable. Conditional score matching uses $\nabla_y\log q(y\mid x)$. Differentiating in $x$ does not remove a normalizer that depends on $x$.
Dropping a boundary term without checking it. The half-line example shows that a normalized density can have a nonzero boundary flux. The oriented contribution includes the outward normal; it is not just the endpoint product.
Assuming an unnormalized formula is automatically a density. The integral defining $Z(x)$ must be finite and positive, and the model must have the appropriate support.
Interpreting the objective's level as a divergence. The score-matching objective may be negative and is not itself the Fisher divergence. Under the boundary condition, it differs from that divergence by a term determined by the true distribution.

When the boundary-free identity holds, matching scores has a clear population interpretation: it minimizes a squared discrepancy between outcome gradients of log densities. If two positive differentiable densities have identical outcome scores on a connected support, their log-density difference is constant there; normalization then forces the densities to agree. This helps explain why score information can characterize a density, while also underscoring that support and boundary assumptions are part of the method rather than technical decoration.

Check Your Understanding

Use the definitions and boundary-aware identity to answer the following.

Why does a conditional normalizing constant $Z(x)$ disappear from the outcome score, and what derivative would fail to remove it?
State the boundary-free conditional score-matching identity and identify the term that is constant with respect to the candidate density.
On $(0,\infty)$, what are the endpoint product and the oriented boundary contribution for the exponential example when $\alpha=\lambda$?
For $h_\theta(y)=\exp(-\theta y^2/2)$, derive the outcome score and its divergence without evaluating the normalizing constant.
Why is pointwise decay of $p(y\mid x)s_q(y,x)$ not by itself a complete boundary condition on an unbounded region?