Source code for cil.optimisation.utilities.StepSizeMethods

#  Copyright 2024 United Kingdom Research and Innovation
#  Copyright 2024 The University of Manchester
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#
# Authors:
# - CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt

from abc import ABC, abstractmethod
import numpy
from numbers import Number
import logging

log = logging.getLogger(__name__)

[docs] class StepSizeRule(ABC): """ Abstract base class for a step size rule. The abstract method, `get_step_size` takes in an algorithm and thus can access all parts of the algorithm (e.g. current iterate, current gradient, objective functions etc) and from this should return a float as a step size. """ def __init__(self): '''Initialises the step size rule ''' pass
[docs] @abstractmethod def get_step_size(self, algorithm): """ Returns -------- the calculated step size:float """ pass
[docs] class ConstantStepSize(StepSizeRule): """ Step-size rule that always returns a constant step-size. Parameters ---------- step_size: float The step-size to be returned with each call. """ def __init__(self, step_size): '''Initialises the constant step size rule Parameters: ------------- step_size : float, the constant step size ''' self.step_size = step_size
[docs] def get_step_size(self, algorithm): """ Returns -------- the calculated step size:float """ return self.step_size
[docs] class ArmijoStepSizeRule(StepSizeRule): r""" Applies the Armijo rule to calculate the step size (step_size). The Armijo rule runs a while loop to find the appropriate step_size by starting from a very large number (`alpha`). The step_size is found by reducing the step size (by a factor `beta`) in an iterative way until a certain criterion is met. To avoid infinite loops, we add a maximum number of times (`max_iterations`) the while loop is run. Parameters ---------- alpha: float, optional, default=1e6 The starting point for the step size iterations beta: float between 0 and 1, optional, default=0.5 The amount the step_size is reduced if the criterion is not met max_iterations: integer, optional, default is numpy.ceil (2 * numpy.log10(alpha) / numpy.log10(2)) The maximum number of iterations to find a suitable step size warmstart: Boolean, default is True If `warmstart = True` the initial step size at each Armijo iteration is the calculated step size from the last iteration. If `warmstart = False` at each Armijo iteration, the initial step size is reset to the original, large `alpha`. In the case of *well-behaved* convex functions, `warmstart = True` is likely to be computationally less expensive. In the case of non-convex functions, or particularly tricky functions, setting `warmstart = False` may be beneficial. Reference ------------ - Algorithm 3.1 in Nocedal, J. and Wright, S.J. eds., 1999. Numerical optimization. New York, NY: Springer New York. https://www.math.uci.edu/~qnie/Publications/NumericalOptimization.pdf) - https://projecteuclid.org/download/pdf_1/euclid.pjm/1102995080 """ def __init__(self, alpha=1e6, beta=0.5, max_iterations=None, warmstart=True): '''Initialises the step size rule ''' self.alpha_orig = alpha if self.alpha_orig is None: # Can be removed when alpha and beta are deprecated in GD self.alpha_orig = 1e6 self.alpha = self.alpha_orig self.beta = beta if self.beta is None: # Can be removed when alpha and beta are deprecated in GD self.beta = 0.5 self.max_iterations = max_iterations if self.max_iterations is None: self.max_iterations = numpy.ceil(2 * numpy.log10(self.alpha_orig) / numpy.log10(2)) self.warmstart=warmstart
[docs] def get_step_size(self, algorithm): """ Applies the Armijo rule to calculate the step size (`step_size`) Returns -------- the calculated step size:float """ k = 0 if not self.warmstart: self.alpha = self.alpha_orig f_x = algorithm.calculate_objective_function_at_point(algorithm.solution) self.x_armijo = algorithm.solution.copy() log.debug("Starting Armijo backtracking with initial step size: %f", self.alpha) while k < self.max_iterations: algorithm.gradient_update.multiply(self.alpha, out=self.x_armijo) algorithm.solution.subtract(self.x_armijo, out=self.x_armijo) f_x_a = algorithm.calculate_objective_function_at_point(self.x_armijo) sqnorm = algorithm.gradient_update.squared_norm() if f_x_a - f_x <= - (self.alpha/2.) * sqnorm: break k += 1. self.alpha *= self.beta log.info("Armijo rule took %d iterations to find step size", k) if k == self.max_iterations: raise ValueError( 'Could not find a proper step_size in {} loops. Consider increasing alpha or max_iterations.'.format(self.max_iterations)) return self.alpha
[docs] class BarzilaiBorweinStepSizeRule(StepSizeRule): r""" Applies the Barzilai- Borwein rule to calculate the step size (step_size). Let :math:`\Delta x=x_k-x_{k-1}` and :math:`\Delta g=g_k-g_{k-1}`. Where :math:`x_k` is the :math:`k` th iterate (current solution after iteration :math:`k` ) and :math:`g_k` is the gradient calculation in the :math:`k` th iterate, found in :code:`algorithm.gradient_update`. A Barzilai-Borwein (BB) iteration is :math:`x_{k+1}=x_k-\alpha_kg_k` where the step size :math:`\alpha _k` is either - :math:`\alpha_k^{LONG}=\frac{\Delta x\cdot\Delta x}{\Delta x\cdot\Delta g}`, or - :math:`\alpha_k^{SHORT}=\frac{\Delta x \cdot\Delta g}{\Delta g \cdot\Delta g}`. Where the operator :math:`\cdot` is the standard inner product between two vectors. This is suitable for use with gradient based iterative methods where the calculated gradient is stored as `algorithm.gradient_update`. Parameters ---------- initial: float, greater than zero The step-size for the first iteration. We recommend something of the order :math:`1/f.L` where :math:`f` is the (differentiable part of) the objective you wish to minimise. mode: One of 'long', 'short' or 'alternate', default is 'short'. This calculates the step-size based on the LONG, SHORT or alternating between the two, starting with short. stabilisation_param: 'auto', float or 'off', default is 'auto' In order to add stability the step-size has an upper limit of :math:`\Delta/\|g_k\|` where by 'default', the `stabilisation_param`, :math:`\Delta` is determined automatically to be the minimium of :math:`\Delta x` from the first 3 iterations. The user can also pass a fixed constant or turn "off" the stabilisation, equivalently passing `np.inf`. Reference --------- - Barzilai, Jonathan; Borwein, Jonathan M. (1988). "Two-Point Step Size Gradient Methods". IMA Journal of Numerical Analysis. 8: 141–148, https://doi.org/10.1093/imanum/8.1.141 - Burdakov, O., Dai, Y. and Huang, N., 2019. STABILIZED BARZILAI-BORWEIN METHOD. Journal of Computational Mathematics, 37(6). https://doi.org/10.4208/jcm.1911-m2019-0171 - https://en.wikipedia.org/wiki/Barzilai-Borwein_method """ def __init__(self, initial, mode='short', stabilisation_param="auto"): '''Initialises the step size rule ''' self.mode=mode if self.mode == 'short': self.is_short = True elif self.mode == 'long' or self.mode == 'alternate': self.is_short = False else: raise ValueError('Mode should be chosen from "long", "short" or "alternate". ') self.store_grad=None self.store_x=None self.initial=initial if stabilisation_param == 'auto': self.adaptive = True stabilisation_param = numpy.inf elif stabilisation_param == "off": self.adaptive = False stabilisation_param = numpy.inf elif ( isinstance(stabilisation_param, Number) and stabilisation_param >=0): self.adaptive = False else: raise TypeError(" The stabilisation_param should be 'auto', a positive number or 'off'") self.stabilisation_param=stabilisation_param
[docs] def get_step_size(self, algorithm): """ Applies the B-B rule to calculate the step size (`step_size`) Returns -------- the calculated step size:float """ #For the first iteration we use an initial step size because the BB step size requires a previous iterate. if self.store_x is None: self.store_x=algorithm.x.copy() # We store the last iterate in order to calculate the BB step size self.store_grad=algorithm.gradient_update.copy()# We store the last gradient in order to calculate the BB step size return self.initial gradient_norm = algorithm.gradient_update.norm() #If the gradient is zero, gradient based algorithms will not update and te step size calculation will divide by zero so we stop iterations. if gradient_norm < 1e-8: raise StopIteration algorithm.x.subtract(self.store_x, out=self.store_x) algorithm.gradient_update.subtract(self.store_grad, out=self.store_grad) if self.is_short: ret = (self.store_x.dot(self.store_grad))/ (self.store_grad.dot(self.store_grad)) else: ret = (self.store_x.dot(self.store_x))/ (self.store_x.dot(self.store_grad)) #This computes the default stabilisation parameter, using the first three iterations if (algorithm.iteration <=3 and self.adaptive): self.stabilisation_param = min(self.stabilisation_param, self.store_x.norm() ) # Computes the step size as the minimum of the ret, above, and :math:`\Delta/\|g_k\|` ignoring any NaN values. ret = numpy.nanmin( numpy.array([ret, self.stabilisation_param/gradient_norm])) # We store the last iterate and gradient in order to calculate the BB step size self.store_x.fill(algorithm.x) self.store_grad.fill(algorithm.gradient_update) if self.mode == "alternate": self.is_short = not self.is_short return ret