diff --git a/pypesto/profile/options.py b/pypesto/profile/options.py index 11bcfb87a..288e3bd6c 100644 --- a/pypesto/profile/options.py +++ b/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 " diff --git a/pypesto/profile/profile.py b/pypesto/profile/profile.py index 9b7934aca..37969a2cd 100644 --- a/pypesto/profile/profile.py +++ b/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, diff --git a/pypesto/profile/walk_along_profile.py b/pypesto/profile/walk_along_profile.py index be85b0523..390ae69bb 100644 --- a/pypesto/profile/walk_along_profile.py +++ b/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] diff --git a/test/profile/test_profile.py b/test/profile/test_profile.py index 5de231925..ae69d153b 100644 --- a/test/profile/test_profile.py +++ b/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))],