diff --git a/Wrappers/Python/cil/framework/block.py b/Wrappers/Python/cil/framework/block.py index e7b40f9b72..ddf84f0d28 100644 --- a/Wrappers/Python/cil/framework/block.py +++ b/Wrappers/Python/cil/framework/block.py @@ -475,7 +475,6 @@ def binary_operations(self, operation, other, *args, **kwargs): # As axpyb cannot return anything we `continue` to skip the rest of the code block continue - else: raise ValueError('Unsupported operation', operation) if out is not None: @@ -600,7 +599,7 @@ def __div__(self, other): # __rdiv__ def __truediv__(self, other): return self.divide(other) - + def __pow__(self, other): return self.power(other) # reverse operand diff --git a/Wrappers/Python/cil/framework/data_container.py b/Wrappers/Python/cil/framework/data_container.py index 2614862f19..f807e2023f 100644 --- a/Wrappers/Python/cil/framework/data_container.py +++ b/Wrappers/Python/cil/framework/data_container.py @@ -567,6 +567,8 @@ def divide(self, other, *args, **kwargs): return _out return self.pixel_wise_binary(numpy.divide, other, *args, **kwargs) + + def power(self, other, *args, **kwargs): return self.pixel_wise_binary(numpy.power, other, *args, **kwargs) @@ -575,6 +577,8 @@ def maximum(self, x2, *args, **kwargs): def minimum(self,x2, out=None, *args, **kwargs): return self.pixel_wise_binary(numpy.minimum, x2=x2, out=out, *args, **kwargs) + + def sapyb(self, a, y, b, out=None, num_threads=NUM_THREADS): diff --git a/Wrappers/Python/cil/optimisation/algorithms/Algorithm.py b/Wrappers/Python/cil/optimisation/algorithms/Algorithm.py index 6a28439422..ddfb42212d 100644 --- a/Wrappers/Python/cil/optimisation/algorithms/Algorithm.py +++ b/Wrappers/Python/cil/optimisation/algorithms/Algorithm.py @@ -51,6 +51,18 @@ def __init__(self, update_objective_interval=1): self.update_objective_interval = update_objective_interval self.iter_string = 'Iter' + def _reset_iteration_state(self): + '''Resets the iteration counter and the objective/iteration history. + + This is intended for internal use by step-size rules that re-run the + algorithm during set-up (e.g. Bayesian optimisation of the step sizes), + so that they do not need to reach into the algorithm's private state. + ''' + self.iteration = -1 + self.__loss = [] + self._iteration = [] + self._total_iterations = 1 + def set_up(self, *args, **kwargs): '''Set up the algorithm''' raise NotImplementedError @@ -185,14 +197,14 @@ def loss(self): @property def update_objective_interval(self): '''gets the update_objective_interval''' - return self.__update_objective_interval + return self._update_objective_interval @update_objective_interval.setter def update_objective_interval(self, value): '''sets the update_objective_interval''' if not isinstance(value, Integral) or value < 0: raise ValueError('interval must be an integer >= 0') - self.__update_objective_interval = value + self._update_objective_interval = value def run(self, iterations=None, callbacks: Optional[List[Callback]] = None, verbose=1): r"""run upto :code:`iterations` with callbacks/logging. diff --git a/Wrappers/Python/cil/optimisation/algorithms/PD3O.py b/Wrappers/Python/cil/optimisation/algorithms/PD3O.py index ac88d328b2..16a6f28d61 100644 --- a/Wrappers/Python/cil/optimisation/algorithms/PD3O.py +++ b/Wrappers/Python/cil/optimisation/algorithms/PD3O.py @@ -125,7 +125,6 @@ class PD3O(Algorithm): Yan, M. A New Primal–Dual Algorithm for Minimizing the Sum of Three Functions with a Linear Operator. J Sci Comput 76, 1698–1717 (2018). https://doi.org/10.1007/s10915-018-0680-3 """ - def __init__(self, f, g, h, operator, delta=None, gamma=None, initial=None, **kwargs): super(PD3O, self).__init__(**kwargs) @@ -139,7 +138,8 @@ def set_up(self, f, g, h, operator, delta=None, gamma=None, initial=None,**kwarg logging.info("{} setting up".format(self.__class__.__name__, )) if isinstance(f, ZeroFunction): - warnings.warn(" If f is the ZeroFunction, then PD3O = PDHG. Please use PDHG instead. Otherwise, select a relatively small parameter gamma ", UserWarning) + warnings.warn( + " If f is the ZeroFunction, then PD3O = PDHG. Please use PDHG instead. Otherwise, select a relatively small parameter gamma ", UserWarning) if gamma is None: gamma = 1.0/operator.norm() @@ -177,7 +177,6 @@ def set_up(self, f, g, h, operator, delta=None, gamma=None, initial=None,**kwarg self.s_old.sapyb(1, self.s, self.delta, out=self.s_old) self.h.proximal_conjugate(self.s_old, self.delta, out=self.s) - def update(self): r""" Performs a single iteration of the PD3O algorithm """ @@ -185,7 +184,6 @@ def update(self): # Following equations 4 in https://link.springer.com/article/10.1007/s10915-018-0680-3 # in this case order of proximal steps we recover the (primal) PDHG, when f=0 - tmp = self.x_old self.x_old = self.x self.x = tmp diff --git a/Wrappers/Python/cil/optimisation/algorithms/PDHG.py b/Wrappers/Python/cil/optimisation/algorithms/PDHG.py index 28abf53a48..39c06328e2 100644 --- a/Wrappers/Python/cil/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/cil/optimisation/algorithms/PDHG.py @@ -20,11 +20,14 @@ from cil.framework import DataContainer, BlockDataContainer from cil.optimisation.algorithms import Algorithm +from cil.optimisation.utilities import StepSizeRule, PDHGStronglyConvexUpdate, PDHGConstantStepSize, PDHGAdaptiveStepSize2013, PDHGAdaptiveStepSize2015 +from cil.optimisation.utilities.StepSizeMethods import _validate_pdhg_step_sizes import warnings import numpy as np from numbers import Number import logging + log = logging.getLogger(__name__) @@ -47,16 +50,10 @@ class PDHG(Algorithm): A convex function with a "simple" proximal. This function must map from the operator domain to the Reals. See below for details. operator : LinearOperator A Linear Operator. - sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (tau*norm(K)**2) if tau is provided - Step size for the dual problem. Needs to obey constraints with tau and operator norm to satisfy convergence guarantees, see below for details. - tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (sigma*norm(K)**2) if sigma is provided - Step size for the primal problem. Needs to obey constraints with sigma and operator's norm to satisfy convergence guarantees, see below for details. + step_size: + Either a PDHG compatible step size rule or a `list` or `tuple` of (tau, sigma) where sigma is the step size for the dual problem and tau is the step size for the primal problem. The step sizes can be either None, scalar or array-objects. If not provided, default values will be set based on the operator norm as described below. initial : `DataContainer`, or `list` or `tuple` of `DataContainer`s, optional, default is a DataContainer of zeros for both primal and dual variables Initial point for the PDHG algorithm. If just one data container is provided, it is used for the primal and the dual variable is initialised as zeros. If a list or tuple is passed, the first element is used for the primal variable and the second one for the dual variable. If either of the two is not provided, it is initialised as a DataContainer of zeros. - gamma_g : positive :obj:`float`, optional, default=None - Strongly convex constant if the function g is strongly convex. Allows primal acceleration of the PDHG algorithm. - gamma_fconj : positive :obj:`float`, optional, default=None - Strongly convex constant if the convex conjugate of f is strongly convex. Allows dual acceleration of the PDHG algorithm. **kwargs: update_objective_interval : :obj:`int`, optional, default=1 @@ -65,20 +62,20 @@ class PDHG(Algorithm): Checks scalar sigma and tau values satisfy convergence criterion and warns if not satisfied. Can be computationally expensive for custom sigma or tau values. theta : Float between 0 and 1, default 1.0 Relaxation parameter for the over-relaxation of the primal variable. - - + gamma_g : positive :obj:`float`, optional, default=None + Note: this is being deprecated. Strongly convex constant if the function g is strongly convex. Allows primal acceleration of the PDHG algorithm. + gamma_fconj : positive :obj:`float`, optional, default=None + Note: this is being deprecated. Strongly convex constant if the convex conjugate of f is strongly convex. Allows dual acceleration of the PDHG algorithm. + sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (tau*norm(K)**2) if tau is provided + Step size for the dual problem. Note: this is being deprecated. In the future, please pass this as part of the `step_size` argument, either as a tuple of (tau,) or using a compatible step size rule. + tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (sigma*norm(K)**2) if sigma is provided + Step size for the primal problem. In the future, please pass this as part of the `step_size` argument, either as a tuple of (tau,) or using a compatible step size rule. Example ------- In our CIL-Demos repository (https://github.com/TomographicImaging/CIL-Demos) you can find examples using the PDHG algorithm for different imaging problems, such as Total Variation denoising, Total Generalised Variation inpainting and Total Variation Tomography reconstruction. More examples can also be found in :cite:`Jorgensen_et_al_2021`, :cite:`Papoutsellis_et_al_2021`. - Note - ---- - - Currently, the strongly convex constants are passed as parameters of PDHG. - In the future, these parameters will be properties of the corresponding functions. - Notes ----- @@ -156,7 +153,7 @@ class PDHG(Algorithm): \sigma = \frac{1}{\tau\|K\|^{2}} - - To monitor the convergence of the algorithm, we compute the primal/dual objectives and the primal-dual gap in :meth:`update_objective`.\ + - To monitor the convergence of the algorithm, we compute the primal/dual objectives and the primal-dual gap in :meth:`objective`.\ The primal objective is @@ -189,73 +186,51 @@ class PDHG(Algorithm): Computing these objectives can be costly, so it is better to compute every some iterations. To do this, use ``update_objective_interval = #number``. - - PDHG algorithm can be accelerated if the functions :math:`f^{*}` and/or :math:`g` are strongly convex. In these cases, the step-sizes :math:`\sigma` and :math:`\tau` are updated using the :meth:`update_step_sizes` method. A function :math:`f` is strongly convex with constant :math:`\gamma>0` if - - .. math:: - - f(x) - \frac{\gamma}{2}\|x\|^{2} \quad\mbox{ is convex. } - - - * For instance the function :math:`\frac{1}{2}\|x\|^{2}_{2}` is :math:`\gamma` strongly convex for :math:`\gamma\in(-\infty,1]`. We say it is 1-strongly convex because it is the largest constant for which :math:`f - \frac{1}{2}\|\cdot\|^{2}` is convex. - - - * The :math:`\|\cdot\|_{1}` norm is not strongly convex. For more information, see `Strongly Convex `_. - - - * If :math:`g` is strongly convex with constant :math:`\gamma` then the step-sizes :math:`\sigma`, :math:`\tau` and :math:`\theta` are updated as: - - - .. math:: - :nowrap: - - \begin{aligned} - - \theta_{n} & = \frac{1}{\sqrt{1 + 2\gamma\tau_{n}}}\\ - \tau_{n+1} & = \theta_{n}\tau_{n}\\ - \sigma_{n+1} & = \frac{\sigma_{n}}{\theta_{n}} - - \end{aligned} - - * If :math:`f^{*}` is strongly convex, we swap :math:`\sigma` with :math:`\tau`. - - Note - ---- - The case where both functions are strongly convex is not available at the moment. """ - def __init__(self, f, g, operator, tau=None, sigma=None, initial=None, gamma_g=None, gamma_fconj=None, **kwargs): + def __init__(self, f, g, operator, step_size=None, initial=None, + **kwargs): """Initialisation of the PDHG algorithm""" + self.initial = initial + self._sigma = kwargs.pop('sigma', None) # To be deprecated + self._tau = kwargs.pop('tau', None) # To be deprecated self._theta = kwargs.pop('theta', 1.0) if self._theta > 1 or self._theta < 0: raise ValueError( "The relaxation parameter theta must be in the range [0,1], passed theta = {}".format(self.theta)) - self._check_convergence = kwargs.pop('check_convergence', True) + if step_size is not None: # To be deprecated + if self._sigma is not None or self._tau is not None: # To be deprecated + raise ValueError("The parameters `sigma` and `tau` are being deprecated in favour of `step_size`. You have passed both. Instead please pass these as part of the `step_size` argument, either as a tuple of (tau, sigma) or using a compatible step size rule.") - super().__init__(**kwargs) + if self._sigma is not None or self._tau is not None: # To be deprecated + warnings.warn("The parameters `sigma` and `tau` are being deprecated. In the future, please pass these as part of the `step_size` argument, either as a tuple of (tau, sigma) or using a compatible step size rule.", category=DeprecationWarning, stacklevel=2) + step_size = (self._tau, self._sigma) - self._tau = None - self._sigma = None + self._gamma_g = kwargs.pop('gamma_g', None) # To be deprecated + self._gamma_fconj = kwargs.pop('gamma_fconj', None) # To be deprecated + if self._gamma_g is not None or self._gamma_fconj is not None: # To be deprecated + warnings.warn("The parameter `gamma_g` is being deprecated. In the future, if you would like to utilise strong convexity you should use the step size method cil.optimisation.utilities.StepSizeMethods.PDHGStronglyConvexUpdate.", category=DeprecationWarning, stacklevel=2) + step_size = PDHGStronglyConvexUpdate(initial_step_size=( + self._tau, self._sigma), gamma_g=self._gamma_g, gamma_fconj=self._gamma_fconj) - # check for gamma_g, gamma_fconj, strongly convex constants - self._gamma_g = None - self._gamma_fconj = None - self.set_gamma_g(gamma_g) - self.set_gamma_fconj(gamma_fconj) + self._check_convergence = kwargs.pop('check_convergence', True) + + super().__init__(**kwargs) - self.set_up(f=f, g=g, operator=operator, tau=tau, - sigma=sigma, initial=initial) + self.set_up(f=f, g=g, operator=operator, + step_size=step_size, initial=initial) @property def tau(self): - """The primal step-size """ + """The primal step-size - Returns the currently being used step size for the primal problem. Note that this can be updated at each iteration if a step size rule is used.""" return self._tau @property def sigma(self): - """The dual step-size """ + """The dual step-size - Returns the currently being used step size for the dual problem. Note that this can be updated at each iteration if a step size rule is used.""" return self._sigma @property @@ -263,61 +238,7 @@ def theta(self): """The relaxation parameter for the over-relaxation of the primal variable """ return self._theta - @property - def gamma_g(self): - """The strongly convex constant for the function g """ - return self._gamma_g - - @property - def gamma_fconj(self): - """The strongly convex constant for the convex conjugate of the function f """ - return self._gamma_fconj - - def set_gamma_g(self, value): - '''Set the value of the strongly convex constant for function `g` - - Parameters - ---------- - value : a positive number or None - ''' - if self.gamma_fconj is not None and value is not None: - raise ValueError("The adaptive update of the PDHG stepsizes in the case where both functions are strongly convex is not implemented at the moment." + - "Currently the strongly convex constant of the convex conjugate of the function f has been specified as ", self.gamma_fconj) - - if isinstance(value, Number): - if value <= 0: - raise ValueError( - "Strongly convex constant is a positive number, {} is passed for the strongly convex function g.".format(value)) - self._gamma_g = value - elif value is None: - pass - else: - raise ValueError( - "Positive float is expected for the strongly convex constant of function g, {} is passed".format(value)) - - def set_gamma_fconj(self, value): - '''Set the value of the strongly convex constant for the convex conjugate of function `f` - - Parameters - ---------- - value : a positive number or None - ''' - if self.gamma_g is not None and value is not None: - raise ValueError("The adaptive update of the PDHG stepsizes in the case where both functions are strongly convex is not implemented at the moment." + - "Currently the strongly convex constant of the function g has been specified as ", self.gamma_g) - - if isinstance(value, Number): - if value <= 0: - raise ValueError( - "Strongly convex constant is positive, {} is passed for the strongly convex conjugate function of f.".format(value)) - self._gamma_fconj = value - elif value is None: - pass - else: - raise ValueError( - "Positive float is expected for the strongly convex constant of the convex conjugate of function f, {} is passed".format(value)) - - def set_up(self, f, g, operator, tau=None, sigma=None, initial=None): + def set_up(self, f, g, operator, step_size=[None, None], initial=None): """Initialisation of the algorithm Parameters @@ -328,14 +249,14 @@ def set_up(self, f, g, operator, tau=None, sigma=None, initial=None): A convex function with a "simple" proximal. operator : LinearOperator A Linear Operator. - sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (tau*norm(K)**2) if tau is provided - Step size for the dual problem. - tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default is 1.0/norm(K) or 1.0/ (sigma*norm(K)**2) if sigma is provided - Step size for the primal problem. initial : `DataContainer`, or `list` or `tuple` of `DataContainer`s, optional, default is a DataContainer of zeros for both primal and dual variables Initial point for the PDHG algorithm. If just one data container is provided, it is used for the primal and the dual variable is initialised as zeros. If a list or tuple is passed, the first element is used for the primal variable and the second one for the dual variable. If either of the two is not provided, it is initialised as a DataContainer of zeros. + step_size: + Either a PDHG compatible step size rule or a `list` or `tuple` of (tau, sigma) where sigma is the step size for the dual problem and tau is the step size for the primal problem. The step sizes can be either scalar or array-objects. If not provided, default values will be set based on the operator norm as described below. + """ + log.info("%s setting up", self.__class__.__name__) # Triplet (f, g, K) @@ -343,10 +264,20 @@ def set_up(self, f, g, operator, tau=None, sigma=None, initial=None): self.g = g self.operator = operator - self.set_step_sizes(sigma=sigma, tau=tau) + if step_size is None: # This line can be removed when sigma and tau deprecated + step_size = (None, None) + if isinstance(step_size, StepSizeRule): + if not hasattr(step_size, 'get_initial_step_size'): + raise ValueError( + "The step-size rule {} does not provide initial primal/dual step sizes " + "and is not compatible with PDHG.".format(type(step_size).__name__)) + self.step_size_rule = step_size + elif isinstance(step_size, (tuple, list)): + self.step_size_rule = PDHGConstantStepSize(step_size=step_size) + else: + raise ValueError("The `step_size` argument must be either None, a PDHG compatible step size rule or a tuple of (tau, sigma) where sigma is the step size for the dual problem and tau is the step size for the primal problem.") + - if self._check_convergence: - self.check_convergence() if isinstance(initial, (tuple, list)): if initial[0] is not None: @@ -370,11 +301,12 @@ def set_up(self, f, g, operator, tau=None, sigma=None, initial=None): self.x_tmp = self.operator.domain_geometry().allocate(0) self.y_tmp = self.operator.range_geometry().allocate(0) - if self.gamma_g is not None: - warnings.warn("Primal Acceleration of PDHG: The function g is assumed to be strongly convex with positive parameter `gamma_g`. You need to be sure that gamma_g = {} is the correct strongly convex constant for g. ".format(self.gamma_g)) + self._tau, self._sigma = self.step_size_rule.get_initial_step_size( + self) + _validate_pdhg_step_sizes(self._tau, self._sigma, self.operator) - if self.gamma_fconj is not None: - warnings.warn("Dual Acceleration of PDHG: The convex conjugate of function f is assumed to be strongly convex with positive parameter `gamma_fconj`. You need to be sure that gamma_fconj = {} is the correct strongly convex constant".format(self.gamma_fconj)) + if self._check_convergence: + self.check_convergence() self.configured = True log.info("%s configured", self.__class__.__name__) @@ -382,7 +314,7 @@ def set_up(self, f, g, operator, tau=None, sigma=None, initial=None): def _update_previous_solution(self): """ Swaps the references to current and previous solution based on the - :func:`~Algorithm.update_previous_solution` of the base class :class:`Algorithm`. + :func:`~Algorithm.previous_solution` of the base class :class:`Algorithm`. """ tmp = self.x_old self.x_old = self.x @@ -392,8 +324,10 @@ def get_output(self): " Returns the current solution. " return self.x_old - def update(self): - """Performs a single iteration of the PDHG algorithm""" + def _pdhg_update(self): + """Applies the primal-dual updates for one PDHG step (see the Notes in the + class docstring for the update equations).""" + # calculate x-bar and store in self.x_tmp self.x_old.sapyb((self.theta + 1.0), self.x, - self.theta, out=self.x_tmp) @@ -412,11 +346,12 @@ def update(self): self.g.proximal(self.x_tmp, self.tau, out=self.x) - # update_previous_solution() called after update by base class - # i.e current solution is now in x_old, previous solution is now in x + def update(self): + """Performs a single iteration of the PDHG algorithm""" + self._pdhg_update() # update the step sizes for special cases - self.update_step_sizes() + self._tau, self._sigma = self.step_size_rule.get_step_size(self) def check_convergence(self): """Check whether convergence criterion for PDHG is satisfied with scalar values of tau and sigma @@ -431,82 +366,23 @@ def check_convergence(self): Li, Y. and Yan, M., 2022. On the improved conditions for some primal-dual algorithms. arXiv preprint arXiv:2201.00139. """ - if isinstance(self.tau, Number) and isinstance(self.sigma, Number): - if self.sigma * self.tau * self.operator.norm()**2 > 4/3: + if isinstance(self.step_size_rule, PDHGConstantStepSize): + if isinstance(self.tau, Number) and isinstance(self.sigma, Number): + if self.sigma * self.tau * self.operator.norm()**2 > 4/3: + warnings.warn( + "Convergence criterion of PDHG for scalar step-sizes is not satisfied.") + return False + return True + else: warnings.warn( - "Convergence criterion of PDHG for scalar step-sizes is not satisfied.") - return False + "Convergence criterion can only be checked for scalar values of tau and sigma, tau={0}, sigma={1}".format(self.tau, self.sigma)) + elif isinstance(self.step_size_rule, PDHGAdaptiveStepSize2013) or isinstance(self.step_size_rule, PDHGAdaptiveStepSize2015): return True - warnings.warn( - "Convergence criterion can only be checked for scalar values of tau and sigma.") - return False - - def set_step_sizes(self, sigma=None, tau=None): - """Sets sigma and tau step-sizes for the PDHG algorithm. The step sizes can be either scalar or array-objects. - - Parameters - ---------- - sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None - Step size for the dual problem. - tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None - Step size for the primal problem. - - The user can set either, both or none. Values passed by the user will be accepted as long as they are positive numbers, - or correct shape array like objects. - """ - # Check acceptable values of the primal-dual step-sizes - if tau is not None: - if isinstance(tau, Number): - if tau <= 0: - raise ValueError( - "The step-sizes of PDHG must be positive, passed tau = {}".format(tau)) - elif tau.shape != self.operator.domain_geometry().shape: - raise ValueError(" The shape of tau = {0} is not the same as the shape of the domain_geometry = {1}".format( - tau.shape, self.operator.domain_geometry().shape)) - - if sigma is not None: - if isinstance(sigma, Number): - if sigma <= 0: - raise ValueError( - "The step-sizes of PDHG are positive, passed sigma = {}".format(sigma)) - elif sigma.shape != self.operator.range_geometry().shape: - raise ValueError(" The shape of sigma = {0} is not the same as the shape of the range_geometry = {1}".format( - sigma.shape, self.operator.range_geometry().shape)) - - # Default sigma and tau step-sizes - if tau is None and sigma is None: - self._sigma = 1/self.operator.norm() - self._tau = 1/self.operator.norm() - elif tau is not None and sigma is not None: - self._sigma = sigma - self._tau = tau - elif sigma is None and isinstance(tau, Number): - self._sigma = 1/(tau*self.operator.norm()**2) - self._tau = tau - elif tau is None and isinstance(sigma, Number): - self._sigma = sigma - self._tau = 1/(self.sigma*self.operator.norm()**2) else: - raise NotImplementedError( - "If using arrays for sigma or tau both must arrays must be provided.") + warnings.warn( + "Convergence checks not currently implemented for this type of step size rule.") + return False - def update_step_sizes(self): - """ - Updates step sizes in the cases of primal or dual acceleration using the strongly convexity property. - The case where both functions are strongly convex is not available at the moment. - """ - # Update sigma and tau based on the strong convexity of G - if self.gamma_g is not None: - self._theta = 1.0 / np.sqrt(1 + 2 * self.gamma_g * self.tau) - self._tau *= self.theta - self._sigma /= self.theta - - # Update sigma and tau based on the strong convexity of F - # Following operations are reversed due to symmetry, sigma --> tau, tau -->sigma - if self.gamma_fconj is not None: - self._theta = 1.0 / np.sqrt(1 + 2 * self.gamma_fconj * self.sigma) - self._sigma *= self.theta - self._tau /= self.theta def update_objective(self): """Evaluates the primal objective, the dual objective and the primal-dual gap.""" diff --git a/Wrappers/Python/cil/optimisation/algorithms/SPDHG.py b/Wrappers/Python/cil/optimisation/algorithms/SPDHG.py index 20f8d11ce1..bf56b81e9f 100644 --- a/Wrappers/Python/cil/optimisation/algorithms/SPDHG.py +++ b/Wrappers/Python/cil/optimisation/algorithms/SPDHG.py @@ -21,7 +21,8 @@ from cil.optimisation.operators import BlockOperator import numpy as np import logging -from cil.optimisation.utilities import Sampler +from cil.optimisation.utilities import Sampler, StepSizeRule, SPDHGConstantStepSize +from cil.optimisation.utilities.StepSizeMethods import _validate_spdhg_step_sizes from numbers import Number import warnings from cil.framework import BlockDataContainer @@ -30,7 +31,7 @@ class SPDHG(Algorithm): - r'''Stochastic Primal Dual Hybrid Gradient (SPDHG) solves separable optimisation problems of the type: + r'''Stochastic Primal Dual Hybrid Gradient (SPDHG) solves separable optimisation problems of the type: .. math:: \min_{x} f(Kx) + g(x) = \min_{x} \sum f_i(K_i x) + g(x) @@ -44,14 +45,10 @@ class SPDHG(Algorithm): A convex function with a "simple" proximal operator : BlockOperator BlockOperator must contain Linear Operators - tau : positive float, optional - Step size parameter for the primal problem. If `None` will be computed by algorithm, see note for details. - sigma : list of positive float, optional - List of Step size parameters for dual problem. If `None` will be computed by algorithm, see note for details. - initial : DataContainer, optional - Initial point for the SPDHG algorithm. The default value is a zero DataContainer in the range of the `operator`. - gamma : float, optional - Parameter controlling the trade-off between the primal and dual step sizes + step_size : tuple of (tau, sigma), optional + A tuple containing the step size parameters for the primal and dual problems. If `None` will be computed by algorithm, see note for details. + initial : `DataContainer`, or `list` or `tuple` of `DataContainer`s, optional, default is a DataContainer of zeros for both primal and dual variables + Initial point for the PDHG algorithm. If just one data container is provided, it is used for the primal and the dual variable is initialised as zeros. If a list or tuple is passed, the first element is used for the primal variable and the second one for the dual variable. If either of the two is not provided, it is initialised as a DataContainer of zeros. sampler: `cil.optimisation.utilities.Sampler`, optional A `Sampler` controllingthe selection of the next index for the SPDHG update. If `None`, a sampler will be created for uniform random sampling with replacement. See notes. @@ -134,17 +131,28 @@ class SPDHG(Algorithm): Physics in Medicine & Biology, Volume 64, Number 22, 2019. https://doi.org/10.1088/1361-6560/ab3d07 ''' - def __init__(self, f=None, g=None, operator=None, tau=None, sigma=None, + def __init__(self, f=None, g=None, operator=None, step_size=None, initial=None, sampler=None, prob_weights=None, **kwargs): - update_objective_interval = kwargs.pop('update_objective_interval', 1) - super(SPDHG, self).__init__( - update_objective_interval=update_objective_interval) - - self.set_up(f=f, g=g, operator=operator, sigma=sigma, tau=tau, + self.initial = initial + self._sigma = kwargs.pop('sigma', None) # To be deprecated + self._tau = kwargs.pop('tau', None) # To be deprecated + + if step_size is not None: # To be deprecated + if self._sigma is not None or self._tau is not None: # To be deprecated + raise ValueError("The parameters `sigma` and `tau` are being deprecated in favour of `step_size`. You have passed both. Instead please pass these as part of the `step_size` argument, either as a tuple of (tau, sigma ) or using a compatible step size rule.") + + if self._sigma is not None or self._tau is not None: # To be deprecated + warnings.warn("The parameters `sigma` and `tau` are being deprecated. In the future, please pass these as part of the `step_size` argument, either as a tuple of (tau, sigma) or using a compatible step size rule.", category=DeprecationWarning, stacklevel=2) + step_size = (self._tau, self._sigma) + + + super(SPDHG, self).__init__(**kwargs) + + self.set_up(f=f, g=g, operator=operator, step_size=step_size, initial=initial, sampler=sampler, prob_weights=prob_weights) - def set_up(self, f, g, operator, sigma=None, tau=None, + def set_up(self, f, g, operator, step_size=None, initial=None, sampler=None, prob_weights=None): '''set-up of the algorithm ''' @@ -178,18 +186,41 @@ def set_up(self, f, g, operator, sigma=None, tau=None, # Set the norms of the operators self._norms = operator.get_norms_as_list() - self.set_step_sizes(sigma=sigma, tau=tau) + if step_size is None: # This line can be removed when sigma and tau deprecated + step_size = (None, None) + if isinstance(step_size, StepSizeRule): + if not hasattr(step_size, 'get_initial_step_size'): + raise ValueError( + "The step-size rule {} does not provide initial primal/dual step sizes " + "and is not compatible with SPDHG.".format(type(step_size).__name__)) + self.step_size_rule = step_size + elif isinstance(step_size, (tuple, list)): + self.step_size_rule = SPDHGConstantStepSize(step_size=step_size) + else: + raise ValueError("The `step_size` argument must be either None, a SPDHG compatible step size rule or a tuple of (tau, sigma) where sigma is the step size for the dual problem and tau is the step size for the primal problem.") + + + if isinstance(initial, (tuple, list)): + if initial[0] is not None: + self.x = initial[0].copy() + else: + self.x = self.operator.domain_geometry().allocate(0) + + + if len(initial) > 1 and initial[1] is not None: + self._y_old = initial[1].copy() + else: + self._y_old = self.operator.range_geometry().allocate(0) - # initialize primal variable - if initial is None: - self.x = self.operator.domain_geometry().allocate(0) else: - self.x = initial.copy() + self._y_old = self.operator.range_geometry().allocate(0) + if initial is None: + self.x = self.operator.domain_geometry().allocate(0) + else: + self.x = initial.copy() self._x_tmp = self.operator.domain_geometry().allocate(0) - - # initialize dual variable to 0 - self._y_old = operator.range_geometry().allocate(0) + # This can be removed once #1863 is fixed if not isinstance(self._y_old, BlockDataContainer): self._y_old = BlockDataContainer(self._y_old) @@ -200,9 +231,14 @@ def set_up(self, f, g, operator, sigma=None, tau=None, # relaxation parameter self._theta = 1 + self._tau, self._sigma = self.step_size_rule.get_initial_step_size( + self) + _validate_spdhg_step_sizes(self._tau, self._sigma, self._ndual_subsets) + self.configured = True logging.info("{} configured".format(self.__class__.__name__, )) + @property def sigma(self): return self._sigma @@ -211,110 +247,7 @@ def sigma(self): def tau(self): return self._tau - def set_step_sizes_from_ratio(self, gamma=1.0, rho=0.99): - r""" Sets gamma, the step-size ratio for the SPDHG algorithm. Currently gamma takes a scalar value. - - The step sizes `sigma` and `tau` are set using the equations: - - .. math:: \sigma_i= \frac{\gamma\rho }{\|K_i\|^2} - - .. math:: \tau = \rho\min_i([ \frac{p_i }{\sigma_i \|K_i\|^2}) - - - Parameters - ---------- - gamma : Positive float - parameter controlling the trade-off between the primal and dual step sizes - rho : Positive float - parameter controlling the size of the product :math:`\sigma\tau` - - - - """ - if isinstance(gamma, Number): - if gamma <= 0: - raise ValueError( - "The step-sizes of SPDHG are positive, gamma should also be positive") - - else: - raise ValueError( - "We currently only support scalar values of gamma") - if isinstance(rho, Number): - if rho <= 0: - raise ValueError( - "The step-sizes of SPDHG are positive, rho should also be positive") - - else: - raise ValueError( - "We currently only support scalar values of gamma") - - self._sigma = [gamma * rho / ni for ni in self._norms] - values = [rho*pi / (si * ni**2) for pi, ni, - si in zip(self._prob_weights, self._norms, self._sigma)] - self._tau = min([value for value in values if value > 1e-8]) - - def set_step_sizes(self, sigma=None, tau=None): - r""" Sets sigma and tau step-sizes for the SPDHG algorithm after the initial set-up. The step sizes can be either scalar or array-objects. - - When setting `sigma` and `tau`, there are 4 possible cases considered by setup function: - - - Case 1: If neither `sigma` or `tau` are provided then `sigma` is set using the formula: - - .. math:: \sigma_i= \frac{0.99}{\|K_i\|^2} - - and `tau` is set as per case 2 - - - Case 2: If `sigma` is provided but not `tau` then `tau` is calculated using the formula - - .. math:: \tau = 0.99\min_i( \frac{p_i}{ (\sigma_i \|K_i\|^2) }) - - - Case 3: If `tau` is provided but not `sigma` then `sigma` is calculated using the formula - - .. math:: \sigma_i= \frac{0.99 p_i}{\tau\|K_i\|^2} - - - Case 4: Both `sigma` and `tau` are provided. - - - Parameters - ---------- - sigma : list of positive float, optional, default= see docstring - List of Step size parameters for dual problem - tau : positive float, optional, default= see docstring - Step size parameter for primal problem - - """ - gamma = 1. - rho = .99 - if sigma is not None: - if len(sigma) == self._ndual_subsets: - if all(isinstance(x, Number) and x > 0 for x in sigma): - pass - else: - raise ValueError( - "Sigma expected to be a positive number.") - - else: - raise ValueError( - "Please pass a list of floats to sigma with the same number of entries as number of operators") - self._sigma = sigma - - elif tau is None: - self._sigma = [gamma * rho / ni for ni in self._norms] - else: - self._sigma = [ - rho*pi / (tau*ni**2) for ni, pi in zip(self._norms, self._prob_weights)] - - if tau is None: - values = [rho*pi / (si * ni**2) for pi, ni, - si in zip(self._prob_weights, self._norms, self._sigma)] - self._tau = min([value for value in values if value > 1e-8]) - - else: - if not (isinstance(tau, Number) and tau > 0): - raise ValueError( - "The step-sizes of SPDHG must be positive, passed tau = {}".format(tau)) - - self._tau = tau + def check_convergence(self): """ Checks whether convergence criterion for SPDHG is satisfied with the current scalar values of tau and sigma @@ -341,11 +274,8 @@ def check_convergence(self): else: raise ValueError( 'Convergence criterion currently can only be checked for scalar values of tau and sigma[i].') - - def update(self): - """ Runs one iteration of SPDHG - - """ + + def _spdhg_update(self, i): # Gradient descent for the primal variable # x_tmp = x - tau * zbar self._zbar.sapyb(self._tau, self.x, -1., out=self._x_tmp) @@ -353,24 +283,23 @@ def update(self): self.g.proximal(self._x_tmp, self._tau, out=self.x) - # Choose subset - i = next(self._sampler) + # Gradient ascent for the dual variable # y_k = y_old[i] + sigma[i] * K[i] x try: - y_k = self.operator[i].direct(self.x) + self.y_k = self.operator[i].direct(self.x) except IndexError: raise IndexError( 'The sampler has outputted an index larger than the number of operators to sample from. Please ensure your sampler samples from {0,1,...,len(operator)-1} only.') - y_k.sapyb(self._sigma[i], self._y_old[i], 1., out=y_k) + self.y_k.sapyb(self._sigma[i], self._y_old[i], 1., out=self.y_k) - y_k = self.f[i].proximal_conjugate(y_k, self._sigma[i]) + self.y_k = self.f[i].proximal_conjugate(self.y_k, self._sigma[i]) # Back-project # x_tmp = K[i]^*(y_k - y_old[i]) - y_k.subtract(self._y_old[i], out=self._y_old[i]) + self.y_k.subtract(self._y_old[i], out=self._y_old[i]) self.operator[i].adjoint(self._y_old[i], out=self._x_tmp) # Update backprojected dual variable and extrapolate @@ -383,9 +312,22 @@ def update(self): self._z.sapyb(1., self._x_tmp, self._theta / self._prob_weights[i], out=self._zbar) - # save previous iteration - self._save_previous_iteration(i, y_k) + + + def update(self): + """ Runs one iteration of SPDHG + """ + # Choose subset + self._index = next(self._sampler) + + self._spdhg_update(self._index) + + # save previous iteration + self._save_previous_iteration(self._index, self.y_k) + + self._tau, self._sigma = self.step_size_rule.get_step_size(self) + def update_objective(self): # p1 = self.f(self.operator.direct(self.x)) + self.g(self.x) p1 = sum(self.f[i](op.direct(self.x)) for i, op in enumerate(self.operator.operators)) diff --git a/Wrappers/Python/cil/optimisation/utilities/StepSizeMethods.py b/Wrappers/Python/cil/optimisation/utilities/StepSizeMethods.py index dc1f980280..712147995b 100644 --- a/Wrappers/Python/cil/optimisation/utilities/StepSizeMethods.py +++ b/Wrappers/Python/cil/optimisation/utilities/StepSizeMethods.py @@ -17,15 +17,152 @@ # - CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt from abc import ABC, abstractmethod -import numpy +from copy import deepcopy from numbers import Number import logging +import numpy as np log = logging.getLogger(__name__) + +def _resolve_pdhg_step_sizes(tau, sigma, operator): + """Validate user-supplied PDHG primal/dual step sizes and fill in defaults. + + Either, both or neither of ``tau``/``sigma`` may be provided. Scalars must be + positive; array-like values must match the operator's domain (``tau``) or + range (``sigma``) geometry. Any missing value is derived from the operator norm. + + Returns + ------- + tuple + ``(tau, sigma)`` with any ``None`` entries replaced by their defaults. + """ + # Check acceptable values of the primal-dual step-sizes + if tau is not None: + if isinstance(tau, Number): + if tau <= 0: + raise ValueError( + "The step-sizes of PDHG must be positive, passed tau = {}".format(tau)) + elif tau.shape != operator.domain_geometry().shape: + raise ValueError(" The shape of tau = {0} is not the same as the shape of the domain_geometry = {1}".format( + tau.shape, operator.domain_geometry().shape)) + + if sigma is not None: + if isinstance(sigma, Number): + if sigma <= 0: + raise ValueError( + "The step-sizes of PDHG are positive, passed sigma = {}".format(sigma)) + elif sigma.shape != operator.range_geometry().shape: + raise ValueError(" The shape of sigma = {0} is not the same as the shape of the range_geometry = {1}".format( + sigma.shape, operator.range_geometry().shape)) + + # Default sigma and tau step-sizes + if tau is None and sigma is None: + sigma = 1.0/operator.norm() + tau = 1.0/operator.norm() + elif tau is not None and sigma is not None: + pass + elif sigma is None and isinstance(tau, Number): + sigma = 1./(tau*operator.norm()**2) + elif tau is None and isinstance(sigma, Number): + tau = 1./(sigma*operator.norm()**2) + else: + raise NotImplementedError( + "If using arrays for sigma or tau both must arrays must be provided.") + return tau, sigma + + +def _validate_pdhg_step_sizes(tau, sigma, operator): + """Validate step sizes produced by a PDHG step-size rule. + + ``tau`` and ``sigma`` must each be a positive scalar or an array-like whose + shape matches the operator's domain (``tau``) or range (``sigma``) geometry. + Raises ``ValueError`` if the step sizes have the wrong type or shape, e.g. + because an incompatible step-size rule was passed to PDHG. + """ + for name, value, shape in ( + ("tau", tau, operator.domain_geometry().shape), + ("sigma", sigma, operator.range_geometry().shape)): + if value is None: + raise ValueError( + "The step-size rule returned {0} = None. This step-size rule is " + "not compatible with PDHG.".format(name)) + if isinstance(value, Number): + if value <= 0: + raise ValueError( + "The step-sizes of PDHG must be positive, got {0} = {1}.".format(name, value)) + elif getattr(value, "shape", None) != shape: + raise ValueError( + "The shape of {0} = {1} is not the same as the expected shape = {2}. " + "This step-size rule may not be compatible with PDHG.".format( + name, getattr(value, "shape", type(value)), shape)) + + +def _validate_spdhg_step_sizes(tau, sigma, n_operators): + """Validate step sizes produced by an SPDHG step-size rule. + + ``tau`` must be a positive scalar and ``sigma`` a list/array of positive + numbers with one entry per operator. Raises ``ValueError`` otherwise, e.g. + because an incompatible step-size rule was passed to SPDHG. + """ + if not (isinstance(tau, Number) and tau > 0): + raise ValueError( + "The primal step-size tau of SPDHG must be a positive number, got tau = {0}. " + "This step-size rule may not be compatible with SPDHG.".format(tau)) + if isinstance(sigma, Number) or len(sigma) != n_operators: + raise ValueError( + "The dual step-size sigma of SPDHG must be a list of {0} positive numbers " + "(one per operator), got sigma = {1}. This step-size rule may not be " + "compatible with SPDHG.".format(n_operators, sigma)) + if not all(isinstance(si, Number) and si > 0 for si in sigma): + raise ValueError( + "The dual step-sizes sigma of SPDHG must all be positive numbers, got sigma = {0}.".format(sigma)) + + +def _spdhg_sigma_from_gamma(gamma, rho, norms): + """Dual SPDHG step sizes from the ratio ``gamma`` and product parameter ``rho``.""" + return [gamma * rho / ni for ni in norms] + + +def _spdhg_tau_from_sigma(sigma, norms, prob_weights, rho): + """Primal SPDHG step size consistent with the dual step sizes ``sigma``. + + Blocks are skipped only when they impose no constraint at all: a zero operator + norm contributes nothing to :math:`\\tau\\sigma_i\\|K_i\\|^2 \\leq \\rho p_i`, and a + block with zero probability is never sampled. Every other block must be honoured, + however small its candidate value -- discarding small candidates would drop + precisely the constraints that bind and return a tau that violates convergence. + """ + values = [rho * pi / (si * ni**2) + for pi, ni, si in zip(prob_weights, norms, sigma) + if ni > 0 and pi > 0] + if not values: + raise ValueError( + "Could not compute an SPDHG primal step size from the given sigma. " + "Every block had either a zero operator norm or a zero probability weight. " + "Check the operator norms and probability weights.") + return min(values) + + 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. + + Notes + ----- + There are two families of step-size rule in CIL: + + * **Gradient-based algorithms** (:class:`~cil.optimisation.algorithms.GD`, + :class:`~cil.optimisation.algorithms.ISTA`, :class:`~cil.optimisation.algorithms.FISTA`) + call :meth:`get_step_size` after each gradient calculation and expect a single scalar + step size to be returned. These rules only need to implement :meth:`get_step_size`. + * **Primal-dual algorithms** (:class:`~cil.optimisation.algorithms.PDHG`, + :class:`~cil.optimisation.algorithms.SPDHG`) additionally require a + ``get_initial_step_size(self, algorithm)`` method which is called once during set-up and + returns the initial ``(tau, sigma)`` step-size pair. Their :meth:`get_step_size` is called + at the end of every iteration and returns the updated ``(tau, sigma)`` pair for the next + iteration. The algorithms check for the presence of ``get_initial_step_size`` (via + ``hasattr``) and raise a clear error if an incompatible, gradient-only, rule is passed. """ def __init__(self): @@ -56,6 +193,7 @@ class ConstantStepSize(StepSizeRule): def __init__(self, step_size): '''Initialises the constant step size rule + Parameters: ------------- step_size : float, the constant step size @@ -90,7 +228,7 @@ class ArmijoStepSizeRule(StepSizeRule): 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)) + max_iterations: integer, optional, default is np.ceil (2 * np.log10(alpha) / np.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`. @@ -103,13 +241,19 @@ def __init__(self, alpha=1e6, beta=0.5, max_iterations=None, warmstart=True): ''' self.alpha_orig = alpha - self.alpha = alpha - self.beta = beta + 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.max_iterations = np.ceil( + 2 * np.log10(self.alpha_orig) / np.log10(2)) - self.warmstart=warmstart + self.warmstart = warmstart def get_step_size(self, algorithm): """ @@ -124,30 +268,35 @@ def get_step_size(self, algorithm): if not self.warmstart: self.alpha = self.alpha_orig - f_x = algorithm.calculate_objective_function_at_point(algorithm.solution) + 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) + 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) + 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 @@ -161,10 +310,11 @@ class BarzilaiBorweinStepSizeRule(StepSizeRule): - :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. + 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 @@ -175,6 +325,8 @@ class BarzilaiBorweinStepSizeRule(StepSizeRule): 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 @@ -188,30 +340,30 @@ def __init__(self, initial, mode='short', stabilisation_param="auto"): '''Initialises the step size rule ''' - self.mode=mode + 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". ') + raise ValueError( + 'Mode should be chosen from "long", "short" or "alternate". ') - self.store_grad=None - self.store_x=None - self.initial=initial + self.store_grad = None + self.store_x = None + self.initial = initial if stabilisation_param == 'auto': self.adaptive = True - stabilisation_param = numpy.inf + stabilisation_param = np.inf elif stabilisation_param == "off": self.adaptive = False - stabilisation_param = numpy.inf - elif ( isinstance(stabilisation_param, Number) and stabilisation_param >=0): + stabilisation_param = np.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 - - + raise TypeError( + " The stabilisation_param should be 'auto', a positive number or 'off'") + self.stabilisation_param = stabilisation_param def get_step_size(self, algorithm): """ @@ -222,37 +374,1451 @@ def get_step_size(self, algorithm): the calculated step size:float """ - #For the first iteration we use an initial step size because the BB step size requires a previous iterate. + # 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 + # We store the last iterate in order to calculate the BB step size + self.store_x = algorithm.x.copy() + # We store the last gradient in order to calculate the BB step size + self.store_grad = algorithm.gradient_update.copy() 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 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) + 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)) + 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)) + 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() ) + # 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])) + ret = np.nanmin(np.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 + self.is_short = not self.is_short return ret + + +class PDHGStronglyConvexUpdate(StepSizeRule): + r'''Updates step sizes (theta, sigma, tau) in the PDHG algorithm in the cases of primal or dual acceleration using the strongly convexity property. + The case where both functions are strongly convex is not available at the moment. + + + The PDHG algorithm can be accelerated if the functions :math:`f^{*}` and/or :math:`g` are strongly convex. In these cases, the step-sizes :math:`\sigma` and :math:`\tau` are updated using the :meth:`update_step_sizes` method. A function :math:`f` is strongly convex with constant :math:`\gamma>0` if + + .. math:: + + f(x) - \frac{\gamma}{2}\|x\|^{2} \quad\mbox{ is convex. } + + + * For instance the function :math:`\frac{1}{2}\|x\|^{2}_{2}` is :math:`\gamma` strongly convex for :math:`\gamma\in(-\infty,1]`. We say it is 1-strongly convex because it is the largest constant for which :math:`f - \frac{1}{2}\|\cdot\|^{2}` is convex. + + + * The :math:`\|\cdot\|_{1}` norm is not strongly convex. For more information, see `Strongly Convex `_. + + + * If :math:`g` is strongly convex with constant :math:`\gamma` then the step-sizes :math:`\sigma`, :math:`\tau` and :math:`\theta` are updated as: + + + .. math:: + :nowrap: + + \begin{aligned} + + \theta_{n} & = \frac{1}{\sqrt{1 + 2\gamma\tau_{n}}}\\ + \tau_{n+1} & = \theta_{n}\tau_{n}\\ + \sigma_{n+1} & = \frac{\sigma_{n}}{\theta_{n}} + + \end{aligned} + + * If :math:`f^{*}` is strongly convex, we swap :math:`\sigma` with :math:`\tau`. + + + + Parameters + ------------- + gamma_g : positive :obj:`float`, optional, default=None + Strongly convex constant if the function g is strongly convex. Allows primal acceleration of the PDHG algorithm. + gamma_fconj : positive :obj:`float`, optional, default=None + Strongly convex constant if the convex conjugate of f is strongly convex. Allows dual acceleration of the PDHG algorithm. + + Note + ---- + This rule is also selected automatically when the (deprecated) ``gamma_g`` or + ``gamma_fconj`` keyword arguments are passed directly to + :class:`~cil.optimisation.algorithms.PDHG`. New code should construct and pass this + rule explicitly via the ``step_size`` argument instead. + ''' + + def __init__(self, initial_step_size=(None, None), gamma_g=None, gamma_fconj=None): + '''Initialises the step size rule''' + + if gamma_g is not None and gamma_fconj is not None: + raise NotImplementedError( + "PDHG strongly convex step size update not implemented for both primal and dual acceleration. Please choose only one of gamma_g or gamma_fconj.") + + if not (gamma_g is None or (isinstance(gamma_g, Number) and gamma_g > 0)): + raise ValueError( + "A positive float is expected for the strongly convex constant of the function g, {} is passed".format(gamma_g)) + self.gamma_g = gamma_g + + if not (gamma_fconj is None or (isinstance(gamma_fconj, Number) and gamma_fconj > 0)): + raise ValueError( + "A positive float is expected for the strongly convex constant of the convex conjugate of function f, {} is passed".format(gamma_fconj)) + self.gamma_fconj = gamma_fconj + + self.initial_step_size = initial_step_size + if len(initial_step_size) != 2: + raise ValueError( + "initial_step_size should be a list or tuple of length two, step_size = {}".format(initial_step_size)) + + def get_initial_step_size(self, algorithm): + """Sets sigma and tau step-sizes for the PDHG algorithm. The step sizes can be either scalar or array-objects. + + Parameters + ---------- + sigma : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None + Step size for the dual problem. + tau : positive :obj:`float`, or `np.ndarray`, `DataContainer`, `BlockDataContainer`, optional, default=None + Step size for the primal problem. + + The user can set either, both or none. Values passed by the user will be accepted as long as they are positive numbers, + or correct shape array like objects. + """ + self.tau = self.initial_step_size[0] + self.sigma = self.initial_step_size[1] + self.tau, self.sigma = _resolve_pdhg_step_sizes( + self.tau, self.sigma, algorithm.operator) + return self.tau, self.sigma + + def get_step_size(self, algorithm): + """ + Applies the PDHG strongly convex step size update to calculate the new primal and dual step sizes + + Returns + -------- + + """ + # Update sigma and tau based on the strong convexity of G + if self.gamma_g is not None: + algorithm._theta = 1.0 / \ + np.sqrt(1 + 2 * self.gamma_g * algorithm._tau) + self.tau *= algorithm._theta + self.sigma /= algorithm._theta + + # Update sigma and tau based on the strong convexity of F + # Following operations are reversed due to symmetry, sigma --> tau, tau -->sigma + if self.gamma_fconj is not None: + algorithm._theta = 1.0 / \ + np.sqrt(1 + 2 * self.gamma_fconj * algorithm._sigma) + self.sigma *= algorithm._theta + self.tau /= algorithm._theta + + return self.tau, self.sigma + + +class PDHGAdaptiveStepSize2013(StepSizeRule): + r"""Adaptively updates the PDHG primal and dual step sizes using the backtracking and residual-balancing method of :cite:`goldstein2013adaptive`. + + This is a PDHG-compatible step-size rule: it provides an initial pair of step sizes via + :meth:`get_initial_step_size` and then, at the end of every PDHG iteration, + :meth:`get_step_size` updates :math:`\tau` and :math:`\sigma` in place based on the + observed behaviour of the algorithm. It combines two mechanisms applied each iteration. + + **1. Backtracking.** Using the primal and dual increments + :math:`\Delta x = x^{n+1} - x^{n}` and :math:`\Delta y = y^{n+1} - y^{n}`, a + backtracking quantity is formed + + .. math:: + + b = \frac{2\,\sigma\tau\,|\langle \Delta y, K\Delta x\rangle|} + {\gamma\,\sigma\,\|\Delta x\|^{2} + \gamma\,\tau\,\|\Delta y\|^{2}}. + + While :math:`b > 1` the step is unstable, so both step sizes are shrunk, + :math:`\tau \leftarrow \tfrac{\beta}{b}\tau`, :math:`\sigma \leftarrow \tfrac{\beta}{b}\sigma`, + the PDHG update is recomputed (via :meth:`~cil.optimisation.algorithms.PDHG._pdhg_update`) + and :math:`b` is re-evaluated. The loop repeats until :math:`b \le 1` or + ``inner_iterations`` is reached. + + **2. Residual balancing.** The primal and dual residuals + + .. math:: + + p = \left\| \tfrac{\Delta x}{\tau} - K^{*}\Delta y \right\|, \qquad + d = \left\| \tfrac{\Delta y}{\sigma} - K\Delta x \right\| + + are compared using the balancing scale :math:`s` and band parameter :math:`\delta > 1`. + If :math:`p < \tfrac{s}{\delta}\, d` the ratio is tilted towards the dual variable + (:math:`\tau \leftarrow (1-\alpha)\tau`, :math:`\sigma \leftarrow \sigma/(1-\alpha)`); + if :math:`s\,\delta\, d < p` it is tilted towards the primal variable + (:math:`\tau \leftarrow \tau/(1-\alpha)`, :math:`\sigma \leftarrow (1-\alpha)\sigma`). + Each time a tilt is applied the adaptation strength decays, :math:`\alpha \leftarrow \eta\,\alpha`, + so the updates become progressively smaller. The product :math:`\sigma\tau` is left + unchanged by the balancing step, preserving the stability guaranteed by the backtracking step. + + Parameters + ------------- + initial_step_size : list of two positive :obj:`float`, optional, default=[10/algorithm.operator.norm(), 10/algorithm.operator.norm()] + Initial values of the primal and dual step sizes used in the adaptive step size method. + initial_alpha : positive :obj:`float`, optional, default=0.95 + Initial value of the adaptation strength :math:`\alpha` controlling the size of the residual-balancing update. + beta : positive :obj:`float`, optional, default=0.95 + The factor :math:`\beta` by which the step sizes are shrunk during backtracking (via :math:`\beta/b`). + gamma : positive :obj:`float`, optional, default=0.9 + The convergence constant :math:`\gamma` appearing in the denominator of the backtracking quantity :math:`b`. + delta : positive :obj:`float`, greater than one, optional, default=1.5 + The band parameter :math:`\delta` setting how far apart the primal and dual residuals are allowed to drift before rebalancing. + s : positive :obj:`float`, optional, default= Norm of the operator A + The balancing scale :math:`s` used to compare the primal and dual residuals. Defaults to the operator norm :math:`\|K\|`. + eta : positive :obj:`float`, optional, default=0.95 + The decay factor :math:`\eta \in (0,1)` applied to :math:`\alpha` each time the step sizes are rebalanced. + inner_iterations : :obj:`int`, optional, default=50 + The maximum number of inner iterations for the backtracking loop. + auto_stop : :obj:`boolean`, optional, default=True + If True, the adaptive step size method automatically stops updating the step sizes when they have not changed for ``auto_stop_patience`` consecutive iterations. + auto_stop_patience : :obj:`int`, optional, default=10 + Number of consecutive iterations with no change to the step sizes after which the adaptive updates are stopped (only used when ``auto_stop=True``). + + + Notes + ----- + This method is memory expensive, requiring the storage of 2 extra image copies and 2 extra data copies. When ``auto_stop=True`` the adaptive updates are switched off once the step sizes have been unchanged for ``auto_stop_patience`` consecutive iterations; the extra images and data are then released, reducing ongoing memory use and computational cost. For a more time expensive, but less memory expensive method, see :class:`PDHGBayesOptimisationStepSize` which does not require the storage of extra images and data. + + See Also + -------- + PDHGAdaptiveStepSize2015 : A related backtracking/residual-balancing rule from :cite:`Goldstein2015`. + + Reference + --------- + Goldstein, T., Li, M., Yuan, X., Esser, E. and Baraniuk, R., 2013. Adaptive primal-dual hybrid gradient methods for saddle-point problems. arXiv preprint arXiv:1305.0546. :cite:`goldstein2013adaptive` + """ + + def __init__(self, initial_step_size=[None, None], initial_alpha=0.95, beta=0.95, gamma=0.9, delta=1.5, s=None, eta=0.95, inner_iterations=50, auto_stop=True, auto_stop_patience=10): + '''Initialises the step size rule''' + self.alpha = initial_alpha + self.eta = eta + self.beta = beta + self.delta = delta + self.s = s + self.gamma = gamma + self.tolerance = 1e-6 + self.p_norm = 100 + self.d_norm = 100 + self.inner_iterations = inner_iterations + self.auto_stop = auto_stop + self.auto_stop_patience = auto_stop_patience + self.count = 0 + + self.y_old = None + self.x_resid = None + self.y_resid = None + + self.adaptive = True + self.initial_step_size = initial_step_size + if len(initial_step_size) != 2: + raise ValueError( + "initial_step_size should be a list or tuple of length two, step_size = {}".format(initial_step_size)) + + def get_initial_step_size(self, algorithm): + tau = self.initial_step_size[0] + sigma = self.initial_step_size[1] + if tau is None: + tau = 10/algorithm.operator.norm() + if sigma is None: + sigma = 10/algorithm.operator.norm() + return tau, sigma + + def get_step_size(self, algorithm): + if self.adaptive: + if self.s is None: + self.s = algorithm.operator.norm() # default balancing scale, ||A|| + if self.y_old is None: + self.y_old = algorithm.operator.range_geometry().allocate(0) # Extra range data 1 + self.x_resid = algorithm.operator.domain_geometry().allocate(0) # Extra image 1 + self.y_resid = algorithm.operator.range_geometry().allocate(0) # Extra range data 2 + # adaptive step sizes only when above tolerance + if self.p_norm > self.tolerance and self.d_norm > self.tolerance: + log.debug('Before adaptive step-size step, tau = {}, sigma = {}'.format( + algorithm._tau, algorithm._sigma)) + b = self._calculate_backtracking(algorithm) + converged = False + for k in range(self.inner_iterations): + if b <= 1: + log.debug('Finished backtracking step, backtracking value b = {}, step sizes are tau = {}, sigma = {}'.format( + b, algorithm._tau, algorithm._sigma)) + converged = True + break + + algorithm._tau *= self.beta/b + algorithm._sigma *= self.beta/b + log.debug(' Backtracking step - multiplying primal and dual step sizes by beta/b = {}, new step sizes are tau = {}, sigma ={}'.format( + self.beta / b, algorithm._tau, algorithm._sigma)) + + algorithm._pdhg_update() + b = self._calculate_backtracking(algorithm) + self.count = 0 + + if not converged: + log.warning('Backtracking step did not converge after {} iterations, backtracking value b = {}, step sizes are tau = {}, sigma = {}'.format( + self.inner_iterations, b, algorithm._tau, algorithm._sigma)) + + self._calculate_pnorm_dnorm(algorithm) + log.debug('Started the rebalancing step with p_norm = {}, d_norm = {}'.format( + self.p_norm, self.d_norm)) + if self.p_norm < (self.s/self.delta)*self.d_norm: + algorithm._tau *= (1 - self.alpha) + algorithm._sigma /= (1 - self.alpha) + self.alpha *= self.eta + self.count = 0 + log.debug('p_norm < (s*delta)*d_norm so rebalancing step sizes, new step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + elif (self.s*self.delta)*self.d_norm < self.p_norm: + algorithm._tau /= (1 - self.alpha) + algorithm._sigma *= (1 - self.alpha) + self.alpha *= self.eta + self.count = 0 + log.debug('(s*delta)*p_norm < d_norm so rebalancing step sizes, new step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + else: + log.debug('No change from the rebalancing step, step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + self.count += 1 + else: + log.debug('No change from the rebalancing step as pnorm and dnorm are below threshold, step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + + self.y_old.fill(algorithm.y) + + if self.count > self.auto_stop_patience and self.auto_stop: + self.adaptive = False + log.debug('Automatic stopping of adaptive step size updates, step sizes have not changed for {} iterations, step sizes are tau = {}, sigma ={}'.format( + self.auto_stop_patience, algorithm._tau, algorithm._sigma)) + del self.x_resid + del self.y_resid + del self.y_old + + return algorithm._tau, algorithm._sigma + + def _calculate_pnorm_dnorm(self, algorithm): + algorithm.operator.adjoint(self.y_resid, out=algorithm.x_tmp) + algorithm.operator.direct(self.x_resid, out=algorithm.y_tmp) + self.x_resid.sapyb((1/algorithm._tau), + algorithm.x_tmp, -1.0, out=algorithm.x_tmp) + self.y_resid.sapyb((1/algorithm._sigma), + algorithm.y_tmp, -1.0, out=algorithm.y_tmp) + self.p_norm = algorithm.x_tmp.norm() + self.d_norm = algorithm.y_tmp.norm() + + def _calculate_backtracking(self, algorithm): + """ Calculates the backtracking parameter b used to update step sizes in the adaptive PDHG algorithm. + Returns + ------- + b : :obj:`float` + Backtracking parameter used to update step sizes in the adaptive PDHG algorithm. + """ + + algorithm.x.sapyb(1.0, algorithm.x_old, -1.0, out=self.x_resid) + x_change_norm = self.x_resid.norm() + algorithm.y.sapyb(1.0, self.y_old, -1.0, out=self.y_resid) + y_change_norm = self.y_resid.norm() + algorithm.operator.direct(self.x_resid, out=algorithm.y_tmp) + cross_term = np.abs(2*algorithm._sigma * + algorithm._tau*self.y_resid.dot(algorithm.y_tmp)) + + denominator = (self.gamma*algorithm._sigma)*x_change_norm**2 + \ + (self.gamma*algorithm._tau)*y_change_norm**2 + if denominator == 0: + # No change in the iterate (x == x_old and y == y_old): nothing to + # backtrack. Return 0 so the caller's `b <= 1` test accepts the step + # (dividing here would give 0/0 = nan and poison the step sizes). + b = 0.0 + else: + b = cross_term/denominator + log.debug('Backtracking value = {}'.format(b)) + return b + + +class PDHGAdaptiveStepSize2015(StepSizeRule): + r"""Adaptively updates the PDHG primal and dual step sizes using the backtracking and residual-balancing method of :cite:`Goldstein2015`. + + This is a PDHG-compatible step-size rule and a variant of :class:`PDHGAdaptiveStepSize2013`. + It provides an initial pair of step sizes via :meth:`get_initial_step_size` and then, at the + end of every PDHG iteration, :meth:`get_step_size` updates :math:`\tau` and :math:`\sigma` + in place. Two mechanisms are applied each iteration. + + **1. Backtracking.** With the increments :math:`\Delta x = x^{n+1} - x^{n}` and + :math:`\Delta y = y^{n+1} - y^{n}`, a backtracking quantity is formed + + .. math:: + + b = c\,\sigma\,\|\Delta x\|^{2} + c\,\tau\,\|\Delta y\|^{2} + - 4\,\sigma\tau\,|\langle \Delta y, K\Delta x\rangle|. + + The step is accepted when :math:`b \ge 0`. Otherwise both step sizes are halved, + :math:`\tau \leftarrow \tfrac{1}{2}\tau`, :math:`\sigma \leftarrow \tfrac{1}{2}\sigma`, the + PDHG update is recomputed (via :meth:`~cil.optimisation.algorithms.PDHG._pdhg_update`) and + :math:`b` is re-evaluated, up to ``inner_iterations`` times. + + **2. Residual balancing.** The primal and dual residuals + + .. math:: + + p = \left\| \tfrac{\Delta x}{\tau} - K^{*}\Delta y \right\|, \qquad + d = \|K\| \left\| \tfrac{\Delta y}{\sigma} - K\Delta x \right\| + + are compared. If :math:`2p < d` the ratio is tilted towards the dual variable + (:math:`\tau \leftarrow (1-\alpha)\tau`, :math:`\sigma \leftarrow \sigma/(1-\alpha)`); + if :math:`2d < p` it is tilted towards the primal variable + (:math:`\tau \leftarrow \tau/(1-\alpha)`, :math:`\sigma \leftarrow (1-\alpha)\sigma`). + Whenever a tilt is applied the adaptation strength decays, :math:`\alpha \leftarrow \eta\,\alpha`. + The balancing step preserves the product :math:`\sigma\tau`. + + Parameters + ------------- + initial_step_size : list of two positive :obj:`float`, optional, default= [10/algorithm.operator.norm(), 10/algorithm.operator.norm()] + Initial values of the primal and dual step sizes used in the adaptive step size method. + initial_alpha : positive :obj:`float`, optional, default=0.95 + Initial value of the adaptation strength :math:`\alpha` controlling the size of the residual-balancing update. + eta : positive :obj:`float`, optional, default=0.95 + The decay factor :math:`\eta \in (0,1)` applied to :math:`\alpha` each time the step sizes are rebalanced. + c : positive :obj:`float`, optional, default=0.9 + The convergence constant :math:`c` appearing in the backtracking quantity :math:`b`. + inner_iterations : :obj:`int`, optional, default=50 + The maximum number of inner iterations for the backtracking loop. + auto_stop : :obj:`boolean`, optional, default=True + If True, the adaptive step size method automatically stops updating the step sizes when they have not changed for ``auto_stop_patience`` consecutive iterations. + auto_stop_patience : :obj:`int`, optional, default=10 + Number of consecutive iterations with no change to the step sizes after which the adaptive updates are stopped (only used when ``auto_stop=True``). + + + Notes + ----- + This method is memory expensive, requiring the storage of 2 extra image copies and 2 extra data copies. When ``auto_stop=True`` the adaptive updates are switched off once the step sizes have been unchanged for ``auto_stop_patience`` consecutive iterations; the extra images and data are then released, reducing ongoing memory use and computational cost. For a more time expensive, but less memory expensive method, see :class:`PDHGBayesOptimisationStepSize` which does not require the storage of extra images and data. + + See Also + -------- + PDHGAdaptiveStepSize2013 : The related backtracking/residual-balancing rule from :cite:`goldstein2013adaptive`. + + Reference + --------- + Goldstein, T., Li, M. and Yuan, X., 2015. Adaptive primal-dual splitting methods for statistical learning and image processing. Advances in Neural Information Processing Systems, 28. :cite:`Goldstein2015` + """ + + def __init__(self, initial_step_size=[None, None], initial_alpha=0.95, eta=0.95, c=0.9, inner_iterations=50, auto_stop=True, auto_stop_patience=10): + '''Initialises the step size rule''' + + self.adaptive = True + self.alpha = initial_alpha + self.eta = eta + self.c = c + self.tolerance = 1e-6 + self.p_norm = 100 + self.d_norm = 100 + self.auto_stop = auto_stop + self.auto_stop_patience = auto_stop_patience + self.count = 0 + self.inner_iterations = inner_iterations + + self.y_old = None + self.x_resid = None + self.y_resid = None + self.initial_step_size = initial_step_size + if len(initial_step_size) != 2: + raise ValueError( + "initial_step_size should be a list or tuple of length two, step_size = {}".format(initial_step_size)) + + def get_initial_step_size(self, algorithm): + tau = self.initial_step_size[0] + sigma = self.initial_step_size[1] + if tau is None: + tau = 10/algorithm.operator.norm() + if sigma is None: + sigma = 10/algorithm.operator.norm() + return tau, sigma + + def get_step_size(self, algorithm): + if self.adaptive: + if self.y_old is None: + self.y_old = algorithm.operator.range_geometry().allocate(0) # Extra range data 1 + self.x_resid = algorithm.operator.domain_geometry().allocate(0) # Extra image 1 + self.y_resid = algorithm.operator.range_geometry().allocate(0) # Extra range data 2 + if self.p_norm > self.tolerance and self.d_norm > self.tolerance: + log.debug('Before adaptive step-size step, tau = {}, sigma = {}'.format( + algorithm._tau, algorithm._sigma)) + + b = self._calculate_backtracking(algorithm) + converged = False + for k in range(self.inner_iterations): + if b >= 0: + log.debug('Finished backtracking step, backtracking value b = {}, step sizes are tau = {}, sigma = {}'.format( + b, algorithm._tau, algorithm._sigma)) + converged = True + break + algorithm._tau *= 0.5 + algorithm._sigma *= 0.5 + log.debug(' Backtracking step - multiplying primal and dual step sizes by 1/2, new step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + + algorithm._pdhg_update() + b = self._calculate_backtracking(algorithm) + self.count = 0 + if not converged: + log.warning('Backtracking step did not converge after {} iterations, backtracking value b = {}, step sizes are tau = {}, sigma = {}'.format( + self.inner_iterations, b, algorithm._tau, algorithm._sigma)) + + + self._calculate_pnorm_dnorm(algorithm) + log.debug('Started the rebalancing step with p_norm = {}, d_norm = {}'.format( + self.p_norm, self.d_norm)) + if 2*self.p_norm < self.d_norm: + algorithm._tau *= (1 - self.alpha) + algorithm._sigma /= (1 - self.alpha) + self.alpha *= self.eta + self.count = 0 + log.debug('2*p_norm < d_norm so rebalancing step sizes, new step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + elif 2*self.d_norm < self.p_norm: + algorithm._tau /= (1 - self.alpha) + algorithm._sigma *= (1 - self.alpha) + self.alpha *= self.eta + self.count = 0 + log.debug('2*d_norm < p_norm so rebalancing step sizes, new step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + else: + log.debug('No change from the rebalancing step, step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + self.count += 1 + else: + log.debug('No change from the rebalancing step as pnorm and dnorm are below threshold, step sizes are tau = {}, sigma ={}'.format( + algorithm._tau, algorithm._sigma)) + self.y_old.fill(algorithm.y) + + if self.count > self.auto_stop_patience and self.auto_stop: + self.adaptive = False + log.debug('Automatic stopping of adaptive step size updates, step sizes have not changed for {} iterations, step sizes are tau = {}, sigma ={}'.format( + self.auto_stop_patience, algorithm._tau, algorithm._sigma)) + del self.y_resid + del self.y_old + del self.x_resid + + return algorithm._tau, algorithm._sigma + + def _calculate_backtracking(self, algorithm): + """ Calculates the backtracking parameter b used to update step sizes in the adaptive PDHG algorithm. + Returns + ------- + b : :obj:`float` + Backtracking parameter used to update step sizes in the adaptive PDHG algorithm. + """ + + algorithm.x.sapyb(1.0, algorithm.x_old, -1.0, out=self.x_resid) + x_change_norm = self.x_resid.norm() + algorithm.y.sapyb(1.0, self.y_old, -1.0, out=self.y_resid) + y_change_norm = self.y_resid.norm() + algorithm.operator.direct(self.x_resid, out=algorithm.y_tmp) + cross_term = np.abs(4*algorithm._sigma*algorithm._tau * + self.y_resid.dot(algorithm.y_tmp)) + b = self.c*algorithm._sigma*x_change_norm**2 + \ + self.c*algorithm._tau*y_change_norm**2 - cross_term + log.debug('Backtracking value = {}'.format(b)) + return b + + def _calculate_pnorm_dnorm(self, algorithm): + """Calculates the primal and dual norms used in the rebalancing step of the adaptive PDHG algorithm. + """ + algorithm.operator.adjoint(self.y_resid, out=algorithm.x_tmp) + algorithm.operator.direct(self.x_resid, out=algorithm.y_tmp) + self.x_resid.sapyb((1/algorithm._tau), + algorithm.x_tmp, -1.0, out=algorithm.x_tmp) + self.y_resid.sapyb((1/algorithm._sigma), + algorithm.y_tmp, -1.0, out=algorithm.y_tmp) + self.p_norm = algorithm.x_tmp.norm() + self.d_norm = algorithm.operator.norm()*algorithm.y_tmp.norm() + + +class _BayesOptimisationStepSizeBase(StepSizeRule): + """Shared implementation for the Bayesian-optimisation step-size rules. + + A Gaussian-process Bayesian optimisation (via ``skopt.gp_minimize``) chooses the + scalar ratio ``gamma`` between the primal and dual step sizes that gives the best + objective after a small number of iterations. The step sizes are then fixed for the + remainder of the algorithm. + + Subclasses implement :meth:`_step_sizes_from_gamma` (the mapping from ``gamma`` to + the ``(tau, sigma)`` step sizes) and may override :meth:`_default_n_iterations`. + """ + + def __init__(self, gamma_bounds=None, n_initial_points=5, n_calls=20, n_iterations=None, seed=None, plot=False): + '''Initialises the step size rule''' + self.gamma_bounds = gamma_bounds + if gamma_bounds is not None: + if len(gamma_bounds) != 2: + raise ValueError( + "gamma_bounds should be a list or tuple of length two, gamma_bounds = {}".format(gamma_bounds)) + if gamma_bounds[0] <= 0 or gamma_bounds[1] <= 0: + raise ValueError( + "gamma_bounds should be positive and strictly greater than zero, gamma_bounds = {}".format(gamma_bounds)) + + self.n_initial_points = n_initial_points + self.n_calls = n_calls + self.n_iterations = n_iterations + self.seed = seed + self.plot = plot + + def _default_n_iterations(self, algorithm): + """Number of iterations per objective evaluation when ``n_iterations`` is None.""" + return 10 + + def _step_sizes_from_gamma(self, algorithm, gamma): + """Map the scalar ratio ``gamma`` to ``(tau, sigma)``. Implemented by subclasses.""" + raise NotImplementedError + + def _algorithm_set_up(self, algorithm, step_size): + """Set up the algorithm with the given step size.""" + raise NotImplementedError + + def get_initial_step_size(self, algorithm): + try: + from skopt import gp_minimize + except ImportError: + raise ImportError( + "skopt is required for the Bayesian optimisation step size rule. Please install scikit-optimize to use this step size rule.") + + if self.n_iterations is None: + self.n_iterations = self._default_n_iterations(algorithm) + log.debug( + "n_iterations not provided, set to {}".format(self.n_iterations)) + + if self.gamma_bounds is None: + ratio = np.sqrt(algorithm.f(algorithm.operator.direct( + 0*algorithm.x)))/algorithm.operator.norm() + log.debug("ratio: {}".format(ratio)) + self.gamma_bounds = (1e-5/ratio, 1e5/ratio) + log.debug( + "gamma_bounds not provided, set to (1e-5/ratio, 1e5/ratio) = {}".format(self.gamma_bounds)) + + log_gamma_bounds = ( + np.log(self.gamma_bounds[0]), np.log(self.gamma_bounds[1])) + update_objective_interval = algorithm.update_objective_interval + + + + def objective_function(log_gamma): + gamma = np.exp(log_gamma[0]) + log.debug( + "Evaluating objective function for gamma = {}".format(gamma)) + # Set the step sizes based on the current gamma + tau, sigma = self._step_sizes_from_gamma(algorithm, gamma) + + self._algorithm_set_up(algorithm, [tau,sigma]) + + try: + with np.errstate(over='raise'): + algorithm.run(self.n_iterations, callbacks=[]) + except FloatingPointError: + log.debug( + "Objective function diverged for gamma = {}".format(gamma)) + return 1e10 # large penalty + + log.debug("Objective function value for gamma = {}: {}".format( + gamma, algorithm.objective[-1])) + if not np.isfinite(algorithm.objective[-1]): + return 1e10 # large penalty + return algorithm.objective[-1] + + + algorithm.update_objective_interval = self.n_iterations-1 + gp_result = gp_minimize(objective_function, [ + log_gamma_bounds], n_random_starts=self.n_initial_points, n_calls=self.n_calls, initial_point_generator="lhs", random_state=self.seed) + + algorithm.update_objective_interval = update_objective_interval + + + # gp_result.x[0] is log(gamma) as the optimisation was over log_gamma_bounds + self.tau, self.sigma = self._step_sizes_from_gamma( + algorithm, np.exp(gp_result.x[0])) + + log.debug("Best gamma found: {}, with objective function value: {}".format( + np.exp(gp_result.x[0]), gp_result.fun)) + log.debug('Initial step sizes are tau = {}, sigma = {}'.format( + self.tau, self.sigma)) + + self._algorithm_set_up(algorithm, [self.tau, self.sigma]) + + if self.plot: + from matplotlib import pyplot as plt + from skopt.plots import plot_convergence, plot_gaussian_process + plot_convergence(gp_result) + plt.show() + plot_gaussian_process(gp_result) + plt.yscale('log') + plt.xlabel('log gamma') + plt.title("Best gamma found: {}, with objective function value: {}".format( + np.exp(gp_result.x[0]), gp_result.fun)) + plt.show() + + return self.tau, self.sigma + + def get_step_size(self, algorithm): + log.debug('Returning step sizes tau = {}, sigma = {}'.format( + self.tau, self.sigma)) + return self.tau, self.sigma + + +class PDHGBayesOptimisationStepSize(_BayesOptimisationStepSizeBase): + r"""The ratio between the primal and dual step sizes (gamma) in the PDHG algorithm is chosen using a guassian process Bayesian optimisation, choosing the gamma that gives the best performance after a small number of iterations. The step sizes are chosen at the beginning of the algorithm and then kept constant throughout the iterations. + Parameters + ------------- + gamma_bounds : list or tuple of length two, optional, the default is an approximation of [1e-5, 1e6]*norm(A)/norm(b) where A is the operator and b is the data, which is a good initial guess for the ratio between the primal and dual step sizes in the PDHG algorithm. + Bounds for the ratio between the primal and dual step sizes (gamma) in the Bayesian optimisation. The gamma that gives the best performance after a small number of iterations is chosen as the ratio between the primal and dual step sizes for the PDHG algorithm. The default bounds are (1e-5, 1e5). + n_initial_points : int, optional, default=5 + Number of initial random evaluations of the objective function in the Bayesian optimisation. + n_calls : int, optional, default=20 + Total number of evaluations of the objective function in the Bayesian optimisation, including the initial random evaluations. + n_iterations : int, optional, default=10 + Number of iterations to run the PDHG algorithm for each evaluation of the objective function in the Bayesian optimisation. The gamma that gives the best performance after this number of iterations is chosen as the ratio between the primal and dual step sizes for the PDHG algorithm. + seed : int, optional, default= None + Random seed for the Bayesian optimisation. This is used to ensure reproducibility of the results. + plot : bool, optional, default=False + If True, plots the convergence of the Bayesian optimisation and the fitted Gaussian process after the step sizes have been chosen. Requires matplotlib and blocks on ``plt.show()``, so it should be left False for headless/non-interactive runs. + + Notes + ----- + The step sizes are parametrised by a single scalar ratio :math:`\gamma` and the operator + norm :math:`\|K\|`, + + .. math:: + + \tau = \frac{1}{\gamma \|K\|}, \qquad \sigma = \frac{\gamma}{\|K\|}, + + so that the product :math:`\sigma\tau = 1/\|K\|^{2}` is fixed and only the primal/dual + balance is optimised. The search is performed over :math:`\log\gamma` (which spans + several orders of magnitude) using ``skopt.gp_minimize``, so `scikit-optimize` must be + installed to use this rule. + + This is a computationally expensive step size rule, as it requires running the PDHG algorithm for a number of iterations for each evaluation of the objective function in the Bayesian optimisation. It is recommended to use this step size rule where you are memory constrained but less time constrained. For the opposite case, where you are more time constrained, not memory constrained, we recommend using the :class:`PDHGAdaptiveStepSize2013` or :class:`PDHGAdaptiveStepSize2015` step size rules, which are adaptive step size rules that update the step sizes at each iteration based on the observed behaviour of the algorithm. + """ + + def __init__(self, gamma_bounds=None, n_initial_points=5, n_calls=20, n_iterations=10, seed=None, plot=False): + '''Initialises the step size rule''' + super().__init__(gamma_bounds=gamma_bounds, n_initial_points=n_initial_points, + n_calls=n_calls, n_iterations=n_iterations, seed=seed, plot=plot) + + def _step_sizes_from_gamma(self, algorithm, gamma): + norm = algorithm.operator.norm() + tau = 1.0 / (gamma * norm) + sigma = 1.0 * gamma / norm + return tau, sigma + + def _algorithm_set_up(self, algorithm, step_size): + """Set up the algorithm with the given step size.""" + algorithm.set_up(initial=algorithm.initial, f=algorithm.f, + g=algorithm.g, operator=algorithm.operator, step_size=step_size) + algorithm._reset_iteration_state() + + + +class PDHGConstantStepSize(StepSizeRule): + r""" + Step-size rule that always returns a constant step-size. + + The user can set either the primal or dual step size, both or none. + + By default, the step sizes :math:`\sigma` and :math:`\tau` are positive scalars and defined as below: + + * If ``sigma`` is ``None`` and ``tau`` is ``None``: + + .. math:: + + \sigma = \frac{1}{\|K\|}, \tau = \frac{1}{\|K\|} + + * If ``tau`` is ``None``: + + .. math:: + + \tau = \frac{1}{\sigma\|K\|^{2}} + + * If ``sigma`` is ``None``: + + .. math:: + + \sigma = \frac{1}{\tau\|K\|^{2}} + + + Parameters + ---------- + step_size : list or tuple of length two, default=[None, None] + Initial values of the primal and dual step sizes. If both are ``None`` they are set to the default values defined above. If one is ``None`` it is calculated based on the other and the norm of the operator. If both are provided, they are used as they are, as long as they are positive numbers. + """ + + def __init__(self, step_size=[None, None]): + '''Initialises the constant step size rule''' + + if len(step_size) != 2: + raise ValueError( + "step_size should be a list or tuple of length two, step_size = {}".format(step_size)) + self.tau = step_size[0] + self.sigma = step_size[1] + + def get_initial_step_size(self, algorithm): + """Sets sigma and tau step-sizes for the PDHG algorithm.""" + self.tau, self.sigma = _resolve_pdhg_step_sizes( + self.tau, self.sigma, algorithm.operator) + return self.tau, self.sigma + + def get_step_size(self, algorithm): + """ + Returns + -------- + the primal and dual step sizes as a tuple ``(tau, sigma)`` + """ + return self.tau, self.sigma + + +class SPDHGConstantStepSize(StepSizeRule): + r"""Step-size rule that always returns a constant step-size for the SPDHG algorithm. + The user can set either the primal or dual step size, both or none. + + When setting `sigma` and `tau`, there are 4 possible cases considered by setup function: + + - Case 1: If neither `sigma` or `tau` are provided then `sigma` is set using the formula: + + .. math:: \sigma_i= \frac{0.99}{\|K_i\|} + + and `tau` is set as per case 2 + + - Case 2: If `sigma` is provided but not `tau` then `tau` is calculated using the formula + + .. math:: \tau = 0.99\min_i( \frac{p_i}{ (\sigma_i \|K_i\|^2) }) + + - Case 3: If `tau` is provided but not `sigma` then `sigma` is calculated using the formula + + .. math:: \sigma_i= \frac{0.99 p_i}{\tau\|K_i\|^2} + + - Case 4: Both `sigma` and `tau` are provided. + + Parameters + ---------- + step_size : list or tuple of length two, default=[None, None] + Initial values of the primal and dual step sizes. If both are ``None`` they are set to the default values defined below. If one is ``None`` it is calculated based on the other and the norm of the operator. If both are provided, they are used as they are, as long as sigma is a list or array of positive numbers of length equal to the number of operators and tau is a positive number. + """ + + def __init__(self, step_size=[None, None]): + '''Initialises the constant step size rule + ''' + + if len(step_size) != 2: + raise ValueError( + "step_size should be a list or tuple of length two, step_size = {}".format(step_size)) + self.tau = step_size[0] + self.sigma = step_size[1] + + + def get_initial_step_size(self, algorithm): + r""" Sets sigma and tau step-sizes for the SPDHG algorithm after the initial set-up. The step sizes can be either scalar or array-objects. + """ + gamma = 1. + rho = .99 + if self.sigma is not None: + + if not isinstance(self.sigma, Number) and len(self.sigma) == algorithm._ndual_subsets: + if all(isinstance(x, Number) and x > 0 for x in self.sigma): + pass + else: + raise ValueError( + "Sigma expected to be a positive number.") + + else: + raise ValueError( + "Please pass a list of floats to sigma with the same number of entries as number of operators") + + elif self.tau is None: + self.sigma = _spdhg_sigma_from_gamma(gamma, rho, algorithm._norms) + else: + self.sigma = [ + rho*pi / (self.tau*ni**2) for ni, pi in zip(algorithm._norms, algorithm._prob_weights)] + + if self.tau is None: + self.tau = _spdhg_tau_from_sigma( + self.sigma, algorithm._norms, algorithm._prob_weights, rho) + + else: + if not (isinstance(self.tau, Number) and self.tau > 0): + raise ValueError( + "The step-sizes of SPDHG must be positive, passed tau = {}".format(self.tau)) + + return self.tau, self.sigma + + def get_step_size(self, algorithm): + """ + Returns + -------- + the primal step size and the list of dual step sizes as a tuple ``(tau, sigma)`` + """ + return self.tau, self.sigma + + +class SPDHGStepSizesFromRatio(StepSizeRule): + r""" Sets gamma, the step-size ratio for the SPDHG algorithm. Currently gamma takes a scalar value. + + The step sizes `sigma` and `tau` are set using the equations: + + .. math:: \sigma_i= \frac{\gamma\rho }{\|K_i\|} + + .. math:: \tau = \rho\min_i\left( \frac{p_i }{\sigma_i \|K_i\|^2}\right) + + where :math:`p_i` is the sampling probability of the :math:`i`-th operator. The dual step + size :math:`\sigma` is therefore a list with one entry per operator, while the primal step + size :math:`\tau` is a scalar. + + Parameters + ---------- + gamma : Positive float + parameter controlling the trade-off between the primal and dual step sizes + rho : Positive float + parameter controlling the size of the product :math:`\sigma\tau` + + + + """ + def __init__(self, gamma, rho): + """Initialises the step size rule""" + self.gamma = gamma + self.rho = rho + + def get_initial_step_size(self, algorithm): + if isinstance(self.gamma, Number): + if self.gamma <= 0: + raise ValueError( + "The step-sizes of SPDHG are positive, gamma should also be positive") + + else: + raise ValueError( + "We currently only support scalar values of gamma") + if isinstance(self.rho, Number): + if self.rho <= 0: + raise ValueError( + "The step-sizes of SPDHG are positive, rho should also be positive") + + else: + raise ValueError( + "We currently only support scalar values of gamma") + + self.sigma = _spdhg_sigma_from_gamma( + self.gamma, self.rho, algorithm._norms) + self.tau = _spdhg_tau_from_sigma( + self.sigma, algorithm._norms, algorithm._prob_weights, self.rho) + + return self.tau, self.sigma + + def get_step_size(self, algorithm): + """ + Returns + -------- + the primal step size and the list of dual step sizes as a tuple ``(tau, sigma)`` + """ + return self.tau, self.sigma + + +class SPDHGBayesOptimisationStepSize(_BayesOptimisationStepSizeBase): + r"""The ratio between the primal and dual step sizes (gamma) in the SPDHG algorithm is chosen using a guassian process Bayesian optimisation, choosing the gamma that gives the best performance after a small number of iterations. The step sizes are chosen at the beginning of the algorithm and then kept constant throughout the iterations. + Parameters + ------------- + gamma_bounds : list or tuple of length two, optional, the default is an approximation of [1e-5, 1e6]*norm(A)/norm(b) where A is the operator and b is the data, which is a good initial guess for the ratio between the primal and dual step sizes in the SPDHG algorithm. + Bounds for the ratio between the primal and dual step sizes (gamma) in the Bayesian optimisation. The gamma that gives the best performance after a small number of iterations is chosen as the ratio between the primal and dual step sizes for the SPDHG algorithm. The default bounds are (1e-5, 1e5). + n_initial_points : int, optional, default=5 + Number of initial random evaluations of the objective function in the Bayesian optimisation. + n_calls : int, optional, default=20 + Total number of evaluations of the objective function in the Bayesian optimisation, including the initial random evaluations. + n_iterations : int, optional, default=None + Number of iterations to run the SPDHG algorithm for each evaluation of the objective function in the Bayesian optimisation. The gamma that gives the best performance after this number of iterations is chosen as the ratio between the primal and dual step sizes for the SPDHG algorithm. If None, set to be 10*number of operators in the SPDHG algorithm. + seed : int, optional, default= None + Random seed for the Bayesian optimisation. This is used to ensure reproducibility of the results. + plot : bool, optional, default=False + If True, plots the convergence of the Bayesian optimisation and the fitted Gaussian process after the step sizes have been chosen. Requires matplotlib and blocks on ``plt.show()``, so it should be left False for headless/non-interactive runs. + + Notes + ----- + The step sizes are parametrised by a single scalar ratio :math:`\gamma` together with + the fixed product parameter :math:`\rho = 0.99`, the operator norms :math:`\|K_i\|` and + the sampling probabilities :math:`p_i`, + + .. math:: + + \sigma_i = \frac{\gamma \rho}{\|K_i\|}, \qquad + \tau = \rho \min_i \frac{p_i}{\sigma_i \|K_i\|^{2}}, + + so only the primal/dual balance is optimised. The search is performed over + :math:`\log\gamma` using ``skopt.gp_minimize``, so `scikit-optimize` must be installed + to use this rule. + + This is a computationally expensive step size rule, as it requires running the SPDHG algorithm for a number of iterations for each evaluation of the objective function in the Bayesian optimisation. It is recommended to use this step size rule where you are memory constrained but less time constrained. + """ + + def __init__(self, gamma_bounds=None, n_initial_points=5, n_calls=20, n_iterations=None, seed=None, plot=False): + '''Initialises the step size rule''' + super().__init__(gamma_bounds=gamma_bounds, n_initial_points=n_initial_points, + n_calls=n_calls, n_iterations=n_iterations, seed=seed, plot=plot) + self.rho = 0.99 + self._pristine_sampler = None + + def _default_n_iterations(self, algorithm): + return 10 * len(algorithm._norms) + + def _step_sizes_from_gamma(self, algorithm, gamma): + sigma = _spdhg_sigma_from_gamma(gamma, self.rho, algorithm._norms) + tau = _spdhg_tau_from_sigma( + sigma, algorithm._norms, algorithm._prob_weights, self.rho) + return tau, sigma + + def _algorithm_set_up(self, algorithm, step_size): + """Set up the algorithm with the given step size. + + The sampler and probability weights are passed back explicitly because + ``SPDHG.set_up`` would otherwise reset them to the uniform default. Each trial + is given an identical copy of the caller's sampler, so that every gamma is + scored on the same sequence of subsets rather than on wherever the previous + trial happened to leave the generator. + """ + if self._pristine_sampler is None: + self._pristine_sampler = deepcopy(algorithm._sampler) + + algorithm.set_up(initial=algorithm.initial, f=algorithm.f, g=algorithm.g, + operator=algorithm.operator, step_size=step_size, + sampler=deepcopy(self._pristine_sampler), + prob_weights=algorithm._prob_weights) + algorithm._reset_iteration_state() + +class _SPDHGAdaptiveStepSizeBase(StepSizeRule): + r"""Shared machinery for the adaptive SPDHG step-size rules of :cite:`chambolle2023stochastic`. + + The two concrete rules (:class:`SPDHGAdaptiveStepSizeBalancing` and + :class:`SPDHGAdaptiveStepSizeAngle`) only *rescale* the step sizes between iterations, + keeping the standard extrapolated SPDHG update. Following equation (2.6) of + :cite:`chambolle2023stochastic`, a single balancing scalar :math:`\gamma` is used at + each iteration, + + .. math:: + + \tau \leftarrow \tau/\gamma, \qquad \sigma_i \leftarrow \gamma\,\sigma_i + \quad (i = 1,\dots,n), + + so the same factor is applied to every dual step size and the products + :math:`\tau\sigma_i` are preserved. Both rules are *balancing only* — there is no + backtracking inner loop (SPDHG has no side-effect-free re-runnable update), so the + step sizes are adjusted once per iteration from quantities measured on the subset that + was sampled that iteration. + + This base class is not meant to be used directly. It owns the buffers required to form + the per-subset primal/dual increments and defers the actual decision to + :meth:`_rebalance`, implemented by each subclass. All state is held by the rule, so the + only cooperation required from SPDHG is that it records the sampled subset index in + ``algorithm._index`` (see :meth:`~cil.optimisation.algorithms.SPDHG.update`). + + Notes + ----- + This rule stores one extra image (the previous primal iterate), a copy of the dual + variable, two domain-sized working buffers and one range-sized working buffer. When + ``auto_stop=True`` the adaptive updates are switched off, and the extra storage + released, once either the step sizes have been unchanged for ``auto_stop_patience`` + consecutive iterations or the adaptation strength :math:`\alpha` has decayed below + ``alpha_tolerance``. The second criterion matters because :math:`\alpha` decays + geometrically on every rebalance while ``count`` is reset by the same event, so a rule + that keeps rebalancing by ever smaller amounts would otherwise never trip the patience + counter. + """ + + def __init__(self, initial_step_size=[None, None], initial_alpha=0.95, eta=0.995, + auto_stop=True, auto_stop_patience=None, alpha_tolerance=1e-3): + '''Initialises the shared adaptive SPDHG step size rule''' + if len(initial_step_size) != 2: + raise ValueError( + "initial_step_size should be a list or tuple of length two, " + "initial_step_size = {}".format(initial_step_size)) + self.initial_step_size = initial_step_size + self.alpha = initial_alpha + self.eta = eta + self.auto_stop = auto_stop + self.auto_stop_patience = auto_stop_patience + self.alpha_tolerance = alpha_tolerance + self.tolerance = 1e-6 + self.count = 0 + self.adaptive = True + # buffers allocated in get_initial_step_size + self.x_prev = None + self.y_prev = None + self.adj_tmp = None + self.v_tmp = None + self.forward_tmp = None + + def _resolve_initial(self, algorithm): + '''Resolve the initial (tau, sigma) pair, filling in any missing value.''' + gamma = 1. + rho = 0.99 + tau = self.initial_step_size[0] + sigma = self.initial_step_size[1] + if sigma is not None and isinstance(sigma, Number): + # a single dual step size is broadcast to every operator + sigma = [sigma] * algorithm._ndual_subsets + if sigma is None and tau is None: + sigma = _spdhg_sigma_from_gamma(gamma, rho, algorithm._norms) + tau = _spdhg_tau_from_sigma( + sigma, algorithm._norms, algorithm._prob_weights, rho) + elif sigma is None: + sigma = [rho * pi / (tau * ni**2) + for ni, pi in zip(algorithm._norms, algorithm._prob_weights)] + elif tau is None: + tau = _spdhg_tau_from_sigma( + sigma, algorithm._norms, algorithm._prob_weights, rho) + if self.auto_stop and self.auto_stop_patience is None: + self.auto_stop_patience = 10 * len(algorithm._norms) + return tau, sigma + + def get_initial_step_size(self, algorithm): + tau, sigma = self._resolve_initial(algorithm) + # buffers owned by the rule; seeded from the initial iterates so the first + # increment (x^0 - x^1, y^0 - y^1) is exact. + self.x_prev = algorithm.x.copy() + self.y_prev = algorithm._y_old.copy() + self.adj_tmp = algorithm.operator.domain_geometry().allocate(0) + self.v_tmp = algorithm.operator.domain_geometry().allocate(0) + # one buffer per subset: element i is the range of operator[i], so subclasses + # must index it with the sampled subset rather than use it whole. + self.forward_tmp = algorithm.operator.range_geometry().allocate(0) + return tau, sigma + + def get_step_size(self, algorithm): + if self.adaptive: + i = algorithm._index + p_i = algorithm._prob_weights[i] + tau = algorithm._tau + sigma = algorithm._sigma + + # primal increment Delta x = x^k - x^{k+1} + algorithm.x.sapyb(-1.0, self.x_prev, 1.0, out=self.x_prev) + # dual increment on the sampled subset Delta y_i = y^k_i - y^{k+1}_i + self.y_prev[i].subtract(algorithm._y_old[i], out=self.y_prev[i]) + # back-projected dual increment A_i^* Delta y_i + algorithm.operator[i].adjoint(self.y_prev[i], out=self.adj_tmp) + # primal residual direction q = Delta x / tau - (1/p_i) A_i^* Delta y_i + self.x_prev.sapyb(1.0 / tau, self.adj_tmp, -1.0 / p_i, out=self.v_tmp) + + # Relative test: the increments scale with the problem, so an absolute + # threshold would silence the rule entirely on rescaled data. + if (self.x_prev.norm() > self.tolerance * algorithm.x.norm() + and self.y_prev[i].norm() + > self.tolerance * algorithm._y_old[i].norm()): + tau_new, sigma_new, changed = self._rebalance( + algorithm, i, p_i, tau, sigma) + if not changed: + self.count += 1 + else: + log.debug('SPDHG adaptive step size: increments below relative ' + 'tolerance, no change to step sizes.') + tau_new, sigma_new = tau, sigma + self.count += 1 + + # refresh the stored iterates for the next call (only subset i changed) + self.x_prev.fill(algorithm.x) + self.y_prev[i].fill(algorithm._y_old[i]) + + if self.auto_stop: + # Two ways the rule can be finished: the step sizes have settled inside + # the band, or alpha has decayed so far that a rebalance would no longer + # move them (both rules rescale by 1 -/+ alpha, which tends to 1). + if self.count > self.auto_stop_patience: + self._stop_adapting( + 'step sizes unchanged for {} iterations'.format( + self.auto_stop_patience)) + elif self.alpha < self.alpha_tolerance: + self._stop_adapting( + 'adaptation strength alpha = {} has decayed below ' + 'alpha_tolerance = {}'.format( + self.alpha, self.alpha_tolerance)) + + return tau_new, sigma_new + + return algorithm._tau, algorithm._sigma + + def _stop_adapting(self, reason): + '''Switch off the adaptive updates and release the working buffers.''' + self.adaptive = False + log.debug('SPDHG adaptive step size: automatically stopping updates, ' + '{}.'.format(reason)) + del self.x_prev, self.y_prev, self.adj_tmp, self.v_tmp, self.forward_tmp + + def _rebalance(self, algorithm, i, p_i, tau, sigma): + """Decide and apply the step-size rescaling for one iteration. + + Returns + ------- + tuple + ``(tau, sigma, changed)`` where ``sigma`` is the list of dual step sizes and + ``changed`` is ``True`` if the step sizes were rescaled this iteration. + """ + raise NotImplementedError + + +class SPDHGAdaptiveStepSizeBalancing(_SPDHGAdaptiveStepSizeBase): + r"""Adaptively balances the SPDHG primal and dual step sizes by tracking the primal and dual progress, following rule (a) of :cite:`chambolle2023stochastic`. + + This is an SPDHG-compatible step-size rule (:class:`SPDHGAdaptiveStepSizeBalancing` + implements A-SPDHG rule (a), Algorithm 3.1 of :cite:`chambolle2023stochastic`). It is + the stochastic, minibatch generalisation of the deterministic adaptive PDHG rule of + :cite:`goldstein2013adaptive` (:class:`PDHGAdaptiveStepSize2013`), to which it reduces + when there is a single operator. + + At the end of each iteration, with sampled subset :math:`i`, sampling probability + :math:`p_i`, primal increment :math:`\Delta x = x^{k} - x^{k+1}` and dual increment + :math:`\Delta y_i = y^{k}_i - y^{k+1}_i`, two :math:`\ell_1` "progress" residuals are + formed (equation (3.6) of :cite:`chambolle2023stochastic`), + + .. math:: + + v = \left\| \frac{\Delta x}{\tau} + - \frac{1}{p_i} A_i^{*}\Delta y_i \right\|_1, \qquad + d = \frac{1}{p_i}\left\| \frac{\Delta y_i}{\sigma_i} + - A_i \Delta x \right\|_1 , + + which track the primal and dual subgradients respectively. In both cases the + :math:`1/p_i` weight makes the quantity measured on the sampled subset an unbiased + estimate of the corresponding full (deterministic) residual: for :math:`v` it + reweights the only nonzero block of :math:`A^{*}\Delta y`, and for :math:`d` it + reweights the sampled block's contribution to :math:`\sum_j \|d_j\|_1`. With a single + operator :math:`p_i = 1` and both reduce to the residuals of + :cite:`goldstein2013adaptive`. They are compared using the balancing scale :math:`s` + and the band parameter :math:`\delta > 1`: + + - if :math:`v > s\,\delta\, d` the primal is making less progress, so the primal step + is boosted and the dual shrunk, + :math:`\tau \leftarrow \tau/(1-\alpha)`, :math:`\sigma_i \leftarrow (1-\alpha)\sigma_i`; + - if :math:`v < s\, d/\delta` the opposite tilt is applied, + :math:`\tau \leftarrow (1-\alpha)\tau`, :math:`\sigma_i \leftarrow \sigma_i/(1-\alpha)`; + - otherwise the step sizes are left unchanged. + + Every dual step size :math:`\sigma_i` is rescaled by the same factor, and each time a + tilt is applied the adaptation strength decays, :math:`\alpha \leftarrow \eta\,\alpha`, + so the changes become progressively smaller (guaranteeing the summable-deviation + condition of :cite:`chambolle2023stochastic`). The products :math:`\tau\sigma_i` are + preserved. + + Parameters + ---------- + initial_step_size : list of length two, optional, default=[None, None] + Initial primal and dual step sizes ``[tau, sigma]``. ``tau`` is a positive scalar + and ``sigma`` is either a positive scalar (broadcast to every operator) or a list + of one positive number per operator. Any entry left as ``None`` is filled in from + the operator norms and sampling probabilities using the standard SPDHG relations + (as in :class:`SPDHGConstantStepSize`). + initial_alpha : positive :obj:`float`, optional, default=0.95 + Initial value of the adaptation strength :math:`\alpha \in (0,1)`. Each rebalance + multiplies or divides :math:`\tau` by :math:`1-\alpha`, so :math:`\alpha` close to + one gives a *large* rescaling (:math:`\alpha = 0.95` changes :math:`\tau` by a + factor of 20) and :math:`\alpha` close to zero a gentle one. + eta : positive :obj:`float`, optional, default=0.95 + The decay factor :math:`\eta \in (0,1)` applied to :math:`\alpha` each time the + step sizes are rebalanced. Since :math:`\alpha` decays only on a rebalance, this + bounds the total drift of :math:`\tau` through + :math:`\sum_k \alpha_k = \alpha_0/(1-\eta)`; values very close to one let + :math:`\alpha` persist for hundreds of rebalances and the step sizes then swing + over several orders of magnitude. + delta : positive :obj:`float`, greater than one, optional, default=1.5 + The band parameter :math:`\delta` setting how far apart the primal and dual residuals are allowed to drift before rebalancing. + s : positive :obj:`float`, optional, default=norm of the operator + The balancing scale :math:`s` used to compare the primal and dual residuals. Defaults to the operator norm :math:`\|A\|`, as recommended in :cite:`chambolle2023stochastic`. + auto_stop : :obj:`bool`, optional, default=True + If True, the adaptive updates stop and the extra storage is released once the step sizes have been unchanged for ``auto_stop_patience`` consecutive iterations. + auto_stop_patience : :obj:`int`, optional, default=``10 * n_operators`` + Number of consecutive iterations with no change to the step sizes after which the adaptive updates are stopped (only used when ``auto_stop=True``). If left as ``None`` it is set to ten times the number of operators, so that the patience scales with the number of iterations needed to sample every subset. + alpha_tolerance : positive :obj:`float`, optional, default=1e-3 + Threshold on the adaptation strength :math:`\alpha` below which the adaptive updates are stopped (only used when ``auto_stop=True``). Once :math:`\alpha` has decayed this far a rebalance changes the step sizes by a relative amount of roughly :math:`\alpha`, so the rule keeps paying for the extra operator applications and storage while making no practical difference. Set to ``0`` to disable this criterion and rely on ``auto_stop_patience`` alone. + + Notes + ----- + The dual step size ``sigma`` is a list with one entry per operator, while ``tau`` is a + scalar; :meth:`get_initial_step_size` and :meth:`get_step_size` return the tuple + ``(tau, sigma)``. Computing the residuals costs one extra forward and one extra adjoint + application per iteration (the overhead discussed in :cite:`chambolle2023stochastic`). + + See Also + -------- + SPDHGAdaptiveStepSizeAngle : The companion angle-alignment rule (rule (b)) from the same paper. + PDHGAdaptiveStepSize2013 : The deterministic PDHG rule this generalises. + + Reference + --------- + Chambolle, A., Delplancke, C., Ehrhardt, M.J., Schönlieb, C.-B. and Tang, J., 2023. Stochastic Primal-Dual Hybrid Gradient Algorithm with Adaptive Step-Sizes. arXiv preprint arXiv:2301.02511. :cite:`chambolle2023stochastic` + """ + + def __init__(self, initial_step_size=[None, None], initial_alpha=0.95, eta=0.95, + delta=1.5, s=None, auto_stop=True, auto_stop_patience=None, + alpha_tolerance=1e-3, global_sigma=True): + '''Initialises the step size rule''' + super().__init__(initial_step_size=initial_step_size, initial_alpha=initial_alpha, + eta=eta, auto_stop=auto_stop, auto_stop_patience=auto_stop_patience, + alpha_tolerance=alpha_tolerance) + self.delta = delta + self.s = s + self.global_sigma = global_sigma + + def _rebalance(self, algorithm, i, p_i, tau, sigma): + if self.s is None: + self.s = algorithm.operator.norm() # default balancing scale, ||A|| + + + v = self.v_tmp.abs().sum() + algorithm.operator[i].direct(self.x_prev, out=self.forward_tmp[i]) + self.y_prev[i].sapyb(1.0 / (p_i * sigma[i]), self.forward_tmp[i], -1.0/p_i, out=self.y_prev[i]) + d = self.y_prev[i].abs().sum() + + log.debug('SPDHG adaptive balancing: v = {}, d = {}, s = {}'.format(v, d, self.s)) + if v > self.s * self.delta * d: + factor = 1 - self.alpha + tau/= factor + if self.global_sigma: + sigma = [s * factor for s in sigma] + else: + sigma[i] *= factor + self.alpha *= self.eta + self.count = 0 + log.debug('SPDHG adaptive balancing: v > s*delta*d, increasing tau and decreasing sigma, tau = {}, sigma = {}'.format(tau, sigma)) + return tau, sigma, True + elif v < self.s * d / self.delta: + factor = 1 - self.alpha + tau*= factor + if self.global_sigma: + sigma = [s / factor for s in sigma] + else: + sigma[i] /= factor + self.alpha *= self.eta + self.count = 0 + log.debug('SPDHG adaptive balancing: v < s*d/delta, decreasing tau and increasing sigma, tau = {}, sigma = {}'.format(tau, sigma)) + return tau, sigma, True + + log.debug('SPDHG adaptive balancing: v and d within band, no change to step sizes.') + return tau, sigma, False + + +class SPDHGAdaptiveStepSizeAngle(_SPDHGAdaptiveStepSizeBase): + r"""Adaptively adjusts the SPDHG primal and dual step sizes using the alignment between successive primal directions, following rule (b) of :cite:`chambolle2023stochastic`. + + This is an SPDHG-compatible step-size rule (A-SPDHG rule (b), Algorithm 3.2 of + :cite:`chambolle2023stochastic`). It is the stochastic extension of the angle-based + adaptive PDHG scheme of Yokota and Hontani, and an alternative to the residual-balancing + :class:`SPDHGAdaptiveStepSizeBalancing`. + + At the end of each iteration, with sampled subset :math:`i`, sampling probability + :math:`p_i`, primal increment :math:`\Delta x = x^{k} - x^{k+1}` and dual increment + :math:`\Delta y_i = y^{k}_i - y^{k+1}_i`, the primal-residual direction is formed + (equation (3.8) of :cite:`chambolle2023stochastic`), + + .. math:: + + q = \frac{\Delta x}{\tau} - \frac{1}{p_i} A_i^{*}\Delta y_i , + + and the cosine of the angle between :math:`\Delta x` and :math:`q` is measured + (equation (3.9)), + + .. math:: + + w = \frac{\langle \Delta x, q\rangle}{\|\Delta x\|_2\,\|q\|_2} . + + The step sizes are then adjusted against the threshold :math:`c` (close to one): + + - if :math:`w \ge c` the directions are well aligned, so the primal step is increased, + :math:`\tau \leftarrow (1+\alpha)\tau`, :math:`\sigma_i \leftarrow \sigma_i/(1+\alpha)`; + - if :math:`w < 0` the directions oppose, so the primal step is decreased, + :math:`\tau \leftarrow \tau/(1+\alpha)`, :math:`\sigma_i \leftarrow (1+\alpha)\sigma_i`; + - otherwise the step sizes are left unchanged. + + Every dual step size :math:`\sigma_i` is rescaled by the same factor, and each time the + step sizes change the adaptation strength decays, :math:`\alpha \leftarrow \eta\,\alpha`. + The products :math:`\tau\sigma_i` are preserved. + + Parameters + ---------- + initial_step_size : list of length two, optional, default=[None, None] + Initial primal and dual step sizes ``[tau, sigma]``. ``tau`` is a positive scalar + and ``sigma`` is either a positive scalar (broadcast to every operator) or a list + of one positive number per operator. Any entry left as ``None`` is filled in from + the operator norms and sampling probabilities using the standard SPDHG relations. + initial_alpha : positive :obj:`float`, optional, default=1.0 + Initial value of the adaptation strength :math:`\alpha`. Each change multiplies or divides :math:`\tau` by :math:`1+\alpha`, so :math:`\alpha_0 = 1` (the value used in the paper) rescales :math:`\tau` by a factor of two. + eta : positive :obj:`float`, optional, default=0.995 + The decay factor :math:`\eta \in (0,1)` applied to :math:`\alpha` each time the step sizes change. Since :math:`\alpha` decays only when the step sizes change, this bounds the total drift of :math:`\tau` through :math:`\sum_k \alpha_k = \alpha_0/(1-\eta)`. + c : :obj:`float`, optional, default=0.999 + The cosine threshold :math:`c` above which the primal step size is increased. Should be close to one in high-dimensional problems, as recommended in :cite:`chambolle2023stochastic`. + auto_stop : :obj:`bool`, optional, default=True + If True, the adaptive updates stop and the extra storage is released once the step sizes have been unchanged for ``auto_stop_patience`` consecutive iterations. + auto_stop_patience : :obj:`int`, optional, default=``10 * n_operators`` + Number of consecutive iterations with no change to the step sizes after which the adaptive updates are stopped (only used when ``auto_stop=True``). If left as ``None`` it is set to ten times the number of operators, so that the patience scales with the number of iterations needed to sample every subset. + alpha_tolerance : positive :obj:`float`, optional, default=1e-3 + Threshold on the adaptation strength :math:`\alpha` below which the adaptive updates are stopped (only used when ``auto_stop=True``). Once :math:`\alpha` has decayed this far a rebalance changes the step sizes by a relative amount of roughly :math:`\alpha`, so the rule keeps paying for the extra operator applications and storage while making no practical difference. Set to ``0`` to disable this criterion and rely on ``auto_stop_patience`` alone. + + Notes + ----- + The dual step size ``sigma`` is a list with one entry per operator, while ``tau`` is a + scalar; :meth:`get_initial_step_size` and :meth:`get_step_size` return the tuple + ``(tau, sigma)``. Computing the direction costs one extra adjoint application per + iteration. + + See Also + -------- + SPDHGAdaptiveStepSizeBalancing : The companion residual-balancing rule (rule (a)) from the same paper. + + Reference + --------- + Chambolle, A., Delplancke, C., Ehrhardt, M.J., Schönlieb, C.-B. and Tang, J., 2023. Stochastic Primal-Dual Hybrid Gradient Algorithm with Adaptive Step-Sizes. arXiv preprint arXiv:2301.02511. :cite:`chambolle2023stochastic` + """ + + def __init__(self, initial_step_size=[None, None], initial_alpha=1.0, eta=0.995, + c=0.999, auto_stop=True, auto_stop_patience=None, alpha_tolerance=1e-3): + '''Initialises the step size rule''' + super().__init__(initial_step_size=initial_step_size, initial_alpha=initial_alpha, + eta=eta, auto_stop=auto_stop, auto_stop_patience=auto_stop_patience, + alpha_tolerance=alpha_tolerance) + self.c = c + + def _rebalance(self, algorithm, i, p_i, tau, sigma): + # cosine of the angle between Delta x (held in x_prev) and the primal + # residual direction q (held in v_tmp), both set by the base class. + w = self.x_prev.dot(self.v_tmp) / \ + (self.x_prev.norm() * self.v_tmp.norm()) + + log.debug('SPDHG adaptive angle: cos(angle) w = {}, index = {}'.format(w, i)) + if w >= self.c: + factor = 1 + self.alpha + tau_new = tau * factor + sigma_new = [si / factor for si in sigma] + self.alpha *= self.eta + self.count = 0 + log.debug('SPDHG adaptive angle: increasing step sizes, tau -> {}, sigma -> {}'.format(tau_new, sigma_new)) + return tau_new, sigma_new, True + elif w < 0: + factor = 1 + self.alpha + tau_new = tau / factor + sigma_new = [si * factor for si in sigma] + self.alpha *= self.eta + self.count = 0 + log.debug('SPDHG adaptive angle: decreasing step sizes, tau -> {}, sigma -> {}'.format(tau_new, sigma_new)) + return tau_new, sigma_new, True + return tau, sigma, False \ No newline at end of file diff --git a/Wrappers/Python/cil/optimisation/utilities/__init__.py b/Wrappers/Python/cil/optimisation/utilities/__init__.py index daa32667a0..5246ec354e 100644 --- a/Wrappers/Python/cil/optimisation/utilities/__init__.py +++ b/Wrappers/Python/cil/optimisation/utilities/__init__.py @@ -19,5 +19,5 @@ from .sampler import Sampler from .sampler import SamplerRandom -from .StepSizeMethods import ConstantStepSize, ArmijoStepSizeRule, StepSizeRule, BarzilaiBorweinStepSizeRule +from .StepSizeMethods import ConstantStepSize, ArmijoStepSizeRule, StepSizeRule, BarzilaiBorweinStepSizeRule, PDHGAdaptiveStepSize2013, PDHGAdaptiveStepSize2015, PDHGStronglyConvexUpdate, PDHGConstantStepSize, PDHGBayesOptimisationStepSize, SPDHGStepSizesFromRatio, SPDHGConstantStepSize, SPDHGBayesOptimisationStepSize, SPDHGAdaptiveStepSizeBalancing, SPDHGAdaptiveStepSizeAngle from .preconditioner import Preconditioner, AdaptiveSensitivity, Sensitivity diff --git a/Wrappers/Python/test/test_algorithm_convergence.py b/Wrappers/Python/test/test_algorithm_convergence.py index a5e5aaa5a8..43ff05adf7 100644 --- a/Wrappers/Python/test/test_algorithm_convergence.py +++ b/Wrappers/Python/test/test_algorithm_convergence.py @@ -1,8 +1,8 @@ from cil.optimisation.algorithms import SPDHG, PDHG, LSQR, FISTA, APGD, GD, PD3O -from cil.optimisation.functions import L2NormSquared, IndicatorBox, BlockFunction, ZeroFunction, KullbackLeibler, OperatorCompositionFunction, LeastSquares, TotalVariation, MixedL21Norm +from cil.optimisation.functions import L2NormSquared, IndicatorBox, BlockFunction, ZeroFunction, KullbackLeibler, OperatorCompositionFunction, LeastSquares, TotalVariation, MixedL21Norm, L1Norm from cil.optimisation.operators import BlockOperator, IdentityOperator, MatrixOperator, GradientOperator -from cil.optimisation.utilities import Sampler, BarzilaiBorweinStepSizeRule, ArmijoStepSizeRule +from cil.optimisation.utilities import Sampler, BarzilaiBorweinStepSizeRule, ArmijoStepSizeRule, PDHGAdaptiveStepSize2013, PDHGAdaptiveStepSize2015, PDHGBayesOptimisationStepSize, SPDHGBayesOptimisationStepSize, SPDHGStepSizesFromRatio, SPDHGAdaptiveStepSizeBalancing, SPDHGAdaptiveStepSizeAngle from cil.framework import AcquisitionGeometry, BlockDataContainer, BlockGeometry, VectorData, ImageGeometry from cil.utilities import dataexample from cil.utilities import noise as applynoise @@ -10,6 +10,7 @@ import numpy as np import unittest from testclass import CCPiTestClass +from utils import has_skopt from scipy.optimize import minimize, rosen from cil.optimisation.functions import Rosenbrock @@ -27,6 +28,8 @@ except ImportError: has_astra = False +import logging +log = logging.getLogger(__name__) class TestSPDHG(CCPiTestClass): @@ -319,3 +322,279 @@ def test_bb_step_size_gd_converge(self): alg = GD(initial=initial, f=f, step_size=ss_rule) alg.run(300, verbose=0) self.assertNumpyArrayAlmostEqual(alg.x.as_array(), x, decimal=4) + +class TestPDHGConvergence(CCPiTestClass): + def test_PDHG_Denoising(self): + # adapted from demo PDHG_TV_Color_Denoising.py in CIL-Demos repository + data = dataexample.PEPPERS.get(size=(256, 256)) + ig = data.geometry + ag = ig + + which_noise = 0 + # Create noisy data. + noises = ['gaussian', 'poisson', 's&p'] + dnoise = noises[which_noise] + + def setup(data, dnoise): + if dnoise == 's&p': + n1 = applynoise.saltnpepper( + data, salt_vs_pepper=0.9, amount=0.2, seed=10) + elif dnoise == 'poisson': + scale = 5 + n1 = applynoise.poisson(data.as_array()/scale, seed=10)*scale + elif dnoise == 'gaussian': + n1 = applynoise.gaussian(data.as_array(), seed=10) + else: + raise ValueError('Unsupported Noise ', noise) + noisy_data = ig.allocate() + noisy_data.fill(n1) + + # Regularisation Parameter depending on the noise distribution + if dnoise == 's&p': + alpha = 0.8 + elif dnoise == 'poisson': + alpha = 1 + elif dnoise == 'gaussian': + alpha = .3 + # fidelity + if dnoise == 's&p': + g = L1Norm(b=noisy_data) + elif dnoise == 'poisson': + g = KullbackLeibler(b=noisy_data) + elif dnoise == 'gaussian': + g = 0.5 * L2NormSquared(b=noisy_data) + return noisy_data, alpha, g + + noisy_data, alpha, g = setup(data, dnoise) + operator = GradientOperator( + ig, correlation=GradientOperator.CORRELATION_SPACE, backend='numpy') + + f1 = alpha * MixedL21Norm() + + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + + # Setup and run the PDHG algorithm + pdhg1 = PDHG(f=f1, g=g, operator=operator, step_size=[tau,sigma]) + self.assertEqual( pdhg1.tau, tau) + self.assertEqual(pdhg1.sigma, sigma) + + pdhg1.update_objective_interval = 200 + pdhg1.run(1000, verbose=0) + + rmse = (pdhg1.get_output() - data).norm() / data.as_array().size + log.info("RMSE %F", rmse) + self.assertLess(rmse, 2e-4) + + which_noise = 1 + noise = noises[which_noise] + noisy_data, alpha, g = setup(data, noise) + operator = GradientOperator( + ig, correlation=GradientOperator.CORRELATION_SPACE, backend='numpy') + + f1 = alpha * MixedL21Norm() + + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + + # Setup and run the PDHG algorithm + pdhg1 = PDHG(f=f1, g=g, operator=operator, step_size=(tau,sigma), + update_objective_interval=200) + pdhg1.run(1000, verbose=0) + + rmse = (pdhg1.get_output() - data).norm() / data.as_array().size + log.info("RMSE %f", rmse) + self.assertLess(rmse, 2e-4) + + which_noise = 2 + noise = noises[which_noise] + noisy_data, alpha, g = setup(data, noise) + operator = GradientOperator( + ig, correlation=GradientOperator.CORRELATION_SPACE, backend='numpy') + + f1 = alpha * MixedL21Norm() + + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) + + # Setup and run the PDHG algorithm + pdhg1 = PDHG(f=f1, g=g, operator=operator, step_size=(tau,sigma) ) + pdhg1.update_objective_interval = 200 + pdhg1.run(1000, verbose=0) + + rmse = (pdhg1.get_output() - data).norm() / data.as_array().size + log.info("RMSE %f", rmse) + self.assertLess(rmse, 2e-4) + + @unittest.skipUnless(has_skopt, "scikit-optimize (skopt) not installed") + def test_PDHG_adaptive_bayes(self): + ig = ImageGeometry(3, 3) + data = ig.allocate(0) + data.fill(np.diag([1, 2, 3])) + ideal = ig.allocate(0) + ideal.fill(np.diag([0.5, 1, 1.5])) + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGBayesOptimisationStepSize( + gamma_bounds=None, n_initial_points=5, n_calls=10, n_iterations=10, seed = 42) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + gamma = np.sqrt(pdhg.sigma / pdhg.tau) + pdhg.run(20, verbose=1) + self.assertAlmostEqual((pdhg.x-ideal).norm(), 0, places=4) + + + + def test_PDHG_adaptive_step_size_2013(self): + ig = ImageGeometry(3, 3) + data = ig.allocate(0) + data.fill(np.diag([1, 2, 3])) + ideal = ig.allocate(0) + ideal.fill(np.diag([0.5, 1, 1.5])) + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013() + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + gamma = np.sqrt(pdhg.sigma / pdhg.tau) + pdhg.run(20, verbose=1) + self.assertAlmostEqual((pdhg.x-ideal).norm(), 0, places=4) + + def test_PDHG_adaptive_step_size_2015(self): + ig = ImageGeometry(3, 3) + data = ig.allocate(0) + data.fill(np.diag([1, 2, 3])) + ideal = ig.allocate(0) + ideal.fill(np.diag([0.5, 1, 1.5])) + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2015() + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + gamma = np.sqrt(pdhg.sigma / pdhg.tau) + pdhg.run(30, verbose=1) + self.assertAlmostEqual((pdhg.x-ideal).norm(), 0, places=4) + + +class TestSPDHGConvergence(CCPiTestClass): + @unittest.skipUnless(has_skopt, "scikit-optimize (skopt) not installed") + def test_SPDHG_adaptive_bayes(self): + subsets = 2 + ig = ImageGeometry(3, 3) + data_ind=ig.allocate(0) + data_ind.fill(np.diag([1, 2, 3])) + ideal_ind = ig.allocate(0) + ideal_ind.fill(np.diag([2/3, 4/3, 2])) + + + self.A = BlockOperator( + *[IdentityOperator(ig) for i in range(subsets)]) + + + # block function + self.F = BlockFunction(*[L2NormSquared(b=data_ind) + for i in range(subsets)]) + self.G = L2NormSquared() + + rule = SPDHGBayesOptimisationStepSize( + gamma_bounds=None, n_initial_points=5, n_calls=10, n_iterations=10, seed = 42) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + spdhg.run(50, verbose=1) + self.assertNumpyArrayAlmostEqual(spdhg.x.as_array() , ideal_ind.as_array(), decimal=4) + + def test_convergence_balancing(self): + # From a poorly balanced start (tau far too large, sigma far too small) the + # residual-balancing rule (rule (a) of Chambolle et al. 2023) should tilt the + # step sizes back, reduce the objective and reach a value comparable to a + # well-tuned constant rule, converging to the same minimiser. + # SPDHG samples subsets at random, so both runs are seeded for reproducibility. + subsets = 2 + ig = ImageGeometry(3, 3) + data_ind=ig.allocate(0) + data_ind.fill(np.diag([1, 2, 3])) + ideal_ind = ig.allocate(0) + ideal_ind.fill(np.diag([2/3, 4/3, 2])) + + self.A = BlockOperator( + *[IdentityOperator(ig) for i in range(subsets)]) + + + # block function + self.F = BlockFunction(*[L2NormSquared(b=data_ind) + for i in range(subsets)]) + self.G = L2NormSquared() + + seed = 7 + np.random.seed(seed) + constant = SPDHG(f=self.F, g=self.G, operator=self.A, + sampler=Sampler.random_with_replacement(subsets, seed=seed)) + constant.run(50) + + rule = SPDHGAdaptiveStepSizeBalancing(initial_step_size=[10.0, 0.001]) + adaptive = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule, + sampler=Sampler.random_with_replacement(subsets, seed=seed)) + adaptive.run(50) + # objective is recorded during the run; the first entry is the initial value + self.assertLess(adaptive.objective[-1], adaptive.objective[0]) + self.assertLess(adaptive.objective[-1], 1.1 * constant.objective[-1] + 1e-3) + self.assertNumpyArrayAlmostEqual(constant.x.as_array() , ideal_ind.as_array(), decimal=4) + # the adaptive run recovers from the bad start but, having spent early iterations + # rebalancing, is a little behind the well-tuned constant rule at 50 iterations. + self.assertNumpyArrayAlmostEqual(adaptive.x.as_array() , ideal_ind.as_array(), decimal=2) + + def test_convergence_angles(self): + # Companion to test_convergence_balancing for the angle/alignment rule (rule (b) + # of Chambolle et al. 2023). This rule *increases* the primal step when successive + # primal directions stay aligned, so the imbalanced start is chosen the other way + # round: tau far too small (sigma filled in from the operator norms). From there it + # should grow tau, reduce the objective and reach a value comparable to a well-tuned + # constant rule, converging to the same minimiser. + # SPDHG samples subsets at random, so both runs are seeded for reproducibility. + subsets = 2 + ig = ImageGeometry(3, 3) + data_ind=ig.allocate(0) + data_ind.fill(np.diag([1, 2, 3])) + ideal_ind = ig.allocate(0) + ideal_ind.fill(np.diag([2/3, 4/3, 2])) + + self.A = BlockOperator( + *[IdentityOperator(ig) for i in range(subsets)]) + + + # block function + self.F = BlockFunction(*[L2NormSquared(b=data_ind) + for i in range(subsets)]) + self.G = L2NormSquared() + + seed = 7 + np.random.seed(seed) + constant = SPDHG(f=self.F, g=self.G, operator=self.A, + sampler=Sampler.random_with_replacement(subsets, seed=seed)) + constant.run(50) + + rule = SPDHGAdaptiveStepSizeAngle(initial_step_size=[0.001, None]) + adaptive = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule, + sampler=Sampler.random_with_replacement(subsets, seed=seed)) + adaptive.run(50) + # objective is recorded during the run; the first entry is the initial value + self.assertLess(adaptive.objective[-1], adaptive.objective[0]) + self.assertLess(adaptive.objective[-1], 1.1 * constant.objective[-1] + 1e-3) + self.assertNumpyArrayAlmostEqual(constant.x.as_array() , ideal_ind.as_array(), decimal=4) + # the adaptive run recovers from the bad start but, having spent early iterations + # rebalancing, is a little behind the well-tuned constant rule at 50 iterations. + self.assertNumpyArrayAlmostEqual(adaptive.x.as_array() , ideal_ind.as_array(), decimal=2) \ No newline at end of file diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py index 1da91a6e3a..e9447ca691 100644 --- a/Wrappers/Python/test/test_algorithms.py +++ b/Wrappers/Python/test/test_algorithms.py @@ -30,7 +30,7 @@ from cil.framework.labels import FillType -from cil.optimisation.utilities import ArmijoStepSizeRule, ConstantStepSize, Sampler, callbacks, Sensitivity, StepSizeRule +from cil.optimisation.utilities import ArmijoStepSizeRule, ConstantStepSize, Sampler, callbacks, Sensitivity, StepSizeRule, SPDHGStepSizesFromRatio, SPDHGConstantStepSize from cil.optimisation.algorithms.APGD import NesterovMomentum, ScalarMomentumCoefficient, ConstantMomentum from cil.optimisation.operators import IdentityOperator, AdjointOperator from cil.optimisation.operators import GradientOperator, BlockOperator, MatrixOperator @@ -840,291 +840,9 @@ def test_init_primal_dual(self): - def test_PDHG_Denoising(self): - # adapted from demo PDHG_TV_Color_Denoising.py in CIL-Demos repository - data = dataexample.PEPPERS.get(size=(256, 256)) - ig = data.geometry - ag = ig - - which_noise = 0 - # Create noisy data. - noises = ['gaussian', 'poisson', 's&p'] - dnoise = noises[which_noise] - - def setup(data, dnoise): - if dnoise == 's&p': - n1 = applynoise.saltnpepper( - data, salt_vs_pepper=0.9, amount=0.2, seed=10) - elif dnoise == 'poisson': - scale = 5 - n1 = applynoise.poisson(data.as_array()/scale, seed=10)*scale - elif dnoise == 'gaussian': - n1 = applynoise.gaussian(data.as_array(), seed=10) - else: - raise ValueError('Unsupported Noise ', noise) - noisy_data = ig.allocate() - noisy_data.fill(n1) - - # Regularisation Parameter depending on the noise distribution - if dnoise == 's&p': - alpha = 0.8 - elif dnoise == 'poisson': - alpha = 1 - elif dnoise == 'gaussian': - alpha = .3 - # fidelity - if dnoise == 's&p': - g = L1Norm(b=noisy_data) - elif dnoise == 'poisson': - g = KullbackLeibler(b=noisy_data) - elif dnoise == 'gaussian': - g = 0.5 * L2NormSquared(b=noisy_data) - return noisy_data, alpha, g - - noisy_data, alpha, g = setup(data, dnoise) - operator = GradientOperator( - ig, correlation=GradientOperator.CORRELATION_SPACE, backend='numpy') - - f1 = alpha * MixedL21Norm() - - # Compute operator Norm - normK = operator.norm() - - # Primal & dual stepsizes - sigma = 1 - tau = 1/(sigma*normK**2) - - # Setup and run the PDHG algorithm - pdhg1 = PDHG(f=f1, g=g, operator=operator, tau=tau, sigma=sigma) - pdhg1.update_objective_interval = 200 - pdhg1.run(1000, verbose=0) - - rmse = (pdhg1.get_output() - data).norm() / data.as_array().size - log.info("RMSE %F", rmse) - self.assertLess(rmse, 2e-4) - - which_noise = 1 - noise = noises[which_noise] - noisy_data, alpha, g = setup(data, noise) - operator = GradientOperator( - ig, correlation=GradientOperator.CORRELATION_SPACE, backend='numpy') - - f1 = alpha * MixedL21Norm() - - # Compute operator Norm - normK = operator.norm() - - # Primal & dual stepsizes - sigma = 1 - tau = 1/(sigma*normK**2) - - # Setup and run the PDHG algorithm - pdhg1 = PDHG(f=f1, g=g, operator=operator, tau=tau, sigma=sigma, - update_objective_interval=200) - - pdhg1.run(1000, verbose=0) - - rmse = (pdhg1.get_output() - data).norm() / data.as_array().size - log.info("RMSE %f", rmse) - self.assertLess(rmse, 2e-4) - - which_noise = 2 - noise = noises[which_noise] - noisy_data, alpha, g = setup(data, noise) - operator = GradientOperator( - ig, correlation=GradientOperator.CORRELATION_SPACE, backend='numpy') - - f1 = alpha * MixedL21Norm() - - # Compute operator Norm - normK = operator.norm() - - # Primal & dual stepsizes - sigma = 1 - tau = 1/(sigma*normK**2) - - # Setup and run the PDHG algorithm - pdhg1 = PDHG(f=f1, g=g, operator=operator, tau=tau, sigma=sigma) - pdhg1.update_objective_interval = 200 - pdhg1.run(1000, verbose=0) - - rmse = (pdhg1.get_output() - data).norm() / data.as_array().size - log.info("RMSE %f", rmse) - self.assertLess(rmse, 2e-4) - - def test_PDHG_step_sizes(self): - ig = ImageGeometry(3, 3) - data = ig.allocate('random', seed=3) - - f = L2NormSquared(b=data) - g = L2NormSquared() - operator = 3*IdentityOperator(ig) - - # check if sigma, tau are None - pdhg = PDHG(f=f, g=g, operator=operator) - self.assertAlmostEqual(pdhg.sigma, 1./operator.norm()) - self.assertAlmostEqual(pdhg.tau, 1./operator.norm()) - - # check if sigma is negative - with self.assertRaises(ValueError): - pdhg = PDHG(f=f, g=g, operator=operator, - sigma=-1) - - # check if tau is negative - with self.assertRaises(ValueError): - pdhg = PDHG(f=f, g=g, operator=operator, tau=-1) - - # check if tau is None - sigma = 3.0 - pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma) - self.assertAlmostEqual(pdhg.sigma, sigma) - self.assertAlmostEqual(pdhg.tau, 1./(sigma * operator.norm()**2)) - - # check if sigma is None - tau = 3.0 - pdhg = PDHG(f=f, g=g, operator=operator, tau=tau) - self.assertAlmostEqual(pdhg.tau, tau) - self.assertAlmostEqual(pdhg.sigma, 1./(tau * operator.norm()**2)) - - # check if sigma/tau are not None - tau = 1.0 - sigma = 1.0 - pdhg = PDHG(f=f, g=g, operator=operator, tau=tau, - sigma=sigma) - self.assertAlmostEqual(pdhg.tau, tau) - self.assertAlmostEqual(pdhg.sigma, sigma) - - # check sigma/tau as arrays, sigma wrong shape - ig1 = ImageGeometry(2, 2) - sigma = ig1.allocate() - with self.assertRaises(ValueError): - pdhg = PDHG(f=f, g=g, operator=operator, - sigma=sigma) - - # check sigma/tau as arrays, tau wrong shape - tau = ig1.allocate() - with self.assertRaises(ValueError): - pdhg = PDHG(f=f, g=g, operator=operator, tau=tau) - - # check sigma not Number or object with correct shape - with self.assertRaises(AttributeError): - pdhg = PDHG(f=f, g=g, operator=operator, - sigma="sigma") - - # check tau not Number or object with correct shape - with self.assertRaises(AttributeError): - pdhg = PDHG(f=f, g=g, operator=operator, - tau="tau") - - # check warning message if condition is not satisfied - sigma = 4/operator.norm() - tau = 1/3 - with self.assertWarnsRegex(UserWarning, "Convergence criterion"): - pdhg = PDHG(f=f, g=g, operator=operator, tau=tau, - sigma=sigma) - - # check no warning message if check convergence is false - sigma = 4/operator.norm() - tau = 1/3 - with warnings.catch_warnings(record=True) as warnings_log: - pdhg = PDHG(f=f, g=g, operator=operator, tau=tau, - sigma=sigma, check_convergence=False) - self.assertEqual(warnings_log, []) - - # check no warning message if condition is satisfied - sigma = 1/operator.norm() - tau = 1/3 - with warnings.catch_warnings(record=True) as warnings_log: - pdhg = PDHG(f=f, g=g, operator=operator, tau=tau, - sigma=sigma) - self.assertEqual(warnings_log, []) - - def test_PDHG_strongly_convex_gamma_g(self): - ig = ImageGeometry(3, 3) - data = ig.allocate('random', seed=3) - - f = L2NormSquared(b=data) - g = L2NormSquared() - operator = IdentityOperator(ig) - - # sigma, tau - sigma = 1.0 - tau = 1.0 - - pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, - gamma_g=0.5) - pdhg.run(1, verbose=0) - self.assertAlmostEqual( - pdhg.theta, 1.0 / np.sqrt(1 + 2 * pdhg.gamma_g * tau)) - self.assertAlmostEqual(pdhg.tau, tau * pdhg.theta) - self.assertAlmostEqual(pdhg.sigma, sigma / pdhg.theta) - pdhg.run(4, verbose=0) - self.assertNotEqual(pdhg.sigma, sigma) - self.assertNotEqual(pdhg.tau, tau) - - # check negative strongly convex constant - with self.assertRaises(ValueError): - pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, - gamma_g=-0.5) - - # check strongly convex constant not a number - with self.assertRaises(ValueError): - pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, - gamma_g="-0.5") - - def test_PDHG_strongly_convex_gamma_fcong(self): - ig = ImageGeometry(3, 3) - data = ig.allocate('random', seed=3) - - f = L2NormSquared(b=data) - g = L2NormSquared() - operator = IdentityOperator(ig) - - # sigma, tau - sigma = 1.0 - tau = 1.0 - - pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, - gamma_fconj=0.5) - pdhg.run(1, verbose=0) - self.assertEqual(pdhg.theta, 1.0 / np.sqrt(1 + - 2 * pdhg.gamma_fconj * sigma)) - self.assertEqual(pdhg.tau, tau / pdhg.theta) - self.assertEqual(pdhg.sigma, sigma * pdhg.theta) - pdhg.run(4, verbose=0) - self.assertNotEqual(pdhg.sigma, sigma) - self.assertNotEqual(pdhg.tau, tau) - - # check negative strongly convex constant - try: - pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, - gamma_fconj=-0.5) - except ValueError as ve: - log.info(str(ve)) - - # check strongly convex constant not a number - try: - pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, - gamma_fconj="-0.5") - except ValueError as ve: - log.info(str(ve)) - - def test_PDHG_strongly_convex_both_fconj_and_g(self): - - ig = ImageGeometry(3, 3) - data = ig.allocate('random', seed=3) - - f = L2NormSquared(b=data) - g = L2NormSquared() - operator = IdentityOperator(ig) - - try: - pdhg = PDHG(f=f, g=g, operator=operator, - gamma_g=0.5, gamma_fconj=0.5) - pdhg.run(verbose=0) - except ValueError as err: - log.info(str(err)) + + def test_pdhg_theta(self): ig = ImageGeometry(3, 3) data = ig.allocate('random', seed=3) @@ -1373,48 +1091,7 @@ def test_SPDHG_defaults_and_setters(self): self.assertListEqual(spdhg._prob_weights, [ 1/self.subsets] * self.subsets) self.assertEqual(spdhg._sampler._type, 'random_with_replacement') - self.assertListEqual( - spdhg.sigma, [rho / ni for ni in spdhg._norms]) - self.assertEqual(spdhg.tau, min([rho*pi / (si * ni**2) for pi, ni, - si in zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) - self.assertNumpyArrayEqual( - spdhg.x.as_array(), self.A.domain_geometry().allocate(0).as_array()) - self.assertEqual(spdhg.update_objective_interval, 1) - - # Test SPDHG setters - "from ratio" - gamma = 3.7 - rho = 5.6 - spdhg.set_step_sizes_from_ratio(gamma, rho) - self.assertListEqual( - spdhg.sigma, [gamma * rho / ni for ni in spdhg._norms]) - self.assertEqual(spdhg.tau, min([pi*rho / (si * ni**2) for pi, ni, - si in zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) - - # Test SPDHG setters - set_step_sizes default values for sigma and tau - gamma = 1. - rho = .99 - spdhg.set_step_sizes() - self.assertListEqual( - spdhg.sigma, [rho / ni for ni in spdhg._norms]) - self.assertEqual(spdhg.tau, min([rho*pi / (si * ni**2) for pi, ni, - si in zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) - - # Test SPDHG setters - set_step_sizes with sigma and tau - spdhg.set_step_sizes(sigma=[1]*self.subsets, tau=100) - self.assertListEqual(spdhg.sigma, [1]*self.subsets) - self.assertEqual(spdhg.tau, 100) - - # Test SPDHG setters - set_step_sizes with sigma - spdhg.set_step_sizes(sigma=[1]*self.subsets, tau=None) - self.assertListEqual(spdhg.sigma, [1]*self.subsets) - self.assertEqual(spdhg.tau, min([(rho*pi / (si * ni**2)) for pi, ni, - si in zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) - - # Test SPDHG setters - set_step_sizes with tau - spdhg.set_step_sizes(sigma=None, tau=100) - self.assertListEqual(spdhg.sigma, [ - gamma * rho*pi / (spdhg.tau*ni**2) for ni, pi in zip(spdhg._norms, spdhg._prob_weights)]) - self.assertEqual(spdhg.tau, 100) + def test_spdhg_non_default_init(self): # Test SPDHG init with non-default values @@ -1452,27 +1129,33 @@ def test_spdhg_set_norms(self): self.assertListEqual(spdhg._norms, [1]*len(self.A2)) def test_spdhg_check_convergence(self): + + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A) - self.assertTrue(spdhg.check_convergence()) gamma = 3.7 rho = 0.9 - spdhg.set_step_sizes_from_ratio(gamma, rho) + rule = SPDHGStepSizesFromRatio(gamma, rho) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) self.assertTrue(spdhg.check_convergence()) gamma = 3.7 rho = 100 - spdhg.set_step_sizes_from_ratio(gamma, rho) + rule = SPDHGStepSizesFromRatio(gamma, rho) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) self.assertFalse(spdhg.check_convergence()) - spdhg.set_step_sizes(sigma=[1]*self.subsets, tau=100) + rule = SPDHGConstantStepSize(step_size=(100, [1]*self.subsets)) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) self.assertFalse(spdhg.check_convergence()) - spdhg.set_step_sizes(sigma=[1]*self.subsets, tau=None) + rule = SPDHGConstantStepSize(step_size=(None, [1]*self.subsets)) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) self.assertTrue(spdhg.check_convergence()) - spdhg.set_step_sizes(sigma=None, tau=100) + rule = SPDHGConstantStepSize(step_size=(100, None)) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) self.assertTrue(spdhg.check_convergence()) @@ -1709,7 +1392,7 @@ def test_compare_with_PDHG(self): sigma = 1./normK tau = 1./normK - pdhg = PDHG(f=F, g=G, operator=K, tau=tau, sigma=sigma, + pdhg = PDHG(f=F, g=G, operator=K, step_size=(tau, sigma), update_objective_interval=10) pdhg.run(500, verbose=0) @@ -1771,8 +1454,7 @@ def test_PD3O_PDHG_denoising_1_iteration(self): G = 0.5 * L2NormSquared(b=self.data) sigma = 1./norm_op tau = 1./norm_op - pdhg = PDHG(f=F, g=G, operator=operator, tau=tau, - sigma=sigma, update_objective_interval=100) + pdhg = PDHG(f=F, g=G, operator=operator, step_size=(tau,sigma), update_objective_interval=100) pdhg.run(1) # setup PD3O denoising (F=ZeroFunction) diff --git a/Wrappers/Python/test/test_stepsizes.py b/Wrappers/Python/test/test_stepsizes.py index 2d774f98d6..94109c507e 100644 --- a/Wrappers/Python/test/test_stepsizes.py +++ b/Wrappers/Python/test/test_stepsizes.py @@ -1,13 +1,21 @@ -from cil.optimisation.algorithms import SIRT, GD, ISTA, FISTA -from cil.optimisation.functions import LeastSquares, IndicatorBox, ZeroFunction +from cil.optimisation.algorithms import SIRT, GD, ISTA, FISTA, PDHG, SPDHG +from cil.optimisation.functions import BlockFunction, LeastSquares, IndicatorBox, ZeroFunction, L2NormSquared from cil.framework import ImageGeometry, VectorGeometry, VectorData -from cil.optimisation.operators import IdentityOperator, MatrixOperator, LinearOperator +from cil.optimisation.operators import BlockOperator, IdentityOperator, MatrixOperator, LinearOperator -from cil.optimisation.utilities import Sensitivity, AdaptiveSensitivity, Preconditioner, ConstantStepSize, ArmijoStepSizeRule, BarzilaiBorweinStepSizeRule +from cil.optimisation.utilities import Sensitivity, AdaptiveSensitivity, Preconditioner, ConstantStepSize, ArmijoStepSizeRule, BarzilaiBorweinStepSizeRule, PDHGStronglyConvexUpdate, PDHGConstantStepSize, PDHGAdaptiveStepSize2013, PDHGAdaptiveStepSize2015, PDHGBayesOptimisationStepSize, StepSizeRule +from cil.optimisation.utilities import SPDHGConstantStepSize, SPDHGBayesOptimisationStepSize, SPDHGStepSizesFromRatio, Sampler +from cil.optimisation.utilities import SPDHGAdaptiveStepSizeBalancing, SPDHGAdaptiveStepSizeAngle import numpy as np +from cil.utilities import dataexample from testclass import CCPiTestClass +from utils import has_skopt +import unittest from unittest.mock import MagicMock +from unittest.mock import Mock +from numbers import Number +import warnings class TestStepSizes(CCPiTestClass): @@ -21,7 +29,7 @@ def test_step_sizes_called(self): step_size_test.get_step_size = MagicMock(return_value=.1) f = LeastSquares(A=A, b=data, c=0.5) alg = GD(initial=ig.allocate('random', seed=10), f=f, step_size=step_size_test, - update_objective_interval=1) + update_objective_interval=1) alg.run(5) @@ -30,24 +38,26 @@ def test_step_sizes_called(self): step_size_test = ConstantStepSize(3) step_size_test.get_step_size = MagicMock(return_value=.1) alg = ISTA(initial=ig.allocate('random', seed=10), f=f, g=IndicatorBox(lower=0), step_size=step_size_test, - update_objective_interval=1) + update_objective_interval=1) alg.run(5) self.assertEqual(len(step_size_test.get_step_size.mock_calls), 5) step_size_test = ConstantStepSize(3) step_size_test.get_step_size = MagicMock(return_value=.1) alg = FISTA(initial=ig.allocate('random', seed=10), f=f, g=IndicatorBox(lower=0), step_size=step_size_test, - update_objective_interval=1) + update_objective_interval=1) alg.run(5) self.assertEqual(len(step_size_test.get_step_size.mock_calls), 5) + class TestStepSizeConstant(CCPiTestClass): def test_constant(self): test_stepsize = ConstantStepSize(0.3) self.assertEqual(test_stepsize.step_size, 0.3) + class TestStepSizeArmijo(CCPiTestClass): - + def setUp(self): self.ig = VectorGeometry(2) self.data = self.ig.allocate('random', seed=3) @@ -55,10 +65,10 @@ def setUp(self): self.A = MatrixOperator(np.diag([1., 1.])) self.f = LeastSquares(self.A, self.data) - def test_armijo_init(self): - test_stepsize = ArmijoStepSizeRule(alpha=1e3, beta=0.4, max_iterations=40, warmstart=False) - self.assertFalse(test_stepsize.warmstart) + test_stepsize = ArmijoStepSizeRule( + alpha=1e3, beta=0.4, max_iterations=40, warmstart=False) + self.assertFalse(test_stepsize.warmstart) self.assertEqual(test_stepsize.alpha_orig, 1e3) self.assertEqual(test_stepsize.beta, 0.4) self.assertEqual(test_stepsize.max_iterations, 40) @@ -71,10 +81,11 @@ def test_armijo_init(self): 2 * np.log10(1e6) / np.log10(2))) def test_armijo_calculation(self): - test_stepsize = ArmijoStepSizeRule(alpha=8, beta=0.5, max_iterations=100, warmstart=False) + test_stepsize = ArmijoStepSizeRule( + alpha=8, beta=0.5, max_iterations=100, warmstart=False) alg = GD(initial=self.ig.allocate(0), f=self.f, - update_objective_interval=1, step_size=test_stepsize) + update_objective_interval=1, step_size=test_stepsize) alg.gradient_update = self.ig.allocate(-1) step_size = test_stepsize.get_step_size(alg) self.assertAlmostEqual(step_size, 4) @@ -88,10 +99,11 @@ def test_armijo_calculation(self): self.assertAlmostEqual(step_size, 2) def test_armijo_ISTA_and_FISTA(self): - test_stepsize = ArmijoStepSizeRule(alpha=8, beta=0.5, max_iterations=100, warmstart=False) + test_stepsize = ArmijoStepSizeRule( + alpha=8, beta=0.5, max_iterations=100, warmstart=False) alg = ISTA(initial=self.ig.allocate(0), f=self.f, g=IndicatorBox(lower=0), - update_objective_interval=1, step_size=test_stepsize) + update_objective_interval=1, step_size=test_stepsize) alg.gradient_update = self.ig.allocate(-1) step_size = test_stepsize.get_step_size(alg) self.assertAlmostEqual(step_size, 4) @@ -105,7 +117,7 @@ def test_armijo_ISTA_and_FISTA(self): self.assertAlmostEqual(step_size, 2) alg = FISTA(initial=self.ig.allocate(0), f=self.f, g=IndicatorBox(lower=0), - update_objective_interval=1, step_size=test_stepsize) + update_objective_interval=1, step_size=test_stepsize) alg.gradient_update = self.ig.allocate(-1) step_size = test_stepsize.get_step_size(alg) self.assertAlmostEqual(step_size, 4) @@ -119,12 +131,12 @@ def test_armijo_ISTA_and_FISTA(self): self.assertAlmostEqual(step_size, 2) def test_warmstart_true(self): - + rule = ArmijoStepSizeRule(warmstart=True, alpha=5000) self.assertTrue(rule.warmstart) self.assertTrue(rule.alpha_orig == 5000) alg = GD(initial=self.ig.allocate(0), f=self.f, - update_objective_interval=1, step_size=rule) + update_objective_interval=1, step_size=rule) alg.update() self.assertFalse(rule.alpha == 5000) @@ -133,10 +145,11 @@ def test_warmstart_false(self): self.assertFalse(rule.warmstart) self.assertTrue(rule.alpha_orig == 5000) alg = GD(initial=self.ig.allocate(0), f=self.f, - update_objective_interval=1, step_size=rule) + update_objective_interval=1, step_size=rule) alg.update() self.assertTrue(rule.alpha_orig == 5000) - self.assertFalse(rule.alpha_orig == rule.alpha) + self.assertFalse(rule.alpha_orig == rule.alpha) + class TestStepSizeBB(CCPiTestClass): def test_bb(self): @@ -149,110 +162,117 @@ def test_bb(self): Aop = MatrixOperator(A) bop = VectorData(b) - ig=Aop.domain + ig = Aop.domain initial = ig.allocate() f = LeastSquares(Aop, b=bop, c=0.5) - - ss_rule=BarzilaiBorweinStepSizeRule(2 ) + + ss_rule = BarzilaiBorweinStepSizeRule(2) self.assertEqual(ss_rule.mode, 'short') self.assertEqual(ss_rule.initial, 2) self.assertEqual(ss_rule.adaptive, True) self.assertEqual(ss_rule.stabilisation_param, np.inf) - - #Check the right errors are raised for incorrect parameters - + + # Check the right errors are raised for incorrect parameters + with self.assertRaises(TypeError): - ss_rule=BarzilaiBorweinStepSizeRule(2,'short',-4, ) + ss_rule = BarzilaiBorweinStepSizeRule(2, 'short', -4, ) with self.assertRaises(TypeError): - ss_rule=BarzilaiBorweinStepSizeRule(2,'long', 'banana', ) + ss_rule = BarzilaiBorweinStepSizeRule(2, 'long', 'banana', ) with self.assertRaises(ValueError): - ss_rule=BarzilaiBorweinStepSizeRule(2, 'banana',3 ) - - - #Check stabilisation parameter unchanged if fixed - ss_rule=BarzilaiBorweinStepSizeRule(2, 'long',3 ) + ss_rule = BarzilaiBorweinStepSizeRule(2, 'banana', 3) + + # Check stabilisation parameter unchanged if fixed + ss_rule = BarzilaiBorweinStepSizeRule(2, 'long', 3) self.assertEqual(ss_rule.mode, 'long') self.assertFalse(ss_rule.adaptive) alg = GD(initial=initial, f=f, step_size=ss_rule) - self.assertEqual(ss_rule.stabilisation_param,3) + self.assertEqual(ss_rule.stabilisation_param, 3) alg.run(2) - self.assertEqual(ss_rule.stabilisation_param,3) - - #Check infinity can be passed - ss_rule=BarzilaiBorweinStepSizeRule(2, 'short',"off" ) + self.assertEqual(ss_rule.stabilisation_param, 3) + + # Check infinity can be passed + ss_rule = BarzilaiBorweinStepSizeRule(2, 'short', "off") self.assertEqual(ss_rule.mode, 'short') self.assertFalse(ss_rule.adaptive) - self.assertEqual(ss_rule.stabilisation_param,np.inf) + self.assertEqual(ss_rule.stabilisation_param, np.inf) alg = GD(initial=initial, f=f, step_size=ss_rule) alg.run(2) - + n = 5 m = 5 A = np.eye(5).astype('float32') - b = (np.array([.5,.5,.5,.5,.5])).astype('float32') + b = (np.array([.5, .5, .5, .5, .5])).astype('float32') Aop = MatrixOperator(A) bop = VectorData(b) - ig=Aop.domain + ig = Aop.domain initial = ig.allocate(0) f = LeastSquares(Aop, b=bop, c=0.5) - ss_rule=BarzilaiBorweinStepSizeRule(0.22, 'long',np.inf ) + ss_rule = BarzilaiBorweinStepSizeRule(0.22, 'long', np.inf) alg = GD(initial=initial, f=f, step_size=ss_rule) self.assertFalse(ss_rule.is_short) - #Check the initial step size was used + # Check the initial step size was used alg.run(1) - self.assertNumpyArrayAlmostEqual( np.array([.11,.11,.11,.11,.11]), alg.x.as_array() ) + self.assertNumpyArrayAlmostEqual( + np.array([.11, .11, .11, .11, .11]), alg.x.as_array()) self.assertFalse(ss_rule.is_short) - #check long + # check long alg.run(1) - x_change= np.array([.11,.11,.11,.11,.11])-np.array([0,0,0,0,0]) - grad_change = -np.array([.39,.39,.39,.39,.39])+np.array([.5,.5,.5,.5,.5]) - step= x_change.dot(x_change)/x_change.dot(grad_change) - self.assertNumpyArrayAlmostEqual( np.array([.11,.11,.11,.11,.11])+step*np.array([.39,.39,.39,.39,.39]), alg.x.as_array() ) + x_change = np.array([.11, .11, .11, .11, .11]) - \ + np.array([0, 0, 0, 0, 0]) + grad_change = -np.array([.39, .39, .39, .39, .39]) + \ + np.array([.5, .5, .5, .5, .5]) + step = x_change.dot(x_change)/x_change.dot(grad_change) + self.assertNumpyArrayAlmostEqual(np.array( + [.11, .11, .11, .11, .11])+step*np.array([.39, .39, .39, .39, .39]), alg.x.as_array()) self.assertFalse(ss_rule.is_short) - - ss_rule=BarzilaiBorweinStepSizeRule(0.22, 'short',np.inf ) + + ss_rule = BarzilaiBorweinStepSizeRule(0.22, 'short', np.inf) alg = GD(initial=initial, f=f, step_size=ss_rule) self.assertTrue(ss_rule.is_short) - #Check the initial step size was used + # Check the initial step size was used alg.run(1) - self.assertNumpyArrayAlmostEqual( np.array([.11,.11,.11,.11,.11]), alg.x.as_array() ) + self.assertNumpyArrayAlmostEqual( + np.array([.11, .11, .11, .11, .11]), alg.x.as_array()) self.assertTrue(ss_rule.is_short) - #check short + # check short alg.run(1) - x_change= np.array([.11,.11,.11,.11,.11])-np.array([0,0,0,0,0]) - grad_change = -np.array([.39,.39,.39,.39,.39])+np.array([.5,.5,.5,.5,.5]) - step= x_change.dot(grad_change)/grad_change.dot(grad_change) - self.assertNumpyArrayAlmostEqual( np.array([.11,.11,.11,.11,.11])+step*np.array([.39,.39,.39,.39,.39]), alg.x.as_array() ) + x_change = np.array([.11, .11, .11, .11, .11]) - \ + np.array([0, 0, 0, 0, 0]) + grad_change = -np.array([.39, .39, .39, .39, .39]) + \ + np.array([.5, .5, .5, .5, .5]) + step = x_change.dot(grad_change)/grad_change.dot(grad_change) + self.assertNumpyArrayAlmostEqual(np.array( + [.11, .11, .11, .11, .11])+step*np.array([.39, .39, .39, .39, .39]), alg.x.as_array()) self.assertTrue(ss_rule.is_short) - - #check stop iteration - ss_rule=BarzilaiBorweinStepSizeRule(1, 'long',np.inf ) + + # check stop iteration + ss_rule = BarzilaiBorweinStepSizeRule(1, 'long', np.inf) alg = GD(initial=initial, f=f, step_size=ss_rule) alg.run(500) self.assertEqual(alg.iteration, 1) - - #check adaptive - ss_rule=BarzilaiBorweinStepSizeRule(0.001, 'long',"auto") + + # check adaptive + ss_rule = BarzilaiBorweinStepSizeRule(0.001, 'long', "auto") alg = GD(initial=initial, f=f, step_size=ss_rule) self.assertEqual(ss_rule.stabilisation_param, np.inf) alg.run(2) self.assertNotEqual(ss_rule.stabilisation_param, np.inf) - - #check stops being adaptive - - ss_rule=BarzilaiBorweinStepSizeRule(0.0000001, 'long',"auto" ) + + # check stops being adaptive + + ss_rule = BarzilaiBorweinStepSizeRule(0.0000001, 'long', "auto") alg = GD(initial=initial, f=f, step_size=ss_rule) self.assertEqual(ss_rule.stabilisation_param, np.inf) alg.run(4) self.assertNotEqual(ss_rule.stabilisation_param, np.inf) - a=ss_rule.stabilisation_param + a = ss_rule.stabilisation_param alg.run(1) self.assertEqual(ss_rule.stabilisation_param, a) - - #Test alternating - ss_rule=BarzilaiBorweinStepSizeRule(0.0000001, 'alternate',"auto" ) + + # Test alternating + ss_rule = BarzilaiBorweinStepSizeRule(0.0000001, 'alternate', "auto") alg = GD(initial=initial, f=f, step_size=ss_rule) self.assertFalse(ss_rule.is_short) alg.run(2) @@ -262,10 +282,1337 @@ def test_bb(self): alg.run(1) self.assertTrue(ss_rule.is_short) + +class TestPDHGConstantStepSize(CCPiTestClass): + + # TODO: remove when deprecated parameters are removed from PDHG + def test_deprecated_parameters(self): + + with self.assertWarns(DeprecationWarning): + pdhg = PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=IdentityOperator(ImageGeometry(2, 2)), + sigma=0.5) + self.assertEqual(pdhg.sigma, 0.5) + self.assertEqual(pdhg.step_size_rule.sigma, 0.5) + self.assertEqual(pdhg.step_size_rule.tau, 1. / + (0.5 * pdhg.operator.norm()**2)) + self.assertTrue(isinstance(pdhg.step_size_rule, PDHGConstantStepSize)) + + with self.assertWarns(DeprecationWarning): + pdhg = PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=IdentityOperator(ImageGeometry(2, 2)), + tau=0.5) + self.assertEqual(pdhg.tau, 0.5) + self.assertEqual(pdhg.step_size_rule.tau, 0.5) + self.assertEqual(pdhg.step_size_rule.sigma, 1. / + (0.5 * pdhg.operator.norm()**2)) + self.assertTrue(isinstance(pdhg.step_size_rule, PDHGConstantStepSize)) + + with self.assertWarns(DeprecationWarning): + pdhg = PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=IdentityOperator(ImageGeometry(2, 2)), + sigma=0.5, tau=0.5) + self.assertEqual(pdhg.sigma, 0.5) + self.assertEqual(pdhg.tau, 0.5) + self.assertTrue(isinstance(pdhg.step_size_rule, PDHGConstantStepSize)) + + def test_PDHG_step_sizes(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = 3*IdentityOperator(ig) + + # check if sigma, tau are None + pdhg = PDHG(f=f, g=g, operator=operator) + self.assertAlmostEqual(pdhg.sigma, 1./operator.norm()) + self.assertAlmostEqual(pdhg.tau, 1./operator.norm()) + self.assertTrue(isinstance(pdhg.step_size_rule, PDHGConstantStepSize)) + + # check if sigma is negative + with self.assertRaises(ValueError): + pdhg = PDHG(f=f, g=g, operator=operator, + step_size=(None, -1)) + + # check if tau is negative + with self.assertRaises(ValueError): + pdhg = PDHG(f=f, g=g, operator=operator, step_size=(-1, None)) + + # check if tau is None + sigma = 3.0 + pdhg = PDHG(f=f, g=g, operator=operator, step_size=(None, sigma)) + self.assertAlmostEqual(pdhg.sigma, sigma) + self.assertAlmostEqual(pdhg.tau, 1./(sigma * operator.norm()**2)) + self.assertTrue(isinstance(pdhg.step_size_rule, PDHGConstantStepSize)) + + # check if sigma is None + tau = 3.0 + pdhg = PDHG(f=f, g=g, operator=operator, step_size=(tau, None)) + self.assertAlmostEqual(pdhg.tau, tau) + self.assertAlmostEqual(pdhg.sigma, 1./(tau * operator.norm()**2)) + self.assertTrue(isinstance(pdhg.step_size_rule, PDHGConstantStepSize)) + + # check if sigma/tau are not None + tau = 1.0 + sigma = 1.0 + pdhg = PDHG(f=f, g=g, operator=operator, step_size=(tau, sigma)) + self.assertAlmostEqual(pdhg.tau, tau) + self.assertAlmostEqual(pdhg.sigma, sigma) + self.assertTrue(isinstance(pdhg.step_size_rule, PDHGConstantStepSize)) + + # check sigma/tau as arrays, sigma wrong shape + ig1 = ImageGeometry(2, 2) + sigma = ig1.allocate() + with self.assertRaises(ValueError): + pdhg = PDHG(f=f, g=g, operator=operator, step_size=(None, sigma)) + + # check sigma/tau as arrays, tau wrong shape + tau = ig1.allocate() + with self.assertRaises(ValueError): + pdhg = PDHG(f=f, g=g, operator=operator, step_size=(tau, None)) + + # check sigma not Number or object with correct shape + with self.assertRaises(AttributeError): + pdhg = PDHG(f=f, g=g, operator=operator, + step_size=("sigma", None)) + + # check tau not Number or object with correct shape + with self.assertRaises(AttributeError): + pdhg = PDHG(f=f, g=g, operator=operator, + step_size=("tau", None)) + + # check warning message if condition is not satisfied + sigma = 4/operator.norm() + tau = 1/3 + with self.assertWarnsRegex(UserWarning, "Convergence criterion"): + pdhg = PDHG(f=f, g=g, operator=operator, step_size=(tau, sigma)) + + # check no warning message if check convergence is false + sigma = 4/operator.norm() + tau = 1/3 + with warnings.catch_warnings(record=True) as warnings_log: + pdhg = PDHG(f=f, g=g, operator=operator, step_size=( + tau, sigma), check_convergence=False) + self.assertEqual(warnings_log, []) + + # check no warning message if condition is satisfied + sigma = 1/operator.norm() + tau = 1/3 + with warnings.catch_warnings(record=True) as warnings_log: + warnings.simplefilter("always") + pdhg = PDHG(f=f, g=g, operator=operator, step_size=[tau, sigma]) + self.assertTrue(pdhg.sigma * pdhg.tau * pdhg.operator.norm()**2 < 4/3) + self.assertTrue(isinstance(pdhg.sigma, Number)) + self.assertTrue(isinstance(pdhg.tau, Number)) + self.assertEqual(warnings_log, []) + + def test_step_size_and_deprecated_sigma_tau_raises_valueerror(self): + # Passing both `step_size` and the deprecated `sigma`/`tau` must raise a + # clear ValueError. Regression: this branch previously passed keyword + # arguments to ValueError, which raised a TypeError instead. + operator = IdentityOperator(ImageGeometry(2, 2)) + with self.assertRaises(ValueError): + PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=operator, + step_size=(0.1, 0.1), sigma=0.5) + with self.assertRaises(ValueError): + PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=operator, + step_size=(0.1, 0.1), tau=0.5) + + def test_incompatible_gd_rule_raises_valueerror(self): + # A GD-style step-size rule returns a single scalar and has no + # get_initial_step_size, so it is not compatible with PDHG and must + # raise a clear ValueError rather than an opaque AttributeError. + operator = IdentityOperator(ImageGeometry(2, 2)) + with self.assertRaises(ValueError): + PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=operator, + step_size=ConstantStepSize(0.1)) + + def test_wrong_shape_rule_raises_valueerror(self): + # A rule that provides get_initial_step_size but returns wrong-shaped + # step sizes (here a list where PDHG expects a scalar/array) must raise + # a clear ValueError from the PDHG step-size validation. + class _BadRule(StepSizeRule): + def get_initial_step_size(self, algorithm): + return 0.1, [0.1, 0.1] + + def get_step_size(self, algorithm): + return 0.1, [0.1, 0.1] + + operator = IdentityOperator(ImageGeometry(2, 2)) + with self.assertRaises(ValueError): + PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=operator, + step_size=_BadRule()) + + +class TestStepSizePDHGStronglyConvex(CCPiTestClass): + + # TODO: remove when deprecated parameters are removed from PDHG + def test_deprecated_parameters(self): + with self.assertWarns(DeprecationWarning): + pdhg = PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=IdentityOperator(ImageGeometry(2, 2)), + gamma_g=0.5) + + with self.assertWarns(DeprecationWarning): + pdhg = PDHG(f=ZeroFunction(), g=ZeroFunction(), operator=IdentityOperator(ImageGeometry(2, 2)), + gamma_fconj=0.5) + self.assertEqual(pdhg.step_size_rule.gamma_fconj, 0.5) + + def test_init_invalid_initial_step_size_length(self): + # initial_step_size must be a length-two list/tuple. A wrong length must + # raise ValueError (regression: the error path previously referenced an + # undefined `step_size` and raised NameError instead). + with self.assertRaises(ValueError): + PDHGStronglyConvexUpdate(initial_step_size=(1.0,), gamma_g=0.5) + with self.assertRaises(ValueError): + PDHGStronglyConvexUpdate( + initial_step_size=(1.0, 2.0, 3.0), gamma_g=0.5) + + def test_PDHG_strongly_convex_gamma_g(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + # sigma, tau + sigma = 1.0 + tau = 1.0 + + step_size_rule = PDHGStronglyConvexUpdate( + initial_step_size=(tau, sigma), gamma_g=0.5) + pdhg = PDHG(f=f, g=g, operator=operator, + step_size=step_size_rule) + pdhg.run(1, verbose=0) + self.assertAlmostEqual( + pdhg.theta, 1.0 / np.sqrt(1 + 2 * step_size_rule.gamma_g * tau)) + self.assertAlmostEqual(pdhg.tau, tau * pdhg.theta) + self.assertAlmostEqual(pdhg.sigma, sigma / pdhg.theta) + pdhg.run(4, verbose=0) + self.assertNotEqual(pdhg.sigma, sigma) + self.assertNotEqual(pdhg.tau, tau) + + # check negative strongly convex constant + with self.assertRaises(ValueError): + pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, + gamma_g=-0.5) + + # check strongly convex constant not a number + with self.assertRaises(ValueError): + pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, + gamma_g="-0.5") + + def test_PDHG_strongly_convex_gamma_fcong(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + # sigma, tau + sigma = 1.0 + tau = 1.0 + step_size_rule = PDHGStronglyConvexUpdate( + initial_step_size=(tau, sigma), gamma_fconj=0.5) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=step_size_rule) + pdhg.run(1, verbose=0) + self.assertEqual(pdhg.theta, 1.0 / np.sqrt(1 + + 2 * step_size_rule .gamma_fconj * sigma)) + self.assertEqual(pdhg.tau, tau / pdhg.theta) + self.assertEqual(pdhg.sigma, sigma * pdhg.theta) + pdhg.run(4, verbose=0) + self.assertNotEqual(pdhg.sigma, sigma) + self.assertNotEqual(pdhg.tau, tau) + + # check negative strongly convex constant + with self.assertRaises(ValueError): + pdhg = PDHG(f=f, g=g, operator=operator, sigma=sigma, tau=tau, + gamma_fconj=-0.5) + + # check strongly convex constant not a number + with self.assertRaises(ValueError): + + step_size_rule = PDHGStronglyConvexUpdate(gamma_fconj="-0.5") + + def test_PDHG_strongly_convex_both_fconj_and_g(self): + + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + with self.assertRaises(NotImplementedError): + pdhg = PDHG(f=f, g=g, operator=operator, + gamma_g=0.5, gamma_fconj=0.5) + pdhg.run(verbose=0) + + +class TestPDHGAdaptive2013(CCPiTestClass): + + def test_init(self): + rule = PDHGAdaptiveStepSize2013(initial_step_size=[1.0, 2.0]) + self.assertEqual(rule.initial_step_size, [1.0, 2.0]) + self.assertAlmostEqual(rule.alpha, 0.95) + self.assertEqual(rule.gamma, 0.9) + self.assertEqual(rule.inner_iterations, 50) + self.assertEqual(rule.tolerance, 1e-06) + self.assertEqual(rule.count, 0) + self.assertEqual(rule.beta, 0.95) + self.assertEqual(rule.delta, 1.5) + + def test_init_invalid(self): + with self.assertRaises(ValueError): + PDHGAdaptiveStepSize2013(initial_step_size=[1.0]) + + def test_no_false_nonconvergence_warning(self): + # Regression: the non-convergence warning used `k == inner_iterations-1`, + # which false-fired when backtracking legitimately converged on the last + # inner iteration. With inner_iterations=1 a converging step (b <= 1) + # breaks at that last index, so no "did not converge" warning should be + # logged. + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], inner_iterations=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + # Force backtracking to report convergence (b <= 1) on the first (and + # only, hence last) inner iteration. + rule._calculate_backtracking = MagicMock(return_value=0.5) + + logger = 'cil.optimisation.utilities.StepSizeMethods' + with self.assertNoLogs(logger, level='WARNING'): + pdhg.run(3) + + def test_initial_step_size_defaults(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013(initial_step_size=[None, None]) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + tau, sigma = rule.get_initial_step_size(pdhg) + self.assertEqual(tau, 10) + self.assertEqual(sigma, 10) + + rule = PDHGAdaptiveStepSize2013(initial_step_size=[3.2, None]) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + tau, sigma = rule.get_initial_step_size(pdhg) + self.assertEqual(tau, 3.2) + self.assertEqual(sigma, 10) + + rule = PDHGAdaptiveStepSize2013(initial_step_size=[None, 3.2]) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + tau, sigma = rule.get_initial_step_size(pdhg) + self.assertEqual(tau, 10) + self.assertEqual(sigma, 3.2) + + def test_backtracking_calculation(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0, gamma=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + rule.x_resid = operator.domain.allocate(0) + rule.y_resid = operator.range.allocate(0) + rule.x_store = operator.domain.allocate(0) + rule.y_old = operator.range.allocate(0) + self.assertEqual(rule.gamma, 1) + self.assertEqual(pdhg.sigma, 1) + self.assertEqual(pdhg.tau, 1) + b = rule._calculate_backtracking(pdhg) + self.assertEqual(rule.x_resid.norm(), 3) + self.assertEqual(rule.y_resid.norm(), 3) + self.assertNumpyArrayAlmostEqual( + rule.x_resid.as_array(), pdhg.y_tmp.as_array()) + self.assertEqual(rule.y_resid.dot(pdhg.y_tmp), 9) + self.assertEqual(b, 1) + + def test_backtracking_no_change_returns_zero(self): + # When the iterate does not change (x == x_old and y == y_old) the + # backtracking denominator is zero; ensure we return 0 (accept) rather + # than 0/0 = nan, which would poison the step sizes. + ig = ImageGeometry(3, 3) + f = L2NormSquared(b=ig.allocate('random', seed=3)) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0, gamma=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(1) + pdhg.x = operator.domain.allocate(1) # x == x_old -> no primal change + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + rule.x_resid = operator.domain.allocate(0) + rule.y_resid = operator.range.allocate(0) + rule.y_old = operator.range.allocate(1) # y == y_old -> no dual change + + b = rule._calculate_backtracking(pdhg) + self.assertTrue(np.isfinite(b)) + self.assertEqual(b, 0.0) + self.assertLessEqual(b, 1) # caller would accept, not backtrack + + def test_backtracking(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0, gamma=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[2.0, 1.5, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 3 + rule.d_norm = 3 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn("Before adaptive step-size step", log_output) + self.assertIn("Backtracking step", log_output) + self.assertIn("Finished backtracking step", log_output) + + self.assertAlmostEqual(pdhg.sigma, 0.95**2 / (2*1.5)) + self.assertAlmostEqual(pdhg.tau, 0.95**2 / (2*1.5)) + self.assertEqual(mock_func.call_count, 3) + + def test_changing_ratio(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, gamma=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[0.5, 0.5, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 1 + rule.d_norm = 4 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + # Assertions + self.assertIn("p_norm < ", log_output) + self.assertAlmostEqual(pdhg.sigma, 1/0.05) + self.assertAlmostEqual(pdhg.tau, 0.05) + self.assertAlmostEqual(rule.alpha, 0.95*0.95) + self.assertEqual(mock_func.call_count, 1) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, gamma=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[0.5, 0.5, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 1 + rule.d_norm = 1 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn("No change", log_output) + self.assertAlmostEqual(pdhg.sigma, 1) + self.assertAlmostEqual(pdhg.tau, 1) + self.assertAlmostEqual(rule.alpha, 0.95) + self.assertEqual(mock_func.call_count, 1) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, gamma=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[0.5, 0.5, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 4 + rule.d_norm = 1 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + print(log_output) + # Assertions + self.assertIn("< d_norm ", log_output) + self.assertAlmostEqual(pdhg.tau, 1/0.05) + self.assertAlmostEqual(pdhg.sigma, 0.05) + self.assertAlmostEqual(rule.alpha, 0.95*0.95) + self.assertEqual(mock_func.call_count, 1) + + def test_inner_iteration_stopping_criterion(self): + + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0, gamma=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(return_value=1.5) + + rule._calculate_backtracking = mock_func + + # Capture logs + with self.assertLogs(level='WARNING') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn("Backtracking step did not converge", log_output) + + self.assertEqual(mock_func.call_count, 51) + + def test_stopping_criterion(self): + + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + + rule = PDHGAdaptiveStepSize2013( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, gamma=1, auto_stop=True) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(return_value=0.5) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 1 + rule.d_norm = 1 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + for i in range(15): + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn('Finished backtracking step', log_output) + self.assertIn("No change", log_output) + self.assertAlmostEqual(pdhg.sigma, 1) + self.assertAlmostEqual(pdhg.tau, 1) + self.assertAlmostEqual(rule.alpha, 0.95) + self.assertEqual(rule.count, 11) + self.assertEqual(rule.adaptive, False) + self.assertEqual(mock_func.call_count, 11) + with self.assertRaises(AttributeError): + rule.x_resid + with self.assertRaises(AttributeError): + rule.y_resid + with self.assertRaises(AttributeError): + rule.y_old + + +class TestPDHGAdaptive2015(CCPiTestClass): + + def test_init(self): + rule = PDHGAdaptiveStepSize2015(initial_step_size=[1.0, 2.0]) + self.assertEqual(rule.initial_step_size, [1.0, 2.0]) + self.assertAlmostEqual(rule.alpha, 0.95) + self.assertEqual(rule.c, 0.9) + self.assertEqual(rule.inner_iterations, 50) + self.assertEqual(rule.tolerance, 1e-06) + self.assertEqual(rule.count, 0) + self.assertEqual(rule.eta, 0.95) + self.assertTrue(rule.adaptive) + self.assertTrue(rule.auto_stop) + + def test_init_invalid(self): + with self.assertRaises(ValueError): + PDHGAdaptiveStepSize2015(initial_step_size=[1.0]) + + def test_initial_step_size_defaults(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2015(initial_step_size=[None, None]) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + tau, sigma = rule.get_initial_step_size(pdhg) + self.assertEqual(tau, 10) + self.assertEqual(sigma, 10) + + rule = PDHGAdaptiveStepSize2015(initial_step_size=[3.2, None]) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + tau, sigma = rule.get_initial_step_size(pdhg) + self.assertEqual(tau, 3.2) + self.assertEqual(sigma, 10) + + rule = PDHGAdaptiveStepSize2015(initial_step_size=[None, 3.2]) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + tau, sigma = rule.get_initial_step_size(pdhg) + self.assertEqual(tau, 10) + self.assertEqual(sigma, 3.2) + + def test_backtracking_calculation(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2015( + initial_step_size=[1.0, 1.0], initial_alpha=0, c=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + rule.x_resid = operator.domain.allocate(0) + rule.y_resid = operator.range.allocate(0) + rule.x_store = operator.domain.allocate(0) + rule.y_old = operator.range.allocate(0) + self.assertEqual(rule.c, 1) + self.assertEqual(pdhg.sigma, 1) + self.assertEqual(pdhg.tau, 1) + b = rule._calculate_backtracking(pdhg) + self.assertEqual(rule.x_resid.norm(), 3) + self.assertEqual(rule.y_resid.norm(), 3) + self.assertNumpyArrayAlmostEqual( + rule.x_resid.as_array(), pdhg.y_tmp.as_array()) + self.assertEqual(rule.y_resid.dot(pdhg.y_tmp), 9) + self.assertEqual(b, -18) + + def test_backtracking(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2015( + initial_step_size=[1.0, 1.0], initial_alpha=0) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[-2, -2, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 3 + rule.d_norm = 3 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn("Before adaptive step-size step", log_output) + self.assertIn("Backtracking step", log_output) + self.assertIn("Finished backtracking step", log_output) + + self.assertAlmostEqual(pdhg.sigma, 0.5**2) + self.assertAlmostEqual(pdhg.tau, 0.5**2) + self.assertEqual(mock_func.call_count, 3) + + def test_changing_ratio(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2015( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, c=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[0.5, 0.5, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 1 + rule.d_norm = 4 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + # Assertions + self.assertIn("< d_norm", log_output) + self.assertAlmostEqual(pdhg.sigma, 1/0.05) + self.assertAlmostEqual(pdhg.tau, 0.05) + self.assertAlmostEqual(rule.alpha, 0.95*0.95) + self.assertEqual(mock_func.call_count, 1) + + rule = PDHGAdaptiveStepSize2015( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, c=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[0.5, 0.5, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 1 + rule.d_norm = 1 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn("No change", log_output) + self.assertAlmostEqual(pdhg.sigma, 1) + self.assertAlmostEqual(pdhg.tau, 1) + self.assertAlmostEqual(rule.alpha, 0.95) + self.assertEqual(mock_func.call_count, 1) + + rule = PDHGAdaptiveStepSize2015( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, c=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(side_effect=[0.5, 0.5, 0.5]) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 4 + rule.d_norm = 1 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + print(log_output) + # Assertions + self.assertIn("< p_norm ", log_output) + self.assertAlmostEqual(pdhg.tau, 1/0.05) + self.assertAlmostEqual(pdhg.sigma, 0.05) + self.assertAlmostEqual(rule.alpha, 0.95*0.95) + self.assertEqual(mock_func.call_count, 1) + + def test_inner_iteration_stopping_criterion(self): + + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGAdaptiveStepSize2015( + initial_step_size=[1.0, 1.0], initial_alpha=0, c=1) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(return_value=-1) + + rule._calculate_backtracking = mock_func + + # Capture logs + with self.assertLogs(level='WARNING') as log: + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn("Backtracking step did not converge", log_output) + + self.assertEqual(mock_func.call_count, 51) + + def test_stopping_criterion(self): + + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + + rule = PDHGAdaptiveStepSize2015( + initial_step_size=[1.0, 1.0], initial_alpha=0.95, c=1, auto_stop=True) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + pdhg.x_old = operator.domain.allocate(0) + pdhg.y_old = operator.range.allocate(0) + pdhg.x = operator.domain.allocate(1) + pdhg.y = operator.range.allocate(1) + pdhg.y_tmp = operator.range.allocate(0) + + mock_func = Mock(return_value=0.5) + + rule._calculate_backtracking = mock_func + + def mock_pnorm_dnorm(algorithm): + rule.x_resid.fill(1) + rule.y_resid.fill(1) + rule.p_norm = 1 + rule.d_norm = 1 + rule._calculate_pnorm_dnorm = mock_pnorm_dnorm + rule.delta = 1 + + # Capture logs + with self.assertLogs(level='DEBUG') as log: + for i in range(15): + rule.get_step_size(pdhg) + + # Combine log messages + log_output = "\n".join(log.output) + + # Assertions + self.assertIn('Finished backtracking step', log_output) + self.assertIn("No change", log_output) + self.assertAlmostEqual(pdhg.sigma, 1) + self.assertAlmostEqual(pdhg.tau, 1) + self.assertAlmostEqual(rule.alpha, 0.95) + self.assertEqual(rule.count, 11) + self.assertEqual(rule.adaptive, False) + self.assertEqual(mock_func.call_count, 11) + with self.assertRaises(AttributeError): + rule.x_resid + with self.assertRaises(AttributeError): + rule.y_resid + with self.assertRaises(AttributeError): + rule.y_old + + +class TestPDHGBayesOpt(CCPiTestClass): + + def test_init(self): + rule = PDHGBayesOptimisationStepSize( + gamma_bounds=[0.1, 2], n_initial_points=5, n_calls=5, n_iterations=10, seed=42) + self.assertEqual(rule.gamma_bounds, [0.1, 2]) + self.assertEqual(rule.n_initial_points, 5) + self.assertEqual(rule.n_calls, 5) + self.assertEqual(rule.n_iterations, 10) + self.assertEqual(rule.seed, 42) + + def test_init_default(self): + rule = PDHGBayesOptimisationStepSize() + self.assertEqual(rule.n_initial_points, 5) + self.assertEqual(rule.n_calls, 20) + self.assertEqual(rule.n_iterations, 10) + self.assertEqual(rule.seed, None) + + def test_init_invalid(self): + with self.assertRaises(ValueError): + PDHGBayesOptimisationStepSize(gamma_bounds=[ + (0.1, 10)], n_initial_points=5) + + with self.assertRaises(ValueError): + PDHGBayesOptimisationStepSize(gamma_bounds=[ + -0, 1.1], n_calls=-5) + + @unittest.skipUnless(has_skopt, "scikit-optimize (skopt) not installed") + def test_get_gamma_bounded(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGBayesOptimisationStepSize( + gamma_bounds=[5, 100], n_initial_points=5, n_calls=5, n_iterations=10, seed=42) + pdhg = PDHG(f=f, g=g, operator=operator, + step_size=rule, update_objective_interval=4) + + gamma = np.sqrt(pdhg.sigma / pdhg.tau) + self.assertTrue(5 <= gamma <= 5.5) + + self.assertEqual(pdhg.iteration, -1) + self.assertEqual(pdhg.update_objective_interval, 4) + self.assertEqual(pdhg.objective, []) + + @unittest.skipUnless(has_skopt, "scikit-optimize (skopt) not installed") + def test_get_gamma(self): + ig = ImageGeometry(3, 3) + data = ig.allocate('random', seed=3) + + f = L2NormSquared(b=data) + g = L2NormSquared() + operator = IdentityOperator(ig) + + rule = PDHGBayesOptimisationStepSize( + gamma_bounds=None, n_initial_points=5, n_calls=10, n_iterations=10, seed=42) + pdhg = PDHG(f=f, g=g, operator=operator, step_size=rule) + gamma = np.sqrt(pdhg.sigma / pdhg.tau) + self.assertAlmostEqual(gamma, 1.7, places=1) + +class TestSPDHGConstantStepSize(CCPiTestClass): + def setUp(self): + self.subsets = 10 + + data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get(size=(16, 16)) + + partitioned_data = data.partition(self.subsets, 'sequential') + self.A = BlockOperator( + *[IdentityOperator(partitioned_data[i].geometry) for i in range(self.subsets)]) + + + # block function + self.F = BlockFunction(*[L2NormSquared(b=partitioned_data[i]) + for i in range(self.subsets)]) + alpha = 0.025 + self.G = alpha * IndicatorBox(lower=0) + def test_init_and_constant_step_size(self): + gamma = 1. + rho = .99 + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A) + self.assertListEqual( + spdhg.sigma, [rho / ni for ni in spdhg._norms]) + self.assertEqual(spdhg.tau, min([rho*pi / (si * ni**2) for pi, ni, + si in zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) + self.assertNumpyArrayEqual( + spdhg.x.as_array(), self.A.domain_geometry().allocate(0).as_array()) + self.assertEqual(spdhg.update_objective_interval, 1) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=(100, [1]*self.subsets)) + self.assertListEqual(spdhg.sigma, [1]*self.subsets) + self.assertEqual(spdhg.tau, 100) + + # Test SPDHG setters - set_step_sizes with sigma + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=(None, [1]*self.subsets)) + self.assertListEqual(spdhg.sigma, [1]*self.subsets) + self.assertEqual(spdhg.tau, min([(rho*pi / (si * ni**2)) for pi, ni, + si in zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) + + # Test SPDHG setters - set_step_sizes with tau + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=(100, None)) + self.assertListEqual(spdhg.sigma, [ + gamma * rho*pi / (spdhg.tau*ni**2) for ni, pi in zip(spdhg._norms, spdhg._prob_weights)]) + self.assertEqual(spdhg.tau, 100) + + def test_init_and_constant_step_size_invalid(self): + with self.assertRaises(ValueError): + SPDHGConstantStepSize(step_size=[1.0]) + with self.assertRaises(ValueError): + SPDHGConstantStepSize(step_size=[1.0, 2.0]).get_initial_step_size(SPDHG(f=self.F, g=self.G, operator=self.A)) + + def test_step_size_and_deprecated_sigma_tau_raises_valueerror(self): + # Passing both `step_size` and the deprecated `sigma`/`tau` must raise a + # clear ValueError. + with self.assertRaises(ValueError): + SPDHG(f=self.F, g=self.G, operator=self.A, + step_size=(100, [1]*self.subsets), sigma=0.5) + with self.assertRaises(ValueError): + SPDHG(f=self.F, g=self.G, operator=self.A, + step_size=(100, [1]*self.subsets), tau=0.5) + + def test_incompatible_gd_rule_raises_valueerror(self): + # A GD-style step-size rule has no get_initial_step_size and is not + # compatible with SPDHG; it must raise a clear ValueError. + with self.assertRaises(ValueError): + SPDHG(f=self.F, g=self.G, operator=self.A, + step_size=ConstantStepSize(0.1)) + + def test_step_sizes_from_ratio(self): + gamma = 3.7 + rho = 5.6 + rule = SPDHGStepSizesFromRatio(gamma,rho) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + self.assertListEqual( + spdhg.sigma, [gamma * rho / ni for ni in spdhg._norms]) + self.assertEqual(spdhg.tau, min([pi*rho / (si * ni**2) for pi, ni, + si in zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) + + +class TestSPDHGBayesStepSize(CCPiTestClass): + def setUp(self): + self.subsets = 4 + + data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get(size=(4, 4)) + + partitioned_data = data.partition(self.subsets, 'sequential') + self.A = BlockOperator( + *[IdentityOperator(partitioned_data[i].geometry) for i in range(self.subsets)]) + + + # block function + self.F = BlockFunction(*[L2NormSquared(b=partitioned_data[i]) + for i in range(self.subsets)]) + alpha = 0.025 + self.G = alpha * IndicatorBox(lower=0) - + def test_init(self): + rule = SPDHGBayesOptimisationStepSize( + gamma_bounds=[0.1, 2], n_initial_points=5, n_calls=10, n_iterations=10, seed=42) + self.assertEqual(rule.gamma_bounds, [0.1, 2]) + self.assertEqual(rule.n_initial_points, 5) + self.assertEqual(rule.n_calls, 10) + self.assertEqual(rule.n_iterations, 10) + self.assertEqual(rule.seed, 42) + + def test_init_default(self): + rule = SPDHGBayesOptimisationStepSize() + self.assertEqual(rule.n_initial_points, 5) + self.assertEqual(rule.n_calls, 20) + self.assertEqual(rule.n_iterations, None) + self.assertEqual(rule.seed, None) + + def test_init_invalid(self): + with self.assertRaises(ValueError): + SPDHGBayesOptimisationStepSize(gamma_bounds=[ + (0.1, 10)], n_initial_points=5) + + with self.assertRaises(ValueError): + SPDHGBayesOptimisationStepSize(gamma_bounds=[ + -0, 1.1], n_calls=-5) + + @unittest.skipUnless(has_skopt, "scikit-optimize (skopt) not installed") + def test_get_gamma_bounded(self): + + rule = SPDHGBayesOptimisationStepSize( + gamma_bounds=[5, 6], n_initial_points=5, n_calls=10, n_iterations=None, seed=42) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, + step_size=rule, update_objective_interval=4) + self.assertEqual(rule.n_iterations, 10*self.subsets) + gamma = min([ pi/(spdhg.tau *ni) for pi, ni in zip(spdhg._prob_weights, spdhg._norms)]) + self.assertTrue(5.8 <= gamma <= 6) + + self.assertEqual(spdhg.iteration, -1) + self.assertEqual(spdhg.update_objective_interval, 4) + self.assertEqual(spdhg.objective, []) + + @unittest.skipUnless(has_skopt, "scikit-optimize (skopt) not installed") + def test_get_gamma(self): + rule = SPDHGBayesOptimisationStepSize( + gamma_bounds=None, n_initial_points=5, n_calls=10, seed=42) + sampler = Sampler.sequential(self.subsets) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule, sampler=sampler) + gamma = min([ pi/(spdhg.tau *ni) for pi, ni in zip(spdhg._prob_weights, spdhg._norms)]) + self.assertAlmostEqual(gamma, 0.4, places=1) + +class TestSPDHGAdaptiveStepSize(CCPiTestClass): + """Tests for the adaptive SPDHG step-size rules based on arXiv:2301.02511.""" + + def setUp(self): + self.subsets = 10 + + data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get(size=(16, 16)) + + partitioned_data = data.partition(self.subsets, 'sequential') + self.A = BlockOperator( + *[IdentityOperator(partitioned_data[i].geometry) for i in range(self.subsets)]) + + self.F = BlockFunction(*[L2NormSquared(b=partitioned_data[i]) + for i in range(self.subsets)]) + alpha = 0.025 + self.G = alpha * IndicatorBox(lower=0) + + + def test_balancing_init(self): + rule = SPDHGAdaptiveStepSizeBalancing( + initial_step_size=[1.0, 2.0], initial_alpha=0.9, eta=0.99, + delta=2.0, s=3.0, auto_stop=False, auto_stop_patience=5) + self.assertEqual(rule.initial_step_size, [1.0, 2.0]) + self.assertEqual(rule.alpha, 0.9) + self.assertEqual(rule.eta, 0.99) + self.assertEqual(rule.delta, 2.0) + self.assertEqual(rule.s, 3.0) + self.assertFalse(rule.auto_stop) + self.assertEqual(rule.auto_stop_patience, 5) + self.assertTrue(rule.adaptive) + self.assertEqual(rule.count, 0) + + def test_balancing_init_defaults(self): + rule = SPDHGAdaptiveStepSizeBalancing() + self.assertEqual(rule.initial_step_size, [None, None]) + self.assertAlmostEqual(rule.alpha, 0.95) + self.assertAlmostEqual(rule.eta, 0.995) + self.assertEqual(rule.delta, 1.5) + self.assertIsNone(rule.s) + self.assertTrue(rule.auto_stop) + self.assertEqual(rule.auto_stop_patience, 10) + + def test_angle_init(self): + rule = SPDHGAdaptiveStepSizeAngle( + initial_step_size=[1.0, 2.0], initial_alpha=0.5, eta=0.9, + c=0.9, auto_stop=False, auto_stop_patience=3) + self.assertEqual(rule.initial_step_size, [1.0, 2.0]) + self.assertEqual(rule.alpha, 0.5) + self.assertEqual(rule.eta, 0.9) + self.assertEqual(rule.c, 0.9) + self.assertFalse(rule.auto_stop) + self.assertEqual(rule.auto_stop_patience, 3) + + def test_angle_init_defaults(self): + rule = SPDHGAdaptiveStepSizeAngle() + self.assertEqual(rule.initial_step_size, [None, None]) + self.assertAlmostEqual(rule.alpha, 1.0) + self.assertAlmostEqual(rule.eta, 0.995) + self.assertEqual(rule.c, 0.999) + self.assertTrue(rule.auto_stop) + + def test_init_invalid(self): + with self.assertRaises(ValueError): + SPDHGAdaptiveStepSizeBalancing(initial_step_size=[1.0]) + with self.assertRaises(ValueError): + SPDHGAdaptiveStepSizeAngle(initial_step_size=[1.0, 2.0, 3.0]) + + + def test_initial_step_size_defaults(self): + rho = 0.99 + for RuleCls in (SPDHGAdaptiveStepSizeBalancing, SPDHGAdaptiveStepSizeAngle): + rule = RuleCls() + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + # tau scalar, sigma list of length n, all positive + self.assertTrue(isinstance(spdhg.tau, Number)) + self.assertEqual(len(spdhg.sigma), self.subsets) + self.assertTrue(all(s > 0 for s in spdhg.sigma)) + # matches the standard SPDHG relations + self.assertListEqual(spdhg.sigma, [rho / ni for ni in spdhg._norms]) + self.assertEqual(spdhg.tau, min([rho*pi / (si * ni**2) for pi, ni, si in + zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) + # per-operator convergence criterion sigma_i * tau * ||A_i||^2 <= p_i + for si, ni, pi in zip(spdhg.sigma, spdhg._norms, spdhg._prob_weights): + self.assertLessEqual(si * spdhg.tau * ni**2, pi + 1e-9) + + def test_initial_step_size_scalar_sigma_broadcast(self): + # a single dual step size is broadcast to a per-operator list + rule = SPDHGAdaptiveStepSizeBalancing(initial_step_size=[3.0, 0.5]) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + self.assertEqual(spdhg.tau, 3.0) + self.assertListEqual(spdhg.sigma, [0.5]*self.subsets) + + def test_initial_step_size_tau_from_sigma_list(self): + # sigma given as a list, tau derived + rho = 0.99 + rule = SPDHGAdaptiveStepSizeAngle(initial_step_size=[None, [0.7]*self.subsets]) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + self.assertListEqual(spdhg.sigma, [0.7]*self.subsets) + self.assertEqual(spdhg.tau, min([rho*pi / (si * ni**2) for pi, ni, si in + zip(spdhg._prob_weights, spdhg._norms, spdhg.sigma)])) + + + + def test_balancing_adapts_step_sizes(self): + # from a deliberately imbalanced ratio, the rule should rescale the step sizes + rule = SPDHGAdaptiveStepSizeBalancing( + initial_step_size=[10.0, 0.001], auto_stop=False) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + start_ratio = spdhg.tau / spdhg.sigma[0] + spdhg.run(20) + end_ratio = spdhg.tau / spdhg.sigma[0] + self.assertNotAlmostEqual(start_ratio, end_ratio) + # step sizes remain valid (positive scalar tau, positive sigma list) + self.assertGreater(spdhg.tau, 0) + self.assertEqual(len(spdhg.sigma), self.subsets) + self.assertTrue(all(s > 0 for s in spdhg.sigma)) + # the adaptation strength decayed as the step sizes were rebalanced + self.assertLess(rule.alpha, 0.95) + + def test_angle_adapts_step_sizes(self): + # a permissive threshold (c=0) makes every iteration rebalance, exercising the + # angle-alignment mechanism (the default c=0.999 is deliberately conservative). + rule = SPDHGAdaptiveStepSizeAngle( + initial_step_size=[10.0, 0.001], c=0.0, auto_stop=False) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + start_ratio = spdhg.tau / spdhg.sigma[0] + spdhg.run(20) + end_ratio = spdhg.tau / spdhg.sigma[0] + self.assertNotAlmostEqual(start_ratio, end_ratio) + self.assertGreater(spdhg.tau, 0) + self.assertEqual(len(spdhg.sigma), self.subsets) + self.assertTrue(all(s > 0 for s in spdhg.sigma)) + + def test_products_preserved_under_rescaling(self): + # rescaling multiplies tau by 1/f and every sigma_i by f, so tau*sigma_i is fixed + for RuleCls in (SPDHGAdaptiveStepSizeBalancing, SPDHGAdaptiveStepSizeAngle): + rule = RuleCls(initial_step_size=[10.0, 0.001], auto_stop=False) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + prod0 = [spdhg.tau * si for si in spdhg.sigma] + spdhg.run(15) + prod1 = [spdhg.tau * si for si in spdhg.sigma] + for a, b in zip(prod0, prod1): + self.assertAlmostEqual(a, b) + + def test_auto_stop_freezes_and_frees(self): + # force the "no change" branch every iteration by making the tolerance gate + # impossible to pass; the rule should then stop after auto_stop_patience. + rule = SPDHGAdaptiveStepSizeBalancing(auto_stop=True, auto_stop_patience=3) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + rule.tolerance = 1e12 + self.assertTrue(rule.adaptive) + spdhg.run(6) + self.assertFalse(rule.adaptive) + # buffers released on auto-stop + self.assertFalse(hasattr(rule, 'x_prev') and rule.x_prev is not None) + # continuing to run after auto-stop keeps returning valid step sizes + tau_frozen = spdhg.tau + sigma_frozen = list(spdhg.sigma) + spdhg.run(3) + self.assertEqual(spdhg.tau, tau_frozen) + self.assertListEqual(list(spdhg.sigma), sigma_frozen) + + def test_auto_stop_disabled_keeps_adapting(self): + rule = SPDHGAdaptiveStepSizeAngle(auto_stop=False, auto_stop_patience=2) + spdhg = SPDHG(f=self.F, g=self.G, operator=self.A, step_size=rule) + spdhg.run(20) + self.assertTrue(rule.adaptive) + self.assertTrue(hasattr(rule, 'x_prev')) + + diff --git a/Wrappers/Python/test/utils.py b/Wrappers/Python/test/utils.py index 0f94d4ad96..d4abde4d1b 100644 --- a/Wrappers/Python/test/utils.py +++ b/Wrappers/Python/test/utils.py @@ -90,6 +90,11 @@ def initialise_tests(): has_cvxpy = True system_state['has_cvxpy']=has_cvxpy +#skopt (scikit-optimize) +module_info = importlib.util.find_spec("skopt") +has_skopt = module_info is not None +system_state['has_skopt'] = has_skopt + #ipp from cil.framework import cilacc diff --git a/docs/source/optimisation.rst b/docs/source/optimisation.rst index 0ec5e4f262..4f60ddbc3d 100644 --- a/docs/source/optimisation.rst +++ b/docs/source/optimisation.rst @@ -149,10 +149,14 @@ Implemented examples are: PDHG ---- .. autoclass:: cil.optimisation.algorithms.PDHG - :members: update, set_step_sizes, update_step_sizes, update_objective + :members: update, update_objective, check_convergence :member-order: bysource :inherited-members: run, update_objective_interval +The primal and dual step sizes of PDHG are controlled by a step-size rule passed via the +``step_size`` argument. See :ref:`PDHG step-size rules` for the available rules, including +adaptive and Bayesian-optimisation options. + LADMM ----- .. autoclass:: cil.optimisation.algorithms.LADMM @@ -192,9 +196,13 @@ Each iteration considers just one index of the sum, potentially reducing computa .. autoclass:: cil.optimisation.algorithms.SPDHG - :members: update, set_step_sizes, set_step_sizes_from_ratio, update_objective + :members: update, update_objective, check_convergence :inherited-members: run, update_objective_interval +The primal step size and the per-operator dual step sizes of SPDHG are controlled by a +step-size rule passed via the ``step_size`` argument. See :ref:`SPDHG step-size rules` for the +available rules. + Approximate gradient methods ---------------------------------- @@ -708,16 +716,18 @@ In each iteration of the :code:`TestAlgo`, the objective :math:`x` is reduced by 15%|███ | 3/20 [00:00<00:00, 11770.73it/s, objective=3.05e-5] -Step size methods +Step size methods ------------------ -A step size method is a class which acts on an algorithm and can be passed to `cil.optimisation.algorithm.GD`, `cil.optimisation.algorithm.ISTA` `cil.optimisation.algorithm.FISTA` and it's method `get_step_size` is called after the calculation of the gradient before the gradient descent step is taken. It outputs a float value to be used as the step-size. +A step size method is a class which acts on an algorithm and can be passed to a selection of the CIL algorithms. Currently in CIL we have a base class: .. autoclass:: cil.optimisation.utilities.StepSizeMethods.StepSizeRule :members: -We also have a number of example classes: +Gradient-based step-size rules +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +For the gradient based algorithms `cil.optimisation.algorithms.GD`, `cil.optimisation.algorithms.ISTA` and `cil.optimisation.algorithms.FISTA`, the method `get_step_size` is called after the calculation of the gradient, before the gradient descent step is taken. It outputs a single float value to be used as the step-size. .. autoclass:: cil.optimisation.utilities.StepSizeMethods.ConstantStepSize :members: @@ -728,6 +738,48 @@ We also have a number of example classes: .. autoclass:: cil.optimisation.utilities.StepSizeMethods.BarzilaiBorweinStepSizeRule :members: +.. _PDHG step-size rules: + +PDHG step-size rules +~~~~~~~~~~~~~~~~~~~~~~ +For `cil.optimisation.algorithms.PDHG`, the rule's `get_initial_step_size` method sets the initial primal and dual step sizes during set-up, and `get_step_size` is called at the end of each iteration to output the primal and dual step sizes for the next iteration. The rules below span three strategies: a fixed step size (`PDHGConstantStepSize`); acceleration exploiting strong convexity (`PDHGStronglyConvexUpdate`); per-iteration backtracking with residual balancing (`PDHGAdaptiveStepSize2013`, `PDHGAdaptiveStepSize2015`); and a one-off Bayesian-optimisation search for the best step-size ratio (`PDHGBayesOptimisationStepSize`). The two adaptive rules are faster but memory-hungry, while the Bayesian-optimisation rule trades compute time for lower memory use. + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.PDHGConstantStepSize + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.PDHGStronglyConvexUpdate + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.PDHGAdaptiveStepSize2013 + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.PDHGAdaptiveStepSize2015 + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.PDHGBayesOptimisationStepSize + :members: + +.. _SPDHG step-size rules: + +SPDHG step-size rules +~~~~~~~~~~~~~~~~~~~~~~~ +For `cil.optimisation.algorithms.SPDHG`, the rule's `get_initial_step_size` method sets the initial step sizes during set-up and `get_step_size` is called at the end of each iteration. Here the primal step size ``tau`` is a scalar while the dual step size ``sigma`` is a list with one entry per operator. The rules below span a fixed step size (`SPDHGConstantStepSize`), a fixed primal/dual ratio (`SPDHGStepSizesFromRatio`), a one-off Bayesian-optimisation search for the best ratio (`SPDHGBayesOptimisationStepSize`), and two per-iteration adaptive rules from :cite:`chambolle2023stochastic` that rescale the step sizes as the algorithm runs — one balancing the primal and dual progress residuals (`SPDHGAdaptiveStepSizeBalancing`) and one exploiting the alignment of successive primal directions (`SPDHGAdaptiveStepSizeAngle`). + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.SPDHGConstantStepSize + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.SPDHGStepSizesFromRatio + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.SPDHGBayesOptimisationStepSize + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.SPDHGAdaptiveStepSizeBalancing + :members: + +.. autoclass:: cil.optimisation.utilities.StepSizeMethods.SPDHGAdaptiveStepSizeAngle + :members: + Preconditioners diff --git a/docs/source/refs.bib b/docs/source/refs.bib index 4ca6c4e543..5db21815bd 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -127,3 +127,26 @@ @article{CombettesValerie doi={10.1137/050626090}, URL={https://doi.org/10.1137/050626090}, eprint={https://doi.org/10.1137/050626090}} +@inproceedings{Goldstein2015, + author = {Goldstein, Tom and Li, Min and Yuan, Xiaoming}, + booktitle = {Advances in Neural Information Processing Systems}, + editor = {C. Cortes and N. Lawrence and D. Lee and M. Sugiyama and R. Garnett}, + pages = {}, + publisher = {Curran Associates, Inc.}, + title = {Adaptive Primal-Dual Splitting Methods for Statistical Learning and Image Processing}, + url = {https://proceedings.neurips.cc/paper_files/paper/2015/file/cd758e8f59dfdf06a852adad277986ca-Paper.pdf}, + volume = {28}, + year = {2015} +} +@article{goldstein2013adaptive, + title={Adaptive primal-dual hybrid gradient methods for saddle-point problems}, + author={Goldstein, Tom and Li, Min and Yuan, Xiaoming and Esser, Ernie and Baraniuk, Richard}, + journal={arXiv preprint arXiv:1305.0546}, + year={2013} +} +@article{chambolle2023stochastic, + title={Stochastic Primal-Dual Hybrid Gradient Algorithm with Adaptive Step-Sizes}, + author={Chambolle, Antonin and Delplancke, Claire and Ehrhardt, Matthias J. and Sch{\"o}nlieb, Carola-Bibiane and Tang, Junqi}, + journal={arXiv preprint arXiv:2301.02511}, + year={2023} +} \ No newline at end of file