# Algorithms for Stochastic Continuous Optimization

## Question

I have a continuous optimization problem of the form $$\max_{x \in \mathbb R^n} f(x),$$ where $$f:\mathbb R^n \rightarrow \mathbb R$$ is mostly smooth and bounded above. The standard approximation algorithms for these kinds of problems, such as basinhopping gradient ascent, take numerical estimates of the derivatives of $$f$$. But I can't compute $$f$$ directly; I can only observe it through noise. Random noise will make the numerical derivatives vary wildly, so I can't use any regular derivative-based optimization algorithm.

What are the best optimization algorithms that avoid this problem? I'm considering simulated annealing, which uses no calculus whatsoever. But I have to imagine that there are differential techniques which are robust to noise in the objective function. Does it help if I know/can control the approximate noise scale?

## EDIT: Context

I didn't initially say I'm working on since it isn't essential to my question, but now I'll add a precise problem statement for those curious. I'm trying to capture as much of Gaussian-distributed population as I can with a fixed number of unit balls. The objective function is the probability mass of the union of the balls. Formally, the problem is $$\max_{p_1,\ldots,p_N \in \mathbb{R}^n} P\left\{1 \geq \min_{n \in [N]} ||x - p_n||_2 \right\}$$ where $$x \sim N(\mu, \Sigma)$$ is a Gaussian-distributed vector.

So in my case, $$f$$ is an aggregation of a Gaussian PDF over a rather ugly domain of integration. My intuition tells me $$f$$ is non-concave, but continuous everywhere and smooth almost everywhere. Yet, I don't think it is possible to compute (or smoothly approximate) this with the efficiency required to support thousands of queries from an optimizer. But a relatively good noisy estimate of $$f$$ can be obtained quickly by simply sampling a large number of points from $$N(\mu, \Sigma)$$ and taking the proportion contained by the balls.

• What an interesting question. Just out of curiosity, do you mind sharing more about the context from which it came? Commented Jun 24, 2021 at 0:38
• @CadeReinberger I just added a section describing the specific optimization problem I'm curious about. Commented Jun 24, 2021 at 1:41
• If I am understanding correctly, you're trying to maximize the coverage you get with $p$ over the values $x$. Probably this is an unsupervised learning algorithm, the closest I know of are things such as k-means. This is just speculation, but I also suspect that adding onto the objective function with an additional regularization term as the shared hyper-volume amongst the unit balls may help. Commented Jun 24, 2021 at 5:09
• k-means is something I had not thought of, but interesting. Right now I'm getting decent-looking results in 2D and 3D with simulated annealing, using a random jumping function which tends to push nearby balls apart. So shared hypervolume isn't punished by the objective function, but still discouraged. Commented Jun 24, 2021 at 14:22

I recommend looking into Simultaneous Perturbation Stochastic Approximation (SPSA).

Very recently, I also happened to be writing an implementation in Python, which you can now view on GitHub. The documentation for the main optimize function goes in-depth about the features it provides, but I'll just provide some of the main ones here.

SPSA doubles as both an optimization and a gradient approximation algorithm. As the name implies, it is capable of handling stochastic problems (you can see the examples.py on my GitHub). It is similar to gradient descent except for the fact that it estimates the gradient using function evaluations instead.

With regards to your concerns on the influence of noise on the gradient approximation, you may use large perturbations to reduce the impact of noise. Usually large perturbations are used initially, and they gradually diminish to provide more accurate average gradient approximations. For more information, you can read the documentation of the optimize function under:

• Parameters.HyperParameters.Perturbation
• OtherFeatures.Momentum
• OtherFeatures.GradientNormalization

I will also suggest keeping your end of the computations to a minimum. The example.py shows an example of polynomial regression, where every function evaluation is a completely random point, with additional noise on top of that. The gradient is estimated from two perturbed coefficient vectors and at two random points over the domain. Despite how bad you might imagine this to be, the plots show that it converges quite well, or at least to an extremely reasonable extent.

• This looks like a cool algorithm, and your documentation is very helpful! If I decide to improve upon my simulated annealing solution, I'll definitely give this a go. It would be fun to benchmark simulated annealing, SPSA, and perhaps some kind of genetic learning algorithm against k-means. Commented Jun 24, 2021 at 14:33
• I haven't gotten around to benchmarking yet, but it's something I probably will. It's worth pointing out that unlike simulated annealing and genetic algorithms, it isn't suitable for combinatorial optimization. This means it's not directly applicable to a problem such as k-means. Commented Jun 24, 2021 at 14:38
• Right, I was still talking about my ball-coverage problem in particular, for which k-means and SPSA might both be successful. But generally, SPSA is still trying to do gradient descent right? So you'd have to combine it with basinhopping or n-restarts to globally optimize nonconvex problems? Commented Jun 24, 2021 at 14:55
• In the documentation under SPSAFeatures.GlobalLocalSearch, you will find information concerning global optimization with SPSA. It performs quite well for me on many of the test functions for optimization compared to gradient descent, especially the ones with lots of local minima. Commented Jun 24, 2021 at 15:49
• I was playing around with your problem and settled on this. The main thing I realized was that your objective function as-is is very bad at detecting changes in values without heavy sampling. By instead considering the similar problem of $$\min_p\mathbb E\left[\min_i\|x-p_i\|_2^2\bigg|x\right],$$ you can get similar results even with only 1 sample per objective function call. Commented Jun 27, 2021 at 0:00