# Optimization by random sampling

Around the internet, I have seen scattered references to the idea of rescaling an objective function and using that as a PDF for the purpose of optimization. (On this site for example: Do optimization techniques map to sampling techniques?) Can someone point me to a place where I can learn more about this method? (Papers, blog posts, lectures, etc.)

The idea, as I've seen it, is to take your objective function $f(x)$ and create a new function $g(x) = e^{kf(x)}$, where $k$ is a very large number for a maximization problem or a very large negative number for a minimization problem. The new function $g(x)$ will then be much higher at the global optimum than anywhere else. If $g(x)$ is then treated as an unnormalized probability density function, most samples drawn from that distribution will be around that optimum.

Things I want to know include, but are not limited to:

1. Which sampling algorithms are effective for these probability functions?
2. Why is this method not used more frequently? (It seems like it could be so effective). In other words, are there arguments against it?
3. Are there any variants of this method that improve efficiency or performance?

It sound like you are interested in studying stochastic optimisation methods. If you re-frame optimisation as a type of sampling problem, you will obtain a stochastic optimisation method, and this latter method will only be advantageous if it provides some improvement over analogous deterministic optimisation methods.

Generally speaking, stochastic optimisation methods of this kind will involve the generation of random values in a way that depends on the function being optimised, and so it is likely to be at least as computationally intensive (and probably more computationally intensive) than corresponding deterministic methods. Stochastic methods are also generally more complicated to understand. Thus, the only advantage that is likely to accrue to (well-constructed) stochastic optimisation methods, is that they generally retain some non-zero probability of "searching" areas of the domain that might be missed by deterministic methods, and so they are arguably more robust, if you are willing to run them for a long time.

Most standard deterministic optimisation methods involve taking iterative steps towards the optima by choosing a deterministic movement direction, and moving in that direction by some deterministic amount (e.g., with steepest ascent we move in the direction of the gradient vector). The direction and length of the steps are usually determined by looking at the gradient of the function, which in practice, is computed by looking at the slope over a small movement increment. By turning the optimisation into a sampling problem, you effectively just move in a random direction by a random amount, but you are still going to want to use information on the slope of the function to determine the stochastic behaviour of this movement. It is likely that the latter method will use the same information as the former, only in a more complex (and thus, more computationally intensive) way. Below I will construct a method using the MH algorithm, based on your description.

Implementation with Metropolis-Hastings algorithm: Suppose you are dealing with a maximisation problem with a distribution over one-dimension, and consider the method of Metropolis-Hasting sampling using Gaussian deviations. This is a well-known method of sampling that is quite robust against nasty behaviour by the sampling density.

To implement the algorithm, we generate a sequence $$\varepsilon_1, \varepsilon_2, \varepsilon_3, ... \sim \text{IID N}(0, 1)$$ of random deviations, which we will later multiply by the "bandwidth" parameter $$\lambda > 0$$ (so this parameter represents the standard deviation of our proposed steps). We also generate a sequence $$U_1,U_2,U_3 ,... \sim \text{U}(0,1)$$ of uniform random variables. We choose a starting value $$x_0$$ arbitrarily and generate the sequence of sample values recursively as:

$$x_{t+1} = \begin{cases} x_t + \lambda \varepsilon_{t+1} & & \text{if } U_{t+1} \leqslant \exp \Big( k \big( f(x_t + \lambda \varepsilon_{t+1}) - f(x_t) \big) \Big), \\[6pt] x_t & & \text{otherwise}. \\[6pt] \end{cases}$$

The MH algorithm has the stationary distribution equal to the target density, so we have the approximate distribution $$X_n \sim \exp( k (f(x))$$ for large $$n$$. Taking a large value for $$k$$ (relative to the bandwidth parameter) will mean that deviations in the direction of the optima will be accepted with high probability, and deviations away from the optima (or overshooting the optima) will be rejected with high probability. Thus, in practice, the sample values should converge near to the maximising point of the density $$f$$. (There would still be cases where this would not occur; e.g., if the density is bimodal, and the algorithm climbs the wrong slope.) If we take a large value of $$k$$ then the distribution $$\exp( k (f(x))$$ is highly concentrated near the maximising value, so in this case the sample mean of the sampled values (discarding some burn-in values) should give us a good estimate of the maximising value in the optimisation.

Now, suppose we counteract the large value of $$k$$ by setting the bandwidth $$\lambda$$ to be small. In that case, we get small values for the deviations, so we have the approximate acceptance probabilities:

$$\exp \Big( k \big( f(x_t + \lambda \varepsilon_{t+1}) - f(x_t) \big) \Big) \approx \exp \Big( k \lambda \cdot \varepsilon_{t+1} \cdot f'(x_t) \Big).$$

We can see that in this case, the acceptance probability is determined by the derivative of the density, the magnitude and direction of the deviation $$\varepsilon_{t+1}$$ and the value $$k \lambda$$.

So, does this algorithm actually perform any better than applicable deterministic optimisation algorithms? Well, we can see that the algorithm requires us to compute the density at each potential step, and if the bandwidth parameter is small then this is tantamount to computing the derivative of the function. So this is quite similar to a form of stochastic gradient ascent, with more computational work than the corresponding deterministic method. The advantage, if any, is that the proposed deviations are random, and all occur with some non-zero probability (albeit vanishing very quickly), so the algorithm has the potential to "escape" regions of high, but non-maximising, density.

• Is there any connection between this approach and the approach where noise is added to the gradient before each step, e.g. $\dot{\theta} = -\nabla L(\theta) + \varepsilon_t$ where $\varepsilon_t \sim N(0, \sigma_t)$? What about to stochastic gradient Langevin dynamics? Commented Sep 11, 2019 at 22:13
• They look like different procedures to me, but it is possible there may be some connections. They are connected at least in the loose sense that involve stochastic iterations where some noise component is added into the standard Newton-Raphson iteration (either in the direction, distance, or gradient).
– Ben
Commented Sep 11, 2019 at 23:06