Optimization as sampling for stochastic functions

Question

Given an input space $X$ and a function $f: X\rightarrow \mathbb R$, we want to find $x^*=argmin_{x\in X} f(x)$. One way is to cast this problem as a sampling, where we define a distribution $p(x)\propto e^{-f(x)}$. The mode of the distribution corresponds to $x^*$. We can draw $N$ samples from $p(x)$ and pick the one that minimizes $f(x)$ as $x^*$. For example, if we use Metropolis-Hastings algorithm as the sampler, then we are doing something similar to simulated annealing.

However, in my problem, $f(x)$ is stochastic, and we want to find the minimizer in expectation, $x^*=argmin_{x\in X} \mathbb E[f(x)]$. I can evaluate $f(x)$ but it is a quite slow procedure, so I would prefer not to e.g. evaluate $f(x)$ 100 times and take the average. In addition, given a specific $y$ from an $f(x)$ evaluation, I don't know its probability mass/density, even up to a constant. Essentially $f(x)$ is just a black-box stochastic procedure that returns a sample after some quite expensive computation.

My question is, can I still use a similar sampling idea for the optimization? A naive way is to pretend that a single $y\sim f(x)$ sample is actually $\mathbb E[f(x)]$, and use that value in the MH-sampler. But I don't know what, if any, distribution is implicitly being sampled.

Another idea is to sample jointly in the $x, y\in X, \mathbb R$ space, but since I can't evaluate the likelihood of $y$, even up to a normalizing constant, under $f(x)$, and running $f(x)$ multiple times is perhaps too expensive, I don't know how to write a sampler with this constraint.

Any ideas are greatly appreciated!

Answer 1

To expand upon the solution which is hinted at in the answer of @Xi'an:

Assume that $f$ is represented as

$$f(x) = \mathbf{E}_{\rho(\xi)} \left[ F(x, \xi) \right]$$

where $\xi$ is some auxiliary source of randomness, and $0 \leqslant F(x, \xi) \leqslant 1$ for all $(x, \xi)$.

One can then develop

\begin{align} \exp(-\beta f(x)) &= \exp \left( -\beta \right) \cdot \exp \left(\beta \left\{1 - f(x) \right\} \right) \\ &= \sum_{n \geqslant 0} \frac{\beta^n e^{-\beta}}{n!} \left\{1 - f(x) \right\}^n \\ &= \mathbf{E}_{N \sim \text{Po}(\beta)} \left[ \left\{1 - f(x) \right\}^N \right] \\ &= \mathbf{E}_{N \sim \text{Po}(\beta)} \left[ \prod_{a = 1}^N \mathbf{E}_{\rho(\xi^a)} \left[ 1 - F \left(x, \xi^a \right) \right] \right]. \end{align}

This implies that if we write down the joint distribution

$$ \Pi \left( x, N, \{ \xi^a \}_{a = 1}^N \right) \propto \frac{\beta^N e^{-\beta}}{N!} \cdot \prod_{a = 1}^N \left\{ \rho(\xi^a) \left[ 1 - F \left(x, \xi^a \right) \right] \right\},$$

then the $x$-marginal is given by $\mu_\beta (x) \propto \exp(-\beta f(x))$.

This enables the application of a Pseudo-Marginal Metropolis-Hastings MCMC algorithm. Consider the proposal

$$Q \left( (x, N, \Xi) \to (x', N', \Xi') \right) = q ( x \to x' ) \cdot \text{Po} ( N' | \beta ) \cdot \prod_{b = 1}^{N'} \rho ( \xi'^b ).$$

Working through the details, one can compute that the Metropolis-Hastings ratio simplifies to

$$r \left( (x, N, \Xi) \to (x', N', \Xi') \right) = \frac{q ( x' \to x )}{q ( x \to x' )} \cdot \frac{ \prod_{b = 1}^{N'} \left[ 1 - F \left(x, \xi'^b \right) \right] }{ \prod_{a = 1}^N \left[ 1 - F \left(x, \xi^a \right) \right]}$$

which can be computed exactly, allowing for a tractable Metropolis-Hastings correction. This means that one can generate a Markov chain with $\Pi \left( x, N, \Xi \right)$ as its invariant measure, and hence the $x$-marginal of the chain will converge to $\mu_\beta$ as desired.

Answer 2

This is a very interesting question for which there is no clear-cut answer. It all depends on the computing budget and the output of a realistic will depend on this computing budget.

My suggestion would be to mix

(i) simulated annealing, that is, simulating from a target like $$h_t(x)\propto e^{-T_t \cdot \mathbb E[f(x)]}\qquad T_t \uparrow \infty$$ where the temperature $T_t$ is slowing increasing with $t$,

(ii) pseudo-marginal Metropolis-Hastings, when the value of the target is replaced with an unbiased estimate at each iteration, and

(iii) debiasing à la Glynn and Rhee, as in Russian roulette estimators, where a converging sequence of biased estimators, $\hat\eta_n$ is turned into a unbiased estimator $$\sum_{n=1}^G \{\eta_{n+1}-\eta_n\}/\mathbb P(G\ge n)$$ $G$ being a integer valued random variable (like a Poisson). This last step involves computing a random number $G$ of realisations of $f(x)$.

An alternative is to use stochastic optimisation, by considering the sequence $(X_n)_n$ such that $$X_{n+1}=X_n-\epsilon_n \nabla f(X_n)\qquad \epsilon_n\downarrow 0$$ where $\nabla f$ denotes a realisation of the gradient of $f$, i.e. $$\mathbb E[\nabla f(X_n] = \nabla \mathbb E[f(X_n]]$$ If this is impossible to obtain, a finite difference approach is the Kiefer-Wolfowitz algorithm $$X_{n+1}=X_n-\epsilon_n \dfrac{f(X_n+\upsilon_n)-f(X_n-\upsilon_n)}{2\upsilon_n}\qquad \epsilon_n,\upsilon_n\downarrow 0$$