Profile: add multistart optimization for profile points

pypesto/profile/options.py

@@ -60,6 +60,14 @@ class ProfileOptions(dict):
     whole_path:
         Whether to profile the whole bounds or only till we get below the
         ratio.
+    profile_n_starts:
+        Number of optimization starts attempted at each profile point.
+    profile_sampling_sigma:
+        Standard deviation used for each free coordinate when sampling
+        additional optimization startpoints from a Gaussian centered at the
+        point suggested by `profile_next_guess`, in the reduced parameter
+        space. The profiled parameter itself is already fixed and is therefore
+        not perturbed.
     step_size_precheck_mode:
         Controls the step-size precheck, which estimates how many profile
         steps the resolved step sizes imply and reports suspiciously small
@@ -83,6 +91,8 @@ def __init__(
         reg_order: int = 4,
         adaptive_target_scaling_factor: float = 1.5,
         whole_path: bool = False,
+        profile_n_starts: int = 6,
+        profile_sampling_sigma: float = 0.01,
         step_size_precheck_mode: str = "warn",
         default_step_size: float | None = None,
         min_step_size: float | None = None,
@@ -130,6 +140,8 @@ def __init__(
         self.reg_order = reg_order
         self.adaptive_target_scaling_factor = adaptive_target_scaling_factor
         self.whole_path = whole_path
+        self.profile_n_starts = profile_n_starts
+        self.profile_sampling_sigma = profile_sampling_sigma
         self.step_size_precheck_mode = step_size_precheck_mode
 
         self.validate()
@@ -210,6 +222,10 @@ def validate_step_size_family(family: str) -> bool:
 
         if self.adaptive_target_scaling_factor < 1:
             raise ValueError("adaptive_target_scaling_factor must be > 1.")
+        if self.profile_n_starts < 1:
+            raise ValueError("profile_n_starts must be >= 1.")
+        if self.profile_sampling_sigma <= 0:
+            raise ValueError("profile_sampling_sigma must be > 0.")
         if self.step_size_precheck_mode not in {"off", "warn", "raise"}:
             raise ValueError(
                 "step_size_precheck_mode must be one of "

pypesto/profile/profile.py

@@ -210,6 +210,34 @@ def create_next_guess(
             *paired_profiles[p_index]
         )
 
+    # Report if profiling found a lower function value than the supplied
+    # optimum so the user notices the profiling is no longer anchored at the
+    # true best known point.
+    better_optima = []
+    for i_par, profiler_result in enumerate(result.profile_result.list[-1]):
+        if profiler_result is None:
+            continue
+        best_fval = float(np.min(profiler_result.fval_path))
+        if best_fval > global_opt or np.isclose(
+            best_fval, global_opt, rtol=1e-10, atol=1e-8
+        ):
+            continue
+        better_optima.append((problem.x_names[i_par], best_fval))
+    if better_optima:
+        par_lines = "\n".join(
+            f"  {name}: best fval = {fval:.6g}" for name, fval in better_optima
+        )
+        logger.warning(
+            "Profiling found lower function values than the supplied "
+            f"optimization result (fval = {global_opt:.6g}) for:\n"
+            f"{par_lines}\n"
+            "This means profiling was started from a suboptimal point, so the "
+            "profile ratios and confidence thresholds are not anchored to the "
+            "true global optimum. Re-optimization of the model is suggested. "
+            "The better profile points can be included as startpoints "
+            "through `x_guesses` of the optimization problem."
+        )
+
     autosave(
         filename=filename,
         result=result,

pypesto/profile/walk_along_profile.py

@@ -14,6 +14,83 @@
 logger = logging.getLogger(__name__)
 
 
+def profile_multistart_optimize(
+    optimizer: Optimizer,
+    problem: Problem,
+    startpoint: np.ndarray,
+    options: ProfileOptions,
+) -> OptimizerResult:
+    """
+    Perform optimization at one profile point using multiple starts.
+
+    Additional starts are sampled in the reduced parameter space from a
+    Gaussian centered at the point suggested by `profile_next_guess`, using
+    `options.profile_sampling_sigma` as the standard deviation for each free
+    coordinate. The original suggested startpoint is always included unchanged
+    as the final start so the helper falls back to the previous single-start
+    behavior if all sampled alternatives fail, raise, or are worse.
+
+    Parameters
+    ----------
+    optimizer:
+        The optimizer to use.
+    problem:
+        The reduced problem with the profiling parameter already fixed.
+    startpoint:
+        Optimization startpoint suggested by `profile_next_guess`, in reduced
+        space.
+    options:
+        Profile options controlling the number of starts and the Gaussian
+        sampling spread around `startpoint`.
+
+    Returns
+    -------
+    The best finite optimization result across all attempted starts, or the
+    result from the original startpoint if all sampled starts fail.
+    """
+    if options.profile_n_starts == 1:
+        return optimizer.minimize(
+            problem=problem,
+            x0=startpoint,
+            id=str(0),
+            optimize_options=OptimizeOptions(allow_failed_starts=True),
+        )
+
+    sampled_startpoints = np.random.normal(
+        loc=startpoint,
+        scale=options.profile_sampling_sigma * np.ones_like(startpoint),
+        size=(options.profile_n_starts - 1, len(startpoint)),
+    )
+    sampled_startpoints = np.clip(sampled_startpoints, problem.lb, problem.ub)
+    startpoints = np.vstack((sampled_startpoints, startpoint[np.newaxis, :]))
+
+    best_optimizer_result = None
+    best_fval = np.inf
+    original_start_result = None
+
+    for i_start, candidate_startpoint in enumerate(startpoints):
+        optimizer_result = optimizer.minimize(
+            problem=problem,
+            x0=candidate_startpoint,
+            id=str(i_start),
+            optimize_options=OptimizeOptions(allow_failed_starts=True),
+        )
+
+        if i_start == len(startpoints) - 1:
+            original_start_result = optimizer_result
+
+        if (
+            np.isfinite(optimizer_result.fval)
+            and optimizer_result.fval < best_fval
+        ):
+            best_fval = optimizer_result.fval
+            best_optimizer_result = optimizer_result
+
+    if best_optimizer_result is not None:
+        return best_optimizer_result
+    return original_start_result
+
+
 def walk_along_profile(
     current_profile: ProfilerResult,
     problem: Problem,
@@ -64,6 +141,8 @@ def walk_along_profile(
     if par_direction not in (-1, 1):
         raise AssertionError("par_direction must be -1 or 1")
 
+    # Warn at most once per non-profiled parameter in each profile half.
+    warned_at_bound: set[int] = set()
     resolved_steps = resolved_steps_by_par[i_par]
 
     # while loop for profiling (will be exited by break command)
@@ -119,13 +198,11 @@ def walk_along_profile(
             startpoint = x_next[problem.x_free_indices]
 
             if startpoint.size > 0:
-                optimizer_result = optimizer.minimize(
+                optimizer_result = profile_multistart_optimize(
+                    optimizer=optimizer,
                     problem=problem,
-                    x0=startpoint,
-                    id=str(0),
-                    optimize_options=OptimizeOptions(
-                        allow_failed_starts=False
-                    ),
+                    startpoint=startpoint,
+                    options=options,
                 )
 
                 if np.isfinite(optimizer_result.fval):
@@ -200,11 +277,11 @@ def walk_along_profile(
             problem.fix_parameters(i_par, x_next[i_par])
             startpoint = x_next[problem.x_free_indices]
 
-            optimizer_result = optimizer.minimize(
+            optimizer_result = profile_multistart_optimize(
+                optimizer=optimizer,
                 problem=problem,
-                x0=startpoint,
-                id=str(0),
-                optimize_options=OptimizeOptions(allow_failed_starts=False),
+                startpoint=startpoint,
+                options=options,
             )
 
             if np.isfinite(optimizer_result.fval):
@@ -271,6 +348,24 @@ def walk_along_profile(
             f"Optimization successful for {problem.x_names[i_par]}={x_next[i_par]:.4f}. "
             f"Start fval {problem.objective(x_next[problem.x_free_indices]):.6f}, end fval {optimizer_result.fval:.6f}."
         )
+
+        for k, j_par in enumerate(problem.x_free_indices):
+            x_j = optimizer_result.x[j_par]
+            lb_j = problem.lb[k]
+            ub_j = problem.ub[k]
+            at_lb = abs(x_j - lb_j) <= 1e-8
+            at_ub = abs(x_j - ub_j) <= 1e-8
+            if (at_lb or at_ub) and j_par not in warned_at_bound:
+                warned_at_bound.add(j_par)
+                bound_val = lb_j if at_lb else ub_j
+                logger.warning(
+                    f"Parameter '{problem.x_names[j_par]}' hit its "
+                    f"{'lower' if at_lb else 'upper'} bound "
+                    f"({bound_val:.4g}) while profiling "
+                    f"'{problem.x_names[i_par]}'. "
+                    "The profile may be constrained near this region."
+                )
+
         if optimizer_result[GRAD] is not None:
             gradnorm = np.linalg.norm(
                 optimizer_result[GRAD][problem.x_free_indices]

test/profile/test_profile.py

@@ -20,6 +20,7 @@
     resolve_profile_step_sizes,
     resolve_profile_step_sizes_for_parameters,
 )
+from pypesto.profile.walk_along_profile import profile_multistart_optimize
 
 from ..conftest import close_fig
 from ..util import rosen_for_sensi
@@ -477,6 +478,8 @@ def test_options_valid():
             "default_step_size_relative": 0.03,
             "max_step_size_relative": 0.02,
         },
+        {"profile_n_starts": 0},
+        {"profile_sampling_sigma": 0},
         {
             "default_step_size_absolute": 0.0,
             "default_step_size_relative": 0.0,
@@ -628,6 +631,70 @@ def test_profile_step_size_precheck_modes(mode, expect_warning, expect_raise):
         assert not precheck_warnings
 
 
+def test_profile_multistart_optimize_uses_best_start(monkeypatch):
+    """Multi-start profiling should tolerate failed starts and keep the best finite result."""
+
+    class DummyOptimizer:
+        def __init__(self):
+            self.calls = []
+
+        def minimize(
+            self,
+            problem,
+            x0=None,
+            id=None,
+            history_options=None,
+            optimize_options=None,
+        ):
+            del problem, history_options
+            self.calls.append(np.array(x0, copy=True))
+            if np.isclose(x0[0], 0.5):
+                if optimize_options.allow_failed_starts:
+                    # Real pyPESTO optimizers back-fill failed tolerated
+                    # starts from history, which leaves them non-finite if no
+                    # useful point was recorded before the exception.
+                    return pypesto.OptimizerResult(
+                        id=id,
+                        x0=np.array(x0, copy=True),
+                        fval=np.inf,
+                        exitflag=-1,
+                        message="sampled start failed",
+                    )
+                raise RuntimeError("sampled start failed")
+            return pypesto.OptimizerResult(
+                id=id,
+                x=np.array(x0, copy=True),
+                fval=float(np.sum(x0**2)),
+            )
+
+    problem = pypesto.Problem(
+        objective=pypesto.Objective(fun=lambda x: x[0] ** 2),
+        lb=np.array([-1.0]),
+        ub=np.array([1.0]),
+    )
+    startpoint = np.array([0.8])
+    options = profile.ProfileOptions(profile_n_starts=3)
+
+    monkeypatch.setattr(
+        np.random,
+        "normal",
+        lambda loc, scale, size: np.array([[0.5], [0.1]]),
+    )
+
+    optimizer = DummyOptimizer()
+    result = profile_multistart_optimize(
+        optimizer=optimizer,
+        problem=problem,
+        startpoint=startpoint,
+        options=options,
+    )
+
+    assert len(optimizer.calls) == options.profile_n_starts
+    assert np.allclose(optimizer.calls[-1], startpoint)
+    assert np.allclose(result.x, np.array([0.1]))
+    assert np.isclose(result.fval, 0.01)
+
+
 @pytest.mark.parametrize(
     "lb,ub",
     [(6 * np.ones(5), 10 * np.ones(5)), (-4 * np.ones(5), 1 * np.ones(5))],