diff --git a/src/tranquilo/acceptance_decision.py b/src/tranquilo/acceptance_decision.py index c5c285cc..62b8a949 100644 --- a/src/tranquilo/acceptance_decision.py +++ b/src/tranquilo/acceptance_decision.py @@ -209,17 +209,29 @@ def accept_classic_line_search( is_accepted = overall_improvement >= min_improvement # ================================================================================== - # Return results + # Process and return results + + if np.isfinite(candidate_fval): + res = _get_acceptance_result( + candidate_x=candidate_x, + candidate_fval=candidate_fval, + candidate_index=candidate_index, + rho=rho, + is_accepted=is_accepted, + old_state=state, + n_evals=1, + ) + else: + res = _get_acceptance_result( + candidate_x=state.x, + candidate_fval=state.fval, + candidate_index=state.index, + rho=-np.inf, + is_accepted=False, + old_state=state, + n_evals=1, + ) - res = _get_acceptance_result( - candidate_x=candidate_x, - candidate_fval=candidate_fval, - candidate_index=candidate_index, - rho=rho, - is_accepted=is_accepted, - old_state=state, - n_evals=1, - ) return res @@ -262,15 +274,26 @@ def _accept_simple( is_accepted = actual_improvement >= min_improvement - res = _get_acceptance_result( - candidate_x=candidate_x, - candidate_fval=candidate_fval, - candidate_index=candidate_index, - rho=rho, - is_accepted=is_accepted, - old_state=state, - n_evals=n_evals, - ) + if np.isfinite(candidate_fval): + res = _get_acceptance_result( + candidate_x=candidate_x, + candidate_fval=candidate_fval, + candidate_index=candidate_index, + rho=rho, + is_accepted=is_accepted, + old_state=state, + n_evals=n_evals, + ) + else: + res = _get_acceptance_result( + candidate_x=state.x, + candidate_fval=state.fval, + candidate_index=state.index, + rho=-np.inf, + is_accepted=False, + old_state=state, + n_evals=n_evals, + ) return res @@ -321,15 +344,26 @@ def accept_noisy( is_accepted = actual_improvement >= min_improvement - res = _get_acceptance_result( - candidate_x=candidate_x, - candidate_fval=candidate_fval, - candidate_index=candidate_index, - rho=rho, - is_accepted=is_accepted, - old_state=state, - n_evals=n_2, - ) + if np.isfinite(candidate_fval): + res = _get_acceptance_result( + candidate_x=candidate_x, + candidate_fval=candidate_fval, + candidate_index=candidate_index, + rho=rho, + is_accepted=is_accepted, + old_state=state, + n_evals=n_2, + ) + else: + res = _get_acceptance_result( + candidate_x=state.x, + candidate_fval=state.fval, + candidate_index=state.index, + rho=-np.inf, + is_accepted=False, + old_state=state, + n_evals=n_2, + ) return res diff --git a/src/tranquilo/process_arguments.py b/src/tranquilo/process_arguments.py index d8c14d04..27553331 100644 --- a/src/tranquilo/process_arguments.py +++ b/src/tranquilo/process_arguments.py @@ -79,6 +79,7 @@ def process_arguments( model_fitter_options=None, cube_subsolver="bntr_fast", sphere_subsolver="gqtpar_fast", + retry_subproblem_with_fallback=True, subsolver_options=None, acceptance_decider=None, acceptance_decider_options=None, @@ -189,6 +190,7 @@ def process_arguments( solve_subproblem = get_subsolver( cube_solver=cube_subsolver, sphere_solver=sphere_subsolver, + retry_with_fallback=retry_subproblem_with_fallback, user_options=subsolver_options, ) diff --git a/src/tranquilo/solve_subproblem.py b/src/tranquilo/solve_subproblem.py index 645417b8..f08ed813 100644 --- a/src/tranquilo/solve_subproblem.py +++ b/src/tranquilo/solve_subproblem.py @@ -14,14 +14,17 @@ gqtpar, ) from tranquilo.subsolvers.gqtpar_fast import gqtpar_fast -from tranquilo.subsolvers.wrapped_subsolvers import ( - slsqp_sphere, - solve_multistart, -) from tranquilo.options import SubsolverOptions +from tranquilo.subsolvers.fallback_subsolvers import ( + robust_cube_solver, + robust_sphere_solver_inscribed_cube, + robust_sphere_solver_norm_constraint, + robust_sphere_solver_reparametrized, + robust_cube_solver_multistart, +) -def get_subsolver(sphere_solver, cube_solver, user_options=None): +def get_subsolver(sphere_solver, cube_solver, retry_with_fallback, user_options=None): """Get an algorithm-function with partialled options. Args: @@ -36,6 +39,8 @@ def get_subsolver(sphere_solver, cube_solver, user_options=None): to be ``lower_bounds`` and ``upper_bounds``. The fourth argument needs to be ``x_candidate``, an initial guess for the solution in the unit space. Moreover, subsolvers can have any number of additional keyword arguments. + retry_with_fallback (bool): Whether to retry solving the subproblem with a + fallback solver if the optimized subsolver raises an exception. user_options (dict): Options for the subproblem solver. The following are supported: - maxiter (int): Maximum number of iterations to perform when solving the @@ -70,13 +75,16 @@ def get_subsolver(sphere_solver, cube_solver, user_options=None): built_in_sphere_solvers = { "gqtpar": gqtpar, "gqtpar_fast": gqtpar_fast, - "slsqp_sphere": slsqp_sphere, + "fallback_reparametrized": robust_sphere_solver_reparametrized, + "fallback_inscribed_cube": robust_sphere_solver_inscribed_cube, + "fallback_norm_constraint": robust_sphere_solver_norm_constraint, } built_in_cube_solvers = { "bntr": bntr, "bntr_fast": bntr_fast, - "multistart": solve_multistart, + "fallback_cube": robust_cube_solver, + "fallback_multistart": robust_cube_solver_multistart, } _sphere_subsolver = get_component( @@ -97,10 +105,29 @@ def get_subsolver(sphere_solver, cube_solver, user_options=None): mandatory_signature=["model", "x_candidate", "lower_bounds", "upper_bounds"], ) + _fallback_sphere_solver = get_component( + name_or_func="fallback_inscribed_cube", + component_name="fallback_sphere_solver", + func_dict=built_in_sphere_solvers, + default_options=SubsolverOptions(), + mandatory_signature=["model", "x_candidate"], + ) + + _fallback_cube_solver = get_component( + name_or_func="fallback_cube", + component_name="fallback_cube_solver", + func_dict=built_in_cube_solvers, + default_options=SubsolverOptions(), + mandatory_signature=["model", "x_candidate"], + ) + solver = partial( _solve_subproblem_template, sphere_solver=_sphere_subsolver, cube_solver=_cube_subsolver, + fallback_sphere_solver=_fallback_sphere_solver, + fallback_cube_solver=_fallback_cube_solver, + retry_with_fallback=retry_with_fallback, ) return solver @@ -111,6 +138,9 @@ def _solve_subproblem_template( trustregion, sphere_solver, cube_solver, + fallback_sphere_solver, + fallback_cube_solver, + retry_with_fallback, ): """Solve the quadratic subproblem. @@ -128,6 +158,8 @@ def _solve_subproblem_template( ``upper_bounds``. The fourth argument needs to be ``x_candidate``, an initial guess for the solution in the unit space. Moreover, subsolvers can have any number of additional keyword arguments. + retry_with_fallback (bool): Whether to retry solving the subproblem with a + fallback solver if the optimized subsolver raises an exception. Returns: @@ -145,16 +177,32 @@ def _solve_subproblem_template( """ old_x_unit = trustregion.map_to_unit(trustregion.center) - solver = sphere_solver if trustregion.shape == "sphere" else cube_solver - - raw_result = solver( - model=model, - x_candidate=old_x_unit, - # bounds can be passed to both solvers because the functions returned by - # `get_component` ignore redundant arguments. - lower_bounds=-np.ones_like(old_x_unit), - upper_bounds=np.ones_like(old_x_unit), - ) + if trustregion.shape == "sphere": + solver = sphere_solver + fallback_solver = fallback_sphere_solver + else: + solver = cube_solver + fallback_solver = fallback_cube_solver + + # try finding a solution using the optimized subsolver. If this raises an + # exception, use a fallback solver if requested, otherwise raise the exception. + try: + raw_result = solver( + model=model, + x_candidate=old_x_unit, + # bounds can be passed to both solvers because the functions returned by + # `get_component` ignore redundant arguments. + lower_bounds=-np.ones_like(old_x_unit), + upper_bounds=np.ones_like(old_x_unit), + ) + except Exception as e: + if not retry_with_fallback: + raise e + + raw_result = fallback_solver( + model=model, + x_candidate=old_x_unit, + ) if trustregion.shape == "cube": raw_result["x"] = np.clip(raw_result["x"], -1.0, 1.0) diff --git a/src/tranquilo/subsolvers/bntr_fast.py b/src/tranquilo/subsolvers/bntr_fast.py index 36c54a45..b02cc196 100644 --- a/src/tranquilo/subsolvers/bntr_fast.py +++ b/src/tranquilo/subsolvers/bntr_fast.py @@ -512,7 +512,7 @@ def _take_preliminary_gradient_descent_step_and_check_for_solution( if not converged: trustregion_radius = min( - max(min_radius, max(trustregion_radius, radius_lower_bound)), max_radius + max(min_radius, trustregion_radius, radius_lower_bound), max_radius ) return ( diff --git a/src/tranquilo/subsolvers/fallback_subsolvers.py b/src/tranquilo/subsolvers/fallback_subsolvers.py new file mode 100644 index 00000000..5192aad8 --- /dev/null +++ b/src/tranquilo/subsolvers/fallback_subsolvers.py @@ -0,0 +1,199 @@ +import numpy as np +from functools import partial +from scipy.optimize import Bounds, NonlinearConstraint, minimize + +from tranquilo.exploration_sample import draw_exploration_sample + + +# ====================================================================================== +# Cube solvers +# ====================================================================================== +def robust_cube_solver(model, x_candidate, radius=1.0): + """Robust cube solver. + + Argument `radius` corresponds to half of the side length of the cube. + + """ + crit, grad = _get_crit_and_grad(model) + + lower_bounds = -radius * np.ones(len(x_candidate)) + upper_bounds = radius * np.ones(len(x_candidate)) + + res = minimize( + crit, + x_candidate, + method="L-BFGS-B", + jac=grad, + bounds=Bounds(lower_bounds, upper_bounds), + options={ + "maxiter": len(x_candidate), + }, + ) + + return { + "x": res.x, + "criterion": res.fun, + "n_iterations": res.nit, + "success": True, + } + + +def robust_cube_solver_multistart(model, x_candidate): + np.random.seed(12345) + start_values = draw_exploration_sample( + x=x_candidate, + lower=-np.ones(len(x_candidate)), + upper=np.ones(len(x_candidate)), + n_samples=100, + sampling_distribution="uniform", + sampling_method="sobol", + seed=1234, + ) + + best_crit = np.inf + accepted_x = None + critvals = [] + + for x in start_values: + res = robust_cube_solver(model, x) + + if res["criterion"] <= best_crit: + accepted_x = res["x"] + best_crit = res["criterion"] + + critvals.append(res["criterion"]) + + return { + "x": accepted_x, + "criterion": best_crit, + "std": np.std(critvals), + "n_iterations": None, + "success": None, + } + + +# ====================================================================================== +# Sphere solvers +# ====================================================================================== + + +def robust_sphere_solver_inscribed_cube(model, x_candidate): + """Robust sphere solver that uses a cube solver in an inscribed cube. + + We let x be in the largest cube that is inscribed inside the unit sphere. Formula + is taken from http://tinyurl.com/4astpuwn. + + This solver cannot find solutions on the hull of the sphere. + + """ + return robust_cube_solver(model, x_candidate, radius=1 / np.sqrt(len(x_candidate))) + + +def robust_sphere_solver_reparametrized(model, x_candidate): + """Robust sphere solver that uses reparametrization. + + We let x be in the cube -1 <= x <= 1, but if the optimizer chooses a point outside + the sphere x is projected onto the sphere inside the criterion function. + + This solver can find solutions on the hull of the sphere. + + """ + + def crit(x): + x_norm = np.linalg.norm(x) + if x_norm <= 1: + x_tilde = x + else: + x_tilde = x / x_norm + return model.predict(x_tilde) + + lower_bounds = -np.ones(len(x_candidate)) + upper_bounds = np.ones(len(x_candidate)) + + res = minimize( + crit, + x_candidate, + method="L-BFGS-B", + bounds=Bounds(lower_bounds, upper_bounds), + options={ + "maxiter": len(x_candidate), + }, + ) + + solution_norm = np.linalg.norm(res.x) + + if solution_norm <= 1: + _minimizer = res.x + else: + _minimizer = res.x / solution_norm + + return { + "x": _minimizer, + "criterion": res.fun, + "n_iterations": res.nit, + "success": True, + } + + +def robust_sphere_solver_norm_constraint(model, x_candidate): + """Robust sphere solver that uses ||x|| <= 1 as a nonlinear constraint. + + This solver can find solutions on the hull of the sphere. + + """ + crit, grad = _get_crit_and_grad(model) + constraint = _get_constraint() + + lower_bounds = -np.ones(len(x_candidate)) + upper_bounds = np.ones(len(x_candidate)) + + res = minimize( + crit, + x_candidate, + method="SLSQP", + bounds=Bounds(lower_bounds, upper_bounds), + jac=grad, + constraints=constraint, + options={"maxiter": 3 * len(x_candidate)}, + ) + + return { + "x": res.x, + "criterion": res.fun, + "success": res.success, + "n_iterations": res.nit, + } + + +# ====================================================================================== +# Criterion, gradient, and spherical constraint +# ====================================================================================== + + +def _get_crit_and_grad(model): + def _crit(x, c, g, h): + return c + x @ g + 0.5 * x @ h @ x + + def _grad(x, g, h): + return g + x @ h + + crit = partial(_crit, c=model.intercept, g=model.linear_terms, h=model.square_terms) + grad = partial(_grad, g=model.linear_terms, h=model.square_terms) + + return crit, grad + + +def _get_constraint(): + def _constr_fun(x): + return x @ x + + def _constr_jac(x): + return 2 * x + + return NonlinearConstraint( + fun=_constr_fun, + lb=-np.inf, + ub=1, + jac=_constr_jac, + keep_feasible=True, + ) diff --git a/src/tranquilo/subsolvers/wrapped_subsolvers.py b/src/tranquilo/subsolvers/wrapped_subsolvers.py deleted file mode 100644 index 1443a5f6..00000000 --- a/src/tranquilo/subsolvers/wrapped_subsolvers.py +++ /dev/null @@ -1,94 +0,0 @@ -from functools import partial - -import numpy as np -from scipy.optimize import Bounds, NonlinearConstraint, minimize - -from tranquilo.exploration_sample import draw_exploration_sample - - -def solve_multistart(model, x_candidate, lower_bounds, upper_bounds): - np.random.seed(12345) - start_values = draw_exploration_sample( - x=x_candidate, - lower=lower_bounds, - upper=upper_bounds, - n_samples=100, - sampling_distribution="uniform", - sampling_method="sobol", - seed=1234, - ) - - def crit(x): - return model.predict(x) - - bounds = Bounds(lower_bounds, upper_bounds) - - best_crit = np.inf - accepted_x = None - critvals = [] - for x in start_values: - res = minimize( - crit, - x, - method="L-BFGS-B", - bounds=bounds, - ) - if res.fun <= best_crit: - accepted_x = res.x - critvals.append(res.fun) - - return { - "x": accepted_x, - "std": np.std(critvals), - "n_iterations": None, - "success": None, - } - - -def slsqp_sphere(model, x_candidate): - crit, grad = get_crit_and_grad(model) - constraints = get_constraints() - - res = minimize( - crit, - x_candidate, - method="slsqp", - jac=grad, - constraints=constraints, - ) - - return { - "x": res.x, - "success": res.success, - "n_iterations": res.nit, - } - - -def get_crit_and_grad(model): - def _crit(x, c, g, h): - return c + x @ g + 0.5 * x @ h @ x - - def _grad(x, g, h): - return g + x @ h - - crit = partial(_crit, c=model.intercept, g=model.linear_terms, h=model.square_terms) - grad = partial(_grad, g=model.linear_terms, h=model.square_terms) - - return crit, grad - - -def get_constraints(): - def _constr_fun(x): - return x @ x - - def _constr_jac(x): - return 2 * x - - constr = NonlinearConstraint( - fun=_constr_fun, - lb=-np.inf, - ub=1, - jac=_constr_jac, - ) - - return (constr,) diff --git a/src/tranquilo/wrap_criterion.py b/src/tranquilo/wrap_criterion.py index 705bf697..13765d54 100644 --- a/src/tranquilo/wrap_criterion.py +++ b/src/tranquilo/wrap_criterion.py @@ -49,13 +49,20 @@ def wrapper_criterion(eval_info): criterion, arguments=arguments, n_cores=effective_n_cores, + error_handling="continue", ) + # The batch evaluator replaces exceptions with their traceback (str) when + # error_handling="continue". We replace these cases with infinity. + raw_evals_with_replaced_traceback = [ + np.inf if isinstance(x, str) else x for x in raw_evals + ] + # replace NaNs but keep infinite values. NaNs would be problematic in many # places, infs are only a problem in model fitting and will be handled there clipped_evals = [ np.nan_to_num(critval, nan=np.inf, posinf=np.inf, neginf=-np.inf) - for critval in raw_evals + for critval in raw_evals_with_replaced_traceback ] history.add_evals( diff --git a/tests/subsolvers/test_fallback_solvers.py b/tests/subsolvers/test_fallback_solvers.py new file mode 100644 index 00000000..68175a35 --- /dev/null +++ b/tests/subsolvers/test_fallback_solvers.py @@ -0,0 +1,58 @@ +import pytest +import numpy as np +from numpy.testing import assert_array_almost_equal as aaae +from tranquilo.models import ScalarModel + +from tranquilo.subsolvers.fallback_subsolvers import ( + robust_cube_solver, + robust_cube_solver_multistart, + robust_sphere_solver_inscribed_cube, + robust_sphere_solver_norm_constraint, + robust_sphere_solver_reparametrized, +) + +FALLBACK_SPHERE_SOLVERS = [ + robust_sphere_solver_inscribed_cube, + robust_sphere_solver_reparametrized, + robust_sphere_solver_norm_constraint, +] + +FALLBACK_CUBE_SOLVERS = [ + robust_cube_solver, + robust_cube_solver_multistart, +] + + +@pytest.fixture +def model(): + """Simple quadratic scalar model. + + - Minimum in cube is at (1, 1) + - Minimum in sphere is at (1/sqrt(2), 1/sqrt(2)) + + """ + return ScalarModel( + intercept=0.0, + linear_terms=-np.ones(2), + square_terms=-np.eye(2), + ) + + +@pytest.mark.parametrize("solver", FALLBACK_SPHERE_SOLVERS) +def test_fallback_sphere_solver(solver, model): + x_candidate = np.zeros(2) + + calculated = solver(model, x_candidate) + expected = np.ones(2) / np.sqrt(2) + + aaae(calculated["x"], expected) + + +@pytest.mark.parametrize("solver", FALLBACK_CUBE_SOLVERS) +def test_fallback_cube_solver(solver, model): + x_candidate = np.zeros(2) + + calculated = solver(model, x_candidate) + expected = np.ones(2) + + aaae(calculated["x"], expected) diff --git a/tests/test_poisedness.py b/tests/test_poisedness.py index ffc466a0..1ec4e89b 100644 --- a/tests/test_poisedness.py +++ b/tests/test_poisedness.py @@ -1,5 +1,6 @@ import numpy as np import pytest +import sys from tranquilo.poisedness import ( _get_minimize_options, _lagrange_poly_matrix, @@ -111,6 +112,7 @@ def evaluate_scalar_model(x, intercept, linear_terms, square_terms): ] +@pytest.mark.skipif(sys.platform == "win32", reason="Test is inaccurate on Windows.") @pytest.mark.parametrize("sample, shape, maxiter, expected", TEST_CASES) def test_improve_poisedness(sample, shape, maxiter, expected): _, got_lambdas = improve_poisedness(sample=sample, shape=shape, maxiter=maxiter) diff --git a/tests/test_rho_noise.py b/tests/test_rho_noise.py index 4cc261d3..70a58f28 100644 --- a/tests/test_rho_noise.py +++ b/tests/test_rho_noise.py @@ -59,7 +59,9 @@ def test_convergence_to_one_if_noise_is_tiny(functype): xs, fvecs, weights=None, region=trustregion, old_model=None ) - subsolver = get_subsolver(sphere_solver="gqtpar", cube_solver="bntr") + subsolver = get_subsolver( + sphere_solver="gqtpar", cube_solver="bntr", retry_with_fallback=False + ) rng = np.random.default_rng(123) diff --git a/tests/test_solve_subproblem.py b/tests/test_solve_subproblem.py index 482c88f3..28fb3da1 100644 --- a/tests/test_solve_subproblem.py +++ b/tests/test_solve_subproblem.py @@ -35,7 +35,9 @@ def test_without_bounds(solver_name): trustregion = Region(center=np.zeros(3), radius=1, bounds=Bounds(None, None)) - solve_subproblem = get_subsolver(sphere_solver=solver_name, cube_solver="bntr") + solve_subproblem = get_subsolver( + sphere_solver=solver_name, cube_solver="bntr", retry_with_fallback=False + ) calculated = solve_subproblem( model=model,