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!

  • 1
    $\begingroup$ It is perhaps worth noting that this bears many similarities with the setting where $f$ is a sum of many terms, i.e. $f(x) = \sum_{i = 1}^N f_i (x)$, and one estimates $f$ by sampling a subset of these terms; this is roughly how stochastic gradient descent works. There are possibly useful ideas in this literature which could be useful. $\endgroup$
    – πr8
    Commented Jun 21, 2020 at 10:45
  • $\begingroup$ Some small clarifying questions: i) is your space $X$ structured in some way? e.g. is it a subset of $\mathbf{R}^d$, or is it discrete, or something else? ii) is your function $f$ structured in some way? e.g. bounded above, smooth, etc. Even quite basic structures like this can be useful in working out which algorithmic approaches are available. $\endgroup$
    – πr8
    Commented Jun 21, 2020 at 10:47
  • $\begingroup$ @πr8 For your questions, in my case $X$ is a discrete combinatorial space, which is why I am thinking about the sampling approach rather than gradient descent. $f$ is bounded above and below, so I can normalize it, and thus the expectation as well, to 0 and 1 very easily. But I am not sure how they help in this specific case. $\endgroup$
    – Y.Z.
    Commented Jun 22, 2020 at 15:03
  • $\begingroup$ A paper from this morning (arxiv.org/abs/2008.00234) addresses roughly your problem, and may be of interest. $\endgroup$
    – πr8
    Commented Aug 4, 2020 at 8:54

2 Answers 2


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.

  • $\begingroup$ Thanks for your explanation. I have two questions. 1. Can you tell me which identity you used to derive the 2nd line in the multi-line equation? 2. As I understand, the last equation defines the transition probability. I can define a transition kernel on $X$, with probability $q(x'\rightarrow x)$ and $q(x\rightarrow x')$, but how should I choose a transition kernel for $N$ and $\Xi$? Are $N$ and $N'$ sampled from $\mathrm{Po}(\beta)$, and $\xi$ implicit when I evaluate the stochastic $F(x, \xi)$? Given that each $F(x, \xi)$ evaluation is time-consuming, can I force $N=1$ throughout? $\endgroup$
    – Y.Z.
    Commented Jun 29, 2020 at 20:11
  • $\begingroup$ The 2nd line in the equation comes from the Taylor expansion of the function $x \mapsto \exp(x)$. The generation of $(N, \Xi)$ are sampled implicitly, as you say. You cannot force $N = 1$ and still get the correct invariant measure; it is important that $N$ can become arbitrarily large. If $F$ is expensive to evaluate, then perhaps there are other structures which you can exploit - can you give more information about $F$? $\endgroup$
    – πr8
    Commented Jun 29, 2020 at 20:56
  • $\begingroup$ I see. Thanks for your reply. In my case, $F$ is the validation accuracy of a deep neural network trained on some data, and $X$ is the hyper-parameter space. So I am essentially trying to find the best hyper-parameter setting that achieves on average best validation accuracy. The stochasticity comes from the random seed and architectural randomness (e.g. dropout). So I don't think there are more insights that can be exploited. $\endgroup$
    – Y.Z.
    Commented Jun 29, 2020 at 21:03
  • $\begingroup$ By the way, I am aware of other techniques such as Bayesian optimization, but I am just wondering if this hyper-parameter search can be approached from this sampling perspective, which is really flexible and does not depend on having a parametric representation of the mapping from the hyper-parameter space to the validation accuracy (as does by e.g. Gaussian process regression). $\endgroup$
    – Y.Z.
    Commented Jun 29, 2020 at 21:04
  • $\begingroup$ Understood. I think in this setting, if you're set on doing something like simulated annealing, it would be a reasonable idea to just do a biased estimate of $\exp(-\beta f(x))$. You will not get the same asymptotic guarantees, but the algorithm may nevertheless be useful. $\endgroup$
    – πr8
    Commented Jun 29, 2020 at 21:27

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$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.