Multivariate Stein Identity by Coordinatewise Integration

Why Integrate One Coordinate at a Time?

For a scalar normal variable, integration by parts relates a centered observation to a derivative of a function. A vector version follows by applying that scalar argument separately to each coordinate. The key is to keep track of which component of a vector-valued function is paired with each coordinate: the resulting derivatives are the matching-coordinate partial derivatives, not a sum over every possible partial derivative.

This tutorial treats \(X\sim N(\mu,\sigma^2 I_p)\), where \(X,\mu\in\mathbb{R}^p\), \(\sigma^2>0\), and \(I_p\) is the \(p\)-dimensional identity matrix. The coordinates are independent normal variables with means \(\mu_i\) and common variance \(\sigma^2\). We derive the identity under explicit conditions that justify integration by parts and the subsequent integration over the other coordinates.

Coordinatewise Derivative and Regularity Conditions

Definition: For a differentiable function \(g=(g_1,\ldots,g_p):\mathbb{R}^p\to\mathbb{R}^p\), define its divergence by \(\operatorname{div}g(x)=\sum_{i=1}^p \frac{\partial g_i(x)}{\partial x_i}\). Each term differentiates component \(g_i\) with respect to the same coordinate \(x_i\).

To use ordinary one-dimensional integration by parts, it is enough to assume that every component \(g_i\) is continuously differentiable and that the following conditions hold for each \(i\):

\(E[|(X_i-\mu_i)g_i(X)|]<\infty\) and \(E[|\partial_i g_i(X)|]<\infty\), where \(\partial_i g_i=\frac{\partial g_i}{\partial x_i}\).
For almost every fixed value of the other coordinates, the boundary terms vanish: \(\lim_{t\to\pm\infty}g_i(x_1,\ldots,x_{i-1},t,x_{i+1},\ldots,x_p)\phi_i(t)=0\), where \(\phi_i\) is the \(N(\mu_i,\sigma^2)\) density.

The integrability assumptions ensure that the relevant expectations exist and that the integrals can be iterated. The boundary condition allows integration by parts along the \(i\)th coordinate. These are sufficient conditions, not a claim that every possible function for which the identity holds must satisfy them.

The Coordinatewise Stein Identity

Theorem: Suppose \(X\sim N(\mu,\sigma^2 I_p)\), with \(\sigma^2>0\), and \(g:\mathbb{R}^p\to\mathbb{R}^p\) satisfies the differentiability, integrability, and boundary conditions above. Then \[ E[(X-\mu)^{\mathsf T}g(X)] =\sigma^2 E[\operatorname{div}g(X)] =\sigma^2 E\left[\sum_{i=1}^p \frac{\partial g_i(X)}{\partial x_i}\right]. \]

To see the mechanism, fix an index \(i\) and hold all coordinates other than \(x_i\) fixed. Write \(x_{-i}\) for those other coordinates. The joint normal density factors as a product of the coordinate densities, so integration over \(x_i\) uses \(\phi_i\), while the other density factors remain unchanged. Since \(\phi_i'(t)=-\frac{t-\mu_i}{\sigma^2}\phi_i(t)\), one-dimensional integration by parts gives

\begin{aligned} \int_{-\infty}^{\infty}(t-\mu_i)g_i(x_{-i},t)\phi_i(t)\,dt &=-\sigma^2\int_{-\infty}^{\infty}g_i(x_{-i},t)\phi_i'(t)\,dt\\ &=\sigma^2\int_{-\infty}^{\infty}\partial_i g_i(x_{-i},t)\phi_i(t)\,dt. \end{aligned}

The boundary term in this integration by parts is zero by assumption. Integrate both sides over the other coordinates. The integrability assumptions permit these iterated integrals, giving \(E[(X_i-\mu_i)g_i(X)]=\sigma^2E[\partial_i g_i(X)]\). Now sum over \(i\). The left side sums to the inner product \((X-\mu)^{\mathsf T}g(X)\), and the right side sums to the defined divergence. This proves the theorem. \(\square\)

Select a coordinate.
Pair \(X_i-\mu_i\) with the \(i\)th component \(g_i(X)\).

Integrate along that coordinate.
Use the derivative of the \(N(\mu_i,\sigma^2)\) density and integrate by parts with respect to \(x_i\).

Integrate over the remaining coordinates and sum.
The individual identities combine into the inner product on the left and the divergence on the right.

Worked Checks of the Identity

Worked Example: An Affine Vector Field

Let \(p=2\), \(X\sim N(\mu,4I_2)\), and set \(Z=X-\mu\). Take \(g(X)=AZ+b\), where \(A=\begin{pmatrix}1&2\\-3&1/2\end{pmatrix}\) and \(b=(2,-1)^{\mathsf T}\). Because \(g\) is affine, its derivatives are constant and its Gaussian expectations are finite.

The divergence is the trace of \(A\): \(\operatorname{div}g=1+\frac12=\frac32\). Thus the right side is \(4\cdot\frac32=6\).

For the left side, \(E[Z]=0\), so \(E[Z^{\mathsf T}b]=0\). The coordinates of \(Z\) are independent, centered, and have variance \(4\), hence \(E[Z_iZ_j]=0\) for \(i\ne j\) and \(E[Z_i^2]=4\). Therefore \(E[Z^{\mathsf T}AZ]=A_{11}E[Z_1^2]+A_{22}E[Z_2^2]=4+2=6\). Both sides equal \(6\).

Worked Example: Cubic Coordinate Functions

Let \(X\sim N(\mu,2I_2)\), write \(Z_i=X_i-\mu_i\), and define \(g_1(X)=2Z_1^3\) and \(g_2(X)=-Z_2^3\). These polynomial components and their derivatives satisfy the integrability and boundary conditions.

The left side is \(E[Z_1g_1(X)+Z_2g_2(X)]=2E[Z_1^4]-E[Z_2^4]\). For a centered normal variable with variance \(2\), the fourth moment is \(3(2^2)=12\). Thus the left side is \(2(12)-12=12\).

The matching-coordinate derivatives are \(\partial_1g_1=6Z_1^2\) and \(\partial_2g_2=-3Z_2^2\). Since \(E[Z_i^2]=2\), the right side is \(2E[6Z_1^2-3Z_2^2]=2(12-6)=12\). The equality holds; the cubic terms also show why the derivative is taken in the coordinate matching the component.

Worked Example: Exponential Functions of Centered Coordinates

Let \(X\sim N(\mu,I_2)\), \(Z=X-\mu\), and \(g_1(X)=2e^{Z_1}\), \(g_2(X)=e^{-2Z_2}\). All required expectations are finite: a normal variable has finite exponential moments for every real exponent, and the derivatives are constant multiples of these same functions.

For \(Z\sim N(0,1)\), completing the square in its density gives \(E[e^{tZ}]=e^{t^2/2}\). Differentiating this explicit expression, or completing the square in the corresponding first-moment integral, gives \(E[Ze^{tZ}]=t e^{t^2/2}\). Consequently the left side is \(2E[Z_1e^{Z_1}]+E[Z_2e^{-2Z_2}] =2e^{1/2}-2e^2\).

The derivatives are \(\partial_1g_1=2e^{Z_1}\) and \(\partial_2g_2=-2e^{-2Z_2}\). Their expected sum is \(2e^{1/2}-2e^2\), equal to the left side because \(\sigma^2=1\).

Worked Example: Cross-Coordinate Dependence

Again let \(X\sim N(\mu,\sigma^2I_2)\), and write \(Z=X-\mu\). Choose \(g_1(X)=Z_2\) and \(g_2(X)=Z_1\). The function components depend on the other coordinate, but their matching-coordinate derivatives are both zero: \(\partial_1g_1=0\) and \(\partial_2g_2=0\).

The left side is \(E[Z_1Z_2+Z_2Z_1]=2E[Z_1Z_2]=0\), since the centered coordinates are independent. The right side is \(\sigma^2E[0+0]=0\). Cross-coordinate dependence is allowed; it simply does not contribute a cross-partial term to this identity.

What the Identity Does—and Does Not—Say

The identity converts an expectation involving centered observations into an expectation involving derivatives of a vector field. Its utility in decision theory is that derivatives can be easier to analyze than the original expectation when studying functions of a multivariate Gaussian observation. At this stage, the result is an integration identity; by itself, it does not assert that any estimator has lower risk or that a particular choice of \(g\) is preferable.

The coordinatewise proof also makes the assumptions visible. The normal density contributes the factor \(1/\sigma^2\) when differentiated, which becomes \(\sigma^2\) after integration by parts. The boundary term must vanish, and the expectations must exist. Omitting those checks can make a formally plausible calculation invalid, especially for functions that grow rapidly or have singularities.

Common Mistakes

Summing all partial derivatives. The divergence here is \(\sum_i\partial_i g_i\), not \(\sum_i\sum_j\partial_jg_i\). For \(g_1=Z_2,\ g_2=Z_1\), the matching derivatives vanish, even though cross-partials do not.
Dropping the variance factor. The identity has \(\sigma^2\) on the derivative side. Differentiating the normal density introduces \(1/\sigma^2\), and integration by parts reverses that factor.
Pairing the wrong components. The left side is \(\sum_i(X_i-\mu_i)g_i(X)\). Each coordinate is paired with its corresponding component before the terms are summed.
Ignoring regularity. A pointwise derivative formula alone does not justify integration by parts over an unbounded interval. Check integrability and the boundary behavior of \(g_i\phi_i\).

Key takeaway: Under suitable regularity conditions, integrate by parts in coordinate \(i\) to obtain \(E[(X_i-\mu_i)g_i(X)]=\sigma^2E[\partial_i g_i(X)]\). Summing these scalar identities yields the multivariate normal identity with the divergence \(\sum_i\partial_i g_i\).

Check Your Understanding

Work through these questions to check both the calculation and its conditions.

For \(g_i(X)=c_i(X_i-\mu_i)\), what are the two sides of the identity?
Why does the boundary condition involve \(g_i(x_{-i},t)\phi_i(t)\), rather than \(g_i(x_{-i},t)\) alone?
For \(p=2\), \(g_1(X)=X_2-\mu_2\) and \(g_2(X)=0\), calculate the divergence and the left side.
Which factor changes if the common coordinate variance is changed from \(\sigma^2\) to \(4\sigma^2\)?
Give one reason why finite expectations and vanishing boundary terms must be checked before applying coordinatewise integration by parts.