2
$\begingroup$

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.)
$\endgroup$
4
  • $\begingroup$ Doesn't the standard Hutchinson trace estimator use a random vector $X$ whose components are i.i.d. Rademacher random variables? I thought the corresponding estimator based on Gaussians is called the Gaussian trace estimator. $\endgroup$ Commented Aug 1, 2022 at 11:25
  • 1
    $\begingroup$ @NawafBou-Rabee, yes, you are correct, Hutchinson's original paper only used Rademacher vectors. The question only relates to the Gaussian trace estimator. $\endgroup$ Commented Aug 1, 2022 at 11:40
  • $\begingroup$ Since $q^*$ has a density relative to a Gaussian reference measure (proportional to $e^{\log(x^T A x)}$), it seems like a preconditioned MCMC method (in the same spirit as elliptical slice sampling) might work well. The convergence of the pMCMC method will, however, depend crucially on the properties of the Hessian of $U(x) = - \log(x^T A x)$. $\endgroup$ Commented Aug 1, 2022 at 12:58
  • 1
    $\begingroup$ This paper will be very interesting to you :) arxiv.org/pdf/2010.09649.pdf $\endgroup$ Commented Aug 16, 2022 at 3:52

0

You must log in to answer this question.