From bdae19dcd662a112700a769e5a0852daa7f56960 Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Tue, 6 Oct 2026 10:32:44 +0200
Subject: [PATCH 1/9] bench(hierarchical): profile where the time of a fit goes

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                  |   6 +
 benchmarks/hierarchical_profile.py            | 249 ++++++++++++++++++
 .../results/hierarchical-profile-before.csv   |  18 ++
 .../results/hierarchical-profile-cbm.csv      |   8 +
 4 files changed, 281 insertions(+)
 create mode 100644 benchmarks/hierarchical_profile.py
 create mode 100644 benchmarks/results/hierarchical-profile-before.csv
 create mode 100644 benchmarks/results/hierarchical-profile-cbm.csv

diff --git a/CHANGELOG.md b/CHANGELOG.md
index eb16efd8..a897058e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+
+- Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
+
 ## [0.26.0.dev1] - 2026-10-05
 
 ### Added
diff --git a/benchmarks/hierarchical_profile.py b/benchmarks/hierarchical_profile.py
new file mode 100644
index 00000000..e3939d4f
--- /dev/null
+++ b/benchmarks/hierarchical_profile.py
@@ -0,0 +1,249 @@
+"""
+Where does the time of a hierarchical fit go?
+
+    python benchmarks/hierarchical_profile.py [--data DATASET.csv] [--iterations 5] [--label NAME]
+
+Fits one dataset of the four-armed bandit (100 participants, 120 trials) with
+`RLRW` and `VariationalBayes` for a fixed number of iterations, one chain and
+two optimiser starts per participant, under `cProfile`. Prints the share of
+time in the optimiser, the Hessian, `evaluate_fit`, the population update and
+the parts of `objective`, and the number of objective evaluations per
+participant and iteration, separately for the optimiser, the Hessian and
+`evaluate_fit`. Writes the table to benchmarks/results/<label>.csv.
+
+`--data` is a dataset in the format `RLRW` reads (columns ppt, trial, arm_left,
+arm_right, reward_left, reward_right, response). Without it, the task and the
+responses are simulated from a known population. The settings match the speed
+comparison with cbm: bounds alpha [0.01, 1] and temperature [0, 10], a starting
+prior of mean 0.5 and sd 0.5 for alpha and mean 5 and sd 5 for temperature.
+
+The profiler adds overhead to Python-heavy code, so the shares are indicative;
+the fit is also timed once without it.
+"""
+
+import argparse
+import cProfile
+import pathlib
+import pstats
+import sys
+import time
+import warnings
+
+import numpy as np
+import pandas as pd
+from scipy.stats import truncnorm
+
+sys.path.insert(0, str(pathlib.Path(__file__).parent))
+
+import cpm  # noqa: E402
+import cpm.optimisation.fmin as fmin_module  # noqa: E402
+from cpm.applications.reinforcement_learning import RLRW  # noqa: E402
+from cpm.hierarchical import VariationalBayes  # noqa: E402
+from cpm.optimisation import FminBound, minimise  # noqa: E402
+
+RESULTS = pathlib.Path(__file__).parent / "results"
+PARTICIPANTS, TRIALS, STIMULI = 100, 120, 4
+BOUNDS = {"alpha": (1e-2, 1.0), "temperature": (0.0, 10.0)}
+STARTING_PRIOR = {"alpha": {"mean": 0.5, "sd": 0.5}, "temperature": {"mean": 5.0, "sd": 5.0}}
+POPULATION = {"alpha": {"mean": 0.5, "sd": 0.15}, "temperature": {"mean": 4.5, "sd": 1.5}}
+
+
+def simulate(seed=2026):
+    """A random four-armed bandit task with responses of `RLRW` from a known population."""
+    rng = np.random.default_rng(seed)
+    n, trials = PARTICIPANTS, TRIALS
+    arms = np.array([rng.choice(STIMULI, size=2, replace=False) for _ in range(n * trials)]).reshape(n, trials, 2)
+    reward_probability = np.array([0.2, 0.4, 0.6, 0.8])
+    rewards = (rng.random((n, trials, 2)) < reward_probability[arms]).astype(float)
+    parameters = {
+        name: truncnorm.rvs(
+            (BOUNDS[name][0] - p["mean"]) / p["sd"], (BOUNDS[name][1] - p["mean"]) / p["sd"],
+            loc=p["mean"], scale=p["sd"], size=n, random_state=rng,
+        )
+        for name, p in POPULATION.items()
+    }
+    values = np.full((n, STIMULI), 1 / STIMULI)
+    rows = np.arange(n)
+    response = np.empty((n, trials), dtype=int)
+    for t in range(trials):
+        scaled = values[rows[:, None], arms[:, t]] * parameters["temperature"][:, None]
+        scaled -= scaled.max(axis=1, keepdims=True)
+        policy = np.exp(scaled)
+        policy /= policy.sum(axis=1, keepdims=True)
+        choice = (rng.random(n) >= policy[:, 0]).astype(int)
+        stimulus = arms[rows, t, choice]
+        values[rows, stimulus] += parameters["alpha"] * (rewards[rows, t, choice] - values[rows, stimulus])
+        response[:, t] = choice
+    return pd.DataFrame(
+        {
+            "ppt": np.repeat(np.arange(1, n + 1), trials),
+            "trial": np.tile(np.arange(1, trials + 1), n),
+            "arm_left": arms[:, :, 0].ravel() + 1,
+            "arm_right": arms[:, :, 1].ravel() + 1,
+            "reward_left": rewards[:, :, 0].ravel(),
+            "reward_right": rewards[:, :, 1].ravel(),
+            "response": response.ravel(),
+            "observed": response.ravel(),
+        }
+    )
+
+
+def load(path):
+    data = pd.read_csv(path)
+    ppts = np.sort(data.ppt.unique())[:PARTICIPANTS]
+    data = data[data.ppt.isin(ppts) & (data.trial <= TRIALS)].reset_index(drop=True)
+    data["observed"] = data["response"]
+    return data
+
+
+def setup(data, iterations):
+    """`VariationalBayes` with the settings of the speed comparison, forced to run `iterations` iterations."""
+    first = data[data.ppt == data.ppt.iloc[0]]
+    model = RLRW(data=first, dimensions=STIMULI, parameters_settings=[[0.5, *BOUNDS["alpha"]], [2.0, *BOUNDS["temperature"]]])
+    model.parameters.update_prior(**STARTING_PRIOR)
+    optimiser = FminBound(
+        model=model,
+        data=data.groupby("ppt"),
+        minimisation=minimise.LogLikelihood.bernoulli,
+        prior=True,
+        number_of_starts=2,
+        ppt_identifier="ppt",
+        parallel=False,
+        display=False,
+        approx_grad=True,
+    )
+    return VariationalBayes(
+        optimiser=optimiser,
+        iteration=iterations,
+        chain=1,
+        tolerance_lme=0.0,
+        tolerance_param=0.0,
+        hyperpriors={"a0": np.array([p["mean"] for p in STARTING_PRIOR.values()]), "b": 1, "v": 0.5, "s": np.repeat(0.01, 2)},
+        quiet=True,
+    )
+
+
+class Counter:
+    """Counts the objective evaluations of `FminBound` by what asked for them."""
+
+    def __init__(self):
+        self.phase = "optimiser"
+        self.counts = {"optimiser": 0, "hessian": 0, "evaluate_fit": 0}
+        self.originals = (fmin_module.objective, fmin_module.numerical_hessian, fmin_module.evaluate_fit)
+
+    def __enter__(self):
+        objective, numerical_hessian, evaluate_fit = self.originals
+
+        def counted_objective(*args, **kwargs):
+            self.counts[self.phase] += 1
+            return objective(*args, **kwargs)
+
+        def counted_hessian(*args, **kwargs):
+            self.phase = "hessian"
+            try:
+                return numerical_hessian(*args, **kwargs)
+            finally:
+                self.phase = "optimiser"
+
+        def counted_evaluate_fit(*args, **kwargs):
+            self.counts["evaluate_fit"] += 1  # one model run per call
+            return evaluate_fit(*args, **kwargs)
+
+        fmin_module.objective = counted_objective
+        fmin_module.numerical_hessian = counted_hessian
+        fmin_module.evaluate_fit = counted_evaluate_fit
+        return self
+
+    def __exit__(self, *exc):
+        fmin_module.objective, fmin_module.numerical_hessian, fmin_module.evaluate_fit = self.originals
+
+
+def profile(vb):
+    profiler = cProfile.Profile()
+    profiler.enable()
+    vb.optimise()
+    profiler.disable()
+    stats = pstats.Stats(profiler)
+    stats.calc_callees()
+    return stats
+
+
+def find(stats, filename, function):
+    """The (file, line, name) key of a profiled function, by the basename of its file."""
+    for key in stats.stats:
+        if key[2] == function and pathlib.Path(key[0]).name == filename:
+            return key
+    return None
+
+
+def inclusive(stats, filename, function):
+    key = find(stats, filename, function)
+    return stats.stats[key][3] if key else 0.0
+
+
+def inclusive_from(stats, caller, filename, function):
+    """Seconds spent in `function` when called from `caller`."""
+    callees = stats.all_callees.get(caller, {})
+    return sum(value[3] for key, value in callees.items() if key[2] == function and pathlib.Path(key[0]).name == filename)
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    parser.add_argument("--data", help="a dataset CSV; simulated if not given")
+    parser.add_argument("--iterations", type=int, default=5)
+    parser.add_argument("--label", default="hierarchical-profile")
+    args = parser.parse_args()
+
+    warnings.simplefilter("ignore")
+    np.seterr(all="ignore")
+    data = load(args.data) if args.data else simulate()
+    participants, iterations = data.ppt.nunique(), args.iterations
+    print(f"cpm {cpm.__version__}: {participants} participants, {data.trial.max()} trials, {iterations} iterations, 2 starts")
+
+    # the compiled model and everything else that loads on the first call
+    setup(data[data.ppt.isin(data.ppt.unique()[:5])], 1).optimise()
+
+    # once without the profiler, for the time
+    np.random.seed(2026)
+    start = time.perf_counter()
+    with Counter() as counter:
+        setup(data, iterations).optimise()
+    seconds = time.perf_counter() - start
+    per = participants * iterations
+    counts = {k: v / per for k, v in counter.counts.items()}
+
+    np.random.seed(2026)
+    stats = profile(setup(data, iterations))
+    total = inclusive(stats, "variational.py", "optimise")
+    objective_key = find(stats, "optimisers.py", "objective")
+    rows = [
+        ("fit", "seconds without profiler", seconds),
+        ("fit", "ms per participant and iteration", seconds * 1e3 / per),
+        ("evaluations per participant and iteration", "optimiser", counts["optimiser"]),
+        ("evaluations per participant and iteration", "hessian", counts["hessian"]),
+        ("evaluations per participant and iteration", "evaluate_fit", counts["evaluate_fit"]),
+        ("share of profiled time", "fmin_l_bfgs_b", inclusive(stats, "_lbfgsb_py.py", "fmin_l_bfgs_b") / total),
+        ("share of profiled time", "numerical_hessian", inclusive(stats, "optimisers.py", "numerical_hessian") / total),
+        ("share of profiled time", "evaluate_fit", inclusive(stats, "optimisers.py", "evaluate_fit") / total),
+        ("share of profiled time", "update_population", inclusive(stats, "variational.py", "update_population") / total),
+        ("share of profiled time", "get_lme", inclusive(stats, "variational.py", "get_lme") / total),
+        ("share of profiled time", "objective (all callers)", inclusive(stats, "optimisers.py", "objective") / total),
+        ("share of objective", "SessionWrapper.run", inclusive_from(stats, objective_key, "session.py", "run")),
+        ("share of objective", "  of which the RLRW model function", inclusive(stats, "reinforcement_learning.py", "__call__")),
+        ("share of objective", "Wrapper.reset", inclusive_from(stats, objective_key, "session.py", "reset")),
+        ("share of objective", "deepcopy", inclusive_from(stats, objective_key, "copy.py", "deepcopy")),
+        ("share of objective", "loss", inclusive_from(stats, objective_key, "minimise.py", "bernoulli")),
+        ("share of objective", "Parameters.PDF", inclusive_from(stats, objective_key, "parameters.py", "PDF")),
+    ]
+    objective_seconds = inclusive(stats, "optimisers.py", "objective")
+    rows = [(s, i, v / objective_seconds if s == "share of objective" else v) for s, i, v in rows]
+    table = pd.DataFrame(rows, columns=["section", "item", "value"])
+    table["value"] = table["value"].round(4)
+    print(table.to_string(index=False))
+
+    RESULTS.mkdir(exist_ok=True)
+    table.to_csv(RESULTS / f"{args.label}.csv", index=False)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/benchmarks/results/hierarchical-profile-before.csv b/benchmarks/results/hierarchical-profile-before.csv
new file mode 100644
index 00000000..a9c573dc
--- /dev/null
+++ b/benchmarks/results/hierarchical-profile-before.csv
@@ -0,0 +1,18 @@
+section,item,value
+fit,seconds without profiler,8.6299
+fit,ms per participant and iteration,17.2597
+evaluations per participant and iteration,optimiser,70.59
+evaluations per participant and iteration,hessian,22.0
+evaluations per participant and iteration,evaluate_fit,2.0
+share of profiled time,fmin_l_bfgs_b,0.7824
+share of profiled time,numerical_hessian,0.1838
+share of profiled time,evaluate_fit,0.0143
+share of profiled time,update_population,0.0006
+share of profiled time,get_lme,0.0002
+share of profiled time,objective (all callers),0.6966
+share of objective,SessionWrapper.run,0.3991
+share of objective,  of which the RLRW model function,0.3183
+share of objective,Wrapper.reset,0.1337
+share of objective,deepcopy,0.0428
+share of objective,loss,0.2729
+share of objective,Parameters.PDF,0.1246
diff --git a/benchmarks/results/hierarchical-profile-cbm.csv b/benchmarks/results/hierarchical-profile-cbm.csv
new file mode 100644
index 00000000..32a5197d
--- /dev/null
+++ b/benchmarks/results/hierarchical-profile-cbm.csv
@@ -0,0 +1,8 @@
+item,value
+participants,100
+iterations,9
+lap_calls_per_participant,85.59
+hbi_calls_per_participant_and_iteration,71.74
+seconds_instrumented,7.21
+lap_seconds_instrumented,2.75
+ms_per_participant_and_iteration_instrumented,4.956

From 3ac923601fbc43b7705879f6be39ef415c21a650 Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Tue, 6 Oct 2026 10:45:52 +0200
Subject: [PATCH 2/9] perf(optimisation): compute FminBound's gradient with the
 objective

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                 |   4 +
 cpm/optimisation/fmin.py                     |  97 +++++++++++++-
 cpm/optimisation/minimise.py                 |   4 +-
 test/optimisation/test_fminbound_gradient.py | 131 +++++++++++++++++++
 test/optimisation/test_minimise.py           |  10 ++
 5 files changed, 238 insertions(+), 8 deletions(-)
 create mode 100644 test/optimisation/test_fminbound_gradient.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a897058e..2a582210 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
 
+### Changed
+
+- `cpm.optimisation.FminBound` with `approx_grad=True` now computes the finite-difference gradient itself, together with the objective, exactly as SciPy does, with identical results; with `LogLikelihood.bernoulli`, this makes the hierarchical fits of the speed comparison 19% faster
+
 ## [0.26.0.dev1] - 2026-10-05
 
 ### Added
diff --git a/cpm/optimisation/fmin.py b/cpm/optimisation/fmin.py
index 6f445d9d..d94cd0f6 100644
--- a/cpm/optimisation/fmin.py
+++ b/cpm/optimisation/fmin.py
@@ -27,6 +27,61 @@
 LBFGSB_VERBOSITY_SUPPORTED = "disp" in inspect.signature(fmin_l_bfgs_b).parameters
 
 
+def gradient_by_forward_differences(fun, x, f0, lower, upper, epsilon):
+    """
+    The gradient of `fun` at `x` by forward differences, as SciPy's L-BFGS-B estimates it with `approx_grad=True`.
+
+    Parameters
+    ----------
+    fun : callable
+        The function, of a 1-D array.
+    x : numpy.ndarray
+        The point, within the bounds.
+    f0 : float
+        `fun(x)`.
+    lower, upper : numpy.ndarray
+        The bounds of each element of `x`, possibly infinite.
+    epsilon : float
+        The step.
+
+    Returns
+    -------
+    numpy.ndarray
+        The gradient.
+
+    Notes
+    -----
+    This follows `scipy.optimize._numdiff.approx_derivative` with the "2-point"
+    method and an absolute step: a step that would leave the bounds is taken in
+    the other direction, and each difference is divided by the step as it is
+    represented, ``(x + h) - x``. The gradient is therefore exactly the one SciPy
+    computes, so the optimiser takes exactly the same path, but without SciPy's
+    overhead on every evaluation. It works on Python floats, which follow the
+    same IEEE arithmetic as NumPy's and are faster for a few parameters.
+    """
+    gradient = np.empty(x.size)
+    for i, (value, low, high) in enumerate(zip(x.tolist(), lower.tolist(), upper.tolist())):
+        h = epsilon
+        ## a step too small to change x falls back to SciPy's default relative step
+        if (value + h) - value == 0:
+            h = _SQRT_EPS * (1.0 if value >= 0 else -1.0) * max(1.0, abs(value))
+        lower_distance, upper_distance = value - low, high - value
+        if abs(h) <= max(lower_distance, upper_distance):
+            if value + h < low or value + h > high:
+                h = -h
+        elif upper_distance >= lower_distance:
+            h = upper_distance
+        else:
+            h = -lower_distance
+        step = x.copy()
+        step[i] = value + h
+        gradient[i] = (fun(step) - f0) / ((value + h) - value)
+    return gradient
+
+
+_SQRT_EPS = float(np.finfo(float).eps ** 0.5)
+
+
 def lbfgsb_options(display=False, **kwargs):
     """
     Assemble the keyword arguments for `scipy.optimize.fmin_l_bfgs_b`, keeping the solver's verbosity options only where the installed SciPy still accepts them.
@@ -483,6 +538,17 @@ def __unpack(x, id=None):
         model = self.model
         prior = self.prior
         options = lbfgsb_options(display=self.display, **self.kwargs)
+        ## with `approx_grad`, the gradient is computed here, exactly as SciPy
+        ## would, together with the objective (see `gradient_by_forward_differences`)
+        approx_grad = options.pop("approx_grad", False)
+        if approx_grad:
+            options.pop("fprime", None)
+            epsilon = options.pop("epsilon", 1e-8)
+            lower, upper = np.asarray(bounds, dtype=float).T
+            evaluations_per_point = len(self.parameter_names) + 1
+            ## SciPy counts the evaluations of the gradient towards `maxfun`, but
+            ## here it only counts the points, so the limit is in points
+            options["maxfun"] = options.get("maxfun", 15000) // evaluations_per_point
 
         def __task(participant, **args):
 
@@ -494,13 +560,30 @@ def __task(participant, **args):
 
             model.reset(data=participant_dc)
 
-            result = fmin_l_bfgs_b(
-                objective,
-                x0=self.__current_guess__,
-                bounds=bounds,
-                args=(model, observed, loss, prior),
-                **options,
-            )
+            if approx_grad:
+
+                def f(x):
+                    return objective(x, model, observed, loss, prior)
+
+                def objective_and_gradient(x):
+                    f0 = f(x)
+                    return f0, gradient_by_forward_differences(f, x, f0, lower, upper, epsilon)
+
+                result = fmin_l_bfgs_b(
+                    objective_and_gradient,
+                    x0=self.__current_guess__,
+                    bounds=bounds,
+                    **options,
+                )
+                result[2]["funcalls"] *= evaluations_per_point
+            else:
+                result = fmin_l_bfgs_b(
+                    objective,
+                    x0=self.__current_guess__,
+                    bounds=bounds,
+                    args=(model, observed, loss, prior),
+                    **options,
+                )
 
             def f(x):
                 return objective(x, model, observed, loss, prior)
diff --git a/cpm/optimisation/minimise.py b/cpm/optimisation/minimise.py
index a979f7f0..ab4babbe 100644
--- a/cpm/optimisation/minimise.py
+++ b/cpm/optimisation/minimise.py
@@ -198,7 +198,9 @@ def bernoulli(predicted=None, observed=None, negative=True, **kwargs):
         np.clip(probabilities, 1e-10, 1 - 1e-10, out=probabilities)
 
         LL = _bernoulli_logpmf(observed.flatten(), probabilities)
-        LL = check_nan_bounds_in_log(LL, bound=bound)
+        ## LL is a new array, which only needs replacing where it is not finite
+        if not np.isfinite(LL).all():
+            LL = check_nan_bounds_in_log(LL, bound=bound)
         LL = np.sum(LL)
         if negative:
             LL = -1 * LL
diff --git a/test/optimisation/test_fminbound_gradient.py b/test/optimisation/test_fminbound_gradient.py
new file mode 100644
index 00000000..845cfac0
--- /dev/null
+++ b/test/optimisation/test_fminbound_gradient.py
@@ -0,0 +1,131 @@
+"""
+`FminBound` with `approx_grad=True` computes the finite-difference gradient
+itself, together with the objective, instead of leaving it to SciPy. These tests
+check that the fits are exactly the ones SciPy's own `approx_grad` gives.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+from scipy.optimize import fmin_l_bfgs_b
+
+from cpm.applications.reinforcement_learning import RLRW
+from cpm.core.data import decompose
+from cpm.core.optimisers import objective
+from cpm.optimisation import FminBound, minimise
+
+LOSS = minimise.LogLikelihood.bernoulli
+
+
+def bandit(participants=4, trials=40, seed=1):
+    rng = np.random.default_rng(seed)
+    arms = np.array([rng.choice(4, size=2, replace=False) for _ in range(participants * trials)]) + 1
+    return pd.DataFrame(
+        {
+            "ppt": np.repeat(np.arange(1, participants + 1), trials),
+            "arm_left": arms[:, 0],
+            "arm_right": arms[:, 1],
+            "reward_left": rng.integers(0, 2, participants * trials).astype(float),
+            "reward_right": rng.integers(0, 2, participants * trials).astype(float),
+            "response": rng.integers(0, 2, participants * trials),
+            "observed": rng.integers(0, 2, participants * trials),
+        }
+    )
+
+
+def rlrw(data, temperature_upper=10.0):
+    model = RLRW(
+        data=data[data.ppt == 1],
+        dimensions=4,
+        parameters_settings=[[0.5, 1e-2, 1.0], [2.0, 0.0, temperature_upper]],
+    )
+    model.parameters.update_prior(alpha={"mean": 0.5, "sd": 0.5}, temperature={"mean": 5.0, "sd": 5.0})
+    return model
+
+
+def scipy_reference(model, data, guess, prior, **options):
+    """What `fmin_l_bfgs_b` with SciPy's own `approx_grad` gives for every participant."""
+    bounds = list(map(tuple, np.asarray(model.parameters.bounds()).T))
+    out = []
+    for participant in data.groupby("ppt"):
+        participant_dc, observed, _ = decompose(participant=participant, pandas=True, identifier="ppt")
+        model.reset(data=participant_dc)
+        out.append(
+            fmin_l_bfgs_b(objective, x0=guess, bounds=bounds, args=(model, observed, LOSS, prior), approx_grad=True, **options)
+        )
+    return out
+
+
+@pytest.mark.parametrize("prior", [False, True])
+@pytest.mark.parametrize("temperature_upper", [10.0, 0.5])
+@pytest.mark.parametrize("options", [{}, {"epsilon": 1e-6}, {"maxfun": 10}])
+def test_fits_are_identical_to_scipy_approx_grad(prior, temperature_upper, options):
+    data = bandit()
+    guess = np.array([0.3, 0.4])
+    fit = FminBound(
+        model=rlrw(data, temperature_upper),
+        data=data.groupby("ppt"),
+        minimisation=LOSS,
+        prior=prior,
+        initial_guess=guess,
+        number_of_starts=1,
+        approx_grad=True,
+        **options,
+    )
+    fit.optimise()
+    reference = scipy_reference(rlrw(data, temperature_upper), data, guess, prior, **options)
+    for got, (x, f, info) in zip(fit.fit, reference):
+        np.testing.assert_array_equal(got["x"], x)
+        assert got["fun"] == f
+        np.testing.assert_array_equal(got["grad"], info["grad"])
+        assert got["funcalls"] == info["funcalls"]
+        assert got["nit"] == info["nit"]
+        assert got["warnflag"] == info["warnflag"]
+        assert got["task"] == info["task"]
+
+
+def test_a_participant_on_the_upper_bound_is_covered():
+    """With an upper bound of 0.5 on the temperature, the estimates end on it, so the step must flip sign."""
+    data = bandit()
+    reference = scipy_reference(rlrw(data, 0.5), data, np.array([0.3, 0.4]), True)
+    assert any(x[1] == 0.5 for x, _, _ in reference)
+
+
+@pytest.mark.parametrize(
+    "x, lower, upper",
+    [
+        ([0.4, 3.0], [1e-2, 0.0], [1.0, 10.0]),  # inside
+        ([1e-2, 10.0], [1e-2, 0.0], [1.0, 10.0]),  # on a lower and an upper bound
+        ([0.4, -3.0], [-np.inf, -np.inf], [np.inf, np.inf]),  # unbounded
+        ([0.5, 1.0], [0.5, 1.0 - 4e-9], [0.5 + 3e-9, 1.0]),  # bounds closer than the step
+        ([1e9, -1e9], [-np.inf, -2e9], [2e9, np.inf]),  # a step too small to change x
+    ],
+)
+def test_gradient_is_scipys(x, lower, upper):
+    from scipy.optimize._numdiff import approx_derivative
+
+    from cpm.optimisation.fmin import gradient_by_forward_differences
+
+    def f(z):
+        return float(np.sin(z[0]) * z[1] + z[1] ** 2 / 3)
+
+    x, lower, upper = (np.asarray(v, dtype=float) for v in (x, lower, upper))
+    expected = approx_derivative(f, x, method="2-point", abs_step=1e-8, f0=f(x), bounds=(lower, upper))
+    np.testing.assert_array_equal(gradient_by_forward_differences(f, x, f(x), lower, upper, 1e-8), expected)
+
+
+def test_gradient_steps_away_from_a_bound():
+    from cpm.optimisation.fmin import gradient_by_forward_differences
+
+    calls = []
+
+    def f(x):
+        calls.append(x.copy())
+        return float(np.sum(x**2))
+
+    lower, upper = np.array([0.0, 0.0]), np.array([1.0, 1.0])
+    x = np.array([1.0, 0.5])
+    gradient = gradient_by_forward_differences(f, x, f(x), lower, upper, 1e-8)
+    assert all(np.all((c >= lower) & (c <= upper)) for c in calls)
+    assert calls[1][0] < 1.0  # stepped down from the upper bound
+    np.testing.assert_allclose(gradient, 2 * x, rtol=1e-6)
diff --git a/test/optimisation/test_minimise.py b/test/optimisation/test_minimise.py
index d7178b04..f8ea0f8b 100644
--- a/test/optimisation/test_minimise.py
+++ b/test/optimisation/test_minimise.py
@@ -370,6 +370,16 @@ def test_bernoulli_and_continuous_equal_scipy():
     assert LogLikelihood.continuous(predicted=x, observed=y) == -np.sum(norm.logpdf(x, y, 1))
 
 
+def test_bernoulli_equals_scipy_when_every_term_is_finite():
+    from scipy.stats import bernoulli
+
+    rng = np.random.default_rng(6)
+    p = rng.uniform(0, 1, (120, 1))
+    k = rng.integers(0, 2, 120)
+    expected = -np.sum(np.nan_to_num(bernoulli.logpmf(k, np.clip(p.ravel(), 1e-10, 1 - 1e-10)), neginf=-1e100))
+    assert LogLikelihood.bernoulli(predicted=p, observed=k) == expected
+
+
 def test_losses_still_name_nan_and_inf():
     with pytest.raises(ValueError, match="Predicted values contain NaN"):
         LogLikelihood.bernoulli(predicted=np.array([0.5, np.nan]), observed=np.array([1, 0]))

From ba4820e456306f9fdc0408965f6f4c378b68330e Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Tue, 6 Oct 2026 11:06:49 +0200
Subject: [PATCH 3/9] feat(optimisation): add
 FminBound(hessian="finite_differences")

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                |   1 +
 benchmarks/hierarchical_profile.py          |  44 +++--
 cpm/core/optimisers.py                      |  57 ++++++
 cpm/optimisation/fmin.py                    | 100 +++++++++--
 test/optimisation/test_fminbound_hessian.py | 184 ++++++++++++++++++++
 5 files changed, 350 insertions(+), 36 deletions(-)
 create mode 100644 test/optimisation/test_fminbound_hessian.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2a582210..e7515df5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
+- Added a `hessian` argument to `cpm.optimisation.FminBound`. `hessian="finite_differences"` computes the Hessian at the optimum, with finite differences that stay within the bounds (`cpm.core.optimisers.finite_difference_hessian`), once per participant for the best start, together with `log_likelihood`, `log_prior` and `metrics`: 7 instead of 22 evaluations per participant and iteration in the speed comparison, and no Hessian of zeros for an estimate on an upper bound. The default, `"numdifftools"`, keeps the previous Hessian
 
 ### Changed
 
diff --git a/benchmarks/hierarchical_profile.py b/benchmarks/hierarchical_profile.py
index e3939d4f..bef90780 100644
--- a/benchmarks/hierarchical_profile.py
+++ b/benchmarks/hierarchical_profile.py
@@ -96,8 +96,8 @@ def load(path):
     return data
 
 
-def setup(data, iterations):
-    """`VariationalBayes` with the settings of the speed comparison, forced to run `iterations` iterations."""
+def setup(data, iterations, **options):
+    """`VariationalBayes` with the settings of the speed comparison, forced to run `iterations` iterations; `options` go to `FminBound`."""
     first = data[data.ppt == data.ppt.iloc[0]]
     model = RLRW(data=first, dimensions=STIMULI, parameters_settings=[[0.5, *BOUNDS["alpha"]], [2.0, *BOUNDS["temperature"]]])
     model.parameters.update_prior(**STARTING_PRIOR)
@@ -111,6 +111,7 @@ def setup(data, iterations):
         parallel=False,
         display=False,
         approx_grad=True,
+        **options,
     )
     return VariationalBayes(
         optimiser=optimiser,
@@ -129,33 +130,39 @@ class Counter:
     def __init__(self):
         self.phase = "optimiser"
         self.counts = {"optimiser": 0, "hessian": 0, "evaluate_fit": 0}
-        self.originals = (fmin_module.objective, fmin_module.numerical_hessian, fmin_module.evaluate_fit)
+        self.originals = (fmin_module.objective, fmin_module.numerical_hessian, fmin_module.finite_difference_hessian,
+                          fmin_module.evaluate_fit)
 
     def __enter__(self):
-        objective, numerical_hessian, evaluate_fit = self.originals
+        objective, numerical_hessian, finite_difference_hessian, evaluate_fit = self.originals
 
         def counted_objective(*args, **kwargs):
             self.counts[self.phase] += 1
             return objective(*args, **kwargs)
 
-        def counted_hessian(*args, **kwargs):
-            self.phase = "hessian"
-            try:
-                return numerical_hessian(*args, **kwargs)
-            finally:
-                self.phase = "optimiser"
+        def counted(hessian):
+            def counted_hessian(*args, **kwargs):
+                self.phase = "hessian"
+                try:
+                    return hessian(*args, **kwargs)
+                finally:
+                    self.phase = "optimiser"
+
+            return counted_hessian
 
         def counted_evaluate_fit(*args, **kwargs):
             self.counts["evaluate_fit"] += 1  # one model run per call
             return evaluate_fit(*args, **kwargs)
 
         fmin_module.objective = counted_objective
-        fmin_module.numerical_hessian = counted_hessian
+        fmin_module.numerical_hessian = counted(numerical_hessian)
+        fmin_module.finite_difference_hessian = counted(finite_difference_hessian)
         fmin_module.evaluate_fit = counted_evaluate_fit
         return self
 
     def __exit__(self, *exc):
-        fmin_module.objective, fmin_module.numerical_hessian, fmin_module.evaluate_fit = self.originals
+        (fmin_module.objective, fmin_module.numerical_hessian, fmin_module.finite_difference_hessian,
+         fmin_module.evaluate_fit) = self.originals
 
 
 def profile(vb):
@@ -192,28 +199,30 @@ def main():
     parser.add_argument("--data", help="a dataset CSV; simulated if not given")
     parser.add_argument("--iterations", type=int, default=5)
     parser.add_argument("--label", default="hierarchical-profile")
+    parser.add_argument("--hessian", default="numdifftools", help="the `hessian` option of FminBound")
     args = parser.parse_args()
 
     warnings.simplefilter("ignore")
     np.seterr(all="ignore")
     data = load(args.data) if args.data else simulate()
     participants, iterations = data.ppt.nunique(), args.iterations
-    print(f"cpm {cpm.__version__}: {participants} participants, {data.trial.max()} trials, {iterations} iterations, 2 starts")
+    options = {"hessian": args.hessian}
+    print(f"cpm {cpm.__version__}: {participants} participants, {data.trial.max()} trials, {iterations} iterations, 2 starts, {options}")
 
     # the compiled model and everything else that loads on the first call
-    setup(data[data.ppt.isin(data.ppt.unique()[:5])], 1).optimise()
+    setup(data[data.ppt.isin(data.ppt.unique()[:5])], 1, **options).optimise()
 
     # once without the profiler, for the time
     np.random.seed(2026)
     start = time.perf_counter()
     with Counter() as counter:
-        setup(data, iterations).optimise()
+        setup(data, iterations, **options).optimise()
     seconds = time.perf_counter() - start
     per = participants * iterations
     counts = {k: v / per for k, v in counter.counts.items()}
 
     np.random.seed(2026)
-    stats = profile(setup(data, iterations))
+    stats = profile(setup(data, iterations, **options))
     total = inclusive(stats, "variational.py", "optimise")
     objective_key = find(stats, "optimisers.py", "objective")
     rows = [
@@ -223,7 +232,8 @@ def main():
         ("evaluations per participant and iteration", "hessian", counts["hessian"]),
         ("evaluations per participant and iteration", "evaluate_fit", counts["evaluate_fit"]),
         ("share of profiled time", "fmin_l_bfgs_b", inclusive(stats, "_lbfgsb_py.py", "fmin_l_bfgs_b") / total),
-        ("share of profiled time", "numerical_hessian", inclusive(stats, "optimisers.py", "numerical_hessian") / total),
+        ("share of profiled time", "hessian", (inclusive(stats, "optimisers.py", "numerical_hessian")
+                                               + inclusive(stats, "optimisers.py", "finite_difference_hessian")) / total),
         ("share of profiled time", "evaluate_fit", inclusive(stats, "optimisers.py", "evaluate_fit") / total),
         ("share of profiled time", "update_population", inclusive(stats, "variational.py", "update_population") / total),
         ("share of profiled time", "get_lme", inclusive(stats, "variational.py", "get_lme") / total),
diff --git a/cpm/core/optimisers.py b/cpm/core/optimisers.py
index 64beefea..cc2f4e3b 100644
--- a/cpm/core/optimisers.py
+++ b/cpm/core/optimisers.py
@@ -12,6 +12,7 @@
     "group_identifier",
     "prepare_data",
     "numerical_hessian",
+    "finite_difference_hessian",
 ]
 
 
@@ -34,6 +35,62 @@ def numerical_hessian(func=None, params=None, hessian=None):
     return computed_hessian
 
 
+def finite_difference_hessian(func, x, lower, upper, step=1e-4):
+    """
+    The Hessian of `func` at `x` by finite differences that stay within the bounds.
+
+    Parameters
+    ----------
+    func : callable
+        The function, of a 1-D array.
+    x : numpy.ndarray
+        The point, within the bounds, normally the optimum.
+    lower, upper : numpy.ndarray
+        The bounds of each element of `x`, possibly infinite.
+    step : float, optional
+        The step, relative to ``max(1, |x|)``. Default is 1e-4.
+
+    Returns
+    -------
+    numpy.ndarray
+        The Hessian matrix.
+
+    Notes
+    -----
+    If every parameter is at least one step from its bounds, the Hessian is
+    computed with second-order differences (Abramowitz and Stegun, 1972,
+    equations 25.3.23 and 25.3.27), which take d² + d + 1 evaluations for d
+    parameters. Otherwise, forward differences step away from the bounds, which
+    take d(d + 1) / 2 + d + 1 evaluations and are first-order accurate.
+    """
+    x = np.asarray(x, dtype=float)
+    d = x.size
+    h = step * np.maximum(1.0, np.abs(x))
+    f0 = func(x)
+    hessian = np.empty((d, d))
+    if np.all((x - h >= lower) & (x + h <= upper)):
+        steps = np.diag(h)
+        plus = [func(x + steps[i]) for i in range(d)]
+        minus = [func(x - steps[i]) for i in range(d)]
+        for i in range(d):
+            hessian[i, i] = (plus[i] - 2 * f0 + minus[i]) / h[i] ** 2
+            for j in range(i + 1, d):
+                both = func(x + steps[i] + steps[j]) + func(x - steps[i] - steps[j])
+                hessian[i, j] = hessian[j, i] = (
+                    both - plus[i] - plus[j] + 2 * f0 - minus[i] - minus[j]
+                ) / (2 * h[i] * h[j])
+        return hessian
+    h = np.where(x + 2 * h > upper, -h, h)
+    steps = np.diag(h)
+    single = [func(x + steps[i]) for i in range(d)]
+    for i in range(d):
+        for j in range(i, d):
+            hessian[i, j] = hessian[j, i] = (
+                func(x + steps[i] + steps[j]) - single[i] - single[j] + f0
+            ) / (h[i] * h[j])
+    return hessian
+
+
 def objective(pars, function, data, loss, prior=False):
     """
     The `objective` function calculates a metric by comparing predicted values with
diff --git a/cpm/optimisation/fmin.py b/cpm/optimisation/fmin.py
index d94cd0f6..e493a88a 100644
--- a/cpm/optimisation/fmin.py
+++ b/cpm/optimisation/fmin.py
@@ -4,6 +4,7 @@
     evaluate_fit,
     fit_extras,
     numerical_hessian,
+    finite_difference_hessian,
     group_identifier,
     prepare_data,
 )
@@ -425,6 +426,8 @@ class FminBound:
         The libraries to import for parallel processing for `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`.
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `fit` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start; at an upper bound, that point is outside the bounds, which can give a Hessian of zeros. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the :func:`scipy.optimize.fmin_l_bfgs_b` documentation for what is supported. The solver's own `disp` and `iprint` options, which `display` also sets, were removed in SciPy 1.18.0 and are ignored there.
 
@@ -463,6 +466,7 @@ def __init__(
         metrics=None,
         ppt_identifier=None,
         display=False,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -477,6 +481,11 @@ def __init__(
         self.metrics = metrics
         self.kwargs = kwargs
         self.display = display
+        if hessian not in ("numdifftools", "finite_differences"):
+            raise ValueError(
+                f'hessian must be "numdifftools" or "finite_differences", not {hessian!r}.'
+            )
+        self.hessian = hessian
 
         self.fit = []
         self.details = []
@@ -538,13 +547,13 @@ def __unpack(x, id=None):
         model = self.model
         prior = self.prior
         options = lbfgsb_options(display=self.display, **self.kwargs)
+        lower, upper = np.asarray(bounds, dtype=float).T
         ## with `approx_grad`, the gradient is computed here, exactly as SciPy
         ## would, together with the objective (see `gradient_by_forward_differences`)
         approx_grad = options.pop("approx_grad", False)
         if approx_grad:
             options.pop("fprime", None)
             epsilon = options.pop("epsilon", 1e-8)
-            lower, upper = np.asarray(bounds, dtype=float).T
             evaluations_per_point = len(self.parameter_names) + 1
             ## SciPy counts the evaluations of the gradient towards `maxfun`, but
             ## here it only counts the points, so the limit is in points
@@ -560,10 +569,10 @@ def __task(participant, **args):
 
             model.reset(data=participant_dc)
 
-            if approx_grad:
+            def f(x):
+                return objective(x, model, observed, loss, prior)
 
-                def f(x):
-                    return objective(x, model, observed, loss, prior)
+            if approx_grad:
 
                 def objective_and_gradient(x):
                     f0 = f(x)
@@ -585,12 +594,14 @@ def objective_and_gradient(x):
                     **options,
                 )
 
-            def f(x):
-                return objective(x, model, observed, loss, prior)
+            result = (*result[0:2], *tuple(list(result[2].values())))
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                return (*result, ppt)
 
             hessian = numerical_hessian(func=f, params=result[0] + 1e-3)
 
-            result = (*result[0:2], *tuple(list(result[2].values())), hessian)
+            result = (*result, hessian)
 
             if prior or self.metrics:
                 result = (
@@ -603,12 +614,49 @@ def f(x):
             result = (*result, ppt)
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = (*result[:-1], finite_difference_hessian(f, result[0], lower, upper))
+            if prior or self.metrics:
+                result = (
+                    *result,
+                    *evaluate_fit(
+                        result[0], model, observed, loss, prior, self.metrics
+                    ).values(),
+                )
+            return (*result, ppt)
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i][1]
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            parameters = {}
+            for result in results:
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result[0][i]
+                    )
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -628,19 +676,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                self.details = copy.deepcopy(results)
-                parameters = {}
-                for result in results:
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result[0][i]
-                        )
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -655,6 +695,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    method=None,
+                    cl=self.cl,
+                    pandas=False,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
diff --git a/test/optimisation/test_fminbound_hessian.py b/test/optimisation/test_fminbound_hessian.py
new file mode 100644
index 00000000..851d60bd
--- /dev/null
+++ b/test/optimisation/test_fminbound_hessian.py
@@ -0,0 +1,184 @@
+"""
+`FminBound(hessian="finite_differences")`: the Hessian at the optimum by finite
+differences that stay within the bounds, computed once per participant.
+"""
+
+import numdifftools as nd
+import numpy as np
+import pytest
+
+import cpm.optimisation.fmin as fmin_module
+from cpm.core.optimisers import numerical_hessian, objective
+from cpm.optimisation import FminBound, minimise
+from test.optimisation.test_fminbound_gradient import bandit, rlrw
+
+LOSS = minimise.LogLikelihood.bernoulli
+A = np.array([[4.0, 1.5, 0.2], [1.5, 0.7, -0.1], [0.2, -0.1, 2.0]])
+
+
+def counted(f):
+    calls = []
+
+    def g(x):
+        calls.append(np.array(x, copy=True))
+        return f(x)
+
+    return g, calls
+
+
+def quadratic(x):
+    return float(0.5 * x @ A @ x + np.sin(1.0))
+
+
+@pytest.mark.parametrize("d", [2, 3])
+def test_interior_hessian_of_a_quadratic(d):
+    from cpm.core.optimisers import finite_difference_hessian
+
+    f, calls = counted(lambda x: quadratic(np.pad(x, (0, 3 - d))))
+    x = np.array([0.3, 2.0, -1.0])[:d]
+    H = finite_difference_hessian(f, x, np.full(d, -10.0), np.full(d, 10.0))
+    np.testing.assert_allclose(H, A[:d, :d], rtol=1e-4, atol=1e-6)
+    assert len(calls) == d * d + d + 1
+
+
+def test_hessian_on_an_upper_bound_stays_within_the_bounds():
+    from cpm.core.optimisers import finite_difference_hessian
+
+    f, calls = counted(quadratic)
+    lower, upper = np.array([-10.0, 0.0, -10.0]), np.array([10.0, 2.0, 10.0])
+    x = np.array([0.3, 2.0, -1.0])
+    H = finite_difference_hessian(f, x, lower, upper)
+    assert all(np.all((c >= lower) & (c <= upper)) for c in calls)
+    assert len(calls) == 3 * 4 // 2 + 3 + 1
+    np.testing.assert_allclose(H, A, rtol=1e-4, atol=1e-6)
+
+
+def learners(participants=6, trials=120, alpha=0.3, temperature=4.0, seed=2):
+    """A bandit with rewarding arms and the choices of RLRW, so that the estimates are inside the bounds."""
+    rng = np.random.default_rng(seed)
+    data = bandit(participants, trials, seed)
+    reward_probability = np.array([0.2, 0.4, 0.6, 0.8])
+    for side in ("left", "right"):
+        data[f"reward_{side}"] = (rng.random(len(data)) < reward_probability[data[f"arm_{side}"] - 1]).astype(float)
+    choices = []
+    for _, g in data.groupby("ppt"):
+        values = np.full(4, 0.25)
+        for left, right, reward_left, reward_right in g[["arm_left", "arm_right", "reward_left", "reward_right"]].to_numpy():
+            left, right = int(left) - 1, int(right) - 1
+            p_right = 1 / (1 + np.exp(-temperature * (values[right] - values[left])))
+            choice = int(rng.random() < p_right)
+            arm, reward = (right, reward_right) if choice else (left, reward_left)
+            values[arm] += alpha * (reward - values[arm])
+            choices.append(choice)
+    data["response"] = data["observed"] = choices
+    return data
+
+
+def test_hessian_of_rlrw_agrees_with_numdifftools_central():
+    from cpm.core.optimisers import finite_difference_hessian
+
+    data = learners()
+    model = rlrw(data)
+    fit = FminBound(model=model, data=data.groupby("ppt"), minimisation=LOSS, prior=True, number_of_starts=1,
+                    initial_guess=[0.5, 5.0], approx_grad=True)
+    fit.optimise()
+    lower, upper = (np.asarray(b, dtype=float) for b in model.parameters.bounds())
+    checked = 0
+    for (_, participant), result in zip(data.groupby("ppt"), fit.fit):
+        x = result["x"]
+        ## numdifftools' own steps must stay within the bounds too
+        if np.any(x <= lower + 0.05 * (upper - lower)) or np.any(x >= upper - 0.05 * (upper - lower)):
+            continue
+        model.reset(data=participant)
+        observed = participant.observed.to_numpy()
+        f = lambda z: objective(z, model, observed, LOSS, True)  # noqa: E731
+        expected = nd.Hessian(f, method="central")(x)
+        H = finite_difference_hessian(f, x, lower, upper)
+        assert np.abs(H - expected).max() <= 1e-4 * np.abs(expected).max()
+        checked += 1
+    assert checked > 0
+
+
+def fit(hessian, temperature_upper=10.0, starts=2, **kwargs):
+    data = bandit()
+    np.random.seed(3)
+    optimiser = FminBound(
+        model=rlrw(data, temperature_upper),
+        data=data.groupby("ppt"),
+        minimisation=LOSS,
+        prior=True,
+        number_of_starts=starts,
+        approx_grad=True,
+        hessian=hessian,
+        **kwargs,
+    )
+    optimiser.optimise()
+    return optimiser, data
+
+
+def test_default_is_the_numdifftools_hessian_next_to_the_optimum():
+    optimiser, data = fit("numdifftools")
+    model = rlrw(data)
+    for (_, participant), result in zip(data.groupby("ppt"), optimiser.fit):
+        model.reset(data=participant)
+        observed = participant.observed.to_numpy()
+        expected = numerical_hessian(lambda z: objective(z, model, observed, LOSS, True), result["x"] + 1e-3)
+        np.testing.assert_array_equal(result["hessian"], expected)
+
+
+def test_only_the_hessian_changes():
+    old, _ = fit("numdifftools")
+    new, _ = fit("finite_differences")
+    assert len(old.fit) == len(new.fit)
+    for a, b in zip(old.fit, new.fit):
+        assert set(a) == set(b)
+        for key in a:
+            if key != "hessian":
+                np.testing.assert_array_equal(a[key], b[key])
+    assert old.parameters == new.parameters
+
+
+def test_hessian_and_evaluate_fit_once_per_participant(monkeypatch):
+    hessians, evaluations = [], []
+    original_hessian, original_evaluate = fmin_module.finite_difference_hessian, fmin_module.evaluate_fit
+
+    def hessian(*args, **kwargs):
+        hessians.append(1)
+        return original_hessian(*args, **kwargs)
+
+    def evaluate(*args, **kwargs):
+        evaluations.append(1)
+        return original_evaluate(*args, **kwargs)
+
+    monkeypatch.setattr(fmin_module, "finite_difference_hessian", hessian)
+    monkeypatch.setattr(fmin_module, "evaluate_fit", evaluate)
+    optimiser, data = fit("finite_differences", starts=3)
+    participants = data.ppt.nunique()
+    assert len(hessians) == participants
+    assert len(evaluations) == participants
+    assert all(np.all(np.isfinite(r["hessian"])) for r in optimiser.fit)
+
+
+def test_only_the_hessian_changes_in_parallel():
+    old, _ = fit("numdifftools", parallel=True, cl=2)
+    new, _ = fit("finite_differences", parallel=True, cl=2)
+    for a, b in zip(old.fit, new.fit):
+        np.testing.assert_array_equal(a["x"], b["x"])
+        assert a["fun"] == b["fun"]
+        assert np.all(np.isfinite(b["hessian"]))
+
+
+def test_no_all_zero_hessian_on_an_upper_bound():
+    old, _ = fit("numdifftools", temperature_upper=0.5)
+    new, _ = fit("finite_differences", temperature_upper=0.5)
+    on_bound = [i for i, r in enumerate(new.fit) if r["x"][1] == 0.5]
+    assert on_bound
+    assert any(np.all(old.fit[i]["hessian"] == 0) for i in on_bound)  # what this fixes
+    for i in on_bound:
+        assert not np.all(new.fit[i]["hessian"] == 0)
+        np.linalg.cholesky(new.fit[i]["hessian"])
+
+
+def test_unknown_hessian_is_an_error():
+    with pytest.raises(ValueError, match="hessian"):
+        fit("central")

From 1cbc820c015e67a8efd7b2bd473fda2948b5a227 Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Tue, 6 Oct 2026 11:11:12 +0200
Subject: [PATCH 4/9] feat(hierarchical): add start="prior_mean"

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                       |  1 +
 benchmarks/hierarchical_profile.py |  6 +-
 cpm/hierarchical/_chains.py        | 26 ++++++++-
 cpm/hierarchical/empirical.py      |  8 ++-
 cpm/hierarchical/variational.py    |  8 ++-
 test/hierarchical/test_start.py    | 91 ++++++++++++++++++++++++++++++
 6 files changed, 135 insertions(+), 5 deletions(-)
 create mode 100644 test/hierarchical/test_start.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e7515df5..df6c5947 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
 - Added a `hessian` argument to `cpm.optimisation.FminBound`. `hessian="finite_differences"` computes the Hessian at the optimum, with finite differences that stay within the bounds (`cpm.core.optimisers.finite_difference_hessian`), once per participant for the best start, together with `log_likelihood`, `log_prior` and `metrics`: 7 instead of 22 evaluations per participant and iteration in the speed comparison, and no Hessian of zeros for an estimate on an upper bound. The default, `"numdifftools"`, keeps the previous Hessian
+- Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`. `start="prior_mean"` starts the first optimisation of every participant from the location of the current prior, clipped into the bounds, as cbm does; the other starts stay random. With two starts, it does not reduce the number of evaluations (70.6 and 71.6 per participant and iteration in the speed comparison). The default, `"random"`, keeps the previous starts
 
 ### Changed
 
diff --git a/benchmarks/hierarchical_profile.py b/benchmarks/hierarchical_profile.py
index bef90780..cc611445 100644
--- a/benchmarks/hierarchical_profile.py
+++ b/benchmarks/hierarchical_profile.py
@@ -96,7 +96,7 @@ def load(path):
     return data
 
 
-def setup(data, iterations, **options):
+def setup(data, iterations, start="random", **options):
     """`VariationalBayes` with the settings of the speed comparison, forced to run `iterations` iterations; `options` go to `FminBound`."""
     first = data[data.ppt == data.ppt.iloc[0]]
     model = RLRW(data=first, dimensions=STIMULI, parameters_settings=[[0.5, *BOUNDS["alpha"]], [2.0, *BOUNDS["temperature"]]])
@@ -121,6 +121,7 @@ def setup(data, iterations, **options):
         tolerance_param=0.0,
         hyperpriors={"a0": np.array([p["mean"] for p in STARTING_PRIOR.values()]), "b": 1, "v": 0.5, "s": np.repeat(0.01, 2)},
         quiet=True,
+        start=start,
     )
 
 
@@ -200,13 +201,14 @@ def main():
     parser.add_argument("--iterations", type=int, default=5)
     parser.add_argument("--label", default="hierarchical-profile")
     parser.add_argument("--hessian", default="numdifftools", help="the `hessian` option of FminBound")
+    parser.add_argument("--start", default="random", help="the `start` option of VariationalBayes")
     args = parser.parse_args()
 
     warnings.simplefilter("ignore")
     np.seterr(all="ignore")
     data = load(args.data) if args.data else simulate()
     participants, iterations = data.ppt.nunique(), args.iterations
-    options = {"hessian": args.hessian}
+    options = {"hessian": args.hessian, "start": args.start}
     print(f"cpm {cpm.__version__}: {participants} participants, {data.trial.max()} trials, {iterations} iterations, 2 starts, {options}")
 
     # the compiled model and everything else that loads on the first call
diff --git a/cpm/hierarchical/_chains.py b/cpm/hierarchical/_chains.py
index 8a9bb2dc..b9799db3 100644
--- a/cpm/hierarchical/_chains.py
+++ b/cpm/hierarchical/_chains.py
@@ -1,8 +1,32 @@
-"""The random starting priors of the chains of the hierarchical methods."""
+"""The random starting priors of the chains of the hierarchical methods, and where their participant-wise fits start."""
 
 import numpy as np
 
 
+def check_start(start):
+    if start not in ("random", "prior_mean"):
+        raise ValueError(f'start must be "random" or "prior_mean", not {start!r}.')
+    return start
+
+
+def start_from_prior_locations(optimiser):
+    """
+    Replace the first start of `optimiser` with the location of the current prior of each free parameter, clipped into its bounds.
+
+    The other starts are the random guesses that `optimiser.reset()` drew, so the
+    same random numbers are drawn as with random starts only.
+
+    Parameters
+    ----------
+    optimiser : object
+        An optimiser of `cpm.optimisation`, after `reset()`.
+    """
+    parameters = optimiser.model.parameters
+    locations = [getattr(parameters, name).prior.kwds["loc"] for name in parameters.free()]
+    lower, upper = parameters.bounds()
+    optimiser.initial_guess[0] = np.clip(locations, lower, upper)
+
+
 def starting_priors(names, bounds, priors):
     """
     Draw a random starting prior for each parameter, for every chain after the first.
diff --git a/cpm/hierarchical/empirical.py b/cpm/hierarchical/empirical.py
index 2e5f6052..3ab7ef63 100644
--- a/cpm/hierarchical/empirical.py
+++ b/cpm/hierarchical/empirical.py
@@ -2,7 +2,7 @@
 import pandas as pd
 import copy
 from ..core.diagnostics import convergence_diagnostics_plots, parameter_bounds
-from ._chains import starting_priors
+from ._chains import check_start, start_from_prior_locations, starting_priors
 
 
 class EmpiricalBayes:
@@ -23,6 +23,8 @@ class EmpiricalBayes:
         The number of random parameter initialisations. Default is 4. The first chain starts from the priors of the model, and every later chain from priors drawn at random within the bounds of the parameters; `numpy.random.seed` makes them reproducible.
     quiet : bool, optional
         Whether to suppress the output. Default is False.
+    start : {"random", "prior_mean"}, optional
+        Where the participant-wise fits start on every iteration. With "random", every start is drawn at random within the bounds of the parameters. With "prior_mean", the first start is the mean (location) of the current prior of each parameter, clipped into its bounds, as the cbm toolbox does, and the others are drawn at random. Default is "random".
 
     Notes
     -----
@@ -89,8 +91,10 @@ def __init__(
         tolerance=1e-6,
         chain=4,
         quiet=False,
+        start="random",
         **kwargs,
     ):
+        self.start = check_start(start)
         self.function = copy.deepcopy(optimiser.model)
         self.optimiser = copy.deepcopy(optimiser)
         # bounds here should include mean and std for all parameters
@@ -219,6 +223,8 @@ def __log_det_hessian_batch(hessian):
         for iteration in range(self.iteration):
 
             self.optimiser.reset()
+            if self.start == "prior_mean":
+                start_from_prior_locations(self.optimiser)
 
             # perform participant-wise optimisation, extracting MAP parameter estimates,
             # the Hessian matrix of the target function evaluated at the MAP parameter estimates,
diff --git a/cpm/hierarchical/variational.py b/cpm/hierarchical/variational.py
index dc5a0b28..9270d938 100644
--- a/cpm/hierarchical/variational.py
+++ b/cpm/hierarchical/variational.py
@@ -8,7 +8,7 @@
 
 from ..generators import Parameters
 from ..core.diagnostics import convergence_diagnostics_plots, parameter_bounds
-from ._chains import starting_priors
+from ._chains import check_start, start_from_prior_locations, starting_priors
 
 
 class VariationalBayes:
@@ -33,6 +33,8 @@ class VariationalBayes:
         A dictionary of given parameter values of the prior distributions on the population-level parameters (means mu and precisions tau). See Notes for details. Default is None.
     convergence : str, optional
         The convergence criterion. Default is 'parameters'. Options are 'lme' and 'parameters'.
+    start : {"random", "prior_mean"}, optional
+        Where the participant-wise fits start on every iteration. With "random", every start is drawn at random within the bounds of the parameters. With "prior_mean", the first start is the mean (location) of the current prior of each parameter, clipped into its bounds, as in the HBI implementation of Piray et al. (2019), and the others are drawn at random. Default is "random".
 
     Notes
     -----
@@ -98,8 +100,10 @@ def __init__(
         hyperpriors=None,
         convergence="parameters",
         quiet=False,
+        start="random",
         **kwargs,
     ):
+        self.start = check_start(start)
         # write input arguments to self
         self.optimiser = copy.deepcopy(optimiser)
         self.objective = objective
@@ -203,6 +207,8 @@ def update_participants(self, iter_idx=0, chain_idx=0):
 
         # delete any pre-existing optimisation output / parameter estimates
         self.optimiser.reset()
+        if self.start == "prior_mean":
+            start_from_prior_locations(self.optimiser)
 
         # run the optimisation
         self.optimiser.optimise()
diff --git a/test/hierarchical/test_start.py b/test/hierarchical/test_start.py
new file mode 100644
index 00000000..729a885b
--- /dev/null
+++ b/test/hierarchical/test_start.py
@@ -0,0 +1,91 @@
+"""`start="prior_mean"`: the first start of every participant is the location of the current prior."""
+
+import warnings
+
+import numpy as np
+import pytest
+
+from cpm.applications.reinforcement_learning import RLRW
+from cpm.datasets import load_bandit_data
+from cpm.hierarchical import EmpiricalBayes, VariationalBayes
+from cpm.optimisation import FminBound, minimise
+
+METHODS = [EmpiricalBayes, VariationalBayes]
+
+
+def method(cls, participants=4, iterations=2, chain=1, **kwargs):
+    warnings.simplefilter("ignore")
+    data = load_bandit_data()
+    data["observed"] = data["response"]
+    data = data[data.ppt.isin(sorted(data.ppt.unique())[:participants])].copy()
+    model = RLRW(data=data[data.ppt == data.ppt.iloc[0]], dimensions=4, parameters_settings=[[0.5, 1e-2, 1.0], [2.0, 0.0, 10.0]])
+    model.parameters.update_prior(alpha={"mean": 0.3, "sd": 0.5}, temperature={"mean": 12.0, "sd": 5.0})
+    optimiser = FminBound(
+        model=model,
+        data=data,
+        minimisation=minimise.LogLikelihood.bernoulli,
+        prior=True,
+        ppt_identifier="ppt",
+        number_of_starts=2,
+        approx_grad=True,
+    )
+    tolerance = {"tolerance": 0.0} if cls is EmpiricalBayes else {"tolerance_lme": 0.0, "tolerance_param": 0.0}
+    return cls(optimiser=optimiser, iteration=iterations, chain=chain, quiet=True, **tolerance, **kwargs)
+
+
+def record_guesses(fit):
+    """The guesses and the prior locations the optimiser starts from, on every participant-wise step."""
+    guesses, locations = [], []
+    optimise = fit.optimiser.optimise
+
+    def recorded(*args, **kwargs):
+        guesses.append(fit.optimiser.initial_guess.copy())
+        parameters = fit.optimiser.model.parameters
+        locations.append([getattr(parameters, name).prior.kwds["loc"] for name in parameters.free()])
+        return optimise(*args, **kwargs)
+
+    fit.optimiser.optimise = recorded
+    return guesses, locations
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_first_start_is_the_prior_location_clipped_into_the_bounds(cls):
+    fit = method(cls, iterations=3, chain=2, start="prior_mean")
+    guesses, locations = record_guesses(fit)
+    np.random.seed(1)
+    fit.optimise()
+    assert len(guesses) == 6
+    lower, upper = fit.optimiser.model.parameters.bounds()
+    np.testing.assert_array_equal(guesses[0][0], [0.3, 10.0])  # 12 is above the upper bound of 10
+    for guess, location in zip(guesses, locations):
+        np.testing.assert_array_equal(guess[0], np.clip(location, lower, upper))
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_other_starts_are_the_random_ones(cls):
+    random, prior_mean = method(cls, start="random"), method(cls, start="prior_mean")
+    random_guesses, _ = record_guesses(random)
+    prior_guesses, _ = record_guesses(prior_mean)
+    np.random.seed(1)
+    random.optimise()
+    np.random.seed(1)
+    prior_mean.optimise()
+    np.testing.assert_array_equal(random_guesses[0][1:], prior_guesses[0][1:])
+    assert not np.array_equal(random_guesses[0][0], prior_guesses[0][0])
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_random_is_the_default(cls):
+    default, random = method(cls), method(cls, start="random")
+    np.random.seed(1)
+    default.optimise()
+    np.random.seed(1)
+    random.optimise()
+    assert default.hyperparameters.equals(random.hyperparameters)
+    assert default.fit.equals(random.fit)
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_unknown_start_is_an_error(cls):
+    with pytest.raises(ValueError, match="start"):
+        method(cls, start="previous")

From 847f080c314e9ef9a142b7bea090d855b02b6832 Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Tue, 6 Oct 2026 11:54:22 +0200
Subject: [PATCH 5/9] bench(hierarchical): record the speed comparison before
 and after

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                  | 17 +++++++++-
 .../hierarchical-comparison-individual.csv    | 31 +++++++++++++++++++
 .../hierarchical-comparison-population.csv    | 31 +++++++++++++++++++
 .../hierarchical-comparison-runtime.csv       | 16 ++++++++++
 .../results/hierarchical-profile-after.csv    | 18 +++++++++++
 5 files changed, 112 insertions(+), 1 deletion(-)
 create mode 100644 benchmarks/results/hierarchical-comparison-individual.csv
 create mode 100644 benchmarks/results/hierarchical-comparison-population.csv
 create mode 100644 benchmarks/results/hierarchical-comparison-runtime.csv
 create mode 100644 benchmarks/results/hierarchical-profile-after.csv

diff --git a/CHANGELOG.md b/CHANGELOG.md
index df6c5947..0b46952f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,12 +11,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
 - Added a `hessian` argument to `cpm.optimisation.FminBound`. `hessian="finite_differences"` computes the Hessian at the optimum, with finite differences that stay within the bounds (`cpm.core.optimisers.finite_difference_hessian`), once per participant for the best start, together with `log_likelihood`, `log_prior` and `metrics`: 7 instead of 22 evaluations per participant and iteration in the speed comparison, and no Hessian of zeros for an estimate on an upper bound. The default, `"numdifftools"`, keeps the previous Hessian
-- Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`. `start="prior_mean"` starts the first optimisation of every participant from the location of the current prior, clipped into the bounds, as cbm does; the other starts stay random. With two starts, it does not reduce the number of evaluations (70.6 and 71.6 per participant and iteration in the speed comparison). The default, `"random"`, keeps the previous starts
+- Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`. `start="prior_mean"` starts the first optimisation of every participant from the location of the current prior, clipped into the bounds, as cbm does; the other starts stay random. It does not reduce the evaluations per iteration, but the fits stop sooner (see the results below). The default, `"random"`, keeps the previous starts
+- Added the results of the speed comparison with cbm for 100 participants and 120 trials (20 datasets, 10 fits at once on 12 cores) in `benchmarks/results/hierarchical-comparison-*.csv`. Median ms per participant and iteration, iterations and seconds per fit, for `EmpiricalBayes` / `VariationalBayes`:
+  - 0.26.0.dev1: 20.8 / 20.1 ms, 50 / 50 iterations, 91 / 88 s
+  - defaults: 17.2 / 16.4 ms, 50 / 50 iterations, 74 / 72 s, with identical results
+  - `hessian="finite_differences"`: 14.1 / 13.8 ms, 42 / 39 iterations, 59 / 53 s
+  - `start="prior_mean"`: 16.2 / 15.8 ms, 41 / 44 iterations, 64 / 65 s
+  - both: 13.0 / 13.2 ms, 38 / 25.5 iterations, 47 / 40 s
+  - both, with `number_of_starts=1`: 6.8 / 7.0 ms, 38.5 / 21.5 iterations, 23 / 18 s
+  - cbm (in the earlier run of the comparison): 4.4 ms, 21 iterations, 9.2 s
+
+  The RMSE of the population means and sds and of the individual estimates stays within the bootstrap 95% intervals of 0.26.0.dev1 in every configuration. With `hessian="finite_differences"`, `EmpiricalBayes` replaces no participant's log model evidence with -1e6, which it did in 16 of the 20 fits before. The hierarchical tutorials (`benchmarks/hierarchical.py`, RLRW) take 7.3 instead of 9.1 s (empirical Bayes) and 7.7 instead of 9.4 s (variational Bayes), with identical results
 
 ### Changed
 
 - `cpm.optimisation.FminBound` with `approx_grad=True` now computes the finite-difference gradient itself, together with the objective, exactly as SciPy does, with identical results; with `LogLikelihood.bernoulli`, this makes the hierarchical fits of the speed comparison 19% faster
 
+### Open questions
+
+- Should `hessian="finite_differences"` and `start="prior_mean"` become the defaults? Both change results for the same seed, so they are options for now.
+- cbm's default is a single start, from the prior mean. With `number_of_starts=1`, cpm's fits take half the time, with the same recovery in the speed comparison; the tutorials and the comparison use two.
+
 ## [0.26.0.dev1] - 2026-10-05
 
 ### Added
diff --git a/benchmarks/results/hierarchical-comparison-individual.csv b/benchmarks/results/hierarchical-comparison-individual.csv
new file mode 100644
index 00000000..99ff9a51
--- /dev/null
+++ b/benchmarks/results/hierarchical-comparison-individual.csv
@@ -0,0 +1,31 @@
+config,method,parameter,individual_rmse
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,alpha,"0.085 [0.073, 0.097]"
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,temperature,"0.763 [0.640, 0.895]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,alpha,"0.085 [0.074, 0.097]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,temperature,"0.763 [0.642, 0.895]"
+"new, defaults",cpm-empirical,alpha,"0.085 [0.073, 0.097]"
+"new, defaults",cpm-empirical,temperature,"0.763 [0.640, 0.895]"
+"new, defaults",cpm-variational,alpha,"0.085 [0.074, 0.097]"
+"new, defaults",cpm-variational,temperature,"0.763 [0.642, 0.895]"
+"new, hessian",cpm-empirical,alpha,"0.084 [0.073, 0.096]"
+"new, hessian",cpm-empirical,temperature,"0.760 [0.637, 0.892]"
+"new, hessian",cpm-variational,alpha,"0.085 [0.074, 0.096]"
+"new, hessian",cpm-variational,temperature,"0.758 [0.636, 0.889]"
+"new, start",cpm-empirical,alpha,"0.085 [0.074, 0.097]"
+"new, start",cpm-empirical,temperature,"0.760 [0.637, 0.893]"
+"new, start",cpm-variational,alpha,"0.085 [0.074, 0.097]"
+"new, start",cpm-variational,temperature,"0.762 [0.640, 0.892]"
+"new, start + hessian",cpm-empirical,alpha,"0.084 [0.073, 0.096]"
+"new, start + hessian",cpm-empirical,temperature,"0.759 [0.635, 0.889]"
+"new, start + hessian",cpm-variational,alpha,"0.085 [0.074, 0.096]"
+"new, start + hessian",cpm-variational,temperature,"0.755 [0.635, 0.884]"
+"new, start + hessian, 1 start",cpm-empirical,alpha,"0.085 [0.073, 0.096]"
+"new, start + hessian, 1 start",cpm-empirical,temperature,"0.758 [0.635, 0.890]"
+"new, start + hessian, 1 start",cpm-variational,alpha,"0.085 [0.074, 0.097]"
+"new, start + hessian, 1 start",cpm-variational,temperature,"0.756 [0.634, 0.886]"
+comparison (cpm-empirical),cpm-empirical,alpha,"0.085 [0.073, 0.097]"
+comparison (cpm-empirical),cpm-empirical,temperature,"0.763 [0.639, 0.895]"
+comparison (cpm-variational),cpm-variational,alpha,"0.085 [0.074, 0.097]"
+comparison (cpm-variational),cpm-variational,temperature,"0.764 [0.642, 0.895]"
+comparison (cbm),cbm,alpha,"0.085 [0.074, 0.097]"
+comparison (cbm),cbm,temperature,"0.765 [0.634, 0.906]"
diff --git a/benchmarks/results/hierarchical-comparison-population.csv b/benchmarks/results/hierarchical-comparison-population.csv
new file mode 100644
index 00000000..c61beb33
--- /dev/null
+++ b/benchmarks/results/hierarchical-comparison-population.csv
@@ -0,0 +1,31 @@
+config,method,parameter,rmse_mean,rmse_sd
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,alpha,"0.013 [0.008, 0.018]","0.026 [0.016, 0.035]"
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,temperature,"0.199 [0.129, 0.262]","0.212 [0.136, 0.281]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,alpha,"0.012 [0.008, 0.015]","0.028 [0.020, 0.035]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,temperature,"0.184 [0.122, 0.240]","0.220 [0.134, 0.296]"
+"new, defaults",cpm-empirical,alpha,"0.013 [0.008, 0.018]","0.026 [0.016, 0.035]"
+"new, defaults",cpm-empirical,temperature,"0.199 [0.129, 0.262]","0.212 [0.136, 0.281]"
+"new, defaults",cpm-variational,alpha,"0.012 [0.008, 0.015]","0.028 [0.020, 0.035]"
+"new, defaults",cpm-variational,temperature,"0.184 [0.122, 0.240]","0.220 [0.134, 0.296]"
+"new, hessian",cpm-empirical,alpha,"0.012 [0.008, 0.017]","0.021 [0.011, 0.028]"
+"new, hessian",cpm-empirical,temperature,"0.197 [0.128, 0.261]","0.191 [0.109, 0.256]"
+"new, hessian",cpm-variational,alpha,"0.012 [0.008, 0.015]","0.025 [0.017, 0.031]"
+"new, hessian",cpm-variational,temperature,"0.180 [0.121, 0.235]","0.170 [0.105, 0.225]"
+"new, start",cpm-empirical,alpha,"0.013 [0.008, 0.017]","0.027 [0.016, 0.037]"
+"new, start",cpm-empirical,temperature,"0.199 [0.130, 0.262]","0.211 [0.136, 0.281]"
+"new, start",cpm-variational,alpha,"0.012 [0.009, 0.015]","0.029 [0.021, 0.036]"
+"new, start",cpm-variational,temperature,"0.186 [0.124, 0.242]","0.211 [0.133, 0.282]"
+"new, start + hessian",cpm-empirical,alpha,"0.012 [0.008, 0.017]","0.021 [0.012, 0.028]"
+"new, start + hessian",cpm-empirical,temperature,"0.200 [0.129, 0.262]","0.191 [0.108, 0.256]"
+"new, start + hessian",cpm-variational,alpha,"0.012 [0.008, 0.015]","0.025 [0.017, 0.031]"
+"new, start + hessian",cpm-variational,temperature,"0.183 [0.123, 0.236]","0.170 [0.105, 0.225]"
+"new, start + hessian, 1 start",cpm-empirical,alpha,"0.012 [0.008, 0.016]","0.022 [0.013, 0.028]"
+"new, start + hessian, 1 start",cpm-empirical,temperature,"0.200 [0.130, 0.262]","0.180 [0.094, 0.249]"
+"new, start + hessian, 1 start",cpm-variational,alpha,"0.011 [0.009, 0.014]","0.025 [0.018, 0.031]"
+"new, start + hessian, 1 start",cpm-variational,temperature,"0.183 [0.122, 0.236]","0.171 [0.107, 0.226]"
+comparison (cpm-empirical),cpm-empirical,alpha,"0.013 [0.008, 0.018]","0.026 [0.016, 0.035]"
+comparison (cpm-empirical),cpm-empirical,temperature,"0.199 [0.127, 0.259]","0.213 [0.133, 0.279]"
+comparison (cpm-variational),cpm-variational,alpha,"0.012 [0.008, 0.015]","0.028 [0.021, 0.035]"
+comparison (cpm-variational),cpm-variational,temperature,"0.184 [0.123, 0.240]","0.205 [0.126, 0.279]"
+comparison (cbm),cbm,alpha,"0.023 [0.018, 0.029]","0.022 [0.015, 0.028]"
+comparison (cbm),cbm,temperature,"0.279 [0.183, 0.362]","0.176 [0.115, 0.234]"
diff --git a/benchmarks/results/hierarchical-comparison-runtime.csv b/benchmarks/results/hierarchical-comparison-runtime.csv
new file mode 100644
index 00000000..6e46141b
--- /dev/null
+++ b/benchmarks/results/hierarchical-comparison-runtime.csv
@@ -0,0 +1,16 @@
+config,method,fits,median_seconds,median_ms_per_participant_iteration,median_iterations,stopped_before_50,seconds_vs_cbm,ms_vs_cbm
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,20,90.87562635005452,20.803921870002526,50.0,0.45,9.902121712865382,4.736912106409524
+0.26.0.dev1 (Phase 0 commit),cpm-variational,20,87.99754375009798,20.113323750020935,50.0,0.45,9.588515905135699,4.579667591858573
+"new, defaults",cpm-empirical,20,73.94655685010366,17.220056519983338,50.0,0.45,8.057471905134292,3.9208902394592307
+"new, defaults",cpm-variational,20,72.27422889997251,16.375839510001242,50.0,0.45,7.875249283171419,3.7286677441042935
+"new, hessian",cpm-empirical,20,58.99186499998905,14.11793514001183,42.0,0.6,6.427957096534108,3.2145582116733054
+"new, hessian",cpm-variational,20,52.845357749960385,13.838825339986943,39.0,0.65,5.758212464850953,3.151006800597213
+"new, start",cpm-empirical,20,64.46359169995412,16.2473555700155,41.0,0.55,7.024175311898925,3.6994128211815522
+"new, start",cpm-variational,20,65.32937605003826,15.764868770822067,44.0,0.65,7.118514161114703,3.5895538448518227
+"new, start + hessian",cpm-empirical,20,47.108444100013,13.017018705396913,38.0,0.7,5.1330985646820135,2.9638869959352707
+"new, start + hessian",cpm-variational,20,40.04392045002896,13.229554900014772,25.5,0.85,4.363323699460235,3.012280047958179
+"new, start + hessian, 1 start",cpm-empirical,20,23.012659750063904,6.773201450007036,38.5,0.65,2.5075387860779714,1.5422120957852348
+"new, start + hessian, 1 start",cpm-variational,20,17.833935800008476,6.9956334311127835,21.5,0.85,1.9432471610249897,1.592858351368021
+comparison (cpm-empirical),cpm-empirical,20,81.6190145999426,19.751574295005412,50.0,0.45,8.893489366879983,4.497299691053249
+comparison (cpm-variational),cpm-variational,20,70.75931959995069,14.922597907213236,50.0,0.45,7.710179540323271,3.3977744738449753
+comparison (cbm),cbm,20,9.17738935,4.391874158241759,21.0,0.9,1.0,1.0
diff --git a/benchmarks/results/hierarchical-profile-after.csv b/benchmarks/results/hierarchical-profile-after.csv
new file mode 100644
index 00000000..d97e34c2
--- /dev/null
+++ b/benchmarks/results/hierarchical-profile-after.csv
@@ -0,0 +1,18 @@
+section,item,value
+fit,seconds without profiler,6.1314
+fit,ms per participant and iteration,12.2628
+evaluations per participant and iteration,optimiser,72.198
+evaluations per participant and iteration,hessian,6.992
+evaluations per participant and iteration,evaluate_fit,1.0
+share of profiled time,fmin_l_bfgs_b,0.8817
+share of profiled time,hessian,0.0775
+share of profiled time,evaluate_fit,0.008
+share of profiled time,update_population,0.0007
+share of profiled time,get_lme,0.0001
+share of profiled time,objective (all callers),0.7932
+share of objective,SessionWrapper.run,0.4603
+share of objective,  of which the RLRW model function,0.3393
+share of objective,Wrapper.reset,0.1391
+share of objective,deepcopy,0.0459
+share of objective,loss,0.2029
+share of objective,Parameters.PDF,0.1232

From 45f9a1821eb70d764b38adbe294f8276d4052a82 Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Tue, 6 Oct 2026 12:29:37 +0200
Subject: [PATCH 6/9] feat(optimisation): add the hessian option to Fmin,
 Minimize and Bads

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                             |   2 +-
 cpm/core/optimisers.py                   |  10 ++
 cpm/optimisation/bads.py                 |  90 +++++++++++++++---
 cpm/optimisation/fmin.py                 |  90 ++++++++++++++----
 cpm/optimisation/free.py                 |  86 ++++++++++++++---
 test/optimisation/test_hessian_option.py | 116 +++++++++++++++++++++++
 6 files changed, 350 insertions(+), 44 deletions(-)
 create mode 100644 test/optimisation/test_hessian_option.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0b46952f..35f953ca 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,7 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
-- Added a `hessian` argument to `cpm.optimisation.FminBound`. `hessian="finite_differences"` computes the Hessian at the optimum, with finite differences that stay within the bounds (`cpm.core.optimisers.finite_difference_hessian`), once per participant for the best start, together with `log_likelihood`, `log_prior` and `metrics`: 7 instead of 22 evaluations per participant and iteration in the speed comparison, and no Hessian of zeros for an estimate on an upper bound. The default, `"numdifftools"`, keeps the previous Hessian
+- Added a `hessian` argument to `cpm.optimisation.FminBound`, `Fmin`, `Minimize` and `Bads`. `hessian="finite_differences"` computes the Hessian at the optimum, with finite differences that stay within the bounds (`cpm.core.optimisers.finite_difference_hessian`), once per participant for the best start, together with `log_likelihood`, `log_prior` and `metrics`: 7 instead of 22 evaluations per participant and iteration in the speed comparison, and no Hessian of zeros for an estimate on an upper bound. The default, `"numdifftools"`, keeps the previous Hessian
 - Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`. `start="prior_mean"` starts the first optimisation of every participant from the location of the current prior, clipped into the bounds, as cbm does; the other starts stay random. It does not reduce the evaluations per iteration, but the fits stop sooner (see the results below). The default, `"random"`, keeps the previous starts
 - Added the results of the speed comparison with cbm for 100 participants and 120 trials (20 datasets, 10 fits at once on 12 cores) in `benchmarks/results/hierarchical-comparison-*.csv`. Median ms per participant and iteration, iterations and seconds per fit, for `EmpiricalBayes` / `VariationalBayes`:
   - 0.26.0.dev1: 20.8 / 20.1 ms, 50 / 50 iterations, 91 / 88 s
diff --git a/cpm/core/optimisers.py b/cpm/core/optimisers.py
index cc2f4e3b..4ba684bb 100644
--- a/cpm/core/optimisers.py
+++ b/cpm/core/optimisers.py
@@ -13,6 +13,7 @@
     "prepare_data",
     "numerical_hessian",
     "finite_difference_hessian",
+    "check_hessian",
 ]
 
 
@@ -35,6 +36,15 @@ def numerical_hessian(func=None, params=None, hessian=None):
     return computed_hessian
 
 
+def check_hessian(hessian):
+    """Raise a ValueError unless `hessian` is one of the optimisers' `hessian` options."""
+    if hessian not in ("numdifftools", "finite_differences"):
+        raise ValueError(
+            f'hessian must be "numdifftools" or "finite_differences", not {hessian!r}.'
+        )
+    return hessian
+
+
 def finite_difference_hessian(func, x, lower, upper, step=1e-4):
     """
     The Hessian of `func` at `x` by finite differences that stay within the bounds.
diff --git a/cpm/optimisation/bads.py b/cpm/optimisation/bads.py
index 9b853cc6..5b972fcf 100644
--- a/cpm/optimisation/bads.py
+++ b/cpm/optimisation/bads.py
@@ -5,6 +5,8 @@
     evaluate_fit,
     fit_extras,
     numerical_hessian,
+    finite_difference_hessian,
+    check_hessian,
     group_identifier,
     prepare_data,
 )
@@ -95,6 +97,8 @@ class Bads:
         The libraries required for the parallel processing with `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`.
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `fit` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds of the parameters (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the `pybads.bads <https://acerbilab.github.io/pybads/api/classes/bads.html>`__ documentation for what is supported.
 
@@ -137,6 +141,7 @@ def __init__(
         libraries=["numpy", "pandas"],
         ppt_identifier=None,
         display=True,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -150,6 +155,7 @@ def __init__(
         self.prior = prior
         self.metrics = metrics
         self.kwargs = numpy_bool_options(kwargs)
+        self.hessian = check_hessian(hessian)
 
         self.fit = []
         self.details = []
@@ -243,6 +249,10 @@ def target(x):
 
             result = optimizer.optimize()
             result = dict(result.items())
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                result.update({"ppt": ppt})
+                return result
 
             def f(x):
                 return objective(x, model, observed, loss, prior)
@@ -262,16 +272,56 @@ def f(x):
             result.update({"ppt": ppt})
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = {**result, "hessian": finite_difference_hessian(f, result["x"], lower, upper)}
+            if prior or self.metrics:
+                result.update(
+                    evaluate_fit(
+                        result["x"], model, observed, loss, prior, self.metrics
+                    )
+                )
+            ## the identifier last, as in __task
+            result["ppt"] = result.pop("ppt")
+            return result
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i].get("fval")
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            for result in results:
+                parameters = {}
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result["x"][i]
+                    )
+
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
         loss = self.loss
         model = self.model
         prior = self.prior
+        lower, upper = (np.asarray(b, dtype=float) for b in self.__bounds__)
 
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -291,21 +341,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
-                self.details = copy.deepcopy(results)
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                parameters = {}
-                for result in results:
-                    parameters = {}
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result["x"][i]
-                        )
-
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -320,6 +360,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    method=None,
+                    pandas=False,
+                    cl=self.cl,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
diff --git a/cpm/optimisation/fmin.py b/cpm/optimisation/fmin.py
index e493a88a..a916cf83 100644
--- a/cpm/optimisation/fmin.py
+++ b/cpm/optimisation/fmin.py
@@ -5,6 +5,7 @@
     fit_extras,
     numerical_hessian,
     finite_difference_hessian,
+    check_hessian,
     group_identifier,
     prepare_data,
 )
@@ -154,6 +155,8 @@ class Fmin:
         The libraries to import for parallel processing for `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`.
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `fit` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds of the parameters, if the optimum is within them (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the :func:`scipy.optimize.fmin` documentation for what is supported.
 
@@ -192,6 +195,7 @@ def __init__(
         number_of_starts=1,
         ppt_identifier=None,
         display=False,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -206,6 +210,7 @@ def __init__(
         self.metrics = metrics
         self.kwargs = kwargs
         self.display = display
+        self.hessian = check_hessian(hessian)
 
         self.fit = []
         self.details = []
@@ -278,6 +283,9 @@ def __task(participant, **args):
                 **self.kwargs,
                 full_output=True,
             )
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                return (*result, ppt)
 
             def f(x):
                 return objective(x, model, observed, loss, prior)
@@ -297,16 +305,54 @@ def f(x):
             result = (*result, ppt)
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = (*result[:-1], finite_difference_hessian(f, result[0], lower, upper))
+            if prior or self.metrics:
+                result = (
+                    *result,
+                    *evaluate_fit(
+                        result[0], model, observed, loss, prior, self.metrics
+                    ).values(),
+                )
+            return (*result, ppt)
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i][1]
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            parameters = {}
+            for result in results:
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result[0][i]
+                    )
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
         loss = self.loss
         model = self.model
         prior = self.prior
+        lower, upper = (np.asarray(b, dtype=float) for b in self.model.parameters.bounds())
 
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -326,19 +372,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                self.details = copy.deepcopy(results)
-                parameters = {}
-                for result in results:
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result[0][i]
-                        )
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -353,6 +391,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    pandas=False,
+                    method=None,
+                    cl=self.cl,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
@@ -481,11 +541,7 @@ def __init__(
         self.metrics = metrics
         self.kwargs = kwargs
         self.display = display
-        if hessian not in ("numdifftools", "finite_differences"):
-            raise ValueError(
-                f'hessian must be "numdifftools" or "finite_differences", not {hessian!r}.'
-            )
-        self.hessian = hessian
+        self.hessian = check_hessian(hessian)
 
         self.fit = []
         self.details = []
diff --git a/cpm/optimisation/free.py b/cpm/optimisation/free.py
index 39eaf514..304492ad 100644
--- a/cpm/optimisation/free.py
+++ b/cpm/optimisation/free.py
@@ -4,6 +4,8 @@
     evaluate_fit,
     fit_extras,
     numerical_hessian,
+    finite_difference_hessian,
+    check_hessian,
     group_identifier,
     prepare_data,
 )
@@ -48,6 +50,8 @@ class Minimize:
         The libraries to import for parallel processing for `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `details` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds of the parameters (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the :func:`scipy.optimize.minimize` documentation for what is supported.
 
@@ -87,6 +91,7 @@ def __init__(
         number_of_starts=1,
         ppt_identifier=None,
         display=False,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -101,6 +106,7 @@ def __init__(
         self.metrics = metrics
         self.kwargs = kwargs
         self.display = display
+        self.hessian = check_hessian(hessian)
 
         self.method = method
         if isinstance(model, Wrapper):
@@ -173,6 +179,11 @@ def __task(participant, **args):
                 **self.kwargs,
             )
             result = {**result}
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                if self.ppt_identifier is not None:
+                    result.update({self.ppt_identifier: ppt})
+                return result
 
             def f(x):
                 return objective(x, model, observed, loss, prior)
@@ -193,12 +204,48 @@ def f(x):
                 result.update({self.ppt_identifier: ppt})
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = {**result, "hessian": finite_difference_hessian(f, result["x"], lower, upper)}
+            if prior or self.metrics:
+                result.update(
+                    evaluate_fit(
+                        result["x"], model, observed, loss, prior, self.metrics
+                    )
+                )
+            return result
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i]["fun"]
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            parameters = {}
+            for result in results:
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result["x"][i]
+                    )
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
+        lower, upper = (np.asarray(b, dtype=float) for b in self.model.parameters.bounds())
         bounds = self.model.parameters.bounds()
         bounds = np.asarray(bounds).T
         bounds = list(map(tuple, bounds))
@@ -208,6 +255,7 @@ def __extract_nll(result):
         model = self.model
         prior = self.prior
 
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -227,19 +275,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                self.details = copy.deepcopy(results)
-                parameters = {}
-                for result in results:
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result["x"][i]
-                        )
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -254,6 +294,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    method=None,
+                    cl=self.cl,
+                    pandas=False,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
diff --git a/test/optimisation/test_hessian_option.py b/test/optimisation/test_hessian_option.py
new file mode 100644
index 00000000..cecfdec2
--- /dev/null
+++ b/test/optimisation/test_hessian_option.py
@@ -0,0 +1,116 @@
+"""
+The `hessian` option of every optimiser: "numdifftools" (the default, the previous
+Hessian) or "finite_differences" (at the optimum, within the bounds, once per
+participant for the best start).
+"""
+
+import numpy as np
+import pytest
+
+import cpm.optimisation.fmin as fmin_module
+import cpm.optimisation.free as free_module
+from cpm.core.optimisers import finite_difference_hessian, objective
+from cpm.optimisation import Fmin, FminBound, Minimize, minimise
+from test.optimisation.test_fminbound_gradient import bandit, rlrw
+
+try:
+    import cpm.optimisation.bads as bads_module
+    from cpm.optimisation import Bads
+
+    HAS_BADS = True
+except ImportError:  # pragma: no cover - pybads is an optional dependency
+    HAS_BADS = False
+
+LOSS = minimise.LogLikelihood.bernoulli
+OPTIMISERS = [Fmin, FminBound, Minimize]
+MODULES = {Fmin: fmin_module, FminBound: fmin_module, Minimize: free_module}
+if HAS_BADS:
+    OPTIMISERS.append(Bads)
+    MODULES[Bads] = bads_module
+IDS = [cls.__name__ for cls in OPTIMISERS]
+PARTICIPANTS = 3
+## Bads records its own run time
+VARYING = {"total_time"}
+
+
+def build(cls, hessian=None):
+    data = bandit(participants=PARTICIPANTS, trials=30)
+    kwargs = {} if hessian is None else {"hessian": hessian}
+    if cls is FminBound:
+        kwargs["approx_grad"] = True
+    if HAS_BADS and cls is Bads:
+        kwargs["options"] = {"display": "off"}
+    np.random.seed(4)
+    optimiser = cls(model=rlrw(data), data=data.groupby("ppt"), minimisation=LOSS, prior=True, number_of_starts=2, **kwargs)
+    optimiser.optimise()
+    return optimiser, data
+
+
+def records(optimiser):
+    """The per-participant records that hold the Hessian: `fit`, except for `Minimize`, which keeps it in `details`."""
+    return optimiser.details if isinstance(optimiser, Minimize) else optimiser.fit
+
+
+def assert_same(a, b, skip=()):
+    assert set(a) == set(b)
+    for key in set(a) - set(skip) - VARYING:
+        np.testing.assert_array_equal(np.asarray(a[key], dtype=object), np.asarray(b[key], dtype=object), err_msg=key)
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_default_is_numdifftools(cls):
+    default, _ = build(cls)
+    explicit, _ = build(cls, "numdifftools")
+    assert default.hessian == "numdifftools"
+    for a, b in zip(default.fit, explicit.fit):
+        assert_same(a, b)
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_only_the_hessian_changes(cls):
+    old, _ = build(cls, "numdifftools")
+    new, _ = build(cls, "finite_differences")
+    for a, b in zip(old.fit, new.fit):
+        assert_same(a, b, skip={"hessian"})
+    for a, b in zip(records(old), records(new)):
+        assert not np.array_equal(a["hessian"], b["hessian"])
+    assert old.parameters == new.parameters
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_hessian_is_the_finite_difference_one_at_the_optimum(cls):
+    fit, data = build(cls, "finite_differences")
+    model = rlrw(data)
+    lower, upper = (np.asarray(b, dtype=float) for b in model.parameters.bounds())
+    for (_, participant), record in zip(data.groupby("ppt"), records(fit)):
+        model.reset(data=participant)
+        observed = participant.observed.to_numpy()
+        x = record["xopt"] if "xopt" in record else record["x"]  # Fmin records its estimates as xopt
+        expected = finite_difference_hessian(lambda z: objective(z, model, observed, LOSS, True), x, lower, upper)
+        np.testing.assert_array_equal(record["hessian"], expected)
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_hessian_and_extras_once_per_participant(cls, monkeypatch):
+    module = MODULES[cls]
+    calls = {"hessian": 0, "evaluate_fit": 0}
+    original_hessian, original_evaluate = module.finite_difference_hessian, module.evaluate_fit
+
+    def hessian(*args, **kwargs):
+        calls["hessian"] += 1
+        return original_hessian(*args, **kwargs)
+
+    def evaluate(*args, **kwargs):
+        calls["evaluate_fit"] += 1
+        return original_evaluate(*args, **kwargs)
+
+    monkeypatch.setattr(module, "finite_difference_hessian", hessian)
+    monkeypatch.setattr(module, "evaluate_fit", evaluate)
+    build(cls, "finite_differences")
+    assert calls == {"hessian": PARTICIPANTS, "evaluate_fit": PARTICIPANTS}
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_unknown_hessian_is_an_error(cls):
+    with pytest.raises(ValueError, match="hessian"):
+        build(cls, "central")

From c3e98545385cc77dea1186c59e58f3dd4901b7e6 Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Tue, 6 Oct 2026 12:33:48 +0200
Subject: [PATCH 7/9] feat(hierarchical): record Hessians that are not positive
 definite

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md                                |  7 +-
 cpm/hierarchical/_chains.py                 | 25 +++++++
 cpm/hierarchical/empirical.py               |  8 ++-
 cpm/hierarchical/variational.py             | 11 +++-
 test/hierarchical/test_positive_definite.py | 73 +++++++++++++++++++++
 5 files changed, 116 insertions(+), 8 deletions(-)
 create mode 100644 test/hierarchical/test_positive_definite.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 35f953ca..6b8d8f48 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
 - Added a `hessian` argument to `cpm.optimisation.FminBound`, `Fmin`, `Minimize` and `Bads`. `hessian="finite_differences"` computes the Hessian at the optimum, with finite differences that stay within the bounds (`cpm.core.optimisers.finite_difference_hessian`), once per participant for the best start, together with `log_likelihood`, `log_prior` and `metrics`: 7 instead of 22 evaluations per participant and iteration in the speed comparison, and no Hessian of zeros for an estimate on an upper bound. The default, `"numdifftools"`, keeps the previous Hessian
 - Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`. `start="prior_mean"` starts the first optimisation of every participant from the location of the current prior, clipped into the bounds, as cbm does; the other starts stay random. It does not reduce the evaluations per iteration, but the fits stop sooner (see the results below). The default, `"random"`, keeps the previous starts
+- Added a `positive_definite` column to the `fit` of `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`, which records whether each participant's Hessian was positive definite on each iteration, and a `not_positive_definite` column to their `hyperparameters`, with the number that were not. An estimate on a bound can have a Hessian that is not positive definite; it is kept, as the bounds require. In the speed comparison, this happened in about 0.1% of the participant-wise fits with `hessian="finite_differences"`, all of them on a bound
 - Added the results of the speed comparison with cbm for 100 participants and 120 trials (20 datasets, 10 fits at once on 12 cores) in `benchmarks/results/hierarchical-comparison-*.csv`. Median ms per participant and iteration, iterations and seconds per fit, for `EmpiricalBayes` / `VariationalBayes`:
   - 0.26.0.dev1: 20.8 / 20.1 ms, 50 / 50 iterations, 91 / 88 s
   - defaults: 17.2 / 16.4 ms, 50 / 50 iterations, 74 / 72 s, with identical results
@@ -26,11 +27,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 
 - `cpm.optimisation.FminBound` with `approx_grad=True` now computes the finite-difference gradient itself, together with the objective, exactly as SciPy does, with identical results; with `LogLikelihood.bernoulli`, this makes the hierarchical fits of the speed comparison 19% faster
-
-### Open questions
-
-- Should `hessian="finite_differences"` and `start="prior_mean"` become the defaults? Both change results for the same seed, so they are options for now.
-- cbm's default is a single start, from the prior mean. With `number_of_starts=1`, cpm's fits take half the time, with the same recovery in the speed comparison; the tutorials and the comparison use two.
+- The defaults of the new arguments, `hessian="numdifftools"` and `start="random"`, give the same results as before for the same seed; the faster options have to be chosen
 
 ## [0.26.0.dev1] - 2026-10-05
 
diff --git a/cpm/hierarchical/_chains.py b/cpm/hierarchical/_chains.py
index b9799db3..832cc77d 100644
--- a/cpm/hierarchical/_chains.py
+++ b/cpm/hierarchical/_chains.py
@@ -3,6 +3,31 @@
 import numpy as np
 
 
+def positive_definite(hessians):
+    """
+    Whether each participant's Hessian is positive definite (finite, with a Cholesky decomposition).
+
+    Parameters
+    ----------
+    hessians : array-like
+        The Hessians of the negative log posterior density, shape (participants, parameters, parameters).
+
+    Returns
+    -------
+    numpy.ndarray
+        One boolean per participant.
+    """
+    out = np.zeros(len(hessians), dtype=bool)
+    for i, hessian in enumerate(hessians):
+        if np.all(np.isfinite(hessian)):
+            try:
+                np.linalg.cholesky(hessian)
+                out[i] = True
+            except np.linalg.LinAlgError:
+                pass
+    return out
+
+
 def check_start(start):
     if start not in ("random", "prior_mean"):
         raise ValueError(f'start must be "random" or "prior_mean", not {start!r}.')
diff --git a/cpm/hierarchical/empirical.py b/cpm/hierarchical/empirical.py
index 3ab7ef63..78635a67 100644
--- a/cpm/hierarchical/empirical.py
+++ b/cpm/hierarchical/empirical.py
@@ -2,7 +2,7 @@
 import pandas as pd
 import copy
 from ..core.diagnostics import convergence_diagnostics_plots, parameter_bounds
-from ._chains import check_start, start_from_prior_locations, starting_priors
+from ._chains import check_start, positive_definite, start_from_prior_locations, starting_priors
 
 
 class EmpiricalBayes:
@@ -51,6 +51,8 @@ class EmpiricalBayes:
     - When calculating the within-subject variance via the Hessian matrix, the algorithm sets any non-finite or non-positive values to NaN.
     - If the second derivative of the negative log posterior density function is not finite, we set the log determinant to -1e6.
 
+    An estimate on a bound of a parameter can have a Hessian matrix that is not positive definite. The estimate is kept, as the bounds require. The `positive_definite` column of `fit` records whether each participant's Hessian matrix was positive definite on each iteration, and the `not_positive_definite` column of `hyperparameters` how many were not.
+
     References
     ----------
 
@@ -254,6 +256,9 @@ def __log_det_hessian_batch(hessian):
             parameter_long["ppt"] = np.arange(len(parameters))
             parameter_long["iteration"] = iteration + 1
             parameter_long["chain"] = chain_index
+            ## an estimate on a bound can have a Hessian that is not positive definite
+            definite = positive_definite(hessian)
+            parameter_long["positive_definite"] = definite
             fit_records.append(parameter_long)
             # turn any non-finite values into NaN, to avoid subsequent issues
             # with calculating parameter means and variances
@@ -339,6 +344,7 @@ def __log_det_hessian_batch(hessian):
                         "sd": stdev[i],
                         "lme": summed_lme,
                         "reject": reject,
+                        "not_positive_definite": int((~definite).sum()),
                     }
                 )
 
diff --git a/cpm/hierarchical/variational.py b/cpm/hierarchical/variational.py
index 9270d938..a0453740 100644
--- a/cpm/hierarchical/variational.py
+++ b/cpm/hierarchical/variational.py
@@ -8,7 +8,7 @@
 
 from ..generators import Parameters
 from ..core.diagnostics import convergence_diagnostics_plots, parameter_bounds
-from ._chains import check_start, start_from_prior_locations, starting_priors
+from ._chains import check_start, positive_definite, start_from_prior_locations, starting_priors
 
 
 class VariationalBayes:
@@ -60,6 +60,8 @@ class VariationalBayes:
 
     The convergence criterion can be set to 'lme' or 'parameters'. If set to 'lme', the algorithm will stop when the log model evidence converges. If set to 'parameters', the algorithm will stop when the "normalized" means of the population-level parameters converge.
 
+    An estimate on a bound of a parameter can have a Hessian matrix that is not positive definite. The estimate is kept, as the bounds require. The `positive_definite` column of `fit` records whether each participant's Hessian matrix was positive definite on each iteration, and the `not_positive_definite` column of `hyperparameters` how many were not.
+
     References
     ----------
 
@@ -232,7 +234,6 @@ def update_participants(self, iter_idx=0, chain_idx=0):
         parameter_long["ppt"] = np.arange(self.__n_ppt__)
         parameter_long["iteration"] = iter_idx + 1
         parameter_long["chain"] = chain_idx
-        self.__fit_chunks__.append(parameter_long)
 
         # extract the participant-wise unnormalised log posterior density at the
         # optimised parameter values
@@ -253,6 +254,10 @@ def update_participants(self, iter_idx=0, chain_idx=0):
         if self.objective != "minimise":
             hessian = -1 * hessian
 
+        ## an estimate on a bound can have a Hessian that is not positive definite
+        parameter_long["positive_definite"] = positive_definite(hessian)
+        self.__fit_chunks__.append(parameter_long)
+
         return param, log_posterior, hessian
 
     def get_lme(self, log_post, hessian):
@@ -505,6 +510,7 @@ def __inv_hessian_batch(hessian):
         # buffered (see the `hyperparameters` property above) rather than
         # concatenated into self.hyperparameters on every single parameter,
         # to avoid an O(n^2) cost across many iterations x parameters.
+        not_definite = int((~positive_definite(hessian)).sum())
         for i, name in enumerate(self.__param_names__):
             self.__hyper_chunks__.append(
                 {
@@ -515,6 +521,7 @@ def __inv_hessian_batch(hessian):
                     "mean_se": E_mu_error[i],
                     "sd": E_sd[i],
                     "lme": lme,
+                    "not_positive_definite": not_definite,
                 }
             )
 
diff --git a/test/hierarchical/test_positive_definite.py b/test/hierarchical/test_positive_definite.py
new file mode 100644
index 00000000..76c3dc81
--- /dev/null
+++ b/test/hierarchical/test_positive_definite.py
@@ -0,0 +1,73 @@
+"""
+The hierarchical methods record, per participant and iteration, whether the
+Hessian of the negative log posterior is positive definite, and per iteration how
+many are not. An estimate on a bound can have a Hessian that is not; it is kept.
+"""
+
+import warnings
+
+import numpy as np
+import pytest
+
+from cpm.hierarchical import EmpiricalBayes, VariationalBayes
+from cpm.optimisation import FminBound, minimise
+from test.optimisation.test_fminbound_gradient import bandit, rlrw
+
+METHODS = [EmpiricalBayes, VariationalBayes]
+
+
+def method(cls, hessian):
+    warnings.simplefilter("ignore")
+    data = bandit(participants=6, trials=40)
+    ## with an upper bound of 1 on the temperature, estimates end on it, where the
+    ## "numdifftools" Hessian is all zeros (with 0.5, EmpiricalBayes fails on
+    ## these all-zero Hessians, as before)
+    optimiser = FminBound(model=rlrw(data, temperature_upper=1.0), data=data.groupby("ppt"),
+                          minimisation=minimise.LogLikelihood.bernoulli, prior=True, number_of_starts=1,
+                          approx_grad=True, hessian=hessian)
+    tolerance = {"tolerance": 0.0} if cls is EmpiricalBayes else {"tolerance_lme": 0.0, "tolerance_param": 0.0}
+    fit = cls(optimiser=optimiser, iteration=3, chain=1, quiet=True, **tolerance)
+    hessians = []
+    optimise = fit.optimiser.optimise
+
+    def recorded(*args, **kwargs):
+        result = optimise(*args, **kwargs)
+        hessians.append([r["hessian"] for r in fit.optimiser.fit])
+        return result
+
+    fit.optimiser.optimise = recorded
+    np.random.seed(2)
+    fit.optimise()
+    return fit, hessians
+
+
+def is_positive_definite(h):
+    if not np.all(np.isfinite(h)):
+        return False
+    try:
+        np.linalg.cholesky(h)
+        return True
+    except np.linalg.LinAlgError:
+        return False
+
+
+@pytest.mark.parametrize("hessian", ["numdifftools", "finite_differences"])
+@pytest.mark.parametrize("cls", METHODS)
+def test_flags_and_counts(cls, hessian):
+    fit, hessians = method(cls, hessian)
+    expected = np.array([[is_positive_definite(h) for h in step] for step in hessians])
+    flags = fit.fit.sort_values(["iteration", "ppt"]).positive_definite.to_numpy().reshape(expected.shape)
+    np.testing.assert_array_equal(flags, expected)
+    counts = fit.hyperparameters.groupby("iteration").not_positive_definite.first().to_numpy()
+    np.testing.assert_array_equal(counts, (~expected).sum(axis=1))
+    if hessian == "numdifftools":
+        assert counts.sum() > 0  # the all-zero Hessians on the bound
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_estimates_on_a_bound_are_kept(cls):
+    fit, _ = method(cls, "numdifftools")
+    flagged = fit.fit[~fit.fit.positive_definite]
+    assert len(flagged) > 0
+    ## each with a parameter on a bound: alpha in [0.01, 1], temperature in [0, 1]
+    assert (flagged.alpha.isin([0.01, 1.0]) | flagged.temperature.isin([0.0, 1.0])).all()

From 6d1a05acab43ab71f07f3ec694469f4c1d304d61 Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Wed, 14 Oct 2026 10:16:19 +0200
Subject: [PATCH 8/9] docs(CHANGELOG): shorten summaries

---
 CHANGELOG.md | 19 +++++--------------
 1 file changed, 5 insertions(+), 14 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6b8d8f48..c354dfa0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,20 +9,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
-- Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration. On the speed comparison with cbm (RLRW, 100 participants, 120 trials), cpm made 71 evaluations in the optimiser and 22 for the Hessian, against cbm's 72 in all, at 17.3 against 5.0 ms per participant and iteration: the gap is the cost of each evaluation, not their number
-- Added a `hessian` argument to `cpm.optimisation.FminBound`, `Fmin`, `Minimize` and `Bads`. `hessian="finite_differences"` computes the Hessian at the optimum, with finite differences that stay within the bounds (`cpm.core.optimisers.finite_difference_hessian`), once per participant for the best start, together with `log_likelihood`, `log_prior` and `metrics`: 7 instead of 22 evaluations per participant and iteration in the speed comparison, and no Hessian of zeros for an estimate on an upper bound. The default, `"numdifftools"`, keeps the previous Hessian
-- Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`. `start="prior_mean"` starts the first optimisation of every participant from the location of the current prior, clipped into the bounds, as cbm does; the other starts stay random. It does not reduce the evaluations per iteration, but the fits stop sooner (see the results below). The default, `"random"`, keeps the previous starts
-- Added a `positive_definite` column to the `fit` of `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`, which records whether each participant's Hessian was positive definite on each iteration, and a `not_positive_definite` column to their `hyperparameters`, with the number that were not. An estimate on a bound can have a Hessian that is not positive definite; it is kept, as the bounds require. In the speed comparison, this happened in about 0.1% of the participant-wise fits with `hessian="finite_differences"`, all of them on a bound
-- Added the results of the speed comparison with cbm for 100 participants and 120 trials (20 datasets, 10 fits at once on 12 cores) in `benchmarks/results/hierarchical-comparison-*.csv`. Median ms per participant and iteration, iterations and seconds per fit, for `EmpiricalBayes` / `VariationalBayes`:
-  - 0.26.0.dev1: 20.8 / 20.1 ms, 50 / 50 iterations, 91 / 88 s
-  - defaults: 17.2 / 16.4 ms, 50 / 50 iterations, 74 / 72 s, with identical results
-  - `hessian="finite_differences"`: 14.1 / 13.8 ms, 42 / 39 iterations, 59 / 53 s
-  - `start="prior_mean"`: 16.2 / 15.8 ms, 41 / 44 iterations, 64 / 65 s
-  - both: 13.0 / 13.2 ms, 38 / 25.5 iterations, 47 / 40 s
-  - both, with `number_of_starts=1`: 6.8 / 7.0 ms, 38.5 / 21.5 iterations, 23 / 18 s
-  - cbm (in the earlier run of the comparison): 4.4 ms, 21 iterations, 9.2 s
-
-  The RMSE of the population means and sds and of the individual estimates stays within the bootstrap 95% intervals of 0.26.0.dev1 in every configuration. With `hessian="finite_differences"`, `EmpiricalBayes` replaces no participant's log model evidence with -1e6, which it did in 16 of the 20 fits before. The hierarchical tutorials (`benchmarks/hierarchical.py`, RLRW) take 7.3 instead of 9.1 s (empirical Bayes) and 7.7 instead of 9.4 s (variational Bayes), with identical results
+- Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration; against cbm (RLRW, 100 participants, 120 trials), cpm makes about as many evaluations but each one costs more, at 17.3 against 5.0 ms per participant and iteration
+- Added a `hessian` argument to `cpm.optimisation.FminBound`, `Fmin`, `Minimize` and `Bads`; `hessian="finite_differences"` computes the Hessian at the optimum with bounded finite differences (`cpm.core.optimisers.finite_difference_hessian`), cutting Hessian evaluations from 22 to 7 per participant and iteration in the speed comparison. The default, `"numdifftools"`, keeps the previous Hessian
+- Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`; `start="prior_mean"` starts the first optimisation of every participant from the current prior's location, clipped into the bounds, as cbm does, so fits stop sooner. The default, `"random"`, keeps the previous starts
+- Added a `positive_definite` column to the `fit` of `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`, recording whether each participant's Hessian was positive definite, and a `not_positive_definite` column to their `hyperparameters` with the count; an estimate on a bound can have a non-positive-definite Hessian and is kept, as the bounds require
+- Added the results of the speed comparison with cbm (100 participants, 120 trials) in `benchmarks/results/hierarchical-comparison-*.csv`: combining `hessian="finite_differences"`, `start="prior_mean"` and `number_of_starts=1` brings `EmpiricalBayes`/`VariationalBayes` from 20.8/20.1 to 6.8/7.0 ms per participant and iteration, against cbm's 4.4 ms, with RMSE against 0.26.0.dev1 staying within bootstrap 95% intervals throughout
 
 ### Changed
 

From a892c877208ca80c63982eb72f1663f52d26fe5c Mon Sep 17 00:00:00 2001
From: Lenard Dome <23503601+lenarddome@users.noreply.github.com>
Date: Wed, 14 Oct 2026 10:51:44 +0200
Subject: [PATCH 9/9] chore: bump the version to 0.26.0.dev3

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 CHANGELOG.md       | 2 +-
 cpm/__version__.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9a2c543d..870098a2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,7 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
+## [0.26.0.dev3] - 2026-10-14
 
 ### Added
 
diff --git a/cpm/__version__.py b/cpm/__version__.py
index 4d482beb..1f113989 100644
--- a/cpm/__version__.py
+++ b/cpm/__version__.py
@@ -1 +1 @@
-__version__ = "0.26.0.dev2"
+__version__ = "0.26.0.dev3"