# Covariance of a linear and quadratic form of a multivariate normal

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.

• You may find some useful identities in the matrix cookbook. In particular, chapter 6.2 and/or 8.2 may have what you need. – Macro Aug 25 '12 at 15:19
• Ok, with your help of your reference I believe I have proven it is 0 if $\boldsymbol{\mu}=\mathbf{0}$ (still need to check it) – Clark Aug 26 '12 at 11:47

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.
• Does Isserlis's theorem also apply in the case when $i, j, k$ are not distinct as would happen in this instance? – Dilip Sarwate Aug 26 '12 at 18:02
• @Dilip, thanks for pointing that out. In that case you don't actually need the theorem, right? For the case where $i=j=k$, $E(y_{i}^3) = 0$ trivially. For the case where $i=k \neq j$, my intuition is that $y_{j}$ is independent of $y_{i}^2$ since the conditional density of $y_{j} | y_{i}^2$ (at a particular point $y_{j} = y$) is symmetric in $y_i$, which has mean zero. – Macro Aug 26 '12 at 18:35
• Hint: $\mathbf y = \mathbf A \mathbf z$ for some $\mathbf z$ a standard normal vector (possibly of smaller dimension than $\mathbf y$). – cardinal Aug 26 '12 at 18:45
• In other words: Set $\mathbf b = \mathbf A^T \mathbf a$ and $\widetilde{\mathbf H} = \mathbf A^T \mathbf H \mathbf A$ and use your proof by dealing with $\mathrm{Cov}( \mathbf b^T \mathbf z, \mathbf z^T \widetilde{\mathbf H}\mathbf z)$. Only the definition of multivariate normality, covariance, independence, and basic properties of the moments of a single iid standard normal random variable are used. – cardinal Aug 26 '12 at 19:47
• Dear Macro, I didn't intend to initiate such a long comment stream and my remarks were only intended to be constructive. Your answer is fine as is; my primary thought was that (as @Dilip alludes to) appealing to the fact that $\mathbb E y_i y_j^2 = 0$ and $\mathbb E y_i y_j y_k = 0$ for zero-mean multivariate normals may appear a bit opaque and this is really where most (or much) of the action actually lies. (It's the only place where any property of normality comes into play.) The comments were intended to show a simple way of bringing out precisely why it was true with no additional work. – cardinal Aug 26 '12 at 23:08

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