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.