Skip to main content
1 of 2

Standard multivariate Normal density function times quadratic form

For a psd real symmetric $d\times d$ matrix $A$ and a function $f: \mathbb{R}^d \to \mathbb{R}$, with $f(x) := x^T A x$ we have that with $p(x) = \mathcal{N}(0_d, I_d)$ (i.e. standard multivariate Normal), the well-known equality

$$\mathbb{E}_{x \sim p(x)}[f(x)] = \textrm{tr}(A).$$

This is Hutchinson's estimator and used in a variety of applications where only matrix-vector products with $A$ are available.

In deriving an importance sampling approach of the form

$$\mathbb{E}_{x \sim p(x)}[f(x)] = \mathbb{E}_{x \sim q(x)}\left[\frac{p(x)}{q(x)} f(x)\right],$$

I am interested in finding the optimal $q(x)$ in terms of variance reduction. This optimal solution is known to be (see Art Owen's book, chapter 9) of the form

$$q^*(x) = \frac{|f(x)|}{\mathbb{E}_{x\sim p}[|f(x)|]} \, p(x),$$

where here the absolute-sign can be dropped due to $f(x) > 0$ for all $x \in \mathbb{R}^d$. Then we have

$$q^*(x) = \frac{x^T A x}{\mathbb{E}_{x\sim p}[x^T A x]} \, p(x),$$

and this simplifies to

$$q^*(x) = \frac{1}{\textrm{tr}(A)} \, (x^T A x) \, p(x),$$

which says that the optimal importance sampling density function is a quadratic form $x^T A x$ multiplied with the multivariate Normal density. While this $q^*$ is not directly useful in practice (we do not know $A$), if we were to use $q^*$ the variance of the Hutchinson estimator becomes zero, as the $x^T A x$ terms cancel.

Questions:

  • Does the $q^*$ distribution have a name or is otherwise known? (There are many works on quadratic forms under normal variates, but here the density function itself is multiplied.)
  • Is there a natural approach for efficient generation from this density? (Elliptical slice sampling seems possible, but perhaps there is a simpler line of attack.)