# Integrating with a multivariate Gaussian

I need to figure out the steps to solve the following integral, where $$Q(\mathbf{w})$$ is a multivariate Gaussian with mean $$\overline{\mathbf{w}}$$ and covariance $$\mathbf{C}$$:

\begin{align}\int Q(\mathbf{w}) \> \frac{1}{2} \mathbf{w}^\mathrm{T}\mathbf{A}\mathbf{w} \> \mathrm{d} \mathbf{w} = \mathrm{Tr}(\mathbf{CA}) + \frac{1}{2} \overline{\mathbf{w}}^\mathrm{T} \mathbf{A} \overline{\mathbf{w}} \end{align}

Essentially, I'm very happy to try and work it out for myself if someone can just point me to relevant textbooks/online info or give a hint - I'm just having trouble getting started.

So far I've tried making a substitution $$\mathbf{u} = \mathbf{w} - \overline{\mathbf{w}}$$ which gets me to:

\begin{align} \int Q(\mathbf{w}) \> \frac{1}{2} \mathbf{w}^\mathrm{T}\mathbf{A}\mathbf{w} \> \mathrm{d} \mathbf{w} &= \int R(\mathbf{u}) \> \frac{1}{2} (\mathbf{u}+ \overline{\mathbf{w}})^\mathrm{T}\mathbf{A}(\mathbf{u}+ \overline{\mathbf{w}}) \> \mathrm{d} \mathbf{u} \\ & = \int R(\mathbf{u}) \{ \frac{1}{2} \mathbf{u}^\mathrm{T} \mathbf{A} \mathbf{u} + \frac{1}{2} \overline{\mathbf{w}}^\mathrm{T}\mathbf{A}\mathbf{u} + \frac{1}{2} \mathbf{u}^\mathrm{T}\mathbf{A}\overline{\mathbf{w}} + \frac{1}{2} \overline{\mathbf{w}}^\mathrm{T} \mathbf{A} \overline{\mathbf{w}} \} \mathrm{d} \mathbf{u} \end{align}

where $$R(\mathbf{u})$$ is a multivariate Gaussian with mean $$\mathbf{0}$$ and covariance $$\mathbf{C}$$. From here, I can see how the second term, $$\frac{1}{2} \overline{\mathbf{w}}^\mathrm{T} \mathbf{A} \overline{\mathbf{w}}$$, might come out of this expression, but not at all how the first term, $$\mathrm{Tr}(\mathbf{CA})$$ might come out.

I also had an idea to try a substitution of the form $$\mathbf{y} = \mathbf{M}^{-1}(\mathbf{w} - \overline{\mathbf{w}})$$, where $$\mathbf{MM}^\mathrm{T} = \mathbf{C}$$.

I've also tried working in Einstein summation notation (as I thought this might make the origin of the trace more clear) but that didn't get me much farther.

Thanks!

Update:

Following a helpful comment by @whuber, I had a look at the one dimensional case in order to gain insight, using the equivalent one dimensional substitution to $$\mathbf{y} = \mathbf{M}^{-1}(\mathbf{w} - \overline{\mathbf{w}})$$.

\begin{align} \int \frac{1}{2} w^2 a \> Q(w)\> \mathrm{d} w &= \int \frac{1}{2} w^2 a \> \frac{\exp(-\frac{1}{2c}(w-\overline{w})^2)}{(2\pi c)^{1/2}} \mathrm{d}w \end{align} Now we make a substitution of $$y = \frac{1}{c^{1/2}} (w - \overline{w})$$. \begin{align} & = \int \frac{1}{2} ( yc^{1/2} + \overline{w})^2 a \frac{\exp(-\frac{1}{2} y^2 )}{(2\pi)^{1/2} } \mathrm{d}y \\ & = \int \frac{1}{2} ( y^2c + 2yc^{1/2}\overline{w} + \overline{w}^2 ) a \> \frac{\exp(-\frac{1}{2} y^2 )}{(2\pi)^{1/2} } \mathrm{d} y \\ \end{align} And we notice that the term with $$2yc^{1/2}\overline{w}$$ is odd, and that $$\overline{w}^2$$ is independent of $$y$$. \begin{align} & =\frac{1}{2} \overline{w}^2 a + \int \frac{1}{2} y^2c a \> \frac{\exp(-\frac{1}{2} y^2 )}{(2\pi)^{1/2} } \mathrm{d} y \\ & =\frac{1}{2} \overline{w}^2 a + \frac{1}{2} ca \frac{1}{(2\pi)^{1/2}} \int y^2\> \exp(-\frac{1}{2} y^2 ) \mathrm{d} y \end{align} Finally, we integrate by parts, using $$u = y$$ and $$\frac{\mathrm{d} v}{\mathrm{d} y} = y \exp(-\frac{1}{2} y^2)$$. \begin{align} \int \frac{1}{2} w^2 a \> Q(w)\> \mathrm{d} w &= \frac{1}{2}\overline{w}^2 a + \frac{1}{2}ca \end{align}

• First principles tell you the answer will be a bilinear function of $C$ and $A$ and a quadratic function of $\bar w,$ so that gets you most of the way there while suggesting what change of variables to employ. Your last idea of a substitution is a good one. If you want to get some sense of what's involved, then first solve the special case where $CA$ is diagonal and $\bar w=0.$ If that's too hard, solve the one-dimensional version. – whuber May 7 at 12:58

The means and covariances already evaluate all the integrals you need, allowing this result to be obtained purely algebraically. It actually has nothing to do with Normal distributions (except insofar as they have finite covariances in the first place).

Let $$X=(X_1,X_2,\ldots,X_n)$$ be a multivariate random variable with mean vector $$\mu = (\mu_1,\mu_2,\ldots,\mu_n)=(E[X_1], E[X_2], \ldots, E[X_n])$$ and covariance matrix $$C = (c_{ij}) = (E[X_iX_j] - E[X_i]E[X_j]) = (E[X_iX_j]-\mu_i\mu_j),$$ which can be rewritten $$E[X_iX_j] = c_{ij} + \mu_i\mu_j.$$

Let $$A=(a_{ij})$$ be any $$n\times n$$ matrix. The definition of matrix multiplication, linearity of expectation, the foregoing formula, and definition of the trace imply

\eqalign{ E\left[X^\prime A X\right] &= E\left[\sum_{i,j} X_i a_{ij} X_j\right] \\ &= \sum_{ij} E[X_iX_j] a_{ij} \\ &= \sum_{ij}(C_{ij} +\mu_i\mu_j) a_{ij} \\ &= \operatorname{Tr}(CA) + \mu^\prime A \mu. }

Multiply everything by $$1/2$$ and let $$X$$ be multivariate normal.

• BTW, using the Einstein summation convention permits you to remove all appearances of "$\sum_{ij}$" but otherwise wouldn't change any of this demonstration. – whuber May 8 at 1:02
• Thanks @whuber - does this mean that the final answer should be $\frac{1}{2} \mathrm{Tr}(\mathbf{CA}) + \frac{1}{2}\overline{\mathbf{w}}^\mathrm{T} \mathbf{A} \overline{\mathbf{w}}$ instead of $\mathrm{Tr}(\mathbf{CA}) + \frac{1}{2}\overline{\mathbf{w}}^\mathrm{T} \mathbf{A} \overline{\mathbf{w}}$? The one dimensional case seems to agree with your answer, but the paper I am looking at (Barber & Bishop, 1998) gives $\mathrm{Tr}(\mathbf{CA}) + \frac{1}{2}\overline{\mathbf{w}}^\mathrm{T} \mathbf{A} \overline{\mathbf{w}}$ as the final result. – Winter May 8 at 8:56
• I believe the factor of $1/2$ applies to both terms. We can check this, as you suggest, by considering the simplest possible case: the one-dimensional standard Normal distribution where $C=(1)$ and $\bar w=(0).$ Letting $A=(1)$ reduces the integral to the expectation of one-half the square, which evaluates to $1/2 = \operatorname{Tr}(CA)/2.$ Thus the factor of $1/2$ must be in there. Assume its absence in your reference is a typographical error. – whuber May 8 at 14:10

The solution posted by whuber gets at this idea but I wanted to make the approach more explicitly use the trace operator, use:

$$\mathbb{E}(u^TAu) = \mathbb{E}(tr(u^TAu)).$$ Note that the quadratic form inside the expectation is a scalar and trace of a scalar is the same scalar. Next use cyclic swap property of the trace operator:

$$\mathbb{E}(tr(u^TAu))= \mathbb{E}(tr(uu^TA)).$$

Now note that $$A$$ is a constant so you can factor that out:

$$\mathbb{E}(tr(uu^TA))=tr(\mathbb{E}(uu^TA)).$$ Once more step: $$tr(\mathbb{E}(uu^TA))= tr(\mathbb{E}(uu^T)A)$$

Next consider that by definition of the covariance matrix of a multivariate Gaussian: $$\mathbb{E}(uu^T) = C.$$

Stitch all of this together and you've got the result.

• This is almost right: you forgot to subtract the means in the definition of $C.$ – whuber May 8 at 3:01
• @whuber I am using the OP's notation where $u=w-\bar{w}$ where $\bar{w}$ is the mean of the original multivariate normal. The mean vectors after subtraction should be 0 vectors. Did I misunderstand your comment? Or maybe the OP's notation? – Lucas Roberts May 8 at 4:21
• I did notice that I had some typos w/the trace operator and once missed inclusion altogether. I think the typos are removed now but let me know if you catch one – Lucas Roberts May 8 at 4:32
• Thank you for the explanation--I hadn't made the connection with the OP's notation for residuals. – whuber May 8 at 5:02