# Optimization as sampling for stochastic functions

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!

• 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.
– πr8
Jun 21, 2020 at 10:45
• 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.
– πr8
Jun 21, 2020 at 10:47
• @π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.
– Y.Z.
Jun 22, 2020 at 15:03
• A paper from this morning (arxiv.org/abs/2008.00234) addresses roughly your problem, and may be of interest.
– πr8
Aug 4, 2020 at 8:54

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.

• 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?
– Y.Z.
Jun 29, 2020 at 20:11
• 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$?
– πr8
Jun 29, 2020 at 20:56
• 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.
– Y.Z.
Jun 29, 2020 at 21:03
• 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).
– Y.Z.
Jun 29, 2020 at 21:04
• 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.
– πr8
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$$