Doubling steps with slice sampler

[segasai]

This is the PR and proposal to add a slice sampling doubling scheme from Neal 2003 (scheme number 4).

The reason for this is this test for example.

from __future__ import (print_function, division)
import numpy as np
import dynesty

nlive = 1000
printing = True  # get_printing()

def loglike(x):
    x1 = (x * 1e7)
    x2 = x1.astype(int)

    rng = np.random.default_rng(x2)

    logl = rng.uniform() * 100 + 1e-9 * (x2 - x1).sum()
    return logl

def prior_transform(x):
    return x

def test_pathology():
    ndim = 4
    rstate = np.random.default_rng(56432)
    SAMPLER = dynesty.DynamicNestedSampler
    kw = dict()
    sampler = SAMPLER(loglike,
                      prior_transform,
                      ndim,
                      nlive=nlive,
                      bound='multi',
                      sample='rslice',
                      rstate=rstate)
    sampler.run_nested(print_progress=printing, **kw)

test_pathology()

When one runs it the code hits the '>10000' expansions and then proceeds extremely slowly. The reason is that this rugged (random) likelihood gets very small 'scale' parameter when sampling, and occasionally when the ellipsoids get much smaller one can end up in the situation when the bounding ellipsoid and scale parameter are equally small -- this leads to the extremely high number of expansions in the slice algorithm. And since the tuning of the scale parameter is only done by a factor of two no more, it takes a long while for the algorithm to get back to normal speed.

The proposed patch just implements the doubling scheme for slice sampling but doesn't enable it yet.

The possible option to enabling it are

• always enable it
• create an option to enable it
• enable it as soon as we hit a 'nexpand> THRESHOLD' once.

I'm inclined to chose the last option, or maybe the first one. The only reason I'm not sure about the last option is additional complexity in tracking the event that we've hit the expansion threshold.

For the example code given above if I don't do anything, the code ~ hangs after a minute of running taking 2 minutes to just do one iteration: 8008it [00:30, 340.39it/s, batch: 0 | bound: 17 | nc: 124 | ncall: 225942 | eff(%): 3.529 | loglstar: -inf < 99.959 < inf8009it [02:58, 1.81s/it, batch: 0 | bound: 18 | nc: 3330173 | ncall: 3556115 | eff(%): 0.225 | loglstar: -inf < 99.959 < inf | logz: 95.337 +/- 0.060 | dlogz: 0.035 > 0.010]/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/sampling.py:378: UserWarning: The slice sample interval was expanded more than 10000 times warnings.warn( The total run takes ~ 7 min. While with the change, the sampling finishes in 56 seconds.

While comparing the number of calls/iteration (before the hanging) I can also see that the doubling scheme requires more function calls ~ 1.5 times more in this case, which would support not making this a default option.

[coveralls]

Pull Request Test Coverage Report for Build 2884197017

• 66 of 69 (95.65%) changed or added relevant lines in 2 files are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage increased (+0.06%) to 90.495%

---

Changes Missing Coverage	Covered Lines	Changed/Added Lines	%
py/dynesty/sampling.py	64	67	95.52%
<!--	Total:	66	69	95.65%	-->

Totals
Change from base Build 2863436648:	0.06%
Covered Lines:	3713
Relevant Lines:	4103

---

💛 - Coveralls

[segasai]

Just ping here @joshspeagle . I will aim to merge this, unless there are objections. I think this is a pure improvement, but some implementation details may be questioned (i.e. the switch from default slice sampling to doubling that happens automatically, rather then through a separate option)

[joshspeagle]

Hi @segasai,

Adding this option looks good and I agree it is a pure improvement. It might be good just to pass a message if/when the switch occurs, but I'm happy with the automatic implementation being chosen based on efficiency.

[segasai]

Hi @segasai,

Adding this option looks good and I agree it is a pure improvement. It might be good just to pass a message if/when the switch occurs, but I'm happy with the automatic implementation being chosen based on efficiency.

Thanks @joshspeagle. Good idea. I've added the message and I will release this as 1.2.4 if all the tests pass.

[segasai]

While checking that the doubling is not coherently enabled, I realized that here https://github.com/joshspeagle/dynesty/blob/b65fbd8f9fba36eea881048b010d1cdd50dc84e0/py/dynesty/sampler.py#L366 the update_proposal is essentially called only every queue_size times and the majority of 'blob's returned by the sampler are just ignored. I didn't quite realize that. That's not great from the tuning point of view especially if queue_size is large (and also if some parts of the posterior are particularly problematic)

I'm not quite certain what's the best to deal with this. Options involve

• not doing anything
• somehow accumulate blobs from multiple calls
• update on every call

None of the options are ideal. One other option is to change update_rwalk, update_rslice etc methods to "accumulate" blobs. I.e. we run the update_method every time but we just collect the number of accepted points etc, but we only change the scale every so often (i.e. when accept+reject or contract+expand is large enough)

[segasai]

Also the current behaviour will lead to results depend on queue_size. I.e the results are deterministic conditional on identical queue_size.

[segasai]

this attempts to resolve the issue https://github.com/joshspeagle/dynesty/commit/54c21fe69f9ae720836bcbc6370e7fd01a8c41b2 by accumulating the acceptance info from all the queued samples and updating using all that information when the queue is empty. That changes the update_X() interface though, which is probably not a problem.