scg: make the Steihaug CG forcing term invariant to the scale of f

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

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

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

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

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]))