PyBADS: Frequently Asked Questions

PyBADS: Frequently Asked Questions#

This FAQ is curated by Luigi Acerbi, and in constant expansion. It is adapted for PyBADS 1.5 from the MATLAB BADS FAQ, with further questions on the Python package.

For a tutorial with detailed examples, see the Jupyter notebook examples.

If you have questions not covered here, please feel free to ask in the lab Discussions forum.

PyBADS is one of the lab’s tools for fitting models to data, which also give the posterior and the model evidence (PyVBMC) and the log-likelihood of models that can only be simulated (PyIBS).

The snippets below use np for NumPy and BADS for the optimizer class:

import numpy as np
from pybads import BADS

Supply your objective function fun, starting point x0, and bounds lb, ub, plb, pub where they appear in a snippet. The bounds are the arguments of BADS:

Variable	`BADS` argument
`lb`	`lower_bounds`
`ub`	`upper_bounds`
`plb`	`plausible_lower_bounds`
`pub`	`plausible_upper_bounds`

Give the starting point and each bound as a one-dimensional NumPy array of D elements, one per variable.

Table of contents#

General
Installing PyBADS
Input arguments (objective function: fun)
Input arguments (domain: x0, lb, ub, plb, pub, non_box_cons)
Output arguments
Noisy objective function
Display
Troubleshooting
Miscellanea

General#

Which kind of problems is PyBADS suited for?#

We recommend PyBADS for problems in which:

the objective function landscape is rough (nonsmooth), typically due to numerical approximations or noise;
the objective function is at least moderately expensive to compute (e.g., more than 0.1 s per function evaluation);
the gradient is unavailable;
the number of input parameters is up to about D = 20.

If your objective function is fully analytical, PyBADS is most likely not suited for your problem (see below).

The performance of BADS on the real model-fitting problems reported in the paper is remarkable. Did you cherry-pick the results?#

No, but we selected projects that we thought BADS would have been suitable for. The benchmark of the paper ran the MATLAB implementation of BADS, and some of its problems have been replicated with PyBADS.

What do I do if PyBADS is not suited for my problem?#

If the objective function is smooth and analytical, we would recommend a gradient-based optimizer instead, such as those of scipy.optimize.minimize (possibly feeding it the analytically calculated gradient). If you can afford tens or even hundreds of thousands of function evaluations, CMA-ES with active covariance adaptation, as implemented in the cma package, is also a valid alternative.

In these cases, you may also consider computing the full posterior, instead of getting only a point estimate via optimization. To this end, you could use Markov Chain Monte Carlo, e.g. via Stan or PyMC. Alternatively, we developed a method to compute approximate posterior distributions, Variational Bayesian Monte Carlo (PyVBMC), which can be used in synergy with PyBADS.

The lab’s page of tools for fitting models to data says which of them a problem needs: PyBADS for the best-fitting parameters, PyVBMC for the posterior and the model evidence, best with up to about 10 parameters, and PyIBS for the log-likelihood of a model that can be simulated but not written down, for data with discrete responses, which PyBADS or PyVBMC then take as a noisy target.

Installing PyBADS#

Where can I download PyBADS?#

Install or upgrade PyBADS with pip:

python -m pip install --upgrade pybads

Or install with Conda:

conda install --channel=conda-forge pybads

See the installation instructions for more details, and the GitHub repository for the source code.

Which external packages does PyBADS require?#

PyBADS 1.5 runs on NumPy 2.0, SciPy 1.13 and matplotlib 3.9 or newer, and builds its Gaussian process models with gpyreg 1.4.0 or newer, the Gaussian process library of our lab. pip and conda install them together with PyBADS. PyBADS does not require MATLAB.

To run the example notebooks you also need Jupyter (see the installation instructions).

Which version of Python do I need?#

PyBADS requires Python 3.10 or newer, and its package on conda-forge Python 3.11 or newer.

How do I know whether a newer version of PyBADS exists?#

Run pybads.check_for_updates(). It asks PyPI for the latest release and, if yours is older, prints the command that updates it: python -m pip install --upgrade pybads, or conda update --channel=conda-forge pybads for conda.

When your release is more than a year old, a run may also remind you to check, without contacting PyPI; options={"show_tips": False} turns the reminder off, with the tips. The check_for_updates API has the details.

Conda installs an older version of PyBADS. Why?#

PyBADS 1.5 requires NumPy 2.0 or newer. When the environment holds NumPy 1.x, or a package that requires it, conda install --channel=conda-forge pybads leaves NumPy as it is and installs the newest PyBADS that fits it, 1.1.0, without a warning. python -m pip install --upgrade pybads installs or upgrades PyBADS and upgrades NumPy when required. Ask conda for the latest release:

conda install --channel=conda-forge "pybads>=1.5"

Conda then upgrades NumPy, or says which package holds it back. A new environment avoids the conflict:

conda create --name pybads --channel=conda-forge python=3.12 pybads

If PyBADS is installed already, conda install keeps it as it is, and conda update --channel=conda-forge pybads updates it. The packages on conda-forge need Python 3.11 or newer: with Python 3.10, conda installs an older release still.

I am having trouble installing PyBADS. Can you help?#

Sure. The PyBADS installation should be pretty straightforward, so tell us in detail which problem you are having in the Discussions forum. Include your operating system, Python version, installation command and the full error message.

Input arguments (objective function: `fun`)#

What is the objective function?#

The objective or target function is the function that you want PyBADS to minimize, a Python callable fun passed as the first argument of BADS.

For a typical model-fitting problem, fun is a function that computes the negative log likelihood of an input parameter vector x, for a given dataset and model.

fun takes a one-dimensional NumPy array of shape (D,), a single point, and returns a finite real scalar: PyBADS evaluates one point per call. For example, to fit the mean and the log standard deviation of normally distributed data:

from scipy.stats import norm

data = np.random.default_rng(0).normal(1.0, 2.0, size=100)


def fun(x):
    mu, log_sigma = x  # x has shape (2,)
    return -np.sum(norm.logpdf(data, loc=mu, scale=np.exp(log_sigma)))


bads = BADS(fun, x0, lb, ub, plb, pub)
optimize_result = bads.optimize()
x_min = optimize_result["x"]
fval = optimize_result["fval"]

A noisy objective that can estimate its own noise returns a pair instead, as explained below.

Why the negative log likelihood?#

By mathematical convention, PyBADS minimizes the objective function, as most other optimization algorithms.

In the typical model-fitting scenario we want to maximize the likelihood. Which is the same as maximizing the log likelihood. Which is the same as minimizing minus the log likelihood, aka the negative log likelihood.

More generally, to maximize a function g, minimize -g, and flip the sign of the returned fval:

bads = BADS(lambda x: -g(x), x0, lb, ub, plb, pub)
optimize_result = bads.optimize()
g_max = -optimize_result["fval"]

My objective function requires additional data/inputs. How do I pass them to PyBADS?#

Suppose that your function takes two inputs, fun(x, data).

The first solution consists of defining a new function

def funwdata(x):
    return fun(x, data)

where data has been defined before in the code. Now you can optimize funwdata, which takes a single input.

Alternatively, you can use functools.partial:

from functools import partial

funwdata = partial(fun, data=data)
bads = BADS(funwdata, x0, lb, ub, plb, pub)

which binds the data argument to the object passed.

Input arguments (domain: `x0`, `lb`, `ub`, `plb`, `pub`, `non_box_cons`)#

How do I choose the starting point `x0`?#

First of all, keep in mind that you should restart PyBADS from different starting points. Probably a minimum of ten, ideally dozens, depending on your problem (see the next question).

We recommend to choose starting points mostly inside the plausible box bounded by plb and pub. Pass x0=None and PyBADS draws the starting point at random inside the plausible box, a draw that the seed of the run decides. Or draw it yourself, for example

rng = np.random.default_rng()
x0 = rng.uniform(plb, pub)

If you think that you would like to also draw points outside plb and pub, then by definition it means that your choice of plb and pub is too narrow (see also below).

How do I run PyBADS from several starting points?#

Create a new BADS object for each run: a BADS object runs a single optimization. For example, with a random starting point and a different seed for each run:

n_runs = 10
results = []
for seed in range(n_runs):
    options = {"random_seed": seed, "display": "off"}
    bads = BADS(fun, None, lb, ub, plb, pub, options=options)
    results.append(bads.optimize())

best = min(results, key=lambda r: r["fval"])
x_best, fval_best = best["x"], best["fval"]

Then compare the runs. If several of them reach nearly the same fval, you can be more confident about the solution; if they are scattered, see this question for a deterministic objective and this one for a noisy one. For a noisy objective, fval is itself an estimate, with standard deviation fsd: runs whose values differ by less than a few fsd cannot be told apart, but evaluating their solutions several more times can tell them apart.

The runs are independent of each other, so you can also run them in parallel processes, for instance with concurrent.futures.ProcessPoolExecutor:

# multistart.py
from concurrent.futures import ProcessPoolExecutor

import numpy as np
from pybads import BADS


def fun(x):
    return np.sum(x**2)


lb, ub = np.full(2, -5.0), np.full(2, 5.0)
plb, pub = np.full(2, -2.0), np.full(2, 2.0)


def run(seed):
    options = {"random_seed": seed, "display": "off"}
    return BADS(fun, None, lb, ub, plb, pub, options=options).optimize()


if __name__ == "__main__":
    with ProcessPoolExecutor() as pool:
        results = list(pool.map(run, range(10)))

Each worker process imports this file, so define the objective, the bounds and run at its top level, as here, and put only the code that starts the runs under if __name__ == "__main__":; a lambda or a function defined inside another function cannot be sent to a worker. From a Jupyter notebook, put these definitions in a .py file and import them.

How do I choose `lb` and `ub`?#

lb and ub are the hard bounds of the optimization. In theory, you could set them to the mathematical limits of your variables. However, using the mathematical limits of a variable is a bad choice for optimization. Instead, we recommend to set them to no wider than their physical or experimental limits.

For example, suppose that you have a parameter sigma that represents the standard deviation (SD) of the movement endpoint of a subject in a task in which people are asked to rapidly touch targets on a screen. Mathematically, sigma, being a SD, could go from 0 to np.inf. However, in this case, it is physically unrealistic that people would have no motor noise. Instead, we set as lb our experimental lower bound, e.g., the resolution of our motion tracker device, or maybe one screen pixel. Similarly, it is physically impossible for people’s pointing error to be larger than, say, the length of their forearms. In fact, we could set as ub the size of the screen.

Importantly, do not set lb to 0 for variables that can only be positive. Choose a small, experimentally meaningful number. Do not pick extremely small numbers such as the machine epsilon, np.finfo(float).eps (about 2e-16), unless they are justified in the context of your problem. A positive lower bound also lets PyBADS work on such a variable in log coordinates.

What if I really have no idea how to choose `lb` and `ub`?#

It is true that occasionally some model parameters might not have an a priori intuitive range of values. One could gain a bit of intuition via preliminary exploration of the function landscape (i.e., manually set some values). Then, you could set bounds to some mid-to-large values, and expand them if needed. We would still not recommend to set incredibly large bounds.

Does PyBADS support (partially) unconstrained optimization?#

Yes and no. You can specify that a variable is (partially) unconstrained by setting its hard bounds to -np.inf or np.inf; passing None for lb or ub leaves every variable unbounded on that side. The plausible bounds of such a variable must then be given, and finite.

However, we encourage users to always set finite, empirically meaningful hard bounds (see above). Infinities are never empirically meaningful, unless perhaps if you are in a black hole.

Can I set `lb = ub` for some variable to fix it to a given value?#

Yes: a variable whose four bounds, lb, ub, plb and pub, are equal is fixed at that value. Where plb and pub are not given they are the hard bounds, so that lb = ub alone fixes the variable. At a fixed variable, give x0 that value, or NaN. For example, to fix the second of three variables at 2:

lb = np.array([-5.0, 2.0, -5.0])
ub = np.array([5.0, 2.0, 5.0])
plb = np.array([-2.0, 2.0, -2.0])
pub = np.array([2.0, 2.0, 2.0])
x0 = np.array([0.0, 2.0, 0.0])

bads = BADS(fun, x0, lb, ub, plb, pub)  # optimizes x[0] and x[2]
optimize_result = bads.optimize()

PyBADS optimizes only the other variables, and the defaults of the options that depend on the number of variables count only those: here max_fun_evals is 500 * 2. PyBADS lists the fixed variables when you create the BADS object. Your objective, non_box_cons and the output function receive points of all the variables, with the fixed ones at their values, so your code needs no change. The result’s x and x0 hold all the variables too, as do the points of precomputed_evaluations, and the indices of options["periodic_vars"] count all of them.

How do I choose `plb` and `pub`?#

plb and pub are the plausible (or reasonable) bounds of the optimization. Set them by thinking of a plausible range in which you would expect to find almost all solutions; as a rule of thumb, you would bet that the minimum lies inside the plausible box they define with probability above 90%. The plausible box naturally represents a good region where to randomly draw starting points for the optimization (see above). The plausible bounds must be finite and satisfy lb <= plb < pub <= ub at every variable that is not fixed.

If you really have no idea about a plausible range, you can set plb and pub equal to lb and ub, or leave them out, in which case PyBADS uses the hard bounds in their place, with a warning that says so; but this should not be the norm.

In the example above (see this question), the plausible bounds for sigma (pointing motor noise) could go from a few pixels for plb to several cm for pub.

The plausible box also sets the scale of the optimization: PyBADS draws its initial design of points inside it, and measures its steps in coordinates in which the box spans [-1, 1] (see the next question).

Does PyBADS rescale or transform my variables?#

Yes, internally. Your objective always receives, and the result always reports, points in the coordinates you gave, but PyBADS runs in coordinates of its own:

it maps the plausible box to [-1, 1] in each variable, so that its steps, and the MeshScale of its display, are relative to the plausible range of each variable;
before that, it takes the logarithm of every variable whose bounds are all positive and whose plausible range spans a factor of 10 or more (pub / plb >= 10), unless the variable is periodic. Its steps along such a variable are then proportional to the variable’s value, which suits scale parameters such as standard deviations, rates or time constants. PyBADS lists these variables when you create the BADS object.

To keep every variable on a linear scale, set options={"nonlinear_scaling": False}. Otherwise, you rarely need to rescale the variables yourself; a parameterization in which the parameters trade off less against each other still makes the problem easier for any optimizer.

How do I prevent PyBADS from evaluating certain inputs or regions of input space?#

If these regions can be identified by coordinate-wise ranges, use lb and ub. Otherwise, use the barrier function non_box_cons (non-box constraints). It takes an array of shape (N, D), one point per row, and returns an array of N values, true (or positive) for each point that is not allowed. PyBADS does not evaluate the objective at such points. For example, to keep the optimization inside the unit ball:

def non_box_cons(X):
    return np.sum(X**2, axis=1) > 1


bads = BADS(fun, x0, lb, ub, plb, pub, non_box_cons=non_box_cons)

The starting point x0 must satisfy the constraints. See Example 2 for a full example.

PyBADS draws its initial points in the plausible box and leaves out those that violate the constraints, so choose a plausible box of which a good part is allowed: in the example above, a box within -1 and 1 rather than a much wider one. A feasible region much thinner than the steps of the optimization, such as a narrow band, can leave PyBADS unable to find allowed points around x0, and the run then ends early: reparameterize such a problem so that its feasible region is wide (the documentation of BADS gives an example). For the same reason, PyBADS does not support equality constraints, such as x[0] + x[1] = 1: write the problem in fewer variables, so that the constraint holds by construction (here, optimize over x[0] alone and set x[1] = 1 - x[0] inside the objective).

Do not have fun(x) return np.inf or np.nan for invalid inputs. PyBADS would simply crash (see this question).

Absolutely do NOT have fun(x) return an arbitrarily large number to enforce PyBADS to avoid certain regions. This strategy may seem innocent enough, but in fact it completely cripples PyBADS by making its models of the objective function nonsensical.

Does PyBADS support integer constraints?#

No, PyBADS does not support integer constraints (that is, variables forced to be integers).

As a simple fix, you could adopt the following hack if you have a single integer variable m:

make the parameter m continuous for the purpose of optimization;
for a given parameter vector, evaluate separately the objective at np.floor(m) and np.ceil(m);
return the linearly interpolated value of the objective;
at the end of the optimization, either return np.round(m), or evaluate your function (several times, if noisy) at both np.floor(m) and np.ceil(m) and pick the best.

For example, with m the variable of index i_int:

i_int = 0  # index of the integer variable


def fun_interp(x):
    m = x[i_int]
    x_lo, x_hi = x.copy(), x.copy()
    x_lo[i_int], x_hi[i_int] = np.floor(m), np.ceil(m)
    if x_lo[i_int] == x_hi[i_int]:
        return fun(x_lo)
    w = m - x_lo[i_int]
    return (1 - w) * fun(x_lo) + w * fun(x_hi)

Note that this will double the cost of each evaluation, so it is worth only if alternative solutions (such as looping over all values of m) would be computationally more expensive. Finally, this approach makes sense only if m is truly an integer (ordered and with an underlying metric), as opposed to merely a categorical variable.

Does PyBADS support periodic variables, such as angles?#

Yes. List the indices of the periodic variables, counted from 0, in the option periodic_vars. The hard bounds lb and ub of a periodic variable are its period, and need to be finite. PyBADS takes lb and ub as the same point and wraps the variable around them, so that a run moves across them as across any other value of the variable, and it models the objective as periodic along the variable. PyBADS lists the periodic variables when you create the BADS object.

The bounds of a periodic variable only mark where its period is cut, not a range of extreme values, so its plausible bounds are usually its hard bounds as well, unlike those of other variables (How do I choose plb and pub?). For example, with an angle in radians as the second variable:

# x[0] is a position, x[1] an angle in radians
lb = np.array([-10.0, -np.pi])
ub = np.array([10.0, np.pi])
plb = np.array([-2.0, -np.pi])
pub = np.array([2.0, np.pi])

options = {"periodic_vars": [1]}
bads = BADS(fun, x0, lb, ub, plb, pub, options=options)
optimize_result = bads.optimize()

The result gives a periodic variable within its period, between lb and ub. A minimum at lb, which is the same point as ub, can come out at lb, just above it or just below ub, so two runs can report nearly the same minimum at opposite ends of the period.

Example 6 runs PyBADS on a function with two periodic variables.

Output arguments#

What does `optimize()` return?#

bads.optimize() returns an OptimizeResult, a dictionary whose entries you read as optimize_result["x"] or as attributes, optimize_result.x. Among them:

x: the solution, a one-dimensional array;
fval: the value of the objective at x, estimated for a noisy objective;
fsd: the standard deviation of that estimate (0 for a deterministic objective);
success, status and message: why the run ended (see below);
func_count and iterations: the number of function evaluations and of iterations;
total_time and overhead.

The documentation of OptimizeResult lists them all.

success is True when the run ended on one of its convergence criteria: the mesh became smaller than options["tol_mesh"] (status 1), or the objective improved by less than options["tol_fun"] over the last options["tol_stall_iters"] iterations (status 2). It is False (status 0) when the run used up its budget of evaluations (options["max_fun_evals"]) or of iterations (options["max_iter"]), or an output function stopped it. message says which. A run that used up its budget can still have found a good solution; if not, give it a larger budget.

How is `fval` computed?#

fval is the (estimated) value of the objective function at x, the returned optimum.

For a deterministic (not-noisy) objective, this is simply fun(x), the lowest value that PyBADS observed. For a noisy objective, x is the point, among those at which the iterations ended, whose value as estimated by PyBADS’s Gaussian process model at the end of the run is lowest after accounting for its uncertainty; it need not be the incumbent at the end of the run. PyBADS normally evaluates fun(x) options["noise_final_samples"] more times (default 10): fval is the average of those evaluations and fsd its standard error. With options["specify_target_noise"], the average and its standard error weight each evaluation by the precision that the objective reports. The evaluations are in optimize_result["yval_vec"] (and the standard deviations that the objective reported in optimize_result["ysd_vec"]). They count towards options["max_fun_evals"].

PyBADS reserves these evaluations from the budget left after initialization, so a small budget can reduce their number. With noise_final_samples=0, or when no budget remains for final samples, yval_vec and ysd_vec are None: fval and fsd keep the estimates available when the run stops, normally from the GP. With only one final sample and no specify_target_noise, PyBADS averages that sample with the observation already recorded at x, and yval_vec contains both. An output function that stops the run at "init" takes no final samples; if samples had been reserved, yval_vec holds the incumbent’s observation alone.

Why do you estimate `fval` by averaging additional function evaluations? Can’t you return the Gaussian process prediction at `x`?#

Glad that you asked. Yes, in theory we could use the Gaussian process (GP) mean prediction at x. However, the GP prediction can occasionally fail, sometimes subtly. While this mismatch is not a major problem during optimization, it could potentially introduce hard-to-detect but substantial biases in fval, which could have catastrophic effects for model selection. For this reason, we chose a more conservative approach for estimating fval.

How is `overhead` defined?#

optimize_result["overhead"] is the fractional overhead, defined as (total running time / total function time - 1). The total running time, optimize_result["total_time"], is the time of bads.optimize(), in seconds; the function time is the time spent evaluating the objective, excluding the automatic noise test at x0, which counts as optimizer time. PyBADS’s own work takes of the order of tens of milliseconds per function evaluation, depending on the problem and the computer.

Typically, you would expect overhead to be (much) smaller than 1 for normal runs of PyBADS. If the overhead is larger than, say, 0.75, your problem affords fast evaluations and it is possible that it would benefit from other algorithms than PyBADS (see above).

For PyBADS test problems and examples, you will find that the reported overhead is astronomical, which is expected since for demonstration purposes we are using simple analytical functions.

Where can I find the internal state and iteration history?#

After bads.optimize(), the BADS object keeps:

bads.optim_state, a dictionary with the state of the optimization;
bads.iteration_history, which records each iteration: among others the point at which the iteration ended, in your coordinates ("x") and in PyBADS’s internal coordinates ("u"), its value ("fval" and "fsd"), the mesh size ("mesh_size"), the number of evaluations so far ("func_count"), and the Gaussian process model ("gp"), which works in the internal coordinates;
bads.function_logger, the log of the evaluated points (see the next question).

These attributes are useful for debugging, but we are not providing explicit support for them. Future versions of PyBADS might change their interface or internal structure.

Is it possible to output or inspect the optimization trajectory?#

Yes. After bads.optimize(), the evaluated points and their observed values are in the function log, whose rows X_flag marks as filled:

log = bads.function_logger
X = log.X_orig[log.X_flag]  # evaluated points, one per row
y = log.Y_orig[log.X_flag].ravel()  # observed function values

Each row holds a point, in the order in which PyBADS first evaluated it. An evaluation at a point that is already in the log adds no row: the second evaluation of x0 that tests the objective for noise, and the final evaluations of a noisy run, which are in optimize_result["yval_vec"] (see above). So the log can have fewer rows than optimize_result["func_count"].

A run given evaluations made before it (the argument precomputed_evaluations of BADS) holds them in the first rows of its log, and does not count them in func_count.

The points at which the iterations ended, and their values, are

x_iter = np.vstack(bads.iteration_history["x"])
f_iter = bads.iteration_history["fval"]

In a noisy run, f_iter holds PyBADS’s estimates of these values as revised at the end of the run, so they can differ from the values shown in the display.

To record the trajectory while the run goes on, use an output function.

Can I monitor or stop a run while it is running?#

Yes, with an output function, passed as options["output_fcn"]. PyBADS calls it as output_fcn(x, optim_state, state) once the initial points have been evaluated (state is "init"), at the end of each iteration ("iter") and when the run ends ("done"). x is the current best point, in your coordinates, and optim_state a copy of the state of the optimization, with for instance optim_state["fval"], optim_state["fsd"] and optim_state["mesh_size"]. When the function returns True, the run stops, and its result says so in message, with success set to False.

For example, to record the point at the end of each iteration, and to stop the run after an hour:

import time

start = time.perf_counter()
trace = []


def output_fcn(x, optim_state, state):
    if state == "iter":
        trace.append((x.copy(), optim_state["fval"]))
    return time.perf_counter() - start > 3600


bads = BADS(fun, x0, lb, ub, plb, pub, options={"output_fcn": output_fcn})
optimize_result = bads.optimize()

Noisy objective function#

What is a noisy objective function?#

A noisy (or stochastic) objective function is an objective that will return different results if evaluated twice at the same point x. A non-noisy objective function is deterministic.

For model fitting, objective functions can be noisy if the log likelihood is evaluated through simulation (e.g., via Monte Carlo methods).

Why are noisy objective functions treated differently?#

For a deterministic objective, we assume that the goal is to minimize f(x). For a noisy objective, we assume that the goal is to minimize the expected value of f(x), also written as E[f(x)]. For this reason, PyBADS will not simply blindly trust whatever f(x) returns, but will do some internal computation to estimate E[f(x)] (effectively, smoothing the observed function values via a Gaussian process).

Incidentally, this means that ideally the function that you provide (and that computes the negative log likelihood) should be an unbiased estimator of the negative log likelihood, such as the one that inverse binomial sampling gives (see How do I estimate the standard deviation of a noisy objective?).

Can I make a noisy objective function deterministic by fixing the noise process?#

Well, technically yes, you could make a noisy objective function deterministic by seeding its random number generator again (e.g., with np.random.default_rng(0)) every time you call it. However, this fix does not really solve the problem, because you are not eliminating the noise in the function observations. In fact, if you do it naively, you might be adding unwanted bias to your fits.

Thus, it is not recommended to ‘fix’ the noise this way (by fixing the random seed at each function call). Instead, let your function be stochastic, and let PyBADS deal with it.

Note that this is different from setting the random seed once at the beginning of an optimization run, for the sake of reproducibility, which is recommended as good practice. The seed of PyBADS, options["random_seed"], governs only the random draws of PyBADS itself, so give your objective a random generator of its own, created once:

rng = np.random.default_rng(12345)  # created once, for the whole run


def fun(x):
    return simulated_nll(x, rng)  # draws new noise at every call

where simulated_nll stands for your simulation-based estimate of the negative log likelihood.

Should I tell PyBADS that my objective is noisy?#

Please do so. Set options={"uncertainty_handling": True} to tell PyBADS that the optimization is noisy.

If you forget about it, PyBADS will determine at initialization whether the provided objective is noisy, by evaluating it twice at x0 (an evaluation that counts towards the budget). Note that this test can fail if a noisy objective happens to return the same value twice, which is rare unless its values take only a few distinct levels. The line Beginning optimization of a STOCHASTIC objective function (or DETERMINISTIC) at the start of the display, and optimize_result["target_type"], say which PyBADS concluded.

Conversely, set options={"uncertainty_handling": False} for a deterministic objective: PyBADS then skips the test, and saves an evaluation. With options={"specify_target_noise": True} (see below), PyBADS treats the objective as noisy without further ado.

Can PyBADS handle any arbitrary amount of noise in the objective?#

No. PyBADS works best if the standard deviation of the objective function, when evaluated in the vicinity of the solution, is small with respect to changes in the objective function itself (that is, there is a good signal-to-noise ratio). In many cases, a standard deviation of order 1 or less should work (this is the default assumption). If you approximately know the magnitude of the noise in the vicinity of the solution, you can help PyBADS by specifying it in advance (set options["noise_size"] = sigma_est, where sigma_est is your estimate of the standard deviation).

If the noise around the solution is too large, PyBADS will perform poorly. In that case, we recommend to increase the precision of your computation of the objective (e.g., by drawing more Monte Carlo samples) such that sigma_est is of order 1 or even lower, as needed by your problem (see also this related question). Note that the noise farther away from the solution can be larger, and this is usually okay.

Does PyBADS assume that the noise is the same for all inputs?#

Yes and no. The Gaussian process (GP) model built by PyBADS is homoskedastic, that is, it assumes constant noise across the input space. However, the GP model is built using only a local set of points, so PyBADS will adapt to local characteristics of the objective function, including amounts of noise that depend on the location.

You can help PyBADS optimize a heteroskedastic objective (i.e., with input-dependent noise) by providing an estimate of the noise at each location, as specified in the questions below.

Should I provide an estimate of the noise associated with each evaluation?#

Yes, you should if you can! This may considerably help PyBADS, especially if the objective is particularly noisy or strongly heteroskedastic. Remember to:

set options["specify_target_noise"] = True;
pass to PyBADS a function fun that returns a tuple (f, sd): the estimate of the objective at x, and an estimate of its standard deviation (SD), a finite positive number.

def fun(x):
    f, sd = estimate_nll(x)  # your estimate and its SD
    return f, sd


bads = BADS(fun, x0, lb, ub, plb, pub, options={"specify_target_noise": True})

Without specify_target_noise, PyBADS refuses such a pair: the objective must return a single number. See Example 4 for further information.

How do I estimate the standard deviation of a noisy objective?#

If you estimate the log-likelihood by inverse binomial sampling (IBS), with PyIBS, one of the lab’s tools for fitting models to data, the algorithm returns the variability of the estimate as second output. Just ensure that the variability is returned as standard deviation (SD) and not as the variance (depending on the implementation, you may have to take the square root of the reported variance). Otherwise, you can also estimate the SD via bootstrap or similar approaches.

Display#

What are the quantities displayed by PyBADS during optimization?#

With options["display"] set to "iter" (the default), PyBADS displays the traces of several optimization quantities:

the Iteration number;
the number of objective function evaluations f-count;
the value of f(x) at the incumbent (current point);
MeshScale, the size of PyBADS’s current steps relative to the plausible box; it shrinks as the run converges, and the run ends when it falls below options["tol_mesh"] (default 1e-6);
the current optimization stage and method, under Method: Initial mesh for the initial design; a Successful search, which improved on the incumbent sufficiently, or an Incremental search, which improved on it by less, with the search method in parentheses (ES-wcm or ES-ell); a Successful poll; or Refine grid, a poll that found no sufficient improvement (it can still move to a slightly better point), after which the mesh shrinks;
additional actions, under Actions, such as the Uncertainty test of the objective’s noise at x0, the evaluation of the Initial points, or re-training the Gaussian process (Train).

If the objective function is noisy, instead of f(x) PyBADS will report the expected value E[f(x)] and its standard deviation SD[f(x)] at the incumbent, both estimated via the current Gaussian process model. In the first lines, before PyBADS has built its model, E[f(x)] is the value observed at the incumbent and SD[f(x)] a placeholder (nan or a nominal value).

For a noisy function, I noticed that the series of displayed `E[f(x)]` values is not monotonically decreasing. Should I worry?#

Well spotted, but nothing to worry about. PyBADS keeps updating the estimate of E[f(x)] at the incumbent, which means that occasionally this value will increase across iterations, and sometimes it will oscillate for a few iterations. Also, if the noise is large, the Gaussian process approximation might occasionally fail, leading to outlier estimates for E[f(x)] (which should then recover in the subsequent iterations). All of this is part of the normal functioning of PyBADS.

For a noisy function, I noticed that the series of displayed `SD[f(x)]` sometimes shows sudden jumps (e.g., from ~1 to ~4). Is that normal?#

First, recall that SD[f(x)] is the estimated posterior standard deviation of the objective function at the incumbent (current best point). This estimate is obtained via the Gaussian process model built by PyBADS every few iterations. When the incumbent changes, or when PyBADS re-trains the Gaussian process model, the uncertainty about the value at the incumbent will also change, sometimes substantially (e.g., if the estimated observation noise parameter has changed, or if the incumbent has moved to a new region). In most cases, such jumps are part of the normal behavior of the algorithm.

Sometimes as `Actions` during optimization I read `Train (failed)`. What does it mean?#

It means that PyBADS was unable to refit the Gaussian process model to the current local training set, usually due to numerical issues. Occasional failures are not reason of concern, in particular at the beginning or towards the end of the optimization. However, if a large number of training attempts are systematically failing, it might mean that PyBADS is having trouble. Sometimes this can be fixed by changing the problem parameterization, or perhaps there are other issues with the model.

How do I silence PyBADS, or send its output elsewhere?#

options["display"] sets how much PyBADS prints: "iter" (the default) prints a line per step, as described above; "final" prints the opening and final messages; "notify" the opening messages alone; "off" nothing but warnings; "full" everything, debug messages included. With "iter" or "full", a run may also print a short tip with a link to the documentation, or, when the installed release is more than a year old, a reminder to check for a newer version; options["show_tips"] = False turns both off.

PyBADS prints through Python’s logging module, with a logger named "BADS", and its warnings, such as the one for unspecified plausible bounds, are log messages too, not Python warnings (the libraries it calls, such as NumPy and gpyreg, can still issue Python warnings of their own). Creating a BADS object calls logging.basicConfig, which sends the messages to the standard output unless your program has configured logging before; a logging.basicConfig call of your own after that takes effect only with force=True. To write PyBADS’s messages to a file instead:

import logging

logger = logging.getLogger("BADS")
logger.addHandler(logging.FileHandler("bads_run.log"))
logger.propagate = False  # not to the console as well

Troubleshooting#

Is there a way to check that PyBADS is running correctly — is it enough that it does not give warnings/errors?#

You can follow the run through its display, and inspect its state while it runs with an output function. And while the fact that no warnings/errors are shown is encouraging, it is not a sufficient condition to guarantee that everything ran correctly. For validation of the results, we recommend usual techniques such as comparing multiple independent runs of the algorithm, and various form of model checking.

PyBADS crashes saying that `The returned function value must be a finite real-valued scalar`. What do I do?#

This ValueError means that your objective function has returned np.inf, np.nan, or something that is not a single real number, such as a complex number or an array of several elements. You should check your code and understand why it returned such a value. A pair (f, sd) gives this error too, unless you set options["specify_target_noise"] = True (see this question).

infs and nans often arise because there are outcomes in your dataset (e.g., responses in a trial) to which the tested model assigns a probability of 0 (usually due to numerical truncation). np.log(0) yields -np.inf, which is then propagated. In these cases, we recommend to make the model more robust, by computing the likelihood directly in log space where you can (for instance with log-density functions such as scipy.stats.norm.logpdf, and scipy.special.logsumexp), or by forcing all outcomes (e.g., the likelihood associated with each trial) to have a minimum non-null probability, such as np.sqrt(np.finfo(float).eps) or some other small value. This should not be necessary if the model already includes a non-zero lapse rate.

nans also arise when you take np.log or np.sqrt of a negative number (of a quantity that should not be negative), which in NumPy gives nan with a RuntimeWarning. You might be setting wrong bounds for your variables, or maybe there are indexing issues.

Note that some optimizers are robust to infs and nans and just keep going, avoiding the problematic region. However, we believe this is dangerous as it might hide deeper issues with the model implementation.

During optimization I received a warning that `The mesh attempted to expand above maximum size too many times`. What does it mean?#

It means that probably your plb or pub bounds are too narrow; try widening them. If these are already as wide as lb and ub, it might be that your hard bounds are too narrow.

If you do not think that this is the case, you can disable this warning by setting options["mesh_overflow_warning"] = np.inf.

I am passing `non_box_cons` to PyBADS, but I get an error that `non_box_cons should be a function that takes an N x D array X`. What am I doing wrong?#

Most likely, your non_box_cons accepts only a single point and returns a single value, whereas you should be sure that it takes an array of shape (N, D), one point per row, and returns an array of N constraint violations.

For example, suppose that your input variables need to be ordered, such that x[0] <= x[1] and x[1] <= x[2]. Then you should set

def non_box_cons(X):
    return (X[:, 0] > X[:, 1]) | (X[:, 1] > X[:, 2])

whereas lambda x: x[0] > x[1] or x[1] > x[2] would yield an error.

I have been running PyBADS with a deterministic objective function from the same starting point, but I get different results each time. Is something wrong?#

Nothing is wrong per se. PyBADS is a stochastic optimizer, so results may differ between different runs, even with the same initial condition x0. If the returned fval varies substantially across runs from the same starting point, it might be a sign that your function landscape is particularly difficult. If fval is similar across runs, but the returned optimum x varies substantially, it is a sign that your function has a plateau or ridge, with trade-offs between parameters.

If you want to have reproducible results (and this advice applies beyond PyBADS), we recommend to fix the random seed of each run to some known quantity (e.g., setting it to the run number), as explained in the next question.

How do I make a run reproducible?#

Pass an integer seed when you create the BADS object:

bads = BADS(fun, x0, lb, ub, plb, pub, options={"random_seed": 42})

Every random draw of PyBADS during the run comes from this seed, including the starting point when x0 is None. With the same seed, objective, inputs and options, a run gives the same result every time on the same computer and installation; another computer, other versions of Python, NumPy, SciPy or gpyreg, or another number of threads for linear algebra can give a different result. On Apple Silicon Macs, two runs with the same seed can end at slightly different points. optimize_result["random_seed"] records the seed.

If you leave random_seed unset, np.random.seed(42) before creating the BADS object also fixes the run.

The seed does not govern your objective: if your objective is noisy, give it a random generator of its own, as explained above.

I have been running PyBADS with a stochastic objective function from different starting points and I get different results each time. What can I do?#

First of all, slightly different results are expected if your objective function is noisy (be sure to have read and understood all the points under the noisy objective function section of the FAQ). So, if you are asking this question, it is because you find wildly different results.

Generally, substantially different results suggest that PyBADS is getting stuck due to excess noise in the objective function with respect to actual improvements of the function in the neighborhood of the current point. For example, even if the expected value of the objective function would have a non-zero gradient, it might be too hard for PyBADS to find a direction of improvement due to low signal-to-noise ratio. For this reason, it is possible for PyBADS to get stuck at very different points in noisy, nearly-flat regions of the input space, with more scattered results for flatter and wider plateaus.

The general solution of this problem, as also mentioned in this question, is to decrease the amount of noise in the objective function (e.g., if you estimate the objective function via Monte Carlo sampling, try increasing the number of samples). While we generally found that a SD of the noise of 1 or less in the vicinity of the solution works for most problems, a particularly difficult (e.g., flat) objective function might need even lower amounts of noise for robust convergence, so YMMV.

On some problems, PyBADS seems to get stuck and stop too early. Is there a way to tune PyBADS to optimize towards a higher precision result or to have it optimize for longer?#

First, check why the run ended, in optimize_result["message"]. If it used up its budget, raise options["max_fun_evals"] (default 500 * D evaluations) or options["max_iter"] (default 200 * D iterations).

Otherwise, there are some options one can modify in PyBADS to make the algorithm poll or search for longer. This might help in some cases (especially for noisy objectives). If you want PyBADS to search for longer at each iteration, you can modify two key options in the options dictionary that you pass to the algorithm:

Set options["complete_poll"] = True (default is False). This will force PyBADS to finish the “poll” step (more info in the paper) instead of skipping it when it thinks that it is not worth continuing.
Change options["search_n_try"]. Be careful that this is an advanced option of PyBADS, and we do not recommend to change it unless you have to. The default value is max(D, floor(3 + D/2)), where D is the number of variables. This quantity represents the number of searches (via local Bayesian optimization) that PyBADS attempts in each round of searches before a poll. You can try and increase it to force PyBADS to search for longer in each iteration.

On some problems, PyBADS seems to find a reasonably good solution, but then it takes a long time to converge, spending many iterations at very small values of `MeshScale`. Is there a way to tune PyBADS to stop earlier once it finds a decent solution?#

This issue is relatively common with noisy objective functions. If it seems that PyBADS is spending too much time to converge at small mesh scales, e.g., dozens and dozens of iterations spent at MeshScale below 0.0001 or so, with little changes to the value of f(x) (or E[f(x)] for noisy objectives) consider modifying options["tol_mesh"].

This option determines a bunch of PyBADS behaviors, including the termination condition when the normalized MeshScale goes below a threshold. The default value is options["tol_mesh"] = 1e-6, which may be exceedingly conservative for stochastic problems. You could try larger values, such as options["tol_mesh"] = 1e-5 or even options["tol_mesh"] = 1e-4. However, use higher values at your own risk, in that you might be terminating a run which is still improving on the solution.

Miscellanea#

This is interesting, but shouldn’t we ideally compute full posterior distributions?#

Yes, we should! However, given the typical class of model-fitting problems PyBADS is designed for (see here), obtaining the full posterior, or even an approximation thereof, can be a challenging task. We developed a method and related toolbox, Variational Bayesian Monte Carlo (PyVBMC), which addresses exactly this problem and is one of the lab’s tools for fitting models to data. Check it out!

Can PyBADS return an approximate posterior, e.g. by computing the Hessian at the optimum?#

First, just to clarify, the Hessian is a matrix of second derivatives, which can be used to build a crude approximation of the posterior via Laplace’s method. The answer to the question is nope, PyBADS cannot return the Hessian or an approximate posterior. The reason is that even if PyBADS tries to build a local Gaussian process approximation of the objective function, this might fail and we cannot trust this approximation at all to represent a valid posterior.

Instead, you should look into Variational Bayesian Monte Carlo (PyVBMC), a method that we developed specifically to compute approximate posterior distributions, one of the lab’s tools for fitting models to data, and that can be used in synergy with PyBADS (see below).

I have run PyBADS on my problem. How do I run PyVBMC?#

PyVBMC, one of the lab’s tools for fitting models to data, computes an approximate posterior distribution over the parameters, and an estimate of the model evidence (see above). Its interface is very similar to the one of PyBADS, and you may only need minor changes to run it on your problem. Note that:

Beware of the sign! The target of PyVBMC is a log density, that is the log likelihood, or the log likelihood plus the log prior, whereas PyBADS minimizes the negative log likelihood.
PyVBMC needs a prior over the parameters, which you pass with prior= or add to the log likelihood.
The plausible bounds of PyVBMC should lie strictly inside its hard bounds, whereas PyBADS takes plausible bounds equal to the hard ones.
PyVBMC does not support variables bounded on one side only, nor fixed variables: give it a function of the other variables that inserts the fixed values (np.insert(x, i, value)), with the bounds of the other variables.
PyVBMC works best with up to about 10 parameters, where PyBADS takes up to about 20.
The solution of PyBADS is a good starting point x0 for PyVBMC.
PyVBMC supports noisy targets too, and works best when the target returns an estimate of its noise.

For example, with a uniform prior over the hard bounds:

from pyvbmc import VBMC
from scipy.stats import uniform


def log_likelihood(x):
    return -fun(x)  # fun is the negative log likelihood minimized by PyBADS


prior = [
    uniform(loc=low, scale=high - low)
    for low, high in zip(np.ravel(lb), np.ravel(ub))
]
vbmc = VBMC(
    log_likelihood, optimize_result["x"], lb, ub, plb, pub, prior=prior
)
vp, results = vbmc.optimize()

See the PyVBMC documentation for installation, priors and examples.

I used BADS in MATLAB. What is different in PyBADS?#

PyBADS implements the same algorithm, with a Python interface:

You create a BADS object and run it, optimize_result = BADS(fun, x0, lb, ub, plb, pub, non_box_cons, options).optimize(), where MATLAB returns [X,FVAL,EXITFLAG,OUTPUT] = bads(...). The result is one dictionary, with x, fval, status (MATLAB’s EXITFLAG), message, func_count and more (see above).
The options are a dictionary. Their names are mostly the MATLAB names in lower case, with underscores (MaxFunEvals becomes max_fun_evals, UncertaintyHandling becomes uncertainty_handling), and their values are Python values: numbers where MATLAB takes a string such as '500*nvars', True or False where MATLAB takes 1, 0, 'on' or 'off', and indices counted from 0 where MATLAB counts from 1 (periodic_vars is [2, 3] where MATLAB’s PeriodicVars is [3 4]). PyBADS raises a ValueError for an option name it does not know and for many values of the wrong kind, such as the string '200*D' for max_iter. The options page lists them all.
The objective receives a one-dimensional array of shape (D,), and non_box_cons an array of shape (N, D). Additional inputs of the objective are bound to it rather than passed to BADS.
A few options of MATLAB BADS, such as plot and restarts, are accepted but have no effect.
Function evaluations made before the run, which MATLAB BADS takes in its option FunValues, are passed to BADS as the argument precomputed_evaluations=(X, y), or (X, y, y_sd) with specify_target_noise. They enter the run’s log and its Gaussian process but not its count of evaluations, and the run starts from x0 and its initial design, as in MATLAB BADS.
Runs of PyBADS and of MATLAB BADS do not match step by step, even from the same starting point: the two draw their random numbers differently, and differ in some numerical details.

Are you planning to port BADS to other languages?#

BADS is currently available as a MATLAB toolbox and as a Python package, PyBADS, both listed with the lab’s other tools for fitting models to data. No other ports are currently planned, but please get in touch if interested. The Python package can also be called from other languages, for example from Julia with PythonCall.jl or from R with reticulate.