Maximizing unknown noisy function

I'm interested in maximizing a function $f(\mathbf \theta)$, where $\theta \in \mathbb R^p$.

The problem is that I don't know the analytic form of the function, or of its derivatives. The only thing that I can do is to evaluate the function point-wise, by plugging in a value $\theta_*$ and get a NOISY estimate $\hat{f}(\theta_*)$ at that point. If I want I can decrease the variability of these estimates, but I have to pay increasing computational costs.

Here is what I have tried so far:

• Stochastic steepest descent with finite differences: it can work but it requires a lot of tuning (ex. gain sequence, scaling factor) and it is often very unstable.

• Simulated annealing: it works and it is reliable, but it requires lots of function evaluations so I found it quite slow.

So I'm asking for suggestions/idea about possible alternative optimization method that can work under these conditions. I'm keeping the problem as general as possible in order to encourage suggestions from research areas different from mine. I must add that I would be very interested in a method that could give me an estimate of the Hessian at convergence. This is because I can use it to estimate the uncertainty of the parameters $\theta$. Otherwise I'll have to use finite differences around the maximum to get an estimate.

-
 If you can't say anything more specific about the noise associated with your function's output, I'm not sure anything more sophisticated than simulated annealing (you'll even have to tune this, to some extent), will be of help. – Aron Ahmadia Oct 9 '12 at 16:47 Unfortunately I don't know much about the random noise associated with each function evaluation. Its distribution is unknown, and it can be a function of $\theta$. On the other hand the noises that affects successive function evaluations are independent. Obviously I'm assuming that the variance of noise is not enormous, otherwise maximization would be impossible. – Jugurtha Oct 9 '12 at 16:52 On the other hand suppose that I know something about the noise distribution, for example that $\hat{f}(\theta_*) \sim N(f(\theta_*),\sigma)$. Would this knowledge help me? – Jugurtha Oct 9 '12 at 17:01 Looks like I stand corrected by Prof. Neumaier :) – Aron Ahmadia Oct 9 '12 at 20:21 Physicists here, i used CMA-ES for optical phase shaping (optimizing the phase of a laser pulse via a pulseshaper), which is quite noisy. – tillsten Oct 11 '12 at 2:27
show 1 more comment

2 Answers

Our Matlab package SnobFit was created precisely for this purpose. No assumption about the distribution of the noise is needed. Moreover, function values can be supplied through text files, thus you can apply it to functions implemented in any system able to write a text file. See
http://www.mat.univie.ac.at/~neum/software/snobfit/

SnobFit had been developed for an application where the function to be optimized didn't even exist, and function values (a measure of manufacturing quality) were obtained by specialized, expensive equipment creating sample products and measuring these by hand, resulting in about 50 function evaluations per day.

-
 Thank you very much for your answer. I've started reading your article regarding the SnobFit package, and I find it really interesting. Also, while reading the introduction to your article, I realized that the problem I'm dealing with (in a statistical context) is pretty frequent in industrial mathematics. There is a vast literature of which I was completely unaware. Actually the approach I was working on is somewhat similar to the quadratic approximation of Powell (2002). – Jugurtha Oct 9 '12 at 19:53 Does snobfit work well with 128 degrees of freedom? Just to know it is worth to try out for my case. – tillsten Oct 11 '12 at 2:28 @tillsten: No methods for noisy problem works well with 128 dof unless you can spend a huge number of function values. You might try our VXQR1, though, which is for not noisy problems, but sometimes handles noisy problems well. – Arnold Neumaier Oct 11 '12 at 8:02 The limit for Snobfit is about 20 variables. if you have more, you need to select by common sense groups of 20 variables which you partially optimize in turn. Or you may let slide some variables simultaneously so that the dimension is reduced. – Arnold Neumaier Oct 11 '12 at 8:04

James Spall's SPSA algorithm (short for Stochastic Perturbation Simulated Annealing, if I recall correctly) has been designed for exactly this kind of problem. He has a couple of papers where he uses it for problems like the one you describe.

-
 I have tried Spall's approach based on a stochastic version of steepest descent and Raphson Newton. I tried Simulated Annealing, but not the version suggested by Spall, I should try it. I'm not really enthusiastic about simulated annealing, because I can't get an estimate of the Hessian at convergence (while, for example, with stochastic Raphson Newton I can get an approximation to the Hessian "for free"). – Jugurtha Oct 11 '12 at 16:12