Skip to main content
changed the title to something more effective
Link

Standard multivariate Normal density function times quadratic form Optimal Monte Carlo Trace Estimator

Source Link

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.)