Principal Hessian Directions From Inverse Moments

Why Look at Response-Weighted Predictor Moments?

A regression direction can be invisible to inverse regression means when the response depends symmetrically on a predictor: for example, the mean of $Z_1^2$ is not a linear function of $Z_1$. A different summary can detect such curvature. It asks whether the response is associated with products of standardized predictors, such as $Z_iZ_j$. The matrix of these associations is the principal Hessian directions (PHD) matrix.

The name reflects a useful connection: under a multivariate normal predictor distribution, this inverse-moment matrix equals an average Hessian of the regression mean. Its eigenspaces can therefore reveal directions along which the mean response bends. The connection is not automatic for arbitrary predictor distributions, and the average curvature can cancel even when the regression function is nonlinear.

Definition and Standardization

Let $X\in\mathbb R^p$ have mean $\mu$ and positive-definite covariance matrix $\Sigma$. Define the standardized predictor $Z=\Sigma^{-1/2}(X-\mu)$, where $\Sigma^{-1/2}$ is the symmetric inverse square root. Then $E[Z]=0$ and $\operatorname{Var}(Z)=I_p$. For the Gaussian identity developed below, assume in addition that $X\sim N(\mu,\Sigma)$, so $Z\sim N(0,I_p)$.

Definition: If $E[|Y|(1+\|Z\|^2)]<\infty$, the population PHD matrix is $$ M_{\mathrm{PHD}}=E\!\left[(Y-E[Y])(ZZ^{\mathsf T}-I_p)\right]. $$ Equivalently, its $(i,j)$ entry is $\operatorname{Cov}(Y,Z_iZ_j)$. The matrix is symmetric. Its eigenvectors are directions in standardized predictor coordinates.

The integrability condition ensures both that $E[Y]$ exists and that every matrix entry is finite. Since $E[ZZ^{\mathsf T}]=I_p$, the same matrix can be written as $E[Y(ZZ^{\mathsf T}-I_p)]$. The centered form makes its covariance interpretation explicit. If an eigenvector in standardized coordinates is $\eta$, its corresponding linear predictor in the original coordinates is $\eta^{\mathsf T}\Sigma^{-1/2}(X-\mu)$, with direction vector $\Sigma^{-1/2}\eta$.

Worked Example: A Quadratic Effect in One Direction

Let $Z\sim N(0,I_2)$, let $\varepsilon$ be independent of $Z$ with $E[\varepsilon]=0$, and set $Y=Z_1^2+\varepsilon$. Assume $E[|\varepsilon|]<\infty$. Then $E[Y]=E[Z_1^2]=1$.

The $(1,1)$ entry is $$ \operatorname{Cov}(Y,Z_1^2)=E[Z_1^4]-E[Z_1^2]^2=3-1=2. $$ Independence of $Z_1$ and $Z_2$, together with $E[Z_2^2]=1$, gives $$ \operatorname{Cov}(Y,Z_2^2)=E[Y]E[Z_2^2]-E[Y]E[Z_2^2]=0. $$ For the off-diagonal entry, $E[YZ_1Z_2]=0$: the term involving $Z_1^3Z_2$ has expectation zero by independence and $E[Z_2]=0$, and the independent error term also has expectation zero. Moreover, $E[Z_1Z_2]=0$. Thus $$ M_{\mathrm{PHD}}= \begin{pmatrix} 2&0\\ 0&0 \end{pmatrix}. $$ Its nonzero eigenvector is the first coordinate direction, which is the direction of the quadratic regression effect. A first-moment summary such as $E[Z\mid Y]$ need not reveal this direction as clearly because the effect is symmetric in $Z_1$.

The Inverse-Moment Identity

The matrix is called an inverse-moment matrix because it can be expressed using the conditional mean of the response given the predictors. Write $m(z)=E[Y\mid Z=z]$, choosing a version of this conditional expectation. Under the integrability condition in the definition, conditioning on $Z$ gives the following identity, whether or not $Z$ is Gaussian:

Identity: If $E[|Y|(1+\|Z\|^2)]<\infty$, then $$ M_{\mathrm{PHD}}=E\!\left[(m(Z)-E[Y])(ZZ^{\mathsf T}-I_p)\right]. $$

Indeed, each entry of $ZZ^{\mathsf T}-I_p$ is a function of $Z$, so the defining property of conditional expectation replaces $Y$ by $m(Z)$ inside its product with that entry. The integrability assumption also gives $E[|m(Z)|(1+\|Z\|^2)]<\infty$, by conditional Jensen's inequality. This identity by itself does not identify PHD with a Hessian; that step uses Gaussian predictors and smoothness.

Gaussian Integration by Parts and the Hessian

Theorem: Suppose $Z\sim N(0,I_p)$, $E[|Y|(1+\|Z\|^2)]<\infty$, and the conditional mean $m:\mathbb R^p\to\mathbb R$ is twice continuously differentiable. Assume $m$ and its first and second partial derivatives have at most polynomial growth. Then $$ M_{\mathrm{PHD}}=E[\nabla^2m(Z)], $$ where $\nabla^2m$ is the matrix of second partial derivatives.

The growth condition is a convenient sufficient condition for the needed Gaussian integrals to be finite and the integration-by-parts boundary terms to vanish. The theorem does not assert that a mean and covariance alone suffice: the standardized predictors must have the stated Gaussian distribution.

Proof. For indices $i,j$, let $\phi(z)=(2\pi)^{-p/2}\exp(-\|z\|^2/2)$ be the standard normal density. Direct differentiation gives $$ \frac{\partial^2\phi(z)}{\partial z_i\,\partial z_j} =(z_i z_j-\delta_{ij})\phi(z), $$ where $\delta_{ij}=1$ if $i=j$ and $0$ otherwise. By the inverse-moment identity and $E[Z_iZ_j]=\delta_{ij}$, $$ (M_{\mathrm{PHD}})_{ij} =E[m(Z)(Z_iZ_j-\delta_{ij})] =\int_{\mathbb R^p}m(z)\frac{\partial^2\phi(z)}{\partial z_i\,\partial z_j}\,dz. $$ Integrate by parts first in coordinate $i$, then in coordinate $j$. The boundary terms vanish because $m$ and its derivatives grow at most polynomially while the Gaussian density and its derivatives decay exponentially. The two integrations each introduce a minus sign, so their product is positive: $$ \int_{\mathbb R^p}m(z)\frac{\partial^2\phi(z)}{\partial z_i\,\partial z_j}\,dz =\int_{\mathbb R^p}\frac{\partial^2m(z)}{\partial z_i\,\partial z_j}\phi(z)\,dz =E\!\left[\frac{\partial^2m(Z)}{\partial z_i\,\partial z_j}\right]. $$ This holds for every $i,j$, proving the matrix identity. $\square$

Worked Example: A Direction After Unequal Predictor Scaling

Suppose $X\sim N(\mu,\operatorname{diag}(4,1))$, and define $Z_1=(X_1-\mu_1)/2$ and $Z_2=X_2-\mu_2$. Because $X$ is multivariate normal with the stated diagonal covariance, $Z\sim N(0,I_2)$ and its coordinates are independent. Let $Y=(Z_1+Z_2)^2+\varepsilon$, where $\varepsilon$ is independent with mean zero and finite first absolute moment.

Here $E[Y]=E[Z_1^2]+2E[Z_1Z_2]+E[Z_2^2]=2$. The diagonal entries are $$ \operatorname{Cov}(Y,Z_1^2)=E[Z_1^4+2Z_1^3Z_2+Z_1^2Z_2^2]-2=3+0+1-2=2, $$ and, by the same moment calculation with the coordinates exchanged, $\operatorname{Cov}(Y,Z_2^2)=2$. For the off-diagonal entry, $$ \operatorname{Cov}(Y,Z_1Z_2) =E[Z_1^3Z_2+2Z_1^2Z_2^2+Z_1Z_2^3]-2E[Z_1Z_2] =0+2+0-0=2. $$ Consequently $M_{\mathrm{PHD}}=\begin{pmatrix}2&2\\2&2\end{pmatrix}$. Its eigenvalues are $4$ and $0$; a unit eigenvector for $4$ is $(1,1)^{\mathsf T}/\sqrt 2$. In the original predictor coordinates the corresponding index is proportional to $$ Z_1+Z_2=\tfrac12(X_1-\mu_1)+(X_2-\mu_2). $$ Thus the original-scale direction is proportional to $(1/2,1)^{\mathsf T}$. The normality assumption is essential to the stated calculation; diagonal covariance alone would not establish independent standard normal coordinates.

When PHD Recovers Regression Directions

The Hessian identity gives a precise recovery condition. Suppose the conditional mean has a $d$-dimensional reduction $m(z)=g(B^{\mathsf T}z)$, where $B$ is a $p\times d$ matrix with linearly independent columns and $g:\mathbb R^d\to\mathbb R$ is twice continuously differentiable. The chain rule gives $$ \nabla^2m(z)=B\,\nabla^2g(B^{\mathsf T}z)\,B^{\mathsf T}. $$ Consequently, for Gaussian $Z$ under the theorem's integrability and smoothness conditions, $$ M_{\mathrm{PHD}}=B A B^{\mathsf T}, \qquad A=E[\nabla^2g(B^{\mathsf T}Z)]. $$

The column space of $M_{\mathrm{PHD}}$ is contained in the column space of $B$. If the $d\times d$ matrix $A$ is nonsingular, then $B$ has full column rank and $B^{\mathsf T}$ maps onto $\mathbb R^d$; it follows that $B A B^{\mathsf T}$ has rank $d$ and its column space equals that of $B$. This is the basic population recovery condition. If $B$ spans the central subspace, PHD then recovers that subspace. More generally, if the conditional distribution of $Y$ given $Z$ depends on $B^{\mathsf T}Z$, the conditional mean has that reduction, but recovery of the full distributional subspace still requires that the mean-curvature matrix $A$ have full rank.

Worked Example: Curvature Cancellation Can Hide a Direction

Let $Z\sim N(0,I_2)$, take $Y=Z_1^4-6Z_1^2+\varepsilon$, and let $\varepsilon$ be independent of $Z$, mean zero, with finite first absolute moment. The conditional mean depends nontrivially on $Z_1$, so the response has a regression direction. But its second derivative in that direction is $$ \frac{d^2}{dz_1^2}(z_1^4-6z_1^2)=12z_1^2-12, $$ whose expectation is $12E[Z_1^2]-12=0$.

The inverse moment gives the same cancellation directly. Since $E[Z_1^2]=1$, $E[Z_1^4]=3$, and $E[Z_1^6]=15$, $$ \operatorname{Cov}(Y,Z_1^2) =E[(Z_1^4-6Z_1^2)(Z_1^2-1)] =15-3-6(3)+6(1)=0. $$ The other diagonal entry is zero because $Z_2$ is independent of $Y$ and $E[Z_2^2-1]=0$. The off-diagonal entry is zero because $Z_2$ is independent and mean zero. Hence $M_{\mathrm{PHD}}=0$, despite the nonconstant regression mean. PHD misses this direction because positive and negative curvature average to zero.

Interpretation and Common Pitfalls

PHD is useful as a curvature-based complement to methods built from inverse regression means. It can identify quadratic effects that are symmetric about the predictor mean, and its population matrix has a direct direction interpretation when the mean-reduction and full-rank curvature conditions hold. Its limitations are just as important: it uses the conditional mean, relies on Gaussian predictors for the Hessian identity, and summarizes curvature by an average that may be singular or zero.

Using covariance as a substitute for normality. Standardizing to covariance $I_p$ does not make $Z$ Gaussian. The Hessian identity requires the Gaussian design assumption (or a separately justified identity for another distribution).
Omitting response integrability. Finiteness of $E[|Y|\|Z\|^2]$ alone need not imply that $E[Y]$ exists. The stated condition $E[|Y|(1+\|Z\|^2)]<\infty$ ensures both centering and finite entries.
Treating a nonzero eigenvector as guaranteed. A regression direction can be present while the average Hessian vanishes. Recovery requires the expected curvature matrix to have adequate rank.
Confusing standardized and original-scale directions. An eigenvector $\eta$ in $Z$-coordinates corresponds to $\Sigma^{-1/2}\eta$ in the original $X$-coordinates.

Key takeaway: For standardized Gaussian predictors, PHD is both a covariance between the response and predictor outer products and the expected Hessian of the conditional mean. Its eigenspace recovers a mean-reduction subspace only when the average curvature spans that subspace.

Check Your Understanding

Use the definitions, identity, and recovery condition to answer the following.

Why is $M_{\mathrm{PHD}}$ a covariance between $Y$ and the entries of $ZZ^{\mathsf T}$ when $E[ZZ^{\mathsf T}]=I_p$?
Which additional distributional assumption is needed to identify PHD with the expected Hessian of the conditional mean?
If $m(z)=g(B^{\mathsf T}z)$, write the expected-curvature matrix whose nonsingularity gives recovery of the column space of $B$.
How can a nonconstant conditional mean produce a zero PHD matrix?
If $\eta$ is a PHD eigenvector for standardized predictors, what is its corresponding direction in the original predictor coordinates?