# Minimizing approximate function

I use a Monte-Carlo algorithm to estimate a target value y=f(x) that depends on a single input parameter x. Because the Monte-Carlo algorithm is stochastic, the target value y fluctuates. I now want to find the value of x that minimizes y. Obviously, I cannot use most of the standard optimization algorithms (like Newton's), because I cannot calculate gradients of y with respect to x.

What optimization algorithms would be useful in this kind of situation?

Here is a simple python code the illustrates the situation

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize

def f(x):
return -np.sin(x + np.pi/2 + 0.5 * np.random.random())

xs = np.linspace(-2, 2)
ys = [f(x) for x in xs]
plt.plot(xs, ys, 'o')

optimize.minimize_scalar(f)


which results in the following 'minimum'

 fun: -0.99996139255316596
nfev: 38
nit: 37  success: True
x: -0.22358360754884515


The function is constructed such that its minimum is at x=0, see the figure produced by the above code: • What optimization algorithms would be useful in this kind of situation? I highlighted it in the questions in case it's difficult to spot. – David Zwicker Jul 31 '16 at 15:32
• But it seems that you have build-in optimization algorithm that works? Moreover, the function itself is strange since it includes random component. Could you tell us a little bit more about your problem? – Tim Jul 31 '16 at 15:36
• It may be possible to directly estimate a noisy derivative estimate. For instance, in this example, -sin(x) plus whatever noise, can be used as a derivative estimate. It may be possible to find such estimates even for very complicated functions which are not known in closed form, but only evaluable by stochastic simulation. – Mark L. Stone Jul 31 '16 at 15:38
• Well the example shows that the minimization does not work very well, since it reports 0.07, where the minimum should be at 0. Obviously, this is caused by the random components, which should illustrate the structure of my problem, where f(x) is a complicated function doing a Monte-Carlo simulation. I'm looking for optimization algorithms that provide a better solution to the problem than what I have so far. – David Zwicker Jul 31 '16 at 15:40
• @DavidZwicker if your f is much more complicated and the noise is not additive then you should change your example to something closer to your problem. I just deleted my answer that was dealing with your example that is actually quite easy to optimize! – Tim Jul 31 '16 at 16:13

It sounds like you're interested in 'stochastic optimization', where the goal is to optimize a stochastic objective function (typically its expected value). Note that some people take this term to include methods for optimizing deterministic functions where the solver uses randomness (not what you want).

These references may be useful:

Hannah (2014). Stochastic Optimization.

Fu et al. (2005). Simulation optimization: A review, new developments, and applications

Amaran et al. (2014). Simulation optimization: A review of algorithms and applications

You may also be interested in Bayesian optimization. In this setting, the objective function can be stochastic, and the goal is to choose parameters that optimize its expected value, given the parameters. As with some other stochastic optimization methods, the objective function can be a black box, meaning you have the ability to evaluate it, but may not have a closed form expression (e.g. it might depend on the result of a simulation or physical experiment). Evaluating it may be very expensive. Bayesian optimization treats evaluations of the objective function as observed data, and uses them to update a probabilistic model of the objective function (e.g. using Gaussian process regression). New evaluation points are chosen in a way that trades off between exploration (sampling from uncertain regions to get a better estimate of the objective function) and exploitation (sampling from regions that are predicted to increase/decrease the objective function).

Brochu et al. (2010). A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning.

Snoek et al. (2012). Practical Bayesian Optimization of Machine Learning Algorithms.

If you do have a closed form expression for the objective function (as in your example), it would make more sense to try to exploit this known structure than to throw it away and treat the function as a black box. For example, in the best case you may be able to derive an expression for the expected value given the parameters, then use standard, deterministic methods to optimize it.

You are trying to optimize a stochastic objective function because of random noise component in the objective function. For a noisy function like yours, You could use an evolutionary optimization methods such as genetic algortithm to solve stochastic optimization problems. Please have a look at this mathworks website that demonstrates how to solve a problem like yours. I would add that in addition to using an evolutionary optimization which is able to find near optimal (it can also easily handle multiple optimal solutions). In order to find a more precise accurate solution, I would club evolutionary optimization with a local optimization. This is called hybrid optimization. The above steps can be easily demonstrated using R. I have specified both global optimization using GA and a hybrid optimization (global-local) optimization.

set.seed(1234)
obj.f <- function(x) {
fx <- sin(x+pi/2+0.5*rnorm(1))
return(fx)
}

#Global optimization only
ga.opt.global <- ga(type = "real-valued", fitness = obj.f, min = -2, max = 2)
summary(ga.opt.global)
plot(ga.opt.global)
# Hybrid Optimziation
ga.opt.hybrid <- ga(type = "real-valued", fitness = obj.f, min = -2, max = 2,optim=TRUE,optimArgs = list(method = "Nelder-Mead"))
summary(ga.opt.hybrid)
plot(ga.opt.hybrid)


Please Note, I'm not sure what your overall objective is, but I have tried to provide an answer from a pure optimization point of you.

As it turns out as evident from your plot any values between -0.5 to 0 is optimal. IF you run the above optimization multiple times you will be able to find almost all the optimal values.

> summary(ga.opt.global)
+-----------------------------------+
|         Genetic Algorithm         |
+-----------------------------------+

GA settings:
Type                  =  real-valued
Population size       =  50
Number of generations =  100
Elitism               =  2
Crossover probability =  0.8
Mutation probability  =  0.1
Search domain =
x1
Min -2
Max  2

GA results:
Iterations             = 100
Fitness function value = 1
Solution =
x1
[1,] 0.4520799
> summary(ga.opt.hybrid)
+-----------------------------------+
|         Genetic Algorithm         |
+-----------------------------------+

GA settings:
Type                  =  real-valued
Population size       =  50
Number of generations =  100
Elitism               =  2
Crossover probability =  0.8
Mutation probability  =  0.1
Search domain =
x1
Min -2
Max  2

GA results:
Iterations             = 100
Fitness function value = 0.9999994
Solution =
x1
[1,] -0.1413327