How to handle randomness in optimizer

[se4u]

Hi, I have prevously worked¹ on a gradient-free optimization algorithm called SPSA²-³, and I have Matlab/mex code⁴ that I can port to python easily. I am interested in benchmarking SPSA against other zero-order optimization algorithms using nevergrad.

I am following the instructions for benchmarking a new optimizer given in adding-your-own-experiments-andor-optimizers-andor-function. My understanding is that I can just add a new SPSA class in nevergrad/optimization/optimizerlib.py and implement

• __init
• _internal_ask
• _internal_provide_recommendation
• _internal_tell

functions and then add SPSA to the optims variable in the right experiment function in the nevergrad/benchmark/experiments.py module and then I should be able to generate graphs like
docs/resources/noise_r400s12_xpresults_namecigar,rotationTrue.png

However, SPSA itself uses an rng in the _internal_ask function. But the optimizer base class does not take any seed in the __init__ function. What will be a good way to make the experiments reproducible in such situation?

[1] Pushpendre Rastogi, Jingyi Zhu, James C. Spall (2016). Efficient implementation of Enhanced Adaptive Simultaneous Perturbation Algorithms. CISS 2016, pdf
[2] https://en.wikipedia.org/wiki/Simultaneous_perturbation_stochastic_approximation
[3] https://www.chessprogramming.org/SPSA
[4] https://github.com/se4u/FASPSA/blob/master/src/SPSA.m

[tomjorquera]

@se4u not one of the authors here, but it seems the lib is relying on random and numpy.random for these matters. In both cases, you can set the seed globally by doing random.seed(1234) or np.random.seed(1234) to make subsequent calls to rngs reproducibles.

You can see how the authors use this here:

nevergrad/nevergrad/benchmark/xpbase.py

Lines 119 to 124 in 7fb29af

	def _run_with_error(self) -> None:
	"""Run an experiment with the provided artificial function and optmizer
	"""
	if self.seed is not None:
	np.random.seed(self.seed)
	random.seed(self.seed)

So I guess if you use either base python or numpy random you should be able to get reproducible experiments the same way.

[jrapin]

Hi, we are happy to welcome new algorithms indeed!

Optimizer implementation

You are correct in your understanding: implementing these 4 functions should make it work. The only one that is absolutely necessary is _internal_ask though.

Randomness

As @tomjorquera has correctly seen (thanks for starting answering ;) ), we rely on seeding during the experiments to make the benchmark reproducible. However, this also means, we avoid seeding inside the optimizer, so that we can run several optimization with different seed. In a nutshell, as long as you do not seed inside the optimizer, everything should work fine. Also, please use np.random over the standard library random if you can, but this is not an absolute requirement (only for consistency)

Experiments

You have several options:

• you update an experiment generator as you mentioned, just adding the name to the optims list should work
• you create another experiment generator, you will be able to run it just as well with the nevergrad.benchmark command.
• you can actually implement a benchmark generator in an external file and use the --import argument of nevergrad.benchmark to import the file at runtime. This is if you do not want to submit the experiment.

My preference go to the second option, so that everybody can use your experiment, and so that it does not modify existing experiments (for reproducibility of the existing experiments)

In any case, I'll try to help you in the process if I am available (Christmas / New year period may slow down my answer speed though :s) and don't hesitate to ask more questions if anything remains unclear!

[teytaud]

Just mentioning that SPSA is very (very) welcome.
Due to complexity associated to the management of seeds in case of multiple runs and averaging, maybe it's better to not seed the optimizers, as previously discussed.

[se4u]

I added SPSA to my fork of nevergrad and compared it on basic, noise and illcond* experiments. I didn't do any hyper-parameter tuning. The hyper-parameter come straight from my earlier SPSA Matlab code. The basic SPSA with averaging never really outperforms TBPSA, but on the sphere function it is competitive with TBPSA on average. The basic first-order SPSA worked badly on illcond problems, so I am not showing those results here. However, there are some more complicated versions of SPSA that try to estimate the hessian and theoretically promise better behavior on ill-conditioned problems.

The results on noise and basic are as follows:

Fig 1a. noise/fight_namesphere	Fig 1b. noise/fight_all

Fig 2a. noise/xpresults_sphere_rotation_false	Fig 2b. noise/xpresults_sphere_rotation_true

Fig 3a. basic/fight_all	Fig 3b. basic/xpresults

The class for SPSA is the following, I hardcoded an rng inside the class so that experiments are reproducible without modifying anything outside the class.

@registry.register
class SPSA(base.Optimizer):
    def __init__(self, dimension: int, budget: Optional[int] = None, num_workers: int = 1) -> None:
        super().__init__(dimension, budget=budget, num_workers=num_workers)
        self.rng = np.random.RandomState(seed=1234)
        self.init = True
        self.idx = 0
        self.delta = self.ym = self.yp = None
        self.t = np.zeros(self.dimension)
        self.avg = np.zeros(self.dimension)
        return

    @staticmethod
    def ck(k) -> float:
        return 1e-1 / (k//2 + 1)**0.101

    @staticmethod
    def ak(k) -> float:
        return 1e-5 / (k//2 + 1 + 10)**0.602

    def _internal_ask(self) -> base.ArrayLike:
        k = self.idx
        if k % 2 == 0:
            if not self.init:
                self.t -= (self.ak(k) * (self.yp - self.ym) / 2 / self.ck(k)) * self.delta
                self.avg += (self.t - self.avg) / (k // 2 + 1)
            self.delta = 2 * self.rng.randint(2, size=self.dimension) - 1
            return self.t - self.ck(k) * self.delta
        return self.t + self.ck(k) * self.delta

    def _internal_tell(self, x: base.ArrayLike, value: float) -> None:
        setattr(self, ('ym' if self.idx % 2 == 0 else 'yp'), value)
        self.idx += 1
        if self.init and self.yp is not None and self.ym is not None:
            self.init = False
        return

    def _internal_provide_recommendation(self) -> base.ArrayLike:
        return self.avg

The code and plots are here - https://github.com/se4u/nevergrad

[jrapin]

I'm happy you were able to obtain results so quickly!
If you have any feedbacks on what could be improved, let us know. It may still take some time to make it completely painless to run, but I hope we will get to it ;)

If you are willing to submit a pull request, you are most welcome to it.
Concerning the random seed generator you use, this is indeed a good option I forgot to mention.
Do you absolutely need to seed it though? Since we seed externally when running the benchmark, we would get reproducibility even without providing a seed to your random generator.

[se4u]

@jrapin I used my own rng because I didn't want to affect the external rng by executing my own code. Also, this makes it easy to test in future what is the effect of the randomness inherent in the SPSA algo. versus the random noise in the problem. Also, this makes the class self-contained. No matter when you initialize it, its behavior will be consistent.

Regarding feedback, I think the code structure is pretty intuitive so that's great. I guess some guidelines about how to submit new algorithms will be useful. For example, I can just submit a pull request after adding the above class, but should I also update the experiments module? Typically anyone will want to compare a new algorithm to the existing ones on some standard problem and if a new experiment function needs to be created for every algorithm that is added then the experiments.py module will become bloated. In my repo I have modified the noise and basic experiments, should I also submit those as part of the PR? Some guidelines about these issues will make it easier to submit PRs.

One more thing is about asynchronous vs synchronous execution. I am not sure exactly what happens in asynchronous execution. SPSA code won't work in an asynchronous manner because of the ask/tell API. SPSA needs to get values of y(θ + δ) and y(θ - δ) and it needs to know the current iteration to chose a step size, so I am maintaining a state in the form of idx which tells me the current iteration index. Therefore I expect that the sequence of calls must be ask, tell, ask, tell, ...

[jrapin]

For testing, I tend to seed outside on a deterministic problem, I'll add some unittests doing just that eventually (we have some internally). Anyway, this is a minor detail, if you prefer it this way I'm fine with it, as long as we dont try to benchmark it on a deterministic problem (in this case, averaging on multiple runs would not make sense...).

I'll try to improve the README indeed, but I think knowing the code too well makes it difficult to realize what is difficult for a newcomer. That is why I am asking feedbacks ;) You are also welcome to improve the README is you want (but I know well that it is not the most interesting of works :D).

Concerning the experiments, you are raising good points that we add not considered so far. @teytaud any thoughts about it?

[jrapin]

As for the sequence of call, if the implementation does not support parallelization, just notify it with a no_parallelization class attribute
See here: https://github.com/facebookresearch/nevergrad/blob/master/nevergrad/optimization/base.py#L53
This is probably among the things that should be explained in the README ;)

[jrapin]

I have kept thinking about how the algorithms should be seeded. The more I think about it, the more I think it should not be seeded:

• for consistency with the current implementation (but this argument is not that important, we could update it
• because we expect stochastic algorithms to be actually stochastic, if we set a seed inside the implementation this assumption is broken.
• to obtain relevant statistics when benchmarking the algorithms on deterministic functions (we will provide some interesting test cases that are deterministic probably before the end of January).
• because we can seed anyway when we need it, this is what happens in the benchmarks, because in this case we do want each independent run to be repeatable.

So, as soon as I have some time (this is not the best time of year for this :D) I will prepare:

• guidelines on how to add an algorithm, as you pointed out this will become necessary. I'll try to detail why it should not be seeded (at least not by default).
• unit tests ensuring that setting the seed outside will make computation repeatable (to make sure the benchmarks will be reproducible). As I mentioned, I already run these tests internally, I'll just need to port them to this repo.

I'll post back here when this is ready ;). Sorry for changing minds about it, it is indeed an important decision you are pointing out here!

[teytaud]

I just wanted to mention that your results are super exciting!

The version of SPSA you are using is supposed to converge at rate simple-regret=O(1/n^alpha) for alpha=2/3 or ... ?

I guess some different noise models with a lot of dissymmetry in terms of variance could be more in favor of SPSA... not sure though.

[se4u]

@teytaud The asymptotic convergence rate using my hyper parameters (a=0.602, γ=0.101) is k^(-(0.602-2*0.101)/2) = 1/k^(1/5). So in your terminology alpha = 1/5.

There are some conditions which basically put restrictions on a, γ.¹ a ∈ [ 0.602, 1] and γ ∈ [a/6, a/2). The fastest rate that can be achieved is when a is highest, i.e. 1, and γ is lowest i.e. 1/6. With those hyperparameters the convergence rate will be alpha = 1/3. But Prof. Spall suggests to use the more conservative values and those are what I used. I think the analysis could probably be simplified and streamlined compared to his approach although I haven't tried it.

[1] Spall, James C. "Multivariate stochastic approximation using a simultaneous perturbation gradient approximation." IEEE transactions on automatic control 37.3 (1992): 332-341. and theorem 7.2 on google books

Also, I wanted to mention that regarding submitting a PR I am waiting for

1. guidelines on how to add an algorithm to the experiments.
2. unit tests ensuring that setting the seed outside will make computation repeatable
3. On my end, I need to rerun things after adding the no_parallelization flag and then remove the hardcoded rng.

[yahayatmab3]

[jrapin]

1. Done in Added guidelines to add algorithms #14. It's still a bit messy though, and the best way to benchmark against other is not clear to me yet. For now you may add a new benchmark in experiments.py (all experiments for which figures are referred somewhere have been moved to frozenexperiments.py to ensure reproducibility. It's probably not the best way to do it but it will do for a time). This is not strictly required though, the figures you have are already good anyway ;) If you have ideas on how to do this better, please open a new issue!
2. Done in Added repeatability test #12: all registered algorithms are automatically tested. You'll need to run the tests locally at least once for them to work though.
3. Up to you ;) (the no_parallelization flag should not change anything to the computation though, it just notifies the benchmark pipeline that it should not be run with more than 1 worker).

[se4u]

I am getting the following two errors during testing. The first one is about my class not supporting parallelization. The second is about reproducibility. I am not sure why the second error is happening.

I get the error regardless of whether I set a fixed seed in my optimizer rng or not

self._rng = np.random.RandomState()
self._rng = np.random.RandomState(1234)

both of the above cause error.

For convenience, I added a WIP PR #16

prastog3@b19 ~/project/nevergrad Wed Jan 02 13:01:49 pr err=1                                        
$ nosetests nevergrad/optimization/test_optimizerlib.py
..........capi_return is NULL                                                                        
Call-back cb_calcfc_in__cobyla__user__routines failed.
............................................E.................S.....capi_return is NULL
Call-back cb_calcfc_in__cobyla__user__routines failed.
.............................................F...............
======================================================================
ERROR: test_optimizers(SPSA) (nevergrad.optimization.test_optimizerlib.OptimizerTests)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/prastog3/.local/lib/python3.6/site-packages/genty-1.3.2-py3.6.egg/genty/genty.py", lin$ 364, in <lambda>
    *dataset
  File "/home/prastog3/project/nevergrad/nevergrad/optimization/test_optimizerlib.py", line 74, in t$st_optimizers
    check_optimizer(optimizer_cls, budget=300, verify_value=verify)
  File "/home/prastog3/project/nevergrad/nevergrad/optimization/test_optimizerlib.py", line 31, in c$eck_optimizer
    optimizer = optimizer_cls(dimension=2, budget=budget, num_workers=num_workers)
  File "/home/prastog3/project/nevergrad/nevergrad/optimization/optimizerlib.py", line 449, in __ini$__
    super().__init__(dimension, budget=budget, num_workers=num_workers)
  File "/home/prastog3/project/nevergrad/nevergrad/optimization/base.py", line 58, in __init__
    raise ValueError(f"{self.__class__.__name__} does not support parallelization")
ValueError: SPSA does not support parallelization

======================================================================
FAIL: test_optimizers_recommendation(SPSA) (nevergrad.optimization.test_optimizerlib.OptimizerTests)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/prastog3/.local/lib/python3.6/site-packages/genty-1.3.2-py3.6.egg/genty/genty.py", lin$ 364, in <lambda>
    *dataset
  File "/home/prastog3/project/nevergrad/nevergrad/optimization/test_optimizerlib.py", line 89, in t$st_optimizers_recommendation
    err_msg="Something has changed, if this is normal, delete "
  File "/home/prastog3/.pyenv/versions/3.6.6/lib/python3.6/site-packages/numpy/testing/_private/util$.py", line 973, in assert_array_almost_equal
    precision=decimal)
  File "/home/prastog3/.pyenv/versions/3.6.6/lib/python3.6/site-packages/numpy/testing/_private/utils
.py", line 789, in assert_array_compare
    raise AssertionError(msg)
AssertionError:
Arrays are not almost equal to 10 decimals
Something has changed, if this is normal, delete /home/prastog3/project/nevergrad/nevergrad/optimizat
ion/recorded_recommendations.csv and rerun to update the values
(mismatch 100.0%)
 x: array([-2.0152900820e-05, -2.0152900820e-05,  7.4004300435e-06,
        2.0152900820e-05])
 y: array([ 1.0604319978e-05, -2.3356762000e-05, -2.3356762000e-05,
        2.3356762000e-05])

----------------------------------------------------------------------
Ran 139 tests in 4.279s

FAILED (SKIP=1, errors=1, failures=1)

[jrapin]

Thanks! Let's follow up on the PR ;)

[jrapin]

Please close the issue if everything is fine for you.
And you're welcome to open new ones if need be ;)

[teytaud]

I should mention that after some more work SPSA turns out to be very good.
A key point is to adapt the scale.
I'd try to do more, but I already got some good results in SPSA which were not visible in the early tests.
https://drive.google.com/file/d/1VOjzHtDuFtdc34SqohvEvZzQMwSdjvsQ/view?usp=sharing
(discussion in the pinned post at https://www.facebook.com/groups/nevergradusers/ )
Thanks to the authors of the SPSA code.

[se4u]

oh wow, really good to hear. I am trying to read the pdf but I must apologize, I wasn't quite able to figure out how to adapt the scale. Also I looked at section 4: "Statistics over all benchmarks" and Section 5: "conclusion" and I didn't see SPSA being competitive. How should I read the results ?