Covariance of a linear and quadratic form of a multivariate normal

Question

Does anyone know of an explicit matrix expression for the covariance of a linear and quadratic form? That is,

$\mathrm{Cov}[\mathbf{a' y},\mathbf{y' Hy}]$ where $\mathbf{y}\sim \mathcal N(\boldsymbol{\mu},\boldsymbol{\Sigma})$.

I'm particularly interested in the case where $\boldsymbol{\mu}=\mathbf{0}$, and I think this simplifies (without the normal assumption) to

$\mathbb E[(\mathbf{a'y})(\mathbf{y'Hy})]$. Since this involves cubic terms it probably isn't going to be simple.

Answer 1

This is straightforward in the case you're interested in (${\boldsymbol \mu} = 0$) without using matrix algebra.

To clarify the notation ${\bf y} = \{y_{1}, ..., y_{n} \}$ is a multivariate normal random vector, ${\bf a} = \{a_{1}, ..., a_{n} \}$ is a row vector, and ${\bf H}$ is an $n \times n$ matrix with entries $\{ h_{jk} \}_{j,k=1}^{n}$. By definition (see e.g. page 3 here) you can re-write this covariance as $$ {\rm cov}({\bf a}'{\bf y}, {\bf y}' {\bf H} {\bf y}) = {\rm cov} \left( \sum_{i=1}^{n} a_i y_i, \sum_{j=1}^{n} \sum_{k=1}^{n} h_{jk} y_{j} y_{k} \right) = \sum_{i,j,k} {\rm cov}( a_i y_i, h_{jk} y_{j} y_{k} ) $$ where the second equality follows from bilinearity of covariance. When ${\boldsymbol \mu} = 0$, each term in the sum is $0$ because $${\rm cov}( a_i y_i, h_{jk} y_{j} y_{k} ) \propto E(y_i y_j y_k) - E(y_i) E(y_k y_j) = 0$$ The second term is zero because $E(y_i) = 0$. The first term is zero because the third order mean-centered moments of a multivariate normal random vector are 0, this can be seen more clearly by looking at the each cases:

• when $i,j,k$ are distinct, then $E(y_i y_j y_k)=0$ by Isserlis' Theorem
• when $i\neq j = k$, then we have $E(y_i y_j y_k) = E(y_i y_{j}^2)$. First we can deduce from here that $E(y_i | y_j=y) = y \cdot \Sigma_{ij}/\Sigma_{jj}$. Therefore, $E(y_{i} y_{j}^2 | y_{j} = y) = y^3 \cdot \Sigma_{ij}/\Sigma_{jj}$. Therefore, by the law of total expectation, $$E(y_i y_{j}^2) = E( E(y_{i} y_{j}^2 | y_{j} = y) ) = E(y^3 \Sigma_{ij}/\Sigma_{jj} ) = E(y^3) \cdot \Sigma_{ij}/\Sigma_{jj} = 0 $$ where $E(y_{j}^3) = 0$ because $y_j$ being symmetrically distributed with mean 0 implies that $y_{j}^3$ is also symmetrically distributed with mean 0.
• when $i=j=k$, $E(y_i y_j y_k) = E(y_{i}^3) = 0$ by the same rationale just given.

Answer 2

This can be proved with multivariate Stein's Lemma. Letting $x=a'y\sim \mathcal{N}(a'\mu,a'\Sigma a)$, we have that $Cov(x,y) = a'\Sigma$. Let $h(y) = y'Hy$. Stein's lemma then tells us $$ Cov(x,h(y)) = Cov(x,y)E\left[\nabla h(y)\right] = a'\Sigma E\left[\left(H + H'\right) y\right] = a'\Sigma \left(H + H'\right) \mu. $$

Some simulations:

library(mvtnorm)

p <- 9
nsim <- 100000
set.seed(1234)
mu <- rnorm(p)
X <- matrix(rnorm(3*p^2),ncol=p)
H <- matrix(rnorm(p^2),ncol=p)
Sigma <- cov(X)
a <- rnorm(p)

Y <- rmvnorm(nsim,mean=mu,sigma=Sigma)
Ya <- Y %*% a
YHY <- rowSums((Y %*% H) * Y)

# empirical
emp <- cov(Ya,YHY)
# theoretical
thr <- a %*% Sigma %*% (H + t(H)) %*% mu
cat('empirical: ',emp,' theoretical: ',thr,' \n')

empirical:  -2.39  theoretical:  -2.385

Answer 3

The covariance between a linear and a quadratic form of a multivariate normal vector is given in Mathai and Provost, page 74, Theorem 3.2d.2.

Let $\mathbf{y} \in \mathbb{R}^p \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})$, $Q = \mathbf{y}'\mathbf{H}\mathbf{y}$ and $P = \mathbf{a}'\mathbf{y}$, with $\mathbf{H} \in \mathbb{R}^{p\times p}$ symmetric and $\mathbf{a} \in \mathbb{R}^p$.

Then, $Cov(Q,P) = 2 \mathbf{\mu}' \mathbf{H} \mathbf{\Sigma} \mathbf{a}$.

In agreement with the other answers, this formula shows that if $\mathbf{\mu} = 0$ the covariance is 0.