From d72cca75d3c28865944b89d127b6c90a0febf070 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Sun, 25 Oct 2026 01:53:21 +0100 Subject: [PATCH] Make the Steihaug CG forcing term invariant to the scale of f (#77) The residual tolerance min(0.5, sqrt(|g|)) * |g| of the scg subproblem solver stays at 0.5 * |g| while |g| >= 0.25. On ill-conditioned problems a single CG iteration (a steepest descent step) often meets it, and the resulting small decrease in f triggered FTOL short of the optimum, depending on the scale of f. Use min(0.5, sqrt(|sg| / |sg0|)) with sg0 the scaled gradient at the first iteration instead. Early steps stay as cheap as before, the tolerance no longer depends on the scale of f and still decreases like sqrt(|sg|) towards the optimum. Only scg is affected. Fixes #77 Co-Authored-By: Claude Opus 5.5 --- fides/minimize.py | 10 ++++++++++ fides/steps.py | 44 ++++++++++++++++++++++++++++++++++++++++-- fides/trust_region.py | 7 ++++++- tests/test_minimize.py | 43 +++++++++++++++++++++++++++++++++++++++++ tests/test_steps.py | 29 ++++++++++++++++++++++++++++ 5 files changed, 130 insertions(+), 3 deletions(-) diff --git a/fides/minimize.py b/fides/minimize.py index ca31023..7d7d020 100644 --- a/fides/minimize.py +++ b/fides/minimize.py @@ -13,6 +13,7 @@ import h5py import numpy as np from numpy.linalg import norm +from scipy import linalg from scipy.sparse import csc_matrix from .constants import ( @@ -338,6 +339,8 @@ def __init__( self.x_min = self.x self.fval_min = self.fval self.grad_min = self.grad + # norm of the scaled gradient at the first iteration + self.sgnorm0: float = np.nan self.hessian_update: HessianApproximation | None = hessian_update self.iterations_since_tr_update: int = 0 @@ -482,6 +485,12 @@ def minimize( v, dv = self.get_affine_scaling() scaling = csc_matrix(np.diag(np.sqrt(np.abs(v)))) + if self.iteration == 1: + # reference for the residual tolerance of the CG subproblem + # solver, scipy's norm does not overflow for huge gradients + self.sgnorm0 = linalg.norm( + scaling.dot(self.grad), check_finite=False + ) theta = max( self.get_option(Options.THETA_MAX), 1 - norm(v * self.grad, np.inf), @@ -502,6 +511,7 @@ def minimize( subspace_dim=self.get_option(Options.SUBSPACE_DIM), stepback_strategy=self.get_option(Options.STEPBACK_STRAT), logger=self.logger, + sgnorm0=self.sgnorm0, ) # sum in the order used by the step classes to compute the diff --git a/fides/steps.py b/fides/steps.py index 59285c4..c4b6a16 100644 --- a/fides/steps.py +++ b/fides/steps.py @@ -372,10 +372,40 @@ class CGStep(Step): """ This class provides the machinery to compute an approximate solution of the trust region subproblem using the conjugate gradients methods + + :ivar sgnorm0: Norm of the rescaled gradient at the first iteration, + reference for the residual tolerance of the CG iteration """ type = 'cg' + def __init__( + self, + x, + sg, + hess, + scaling, + g_dscaling, + delta, + theta, + ub, + lb, + logger, + sgnorm0: float | None = None, + ): + """ + :param sgnorm0: + Norm of the rescaled gradient at the first iteration, defaults + to the norm of `sg`, i.e., the step is treated as the first step + """ + super().__init__( + x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb, logger + ) + # scipy's norm does not overflow for entries beyond ~1e154 + self.sgnorm0 = ( + linalg.norm(sg, check_finite=False) if sgnorm0 is None else sgnorm0 + ) + def calculate(self): # the subproblem solution is invariant to a joint rescaling of shess # and sg. Rescale such that the largest gradient entry is in @@ -384,8 +414,18 @@ def calculate(self): # the rescaling exact. scale = np.ldexp(1.0, np.frexp(np.max(np.abs(self.sg)))[1]) nsg = norm(self.sg / scale) - # residual tolerance min(0.5, sqrt(|sg|)) * |sg| in rescaled units - eps = min(0.5, np.sqrt(scale * nsg)) * nsg + # residual tolerance eta * |sg|. The textbook forcing term + # eta = min(0.5, sqrt(|sg|)) [NocedalWright2006] depends on the scale + # of the objective function and stays at 0.5 while |sg| >= 0.25, + # where a single CG iteration (a steepest descent step) often meets + # it on ill-conditioned problems. Relative to the gradient norm at + # the first iteration, eta is invariant to that scale, keeps early + # steps cheap and still decreases like sqrt(|sg|) towards the + # optimum. + eta = 0.5 + if 0 < scale * nsg < self.sgnorm0: + eta = min(eta, np.sqrt(scale * nsg / self.sgnorm0)) + eps = eta * nsg ss = self.conj_grad(self.shess / scale, self.sg / scale, eps) # represent the whole CG step as 1D subspace (with ss0 = 0), such # that step back, reflection and truncation act on the step that is diff --git a/fides/trust_region.py b/fides/trust_region.py index f4af887..4b0bc25 100644 --- a/fides/trust_region.py +++ b/fides/trust_region.py @@ -7,6 +7,7 @@ """ import logging +from functools import partial import numpy as np from scipy.sparse import csc_matrix @@ -37,6 +38,7 @@ def trust_region( subspace_dim: SubSpaceDim, stepback_strategy: StepBackStrategy, logger: logging.Logger, + sgnorm0: float | None = None, ) -> Step: """ Compute a step according to the solution of the trust-region subproblem. @@ -72,6 +74,9 @@ def trust_region( optimization boundary. :param logger: logging.Logger instance to be used for logging + :param sgnorm0: + Norm of the scaled gradient at the first iteration, reference for + the residual tolerance of the conjugate gradient subproblem solver :return: s: proposed step, @@ -83,7 +88,7 @@ def trust_region( step_options = { SubSpaceDim.TWO: TRStep2D, SubSpaceDim.FULL: TRStepFull, - SubSpaceDim.STEIHAUG: TRStepSteihaug, + SubSpaceDim.STEIHAUG: partial(TRStepSteihaug, sgnorm0=sgnorm0), } tr_step = step_options[subspace_dim]( x, sg, hess, scaling, g_dscaling, delta, theta, ub, lb, logger diff --git a/tests/test_minimize.py b/tests/test_minimize.py index 9796530..9e631d5 100644 --- a/tests/test_minimize.py +++ b/tests/test_minimize.py @@ -951,6 +951,49 @@ def timeout(signum, frame): assert np.allclose(opt.x, -50) +@pytest.mark.parametrize('fscale', [1e-4, 1.0, 1e4]) +def test_steihaug_ill_conditioned_no_premature_ftol(fscale): + """ + On this ill-conditioned quadratic (condition number 6.5e7), the forcing + term min(0.5, sqrt(|g|)) of the Steihaug CG accepted single CG + iterations (steepest descent steps) and FTOL was declared up to 0.8 % + above the optimum, depending on the scale of the objective function + (fides-dev/fides#77). + """ + hess = np.array( + [ + [15106736.52931055, 21141502.501449488], + [21141502.50144949, 29587009.60209251], + ] + ) + c = np.array([9.66093597752945, -8.568774454162156]) + lb = np.array([-0.9885675300285988, -1.8409960075420684]) + ub = np.array([2.8159846577097527, 1.6686554904476258]) + # bounded least squares solution (scipy bvls) + fopt = 46.03711167281389 + + def fun(x): + d = x - c + g = hess @ d + return fscale * 0.5 * d @ g, fscale * g, fscale * hess + + opt = Optimizer( + fun, + ub=ub, + lb=lb, + verbose=logging.WARNING, + options={ + Options.SUBSPACE_DIM: SubSpaceDim.STEIHAUG, + # absolute tolerances would make convergence depend on fscale + Options.FATOL: 0, + Options.GATOL: 0, + }, + ) + fval, *_ = opt.minimize(lb.copy()) + assert opt.exitflag == fides.ExitFlag.FTOL + assert fval / fscale == pytest.approx(fopt, rel=1e-7) + + def test_gtol_not_triggered_by_rejected_step(): """ A rejected step must not trigger gradient based convergence: the diff --git a/tests/test_steps.py b/tests/test_steps.py index 68f0832..a3f433c 100644 --- a/tests/test_steps.py +++ b/tests/test_steps.py @@ -562,3 +562,32 @@ def test_steihaug_non_finite_curvature(g): step.calculate() assert np.all(np.isfinite(step.ss)) assert norm(step.ss) <= 1.0 + + +def test_steihaug_tolerance_invariant_to_objective_scale(): + """ + The residual tolerance of the CG iteration must not depend on the scale + of the objective function (fides-dev/fides#77). With the forcing term + min(0.5, sqrt(|g|)), a single CG iteration (a steepest descent step) was + accepted for large gradients, whereas small gradients yielded the + Newton step. + """ + hess = np.diag([1e6, 1.0]) + # one steepest descent step reduces the residual to ~0.1 * |g| + g = hess @ np.array([1.0, 1e5]) + x = np.zeros(2) + lb, ub = np.full(2, -np.inf), np.full(2, np.inf) + steps = [] + for scale in 2.0 ** np.array([-40, 0, 40]): + # gradient decreased by 1e4 since the first iteration, eta = 0.01 + step = TRStepSteihaug( + *step_args(x, scale * g, scale * hess, lb, ub, delta=1e6), + sgnorm0=scale * 1e4 * norm(g), + ) + step.calculate() + steps.append(step.ss) + # power of two scales keep the computation exact + for ss in steps: + assert np.array_equal(ss, steps[0]) + # interior Newton step + assert steps[0] == pytest.approx(-np.array([1.0, 1e5]))