# Optimizing noisy objective

**URL:** <https://discourse.julialang.org/t/optimizing-noisy-objective/25260>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [June 13, 2019, 5:39pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260 "2019-06-13T17:39:18Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [June 14, 2019, 2:28pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/21 "2019-06-14T14:28:29Z")

</div>

I don’t know if any Julia package implements it, but SPSA optimization is specifically made for this kind of thing. It gives you pseudo-gradients in O(1) function eval.

[https://www.jhuapl.edu/SPSA/](https://www.jhuapl.edu/SPSA/)

Easy to code up yourself if needed.

---

<div class="post-metadata">

**Author:** ![jonathanBieler](https://avatars.discourse-cdn.com/v4/letter/j/82dd89/32.png) [@jonathanBieler](https://discourse.julialang.org/u/jonathanBieler)\
**Post date:** [June 14, 2019, 3:21pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/22 "2019-06-14T15:21:37Z")

</div>

My version of [CMAES](https://github.com/jonathanBieler/CMAES.jl) gets down to 1e-6 on the abs2 case, but it take about 4000 function evaluations to get there.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 14, 2019, 4:09pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/23 "2019-06-14T16:09:43Z")

</div>

I did not check out your specific CMAES implementation, but in my experience it is usually rather inefficient.

Currently the most successful solution is NLopt’s `:GN_CRS2_LM` suggested by @aaowens above. Simulated annealing works fine on the MWE but stalls on the real problems (I had to kill it in after hour, CRS2 gets me to a reasonable optimum in a minute).

I am now thinking about making this problem differentiable, but it will be tricky. CRS2 consistently finds optima in about 500-1000 iterations, which is fine for me at the moment.

---

<div class="post-metadata">

**Author:** ![aharoun](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aharoun/32/6887_2.png) [@aharoun](https://discourse.julialang.org/u/aharoun)\
**Post date:** [June 14, 2019, 5:00pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/24 "2019-06-14T17:00:14Z")

</div>

Regarding the objective function: normalizing it with the covariance (like GMM) or using percentage difference between model and data may help.

I don’t have much experience with discrete choice models but I wonder why fixing the seed of the random number generator doesn’t kill the stochasticity.

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [June 15, 2019, 4:36am UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/25 "2019-06-15T04:36:53Z")

</div>

We wrapped a python based algorithm for gradient free nonlinear least squares a few months ago, see [https://github.com/QuantEcon/DFOLS.jl](https://github.com/QuantEcon/DFOLS.jl)

It turned out not to be useful for our problem, but if you think it matches your problem, give it a shot.

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [June 15, 2019, 6:11am UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/26 "2019-06-15T06:11:27Z")

</div>

The behavior of SA is very much dependent on the rt control. You might try setting it to rt=0.5, if you’re inclined to experiment more. SA has the advantage that it is very robust, if rt is set high enough to avoid focusing on a local min. But what is “high enough” is problem-specific. Too high, and it’s slow, too low, and you may skip over the global min.

Given that this is an indirect inference problem, I don’t understand from the code where the number of simulation replications comes in. With more replications, the objective function becomes smoother and easier to optimize. Also, as you know, fixing the random seed so that the simulations use common random numbers across different objective function evaluations is needed to get convergence of parameters. I don’t see that in the MWE.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 15, 2019, 6:20am UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/27 "2019-06-15T06:20:29Z")

</div>

> [@mcreel](#):
>
> fixing the random seed so that the simulations use common random numbers across different objective function evaluations is needed to get convergence of parameters. I don’t see that in the MWE.

As I mentioned multiple times above, there are discrete outcomes (and branches) in the simulation. Common random numbers or seeding the RNG does not help with continuity at all, unfortunately.

If necessary, I can make an MWE that makes this feature prominent, but it is equally useful just to pretend that we can’t use CRN.

> [@mcreel](#):
>
> I don’t understand from the code where the number of simulation replications comes in. With more replications, the objective function becomes smoother and easier to optimize.

Sure, but essentially we are trading runtime for O(1/\sqrt{n}) noise. The question is eg whether it is better to run 16 repications for 1/4 of the noise, _or_ just use an algorithm that can make use of 16x samples as is. CRS2 seems to be very good at this.

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [June 15, 2019, 7:07am UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/28 "2019-06-15T07:07:22Z")

</div>

Common random numbers won’t help with continuity, but they will keep the global minimum at the same point in the parameter space across different trial parameter values. With random draws that change at every function evaluation, the global minimum moves at every iteration. This is the phenomenon known as “chatter”, and Daniel McFadden’s papers discussed the need to avoid it.

Smoothing the objective function is a possibility. I think that’s a secondary issue, though.

Increasing the number of simulations affects the asymptotic variance of the parameter estimates, lowering it. You can use antithetic random draws, and other tricks. This is unrelated to the number of iterations that the minimization algorithm uses. Being able to try many different parameter values for an objective function that is less informative than is possible will not lead to the best possible parameter estimates. Trying to make the objective smoother (more simulations for each parameter value) and stable (common random numbers across trial parameter value) should help, I would think.

---

<div class="post-metadata">

**Author:** ![robsmith11](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robsmith11/32/29641_2.png) [@robsmith11](https://discourse.julialang.org/u/robsmith11)\
**Post date:** [June 15, 2019, 4:13pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/29 "2019-06-15T16:13:56Z")

</div>

I don’t think CRN are appropriate here. The goal is to optimize the expected value of the noisy objective function, not to find the minimum value out of all the evaluations (since the location of that point will be heavily influenced by the idiosyncratic noise), so it shouldn’t be a problem if the global minimum moves around as a result of the noise.

As the noise gets higher, I think it actually becomes increasingly important to explicitly model the noise with samples using the same parameter values (but different RNG seeds). Otherwise, how would one distinguish between variation in the objective caused by idiosyncratic noise vs variation caused by change in parameter values (what one is actually trying to optimize over)?

A couple references on the value of using replications with Bayesian Optimization:

> **[hetGP\_vignette.pdf](https://cran.r-project.org/web/packages/hetGP/vignettes/hetGP_vignette.pdf)**
>
> 970.77 KB

> **[Constrained Bayesian Optimization with Noisy Experiments](https://arxiv.org/abs/1706.07094)**
>
> Randomized experiments are the gold standard for evaluating the effects of changes to real-world systems. Data in these tests may be difficult to collect and outcomes may have high variance, resulting in potentially large measurement error. Bayesian...

If using CRN resulted in noise that were highly correlated between points nearby in parameter space, I could see some value in that approach, but that’s not the case at least for my simulations, which also have discrete actions and a high amount of branching.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 16, 2019, 5:17am UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/30 "2019-06-16T05:17:11Z")

</div>

> [@robsmith11](#):
>
> As the noise gets higher, I think it actually becomes increasingly important to explicitly model the noise with samples using the same parameter values (but different RNG seeds).

I agree with you and Bayesian optimization was the first thing I started with. But in practice, it is very slow and makes no steps most of the time once _somewhat_ near the optimum. Maybe I am tuning it wrong.

That said, I don’t think my noise is that large.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [June 16, 2019, 12:48pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/31 "2019-06-16T12:48:06Z")

</div>

Did you try the SPSA method I linked above? It works very well for this exact kind of problem.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 16, 2019, 12:57pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/32 "2019-06-16T12:57:12Z")

</div>

Not yet, it is on my TODO list. Is there a recommended Julia implementation?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [June 16, 2019, 1:30pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/33 "2019-06-16T13:30:17Z")

</div>

not that I know of, but if you point me to a howto on registering a package I’d be happy to create one, it’s like 20-50 likes of code and is v useful.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 16, 2019, 2:28pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/34 "2019-06-16T14:28:49Z")

</div>

Thanks, I want to understand it so I think I will just implement it.

---

<div class="post-metadata">

**Author:** ![asprionj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asprionj/32/6856_2.png) [@asprionj](https://discourse.julialang.org/u/asprionj)\
**Post date:** [June 16, 2019, 8:58pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/35 "2019-06-16T20:58:04Z")

</div>

Just some ideas:

- Particle-Swarm Optimisation: Had some good results on various “black-box” problems with this, but probably will require too many iterations. (Also can’t point at any reference implementation right now…)
- DIRECT (e.g. [https://github.com/npinto/direct](https://github.com/npinto/direct))
- Some sort of surrogate modelling ([Wikipedia](https://en.wikipedia.org/wiki/Surrogate_model))?
- Sobol sampling and Gaussian Process Model (GPM) as surrogate: Did this once for a parameter optimisation that used a very expensive simulation (hours per run) as evaluations. We just sampled the parameter space using N Sobol points (we could obtain over night), trained a GPM and figured out if the model accuracy was already stagnating (by cross-validation). Went on with sampling (overnight) until the model accuracy was not improving anymore, then applied a (very expensive) global optimisation algorithm to the GPM which evaluated lightyears faster than the original simulation… don’t know whether this could be successfully applied to such high-dimensional problems though. Maybe sparse GPM’s or a simpler type of surrogate model?
- Polynomial Chaos Expansion (PCE, as surrogate model): according to a colleague of mine, PCE yields reliable (accurate) models with even less samples than e.g. GPM’s. They use it for surrogate modelling of astronomic processes… Just some pointers I found by a quick google:
  - [Polynomial chaos expansion for surrogate modelling: Theory and software](https://www.researchgate.net/publication/327594993_Polynomial_chaos_expansion_for_surrogate_modelling_Theory_and_software)
  - [Data-driven polynomial chaos expansion for machine learning regression](https://www.ethz.ch/content/dam/ethz/special-interest/baug/ibk/risk-safety-and-uncertainty-dam/publications/reports/RSUQ-2018-005.pdf)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 16, 2019, 11:47pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/36 "2019-06-16T23:47:33Z")

</div>

> [@asprionj](#):
>
> - Some sort of surrogate modelling ([Wikipedia](https://en.wikipedia.org/wiki/Surrogate_model))?
> - Sobol sampling and Gaussian Process Model (GPM) as surrogate: Did this once for a parameter optimisation that used a very expensive simulation (hours per run) as evaluations. We just sampled the parameter space using N Sobol points (we could obtain over night), trained a GPM and figured out if the model accuracy was already stagnating (by cross-validation). Went on with sampling (overnight) until the model accuracy was not improving anymore, then applied a (very expensive) global optimisation algorithm to the GPM which evaluated lightyears faster than the original simulation… don’t know whether this could be successfully applied to such high-dimensional problems though. Maybe sparse GPM’s or a simpler type of surrogate model?

Note that we are building out [GitHub - SciML/Surrogates.jl: Surrogate modeling and optimization for scientific machine learning (SciML)](https://github.com/JuliaDiffEq/Surrogates.jl) for this as part of JSoC.

---

<div class="post-metadata">

**Author:** ![tcovert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tcovert/32/3855_2.png) [@tcovert](https://discourse.julialang.org/u/tcovert)\
**Post date:** [June 17, 2019, 3:13pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/37 "2019-06-17T15:13:58Z")

</div>

Just want to echo @mcreel’s points: if you use different random number draws for each function evaluation, whatever optimizer you pick is going to have a really hard time convincing itself that its found even a local optima. Further, regardless of whether the function that generates outcomes at the “true” parameter value is smooth in the realizations of the random draws, what matters for asymptotics is whether its expectation is smooth, which it probably (???) is.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 17, 2019, 3:45pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/38 "2019-06-17T15:45:07Z")

</div>

> [@tcovert](#):
>
> if you use different random number draws for each function evaluation, whatever optimizer you pick is going to have a really hard time convincing itself that its found even a local optima

I experimented with this (setting the random seed before each simulation), and it makes no difference.

The simulation has about 10⁶ branches, each affecting the number of random draws needed, so practically it is IID even with a small difference in inputs.

I really thought that this is the textbook use case for Bayesian optimization, and was surprised to see that CRS2 dominates it.

---

<div class="post-metadata">

**Author:** ![tcovert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tcovert/32/3855_2.png) [@tcovert](https://discourse.julialang.org/u/tcovert)\
**Post date:** [June 17, 2019, 4:44pm UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/39 "2019-06-17T16:44:07Z")

</div>

interesting. I think if you were to stick to the textbook mcfadden approach here, you’d really need to pre-draw the entirety of every draw that you could need, for every identical path of branching. This would guarantee that if two different parameter values result in the same branching sequence for some part of your data, the contribution of that part of your data to the objective function would need to be the same. maybe your code is written in such a way that fixing the seed does this, but if different parameter values result in the same branches happening at different points of time during the simulation, maybe not…

---

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [June 18, 2019, 12:36am UTC](https://discourse.julialang.org/t/optimizing-noisy-objective/25260/40 "2019-06-18T00:36:35Z")

</div>

Random thought: would a splittable PRNG help here? For example, it looks like [SplitMix](https://dl.acm.org/citation.cfm?doid=2660193.2660195) might allow you to have the same random numbers on a given branch of the simulation.

[Previous page](https://discourse.julialang.org/t/optimizing-noisy-objective/25260.md?page=1)

[Next page](https://discourse.julialang.org/t/optimizing-noisy-objective/25260.md?page=3)