JohannesBuchner / UltraNest

Fit and compare complex models reliably and rapidly. Advanced nested sampling.
https://johannesbuchner.github.io/UltraNest/
Other
142 stars 30 forks source link

Eggbox problem intractable for Ultranest #92

Closed rruizdeaustri closed 10 months ago

rruizdeaustri commented 1 year ago

Description

I'm trying to estimate the posterior for the eggbox problem.

For the 2D case the number of likelihood evaluations I need to converge is of order O(10k), both using the reactive and the NestedSampler options.

If I increase the number of dimensions, this blows up by orders of magnitude. For instance, the case 5D is intractable since I need at least order 100 M likelihoods !!

For multinest, polychord and dynesty one needs
about 1 M likelihood evaluations to reach the convergence of the algorithm.

There must be some flags that I'm missing.

This is the code:

def main(args):

def loglike(z):
    chi = (cos(z / 2.)).prod(axis=1)
    return (2. + chi)**5

def transform(x):
    return x * 10 * pi

import string
paramnames = list(string.ascii_lowercase)[:args.x_dim]

args.reactive = False
if args.reactive:
    from ultranest import ReactiveNestedSampler
    sampler = ReactiveNestedSampler(paramnames, loglike, transform=transform,
        log_dir=args.log_dir, resume='overwrite',
        draw_multiple=False, vectorized=True,
    )
    sampler.run(log_interval=20,
        max_num_improvement_loops=10, min_num_live_points=args.num_live_points,)
    sampler.print_results()
else:
    from ultranest import NestedSampler
    sampler = NestedSampler(paramnames, loglike, transform=transform,
        num_live_points=args.num_live_points, vectorized=True,
        log_dir=args.log_dir, resume='overwrite')
        #log_dir=None)
    sampler.run(log_interval=20)
    sampler.print_results()

if name == 'main': parser = argparse.ArgumentParser()

parser.add_argument('--x_dim', type=int, default=2,
                    help="Dimensionality")
parser.add_argument("--num_live_points", type=int, default=1000)
parser.add_argument('--noise', type=float, default=-1)
parser.add_argument('--log_dir', type=str, default='logs/eggbox')
parser.add_argument('--reactive', action='store_false')

args = parser.parse_args()
main(args)

Thanks !!

Roberto

JohannesBuchner commented 1 year ago

You can try setting max_num_improvement_loops=0 to prevent too much refinement.

JohannesBuchner commented 1 year ago

Maybe first analytically calculate the number of posterior modes for the 5D problem to appreciate the difficulty of this problem. Other samplers giving answers quickly does not mean their answers are reliable, which is something you may care about.

JohannesBuchner commented 10 months ago

It may be 486 modes in 5 dimensions.

Since it is a product over the parameters, and each axis has a cos(x/2) function which has three +1 peaks, I get 3^d positive interferences and then as many for negative interferences (-1 * -1 is also a peak).

This can be solved by convergence assuming the modes are similar, but if you set cluster_num_live_points (default 40), then it will want to add more live points for each newly found cluster.

Closing this for now, please reopen if there is something more to discuss / resolve.