From 830c407fd7df0a0c40711103112116df423ea9fc Mon Sep 17 00:00:00 2001 From: timmens Date: Wed, 30 Sep 2026 11:33:13 +0200 Subject: [PATCH] Fix gradient norm scaling in pounders convergence checks (#657) The residual model is fitted in coordinates scaled by the trust-region radius, so main_model.linear_terms already equals delta times the gradient. The convergence checks multiplied this by delta again, so the gtol_abs, gtol_rel and gtol_scaled tests compared delta**2 * ||grad|| against the tolerances. Once the trust region shrank, runs could report gradient convergence at points far from a minimum. Divide the model gradient by delta instead, so the checks use the gradient in the original coordinates, as in Stefan Wild's reference POUNDERs implementation (IBCDFO formquad.m returns G / delta). Remove the Linux and Windows skips on test_bntr, assert the final criterion value, and add a regression test with jittered start values. --- src/optimagic/optimizers/pounders.py | 8 +-- .../optimizers/test_pounders_integration.py | 64 +++++++++++++++---- 2 files changed, 56 insertions(+), 16 deletions(-) diff --git a/src/optimagic/optimizers/pounders.py b/src/optimagic/optimizers/pounders.py index 2fa2650b2..64b8530a4 100644 --- a/src/optimagic/optimizers/pounders.py +++ b/src/optimagic/optimizers/pounders.py @@ -318,8 +318,9 @@ def internal_solve_pounders( ) x_accepted = history.get_best_x() - gradient_norm_initial = np.linalg.norm(main_model.linear_terms) - gradient_norm_initial *= delta + # The model is fitted in coordinates scaled by the trust-region radius, so its + # linear terms equal delta times the gradient in the original coordinates. + gradient_norm_initial = np.linalg.norm(main_model.linear_terms) / delta valid = True n_modelpoints = n + 1 @@ -548,8 +549,7 @@ def internal_solve_pounders( main_model = create_main_from_residual_model(residual_model) - gradient_norm = np.linalg.norm(main_model.linear_terms) - gradient_norm *= delta + gradient_norm = np.linalg.norm(main_model.linear_terms) / delta ( last_model_indices, diff --git a/tests/optimagic/optimizers/test_pounders_integration.py b/tests/optimagic/optimizers/test_pounders_integration.py index 72b2750c6..b3a73695c 100644 --- a/tests/optimagic/optimizers/test_pounders_integration.py +++ b/tests/optimagic/optimizers/test_pounders_integration.py @@ -1,6 +1,5 @@ """Test suite for the internal pounders interface.""" -import sys from functools import partial from itertools import product @@ -90,13 +89,11 @@ def trustregion_subproblem_options(): TEST_CASES = universal_tests + specific_tests -@pytest.mark.skipif(sys.platform == "win32", reason="Not accurate on Windows.") -@pytest.mark.skipif( - sys.platform == "linux" and sys.version_info[:2] >= (3, 10), - reason="Not accurate on Linux with Python 3.10 or higher.", -) -@pytest.mark.parametrize("start_vec, conjugate_gradient_method_sub", TEST_CASES) -def test_bntr( +X_EXPECTED = np.array([0.1902789114691, 0.006131410288292, 0.01053088353832]) +CRITERION_EXPECTED = 2384.477 + + +def _solve_bntr( start_vec, conjugate_gradient_method_sub, criterion, @@ -136,9 +133,53 @@ def batch_fun(x_list, n_cores): batch_fun=batch_fun, **pounders_options, ) + return result + + +@pytest.mark.parametrize("start_vec, conjugate_gradient_method_sub", TEST_CASES) +def test_bntr( + start_vec, + conjugate_gradient_method_sub, + criterion, + pounders_options, + trustregion_subproblem_options, +): + result = _solve_bntr( + start_vec, + conjugate_gradient_method_sub, + criterion, + pounders_options, + trustregion_subproblem_options, + ) + + aaae(result.x, X_EXPECTED, decimal=3) + assert np.sum(criterion(result.x) ** 2) == pytest.approx( + CRITERION_EXPECTED, abs=1e-2 + ) + + +@pytest.mark.slow +def test_bntr_no_false_convergence_under_start_value_jitter( + criterion, pounders_options, trustregion_subproblem_options +): + """Tiny perturbations of x0 must not lead to convergence at a non-minimum. + + The path from this start vector is chaotic, so relative perturbations of order + 1e-15 lead to different iterates. With a mis-scaled gradient norm some of these + runs used to report convergence far away from the minimum (see #657). + + """ + rng = np.random.default_rng(0) + start_vec = np.array([1e-6, 1e-6, 1e-6]) + criterion_values = [] + for _ in range(10): + x0 = start_vec * (1 + 1e-15 * rng.standard_normal(3)) + result = _solve_bntr( + x0, "cg", criterion, pounders_options, trustregion_subproblem_options + ) + criterion_values.append(np.sum(criterion(result.x) ** 2)) - x_expected = np.array([0.1902789114691, 0.006131410288292, 0.01053088353832]) - aaae(result.x, x_expected, decimal=3) + np.testing.assert_allclose(criterion_values, CRITERION_EXPECTED, atol=1e-2) @pytest.mark.parametrize("start_vec", [(np.array([0.15, 0.008, 0.01]))]) @@ -177,5 +218,4 @@ def batch_fun(x_list, n_cores): **pounders_options, ) - x_expected = np.array([0.1902789114691, 0.006131410288292, 0.01053088353832]) - aaae(result.x, x_expected, decimal=4) + aaae(result.x, X_EXPECTED, decimal=4)