Tutorials › Statistical Decision Theory › Gaussian Score Functions and Integration by Parts

Stein Identity and Risk Improvement · Tutorial 415 of 982

Gaussian Score Functions and Integration by Parts

Learn to differentiate a multivariate normal log density and use its score in integration-by-parts calculations for Gaussian expectations.

Advanced 10 min read

What You'll Learn

  • Derive the score vector for a multivariate normal density with positive-definite covariance.
  • Distinguish the score from the centered observation and track the covariance scaling.
  • Prove an integration-by-parts relation for a spherical multivariate normal and suitable test functions.
  • Apply the relation to polynomial and trigonometric test functions.
  • Recognize why boundary behavior and integrability matter.

Why the Gaussian Score Matters

A density’s score records how quickly its log density changes as the observation moves. For a Gaussian model, that change has a particularly simple form: it points from the observation back toward the mean. Integration by parts turns this density derivative into a derivative of a test function inside an expectation. The result is a useful calculation tool, provided the test function and the Gaussian tails make the boundary terms vanish.

We first derive the score for a general nonsingular multivariate normal distribution. We then use integration by parts to obtain an expectation relation in the spherical case, where the covariance is \(\sigma^2 I_d\). Keeping these steps separate makes the role of the covariance matrix and the boundary conditions explicit.

The Multivariate Normal Score

Definition: Let \(X\) be a \(d\)-dimensional random vector with distribution \(N_d(\mu,\Sigma)\), where \(\mu\in\mathbb{R}^d\) and \(\Sigma\) is a symmetric positive-definite \(d\times d\) matrix. Its density at \(x\in\mathbb{R}^d\) is $$ f_{\mu,\Sigma}(x)=(2\pi)^{-d/2}|\Sigma|^{-1/2} \exp\left\{-\frac{1}{2}(x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu)\right\}. $$ The score with respect to the observation \(x\) is the gradient of the log density, \(\nabla_x\log f_{\mu,\Sigma}(x)\).

Here, \(\nabla_x\) is the column vector of partial derivatives with respect to the coordinates of \(x\). The score is a function of the observed value \(x\), with \(\mu\) and \(\Sigma\) held fixed. It is not the score with respect to a model parameter.

Result: For \(X\sim N_d(\mu,\Sigma)\), the score with respect to \(x\) is $$ \nabla_x\log f_{\mu,\Sigma}(x)=-\Sigma^{-1}(x-\mu). $$

To derive it, take the logarithm of the density:

$$ \log f_{\mu,\Sigma}(x) =-\frac d2\log(2\pi)-\frac12\log|\Sigma| -\frac12(x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu). $$

Only the last term depends on \(x\). Since \(\Sigma^{-1}\) is symmetric, differentiating the quadratic form gives \(\nabla_x\{(x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu)\}=2\Sigma^{-1}(x-\mu)\). Therefore,

$$ \nabla_x\log f_{\mu,\Sigma}(x) =-\frac12\left\{2\Sigma^{-1}(x-\mu)\right\} =-\Sigma^{-1}(x-\mu). $$

The minus sign reflects the fact that the density decreases as one moves away from its center in the direction measured by \(\Sigma^{-1}\). Equivalently, the centered observation can be written as \(x-\mu=-\Sigma\nabla_x\log f_{\mu,\Sigma}(x)\). This algebraic form is one reason the score is useful in expectation calculations.

Worked Example: A Diagonal Covariance Score

Let \(X\sim N_2(\mu,\Sigma)\), where \(\mu=(1,-2)^\mathsf{T}\) and \(\Sigma=\begin{pmatrix}4&0\\0&9\end{pmatrix}\). Find the score at \(x=(3,1)^\mathsf{T}\).

The centered observation and inverse covariance are $$ x-\mu=\begin{pmatrix}2\\3\end{pmatrix}, \qquad \Sigma^{-1}=\begin{pmatrix}1/4&0\\0&1/9\end{pmatrix}. $$

Thus $$ \nabla_x\log f_{\mu,\Sigma}(x) =-\begin{pmatrix}1/4&0\\0&1/9\end{pmatrix} \begin{pmatrix}2\\3\end{pmatrix} =\begin{pmatrix}-1/2\\-1/3\end{pmatrix}. $$ The coordinates are scaled differently because their variances differ. In particular, the score is not simply \(-(x-\mu)\).

Integration by Parts in the Spherical Case

For the expectation relation below, assume \(X\sim N_d(\mu,\sigma^2 I_d)\), with \(\sigma^2>0\). Its density is

$$ \phi_{\mu,\sigma^2}(x) =(2\pi\sigma^2)^{-d/2} \exp\left\{-\frac{\|x-\mu\|^2}{2\sigma^2}\right\}, $$

and its score is \(\nabla_x\log\phi_{\mu,\sigma^2}(x)=-(x-\mu)/\sigma^2\). We use the following sufficient conditions on a test function \(h:\mathbb{R}^d\to\mathbb{R}\): \(h\) is continuously differentiable, and there are finite constants \(C\) and \(m\geq 0\) such that \(|h(x)|\) and every \(|\partial h(x)/\partial x_i|\) are at most \(C(1+\|x\|^m)\). These conditions ensure the expectations below exist and that the boundary terms vanish. They are sufficient conditions, not a claim that no weaker assumptions can work.

Theorem: If \(X\sim N_d(\mu,\sigma^2 I_d)\) and \(h\) satisfies the conditions just stated, then for each \(i=1,\ldots,d\), $$ E[(X_i-\mu_i)h(X)] =\sigma^2 E\left[\frac{\partial h(X)}{\partial x_i}\right]. $$

Proof. Fix a coordinate \(i\), and hold the other \(d-1\) coordinates fixed. The one-dimensional Gaussian density factor in coordinate \(i\), as a function of \(x_i\), satisfies $$ \frac{\partial}{\partial x_i}\phi_{\mu,\sigma^2}(x) =-\frac{x_i-\mu_i}{\sigma^2}\phi_{\mu,\sigma^2}(x). $$ On the interval \([-R,R]\) in coordinate \(x_i\), integration by parts gives

$$ \int_{-R}^{R}h(x)\frac{\partial\phi_{\mu,\sigma^2}(x)}{\partial x_i}\,dx_i =\left[h(x)\phi_{\mu,\sigma^2}(x)\right]_{-R}^{R} -\int_{-R}^{R}\frac{\partial h(x)}{\partial x_i} \phi_{\mu,\sigma^2}(x)\,dx_i. $$

As \(R\to\infty\), the boundary term tends to zero: for fixed values of the other coordinates, \(h(x)\) grows at most polynomially in \(|x_i|\), while the density decays as a Gaussian exponential. The two integrands after substituting the density derivative are absolutely integrable over \(\mathbb{R}^d\). Indeed, they are bounded in absolute value by a polynomial in \(\|x\|\) times a Gaussian density, whose integral is finite. Fubini’s theorem therefore allows integration over the remaining coordinates and passage to the full-space integrals. Substituting the derivative of the density yields $$ -\frac{1}{\sigma^2}E[(X_i-\mu_i)h(X)] =-E\left[\frac{\partial h(X)}{\partial x_i}\right]. $$ Multiplying by \(-\sigma^2\) proves the claimed relation. \(\square\)

Worked Example: Recovering a Gaussian Variance

Let \(X\sim N_2((1,-2)^\mathsf{T},4I_2)\), and choose \(h(x)=x_1\). Then \(\partial h/\partial x_1=1\), so the theorem gives $$ E[(X_1-1)X_1]=4E[1]=4. $$

The same result follows by expanding \(X_1=(X_1-1)+1\): $$ E[(X_1-1)X_1] =E[(X_1-1)^2]+E[X_1-1] =4+0=4. $$ This checks both the variance factor and the centering in the formula.

Worked Example: A Quadratic Test Function

Under the same distribution, take \(h(x)=x_1^2\). The derivative is \(2x_1\), and \(E[X_1]=1\). Applying the relation, $$ E[(X_1-1)X_1^2] =4E[2X_1]=8. $$

To check directly, write \(X_1=1+Z\), where \(Z\sim N(0,4)\). Then $$ E[(X_1-1)X_1^2] =E[Z(1+Z)^2] =E[Z+2Z^2+Z^3] =0+2(4)+0=8. $$ The odd centered moments vanish by symmetry, and the result agrees with the integration-by-parts calculation.

Worked Example: A Trigonometric Test Function

Let \(X\sim N_d(\mu,\sigma^2I_d)\), and let \(t\in\mathbb{R}^d\). Choose \(h(x)=\sin(t^\mathsf{T}x)\), a bounded function with bounded derivatives. Its \(i\)th partial derivative is \(t_i\cos(t^\mathsf{T}x)\), so $$ E[(X_i-\mu_i)\sin(t^\mathsf{T}X)] =\sigma^2t_iE[\cos(t^\mathsf{T}X)]. $$

The scalar \(t^\mathsf{T}X\) is normal with mean \(t^\mathsf{T}\mu\) and variance \(\sigma^2\|t\|^2\). Its normal characteristic function gives \(E[\cos(t^\mathsf{T}X)]=\exp\{-\sigma^2\|t\|^2/2\}\cos(t^\mathsf{T}\mu)\). Hence $$ E[(X_i-\mu_i)\sin(t^\mathsf{T}X)] =\sigma^2t_i\exp\left\{-\frac{\sigma^2\|t\|^2}{2}\right\} \cos(t^\mathsf{T}\mu). $$ The derivative factor \(t_i\) shows how the result responds to changes in the test function along coordinate \(i\).

Why the Conditions and Scaling Matter

The integration-by-parts step is not merely a formal manipulation. Its boundary term must vanish, and the integrals must exist. Gaussian tails handle any test function and first derivative of at most polynomial growth, but they do not control every rapidly growing function. For example, in one dimension, if \(X\sim N(0,1)\) and \(h(x)=\exp(x^2/2)\), then \(h(x)\phi_{0,1}(x)\) is a positive constant. Its integral over the real line is infinite, so the expectation in the proposed formula is not defined. The theorem cannot be applied.

A second common error is to omit the variance factor. The derivative of a \(N(\mu,\sigma^2)\) density is \(-(x-\mu)\phi(x)/\sigma^2\), not \(-(x-\mu)\phi(x)\). Consequently, integration by parts produces \(\sigma^2\) multiplying the derivative expectation. When \(\sigma^2=1\), that factor happens to disappear, which can conceal the error.

The score derivation works for any positive-definite covariance matrix, while the expectation relation proved here assumes spherical covariance \(\sigma^2 I_d\). These are distinct claims: the density’s score identifies the direction and scaling of local change, and the integration-by-parts calculation additionally depends on a valid test function and vanishing boundary terms. Making these distinctions explicit helps prevent an algebraic score formula from being used as an expectation identity without checking its conditions.

Check Your Understanding

Use the density derivative and the integration-by-parts argument to answer the following.

  1. For \(N_d(\mu,\Sigma)\), what is the score with respect to \(x\), and why does positive definiteness of \(\Sigma\) matter?
  2. For \(X\sim N_d(\mu,\sigma^2I_d)\), apply the result with \(h(x)=x_j-\mu_j\) to obtain \(E[(X_i-\mu_i)(X_j-\mu_j)]\).
  3. What boundary term appears when integrating \(h(x)\,\partial\phi(x)/\partial x_i\) over a finite interval in coordinate \(i\)?
  4. Why is a polynomial-growth condition on \(h\) and its first derivatives sufficient for the Gaussian integrals used in the proof to be finite?
  5. Identify the missing factor in the incorrect claim \(E[(X_i-\mu_i)h(X)]=E[\partial h(X)/\partial x_i]\) when \(X\sim N_d(\mu,\sigma^2I_d)\).
Clinical Biostats