Optimisation framework#
This package allows rapid prototyping of optimisationbased reconstruction problems, i.e. defining and solving different optimization problems to enforce different properties on the reconstructed image.
Firstly, it provides an objectoriented framework for defining mathematical operators and functions as well a collection of useful example operators and functions. Both smooth and nonsmooth functions can be used.
Further, it provides a number of highlevel generic implementations of optimisation algorithms to solve generically formulated optimisation problems constructed from operator and function objects.
The fundamental components are:
Operator
: A class specifying a (currently linear) operator.Function
: A class specifying mathematical functions such as a least squares data fidelity.Algorithm
: Implementation of an iterative optimisation algorithm to solve a particular generic optimisation problem. Algorithms are iterable Python object which can be run in a for loop. Can be stopped and warm restarted.
Algorithms (Deterministic)#
A number of generic algorithm implementations are provided including Gradient Descent (GD), Conjugate Gradient Least Squares (CGLS), Simultaneous Iterative Reconstruction Technique (SIRT), Primal Dual Hybrid Gradient (PDHG), Iterative Shrinkage Thresholding Algorithm (ISTA), and Fast Iterative Shrinkage Thresholding Algorithm (FISTA).
An algorithm is designed for a particular generic optimisation problem accepts and number of
instances of Function
derived classes and/or Operator
derived classes as input to
define a specific instance of the generic optimisation problem to be solved.
They are iterable objects which can be run in a for loop.
The user can provide a stopping criterion different than the default max_iteration.
New algorithms can be easily created by extending the Algorithm
class.
The user is required to implement only 4 methods: set_up, __init__, update and update_objective.
set_up
and__init__
are used to configure the algorithmupdate
is the actual iteration updating the solutionupdate_objective
defines how the objective is calculated.
For example, the implementation of the update of the Gradient Descent algorithm to minimise a Function will only be:
def update(self):
self.x += self.rate * self.objective_function.gradient(self.x)
def update_objective(self):
self.loss.append(self.objective_function(self.x))
The Algorithm
provides the infrastructure to continue iteration, to access the values of the
objective function in subsequent iterations, the time for each iteration, and to provide a nice
print to screen of the status of the optimisation.
Base class#
 class cil.optimisation.algorithms.Algorithm(update_objective_interval=1, max_iteration=None, log_file=None)[source]#
Base class for iterative algorithms
provides the minimal infrastructure.
Algorithms are iterables so can be easily run in a for loop. They will stop as soon as the stop criterion is met. The user is required to implement the
set_up
,__init__
,update
and andupdate_objective
methodsA courtesy method
run
is available to runn
iterations. The method accepts acallbacks
list of callables, each of which receive the current Algorithm object (which in turn contains the iteration number and the actual objective value) and can be used to trigger print to screens and other user interactions. Therun
method will stop when the stopping criterion is met or StopIteration is raised. should_stop()[source]#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
 is_provably_convergent()[source]#
Check if the algorithm is convergent based on the provable convergence criterion.
 get_last_loss(return_all=False)[source]#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)[source]#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
GD#
 class cil.optimisation.algorithms.GD(initial=None, objective_function=None, step_size=None, rtol=1e05, atol=1e08, preconditioner=None, **kwargs)[source]#
Gradient Descent algorithm
 Parameters:
initial (DataContainer (e.g. ImageData)) – The initial point for the optimisation
objective_function (CIL function (
Function()
. ) with a defined gradient method) – The function to be minimised.step_size (positive real float or subclass of
StepSizeRule()
, default = None) – If you pass a float this will be used as a constant step size. If left as None and do not pass a step_size_rule then the Armijio rule will be used to perform backtracking to choose a step size at each iteration. If a child class ofcil.optimisation.utilities.StepSizeRule()
’ is passed then it’s method get_step_size is called for each update.preconditioner (class with a apply method or a function that takes an initialised CIL function as an argument and modifies a provided gradient.) – This could be a custom preconditioner or one provided in
preconditioner()
. If None is passed then self.gradient_update will remain unmodified.rtol (positive float, default 1e5) – optional parameter defining the relative tolerance comparing the current objective function to 0, default 1e5, see numpy.isclose
atol (positive float, default 1e8) – optional parameter defining the absolute tolerance comparing the current objective function to 0, default 1e8, see numpy.isclose
 set_up(initial, objective_function, step_size, preconditioner)[source]#
initialisation of the algorithm
 Parameters:
initial (DataContainer (e.g. ImageData)) – The initial point for the optimisation
objective_function (CIL function with a defined gradient) – The function to be minimised.
step_size (positive real float or subclass of
StepSizeRule()
, default = None) – If you pass a float this will be used as a constant step size. If left as None and do not pass a step_size_rule then the Armijio rule will be used to perform backtracking to choose a step size at each iteration. If a child class ofcil.optimisation.utilities.StepSizeRule()
’ is passed then it’s method get_step_size is called for each update.preconditioner (class with a apply method or a function that takes an initialised CIL function as an argument and modifies a provided gradient.) – This could be a custom preconditioner or one provided in
preconditioner()
. If None is passed then self.gradient_update will remain unmodified.
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_output()#
Returns the current solution.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
CGLS#
 class cil.optimisation.algorithms.CGLS(initial=None, operator=None, data=None, tolerance=1e06, **kwargs)[source]#
Conjugate Gradient Least Squares algorithm
Problem:
\[\min  A x  b ^2_2\]Parameters :
:parameter operator : Linear operator for the inverse problem :parameter initial : Initial guess ( Default initial = 0) :parameter data : Acquired data to reconstruct :parameter tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm
 set_up(initial, operator, data, tolerance=1e06)[source]#
initialisation of the algorithm
 Parameters:
operator – Linear operator for the inverse problem
initial – Initial guess ( Default initial = 0)
data – Acquired data to reconstruct
tolerance – Tolerance/ Stopping Criterion to end CGLS algorithm
 should_stop()[source]#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_output()#
Returns the current solution.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
SIRT#
 class cil.optimisation.algorithms.SIRT(initial, operator, data, lower=None, upper=None, constraint=None, **kwargs)[source]#
Simultaneous Iterative Reconstruction Technique, see [8].
Simultaneous Iterative Reconstruction Technique (SIRT) solves the following problem
\[A x = b\]The SIRT algorithm is
\[x^{k+1} = \mathrm{proj}_{C}( x^{k} + \omega * D ( A^{T} ( M * (b  Ax^{k}) ) ) ),\]where, \(M = \frac{1}{A*\mathbb{1}}\), \(D = \frac{1}{A^{T}\mathbb{1}}\), \(\mathbb{1}\) is a
DataContainer
of ones, \(\mathrm{prox}_{C}\) is the projection over a set \(C\), and \(\omega\) is the relaxation parameter. Parameters:
initial (DataContainer, default = None) – Starting point of the algorithm, default value = Zero DataContainer
operator (LinearOperator) – The operator A.
data (DataContainer) – The data b.
lower (
float
, default = None) – Lower bound constraintupper (
float
, default = None) – Upper bound constraintconstraint (Function, default = None) – A function with
proximal
method, e.g.,IndicatorBox
function andIndicatorBox.proximal()
, orTotalVariation
function andTotalVariation.proximal()
.kwargs – Keyword arguments used from the base class
Algorithm
.
Note
If
constraint
is not passed,lower
andupper
are used to create anIndicatorBox
and apply itsproximal
.If
constraint
is passed,proximal
method is required to be implemented.Note
The preconditioning arrays (weights)
M
andD
used in SIRT are defined as\[M = \frac{1}{A*\mathbb{1}} = \frac{1}{\sum_{j}a_{i,j}}\]\[D = \frac{1}{A*\mathbb{1}} = \frac{1}{\sum_{i}a_{i,j}}\]Examples
\[\underset{x}{\mathrm{argmin}} \frac{1}{2}\ x  d\^{2}\]>>> sirt = SIRT(initial = ig.allocate(0), operator = A, data = d, max_iteration = 5)
 set_up(initial, operator, data, lower=None, upper=None, constraint=None)[source]#
Initialisation of the algorithm
 set_relaxation_parameter(value=1.0)[source]#
Set the relaxation parameter \(\omega\)
 Parameters:
value (float) – The relaxation parameter to be applied to the update. Must be between 0 and 2 to guarantee asymptotic convergence.
 update()[source]#
Performs a single iteration of the SIRT algorithm
\[x^{k+1} = \mathrm{proj}_{C}( x^{k} + \omega * D ( A^{T} ( M * (b  Ax) ) ) )\]
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_output()#
Returns the current solution.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
 should_stop()#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
ISTA#
 class cil.optimisation.algorithms.ISTA(initial, f, g, step_size=None, preconditioner=None, **kwargs)[source]#
Iterative ShrinkageThresholding Algorithm, see [1], [2].
Iterative ShrinkageThresholding Algorithm (ISTA)
\[x^{k+1} = \mathrm{prox}_{\alpha^{k} g}(x^{k}  \alpha^{k}\nabla f(x^{k}))\]is used to solve
\[\min_{x} f(x) + g(x)\]where \(f\) is differentiable, \(g\) has a simple proximal operator and \(\alpha^{k}\) is the
step_size
per iteration.Note
For a constant step size, i.e., \(a^{k}=a\) for \(k\geq1\), convergence of ISTA is guaranteed if
\[\alpha\in(0, \frac{2}{L}),\]where \(L\) is the Lipschitz constant of \(f\), see [4].
 Parameters:
initial (DataContainer) – Initial guess of ISTA.
f (Function) – Differentiable function. If None is passed, the algorithm will use the ZeroFunction.
g (Function or None) – Convex function with simple proximal operator. If None is passed, the algorithm will use the ZeroFunction.
step_size (positive
float
or child class ofcil.optimisation.utilities.StepSizeRule()
’, default = None) – Step size for the gradient step of ISTA. If a float is passed, this is used as a constant step size. If a child class of
cil.optimisation.utilities.StepSizeRule()
’ is passed then it’s method get_step_size is called for each update. The default
step_size
is a constant \(\frac{0.99*2}{L}\) or 1 if f=None. preconditioner: class with a apply method or a function that takes an initialised CIL function as an argument and modifies a provided gradient.
This could be a custom preconditioner or one provided in
preconditoner()
. If None is passed then self.gradient_update will remain unmodified.
 Step size for the gradient step of ISTA. If a float is passed, this is used as a constant step size. If a child class of
 kwargs: Keyword arguments
Arguments from the base class
Algorithm
.
Note
If the function g is set to None or to the ZeroFunction then the ISTA algorithm is equivalent to Gradient Descent.
If the function f is set to None or to the ZeroFunction then the ISTA algorithm is equivalent to a Proximal Point Algorithm.
Examples
\[\underset{x}{\mathrm{argmin}}\A x  b\^{2}_{2}\]>>> f = LeastSquares(A, b=b, c=0.5) >>> g = ZeroFunction() >>> ig = Aop.domain >>> ista = ISTA(initial = ig.allocate(), f = f, g = g, max_iteration=10) >>> ista.run()
 __init__(initial, f, g, step_size=None, preconditioner=None, **kwargs)[source]#
Set the minimal number of parameters:
 Parameters:
update_objective_interval (int, optional, default 1) – the interval every which we would save the current objective. 1 means every iteration, 2 every 2 iteration and so forth. This is by default 1 and should be increased when evaluating the objective is computationally expensive.
 update()[source]#
Performs a single iteration of ISTA
\[x_{k+1} = \mathrm{prox}_{\alpha g}(x_{k}  \alpha\nabla f(x_{k}))\]
 __delattr__(name, /)#
Implement delattr(self, name).
 __dir__()#
Default dir() implementation.
 __eq__(value, /)#
Return self==value.
 __format__(format_spec, /)#
Default object formatter.
Return str(self) if format_spec is empty. Raise TypeError otherwise.
 __ge__(value, /)#
Return self>=value.
 __getattribute__(name, /)#
Return getattr(self, name).
 __getstate__()#
Helper for pickle.
 __gt__(value, /)#
Return self>value.
 __hash__()#
Return hash(self).
 __init_subclass__()#
This method is called when a class is subclassed.
The default implementation does nothing. It may be overridden to extend subclasses.
 __iter__()#
Algorithm is an iterable
 __le__(value, /)#
Return self<=value.
 __lt__(value, /)#
Return self<value.
 __ne__(value, /)#
Return self!=value.
 __new__(**kwargs)#
 __next__()#
Algorithm is an iterable
This method triggers
update()
andupdate_objective()
 __reduce__()#
Helper for pickle.
 __reduce_ex__(protocol, /)#
Helper for pickle.
 __repr__()#
Return repr(self).
 __setattr__(name, value, /)#
Implement setattr(self, name, value).
 __sizeof__()#
Size of object in memory, in bytes.
 __str__()#
Return str(self).
 __subclasshook__()#
Abstract classes can override this to customize issubclass().
This is invoked early on by abc.ABCMeta.__subclasscheck__(). It should return True, False or NotImplemented. If it returns NotImplemented, the normal algorithm is used. Otherwise, it overrides the normal algorithm (and the outcome is cached).
 __weakref__#
list of weak references to the object
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
 should_stop()#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
FISTA#
 class cil.optimisation.algorithms.FISTA(initial, f, g, step_size=None, preconditioner=None, **kwargs)[source]#
Fast Iterative ShrinkageThresholding Algorithm, see [1], [2].
Fast Iterative ShrinkageThresholding Algorithm (FISTA)
\[\begin{split}\begin{cases} y_{k} = x_{k}  \alpha\nabla f(x_{k}) \\ x_{k+1} = \mathrm{prox}_{\alpha g}(y_{k})\\ t_{k+1} = \frac{1+\sqrt{1+ 4t_{k}^{2}}}{2}\\ y_{k+1} = x_{k} + \frac{t_{k}1}{t_{k1}}(x_{k}  x_{k1}) \end{cases}\end{split}\]is used to solve
\[\min_{x} f(x) + g(x)\]where \(f\) is differentiable, \(g\) has a simple proximal operator and \(\alpha^{k}\) is the
step_size
per iteration. Parameters:
initial (DataContainer) – Starting point of the algorithm
f (Function) – Differentiable function. If None is passed, the algorithm will use the ZeroFunction.
g (Function or None) – Convex function with simple proximal operator. If None is passed, the algorithm will use the ZeroFunction.
step_size (positive
float
or child class ofcil.optimisation.utilities.StepSizeRule()
’, default = None) – Step size for the gradient step of ISTA. If a float is passed, this is used as a constant step size. If a child class ofcil.optimisation.utilities.StepSizeRule()
’ is passed then it’s method get_step_size is called for each update. The defaultstep_size
is a constant \(\frac{1}{L}\) or 1 if f=None.preconditioner (class with a apply method or a function that takes an initialised CIL function as an argument and modifies a provided gradient.) – This could be a custom preconditioner or one provided in
preconditoner()
. If None is passed then self.gradient_update will remain unmodified.kwargs (Keyword arguments) – Arguments from the base class
Algorithm
.
Note
If the function g is set to None or to the ZeroFunction then the FISTA algorithm is equivalent to Accelerated Gradient Descent by Nesterov ([9] algorithm 2.2.9).
If the function f is set to None or to the ZeroFunction then the FISTA algorithm is equivalent to Guler’s First Accelerated Proximal Point Method ([6] sec 2).
Examples
\[\underset{x}{\mathrm{argmin}}\A x  b\^{2}_{2}\]>>> f = LeastSquares(A, b=b, c=0.5) >>> g = ZeroFunction() >>> ig = Aop.domain >>> fista = FISTA(initial = ig.allocate(), f = f, g = g, max_iteration=10) >>> fista.run()
 __init__(initial, f, g, step_size=None, preconditioner=None, **kwargs)[source]#
Set the minimal number of parameters:
 Parameters:
update_objective_interval (int, optional, default 1) – the interval every which we would save the current objective. 1 means every iteration, 2 every 2 iteration and so forth. This is by default 1 and should be increased when evaluating the objective is computationally expensive.
 update()[source]#
Performs a single iteration of FISTA
\[\begin{split}\begin{cases} x_{k+1} = \mathrm{prox}_{\alpha g}(y_{k}  \alpha\nabla f(y_{k}))\\ t_{k+1} = \frac{1+\sqrt{1+ 4t_{k}^{2}}}{2}\\ y_{k+1} = x_{k} + \frac{t_{k}1}{t_{k1}}(x_{k}  x_{k1}) \end{cases}\end{split}\]
 __delattr__(name, /)#
Implement delattr(self, name).
 __dir__()#
Default dir() implementation.
 __eq__(value, /)#
Return self==value.
 __format__(format_spec, /)#
Default object formatter.
Return str(self) if format_spec is empty. Raise TypeError otherwise.
 __ge__(value, /)#
Return self>=value.
 __getattribute__(name, /)#
Return getattr(self, name).
 __getstate__()#
Helper for pickle.
 __gt__(value, /)#
Return self>value.
 __hash__()#
Return hash(self).
 __init_subclass__()#
This method is called when a class is subclassed.
The default implementation does nothing. It may be overridden to extend subclasses.
 __iter__()#
Algorithm is an iterable
 __le__(value, /)#
Return self<=value.
 __lt__(value, /)#
Return self<value.
 __ne__(value, /)#
Return self!=value.
 __new__(**kwargs)#
 __next__()#
Algorithm is an iterable
This method triggers
update()
andupdate_objective()
 __reduce__()#
Helper for pickle.
 __reduce_ex__(protocol, /)#
Helper for pickle.
 __repr__()#
Return repr(self).
 __setattr__(name, value, /)#
Implement setattr(self, name, value).
 __sizeof__()#
Size of object in memory, in bytes.
 __str__()#
Return str(self).
 __subclasshook__()#
Abstract classes can override this to customize issubclass().
This is invoked early on by abc.ABCMeta.__subclasscheck__(). It should return True, False or NotImplemented. If it returns NotImplemented, the normal algorithm is used. Otherwise, it overrides the normal algorithm (and the outcome is cached).
 __weakref__#
list of weak references to the object
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_output()#
Returns the current solution.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
 set_up(initial, f, g, step_size, preconditioner, **kwargs)#
Set up of the algorithm
 should_stop()#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
 update_objective()#
Updates the objective
\[f(x) + g(x)\]
PDHG#
 class cil.optimisation.algorithms.PDHG(f, g, operator, tau=None, sigma=None, initial=None, gamma_g=None, gamma_fconj=None, **kwargs)[source]#
Primal Dual Hybrid Gradient (PDHG) algorithm, see [3], [5].
 Parameters:
f (Function) – A convex function with a “simple” proximal method of its conjugate.
g (Function) – A convex function with a “simple” proximal.
operator (LinearOperator) – A Linear Operator.
sigma (positive
float
, or np.ndarray, DataContainer, BlockDataContainer, optional, default=None) – Step size for the dual problem.tau (positive
float
, or np.ndarray, DataContainer, BlockDataContainer, optional, default=None) – Step size for the primal problem.initial (DataContainer, optional, default=None) – Initial point for the PDHG algorithm.
gamma_g (positive
float
, optional, default=None) – Strongly convex constant if the function g is strongly convex. Allows primal acceleration of the PDHG algorithm.gamma_fconj (positive
float
, optional, default=None) – Strongly convex constant if the convex conjugate of f is strongly convex. Allows dual acceleration of the PDHG algorithm.**kwargs –
Keyward arguments used from the base class
Algorithm
. max_iteration
int
, optional, default=0 Maximum number of iterations of the PDHG.
 update_objective_interval
int
, optional, default=1 Evaluates objectives, e.g., primal/dual/primaldual gap every
update_objective_interval
. check_convergence
boolean
, default=True Checks scalar sigma and tau values satisfy convergence criterion
 max_iteration
Example
In our CILDemos repositoryyou can find examples using the PDHG algorithm for different imaging problems, such as Total Variation denoising, Total Generalised Variation inpaintingand Total Variation Tomography reconstruction. More examples can also be found in [7], [10].
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
A firstorder primaldual algorithm for convex optimization problems with known saddlepoint structure with applications in imaging.
The general problem considered in the PDHG algorithm is the generic saddlepoint problem
\[\min_{x\in X}\max_{y\in Y} \langle Kx, y \rangle + g(x)  f^{*}(x)\]where \(f\) and \(g\) are convex functions with “simple” proximal operators.
\(X\) and \(Y\) are two two finitedimensional vector spaces with an inner product and representing the domain of \(g\) and \(f^{*}\), the convex conjugate of \(f\), respectively.
The operator \(K\) is a continuous linear operator with operator norm defined as
\[\K\ = \max\{ \Kx\ : x\in X, \x\\leq1\}\]The saddle point problem is decomposed into the primal problem:
\[\min_{x\in X} f(Kx) + g(x),\]and its corresponding dual problem
\[\max_{y\in Y}  g^{*}(K^{*}y)  f^{*}(y).\]The PDHG algorithm consists of three steps:
gradient ascent step for the dual problem,
gradient descent step for the primal problem and
an overrelaxation of the primal variable.
\[y^{n+1} = \mathrm{prox}_{\sigma f^{*}}(y^{n} + \sigma K \bar{x}^{n})\]\[x^{n+1} = \mathrm{prox}_{\tau g}(x^{n}  \tau K^{*}y^{n+1})\]\[\bar{x}^{n+1} = x^{n+1} + \theta (x^{n+1}  x^{n})\]Notes
Convergence is guaranteed if \(\theta\) = 1.0, the operator norm \(\K\\), the dual step size \(\sigma\) and the primal step size \(\tau\), satisfy the following inequality:
\[\tau \sigma \K\^2 < 1\]By default, the step sizes \(\sigma\) and \(\tau\) are positive scalars and defined as below:
If
sigma
isNone
andtau
isNone
:
\[\sigma = \frac{1}{\K\}, \tau = \frac{1}{\K\}\]If
tau
isNone
:
\[\tau = \frac{1}{\sigma\K\^{2}}\]If
sigma
isNone
:
\[\sigma = \frac{1}{\tau\K\^{2}}\]To monitor the convergence of the algorithm, we compute the primal/dual objectives and the primaldual gap in
update_objective()
.The primal objective is
\[f(Kx) + g(x)\]and the dual objective is
\[ g^{*}(K^{*}y)  f^{*}(y)\]The primaldual gap (or duality gap) is
\[f(Kx) + g(x) + g^{*}(K^{*}y) + f^{*}(y)\]and measures how close is the primaldual pair (x,y) to the primaldual solution. It is always nonnegative and is used to monitor convergence of the PDHG algorithm. For more information, see Duality Gap.
Note
The primal objective is printed if verbose=1,
pdhg.run(verbose=1)
.All the objectives are printed if verbose=2,
pdhg.run(verbose=2)
.
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 \(f^{*}\) and/or \(g\) are strongly convex. In these cases, the stepsizes \(\sigma\) and \(\tau\) are updated using the
update_step_sizes()
method. A function \(f\) is strongly convex with constant \(\gamma>0\) if\[f(x)  \frac{\gamma}{2}\x\^{2} \quad\mbox{ is convex. }\]For instance the function \(\frac{1}{2}\x\^{2}_{2}\) is \(\gamma\) strongly convex for \(\gamma\in(\infty,1]\). We say it is 1strongly convex because it is the largest constant for which \(f  \frac{1}{2}\\cdot\^{2}\) is convex.
The \(\\cdot\_{1}\) norm is not strongly convex. For more information, see Strongly Convex.
If \(g\) is strongly convex with constant \(\gamma\) then the stepsizes \(\sigma\), \(\tau\) and \(\theta\) are updated as:
\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 \(f^{*}\) is strongly convex, we swap \(\sigma\) with \(\tau\).
Note
The case where both functions are strongly convex is not available at the moment.
Todo
Implement acceleration of PDHG when both functions are strongly convex.
 set_gamma_g(value)[source]#
Set the value of the strongly convex constant for function g
 Parameters:
value (a positive number or None)
 set_gamma_fconj(value)[source]#
Set the value of the strongly convex constant for the convex conjugate of function f
 Parameters:
value (a positive number or None)
 set_up(f, g, operator, tau=None, sigma=None, initial=None, **kwargs)[source]#
Initialisation of the algorithm
 Parameters:
f (Function) – A convex function with a “simple” proximal method of its conjugate.
g (Function) – A convex function with a “simple” proximal.
operator (LinearOperator) – A Linear Operator.
sigma (positive
float
, or np.ndarray, DataContainer, BlockDataContainer, optional, default=None) – Step size for the dual problem.tau (positive
float
, or np.ndarray, DataContainer, BlockDataContainer, optional, default=None) – Step size for the primal problem.initial (DataContainer, optional, default=None) – Initial point for the PDHG algorithm.
theta (Relaxation parameter, Number, default 1.0)
 check_convergence()[source]#
Check whether convergence criterion for PDHG is satisfied with scalar values of tau and sigma
 Returns:
True if convergence criterion is satisfied. False if not satisfied or convergence is unknown.
 Return type:
Boolean
 set_step_sizes(sigma=None, tau=None)[source]#
Sets sigma and tau stepsizes for the PDHG algorithm. The step sizes can be either scalar or arrayobjects.
 Parameters:
sigma (positive
float
, or np.ndarray, DataContainer, BlockDataContainer, optional, default=None) – Step size for the dual problem.tau (positive
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.
 update_step_sizes()[source]#
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_objective()[source]#
Evaluates the primal objective, the dual objective and the primaldual gap.
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
 should_stop()#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
LADMM#
 class cil.optimisation.algorithms.LADMM(f=None, g=None, operator=None, tau=None, sigma=1.0, initial=None, **kwargs)[source]#
LADMM is the Linearized Alternating Direction Method of Multipliers (LADMM)
General form of ADMM : min_{x} f(x) + g(y), subject to Ax + By = b
Case: A = Id, B = K, b = 0 ==> min_x f(Kx) + g(x)
The quadratic term in the augmented Lagrangian is linearized for the xupdate.
Main algorithmic difference is that in ADMM we compute two proximal subproblems, where in the PDHG a proximal and proximal conjugate.
Reference (Section 8) : https://link.springer.com/content/pdf/10.1007/s1010701813211.pdf
x^{k} = prox_{ au f } (x^{k1}  tau/sigma A^{T}(Ax^{k1}  z^{k1} + u^{k1} )
z^{k} = prox_{sigma g} (Ax^{k} + u^{k1})
u^{k} = u^{k1} + Ax^{k}  z^{k}
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_output()#
Returns the current solution.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 property objective#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
 should_stop()#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
Algorithms (Stochastic)#
Consider optimisation problems that take the form of a separable sum:
where \(n\) is the number of functions. Where there is a large number of \(f_i\) or their gradients are expensive to calculate, stochastic optimisation methods could prove more efficient. There is a growing range of Stochastic optimisation algorithms available with potential benefits of faster convergence in number of iterations or in computational cost. This is an area of continued development for CIL and, depending on the properties of the \(f_i\) and the regulariser \(g\), there is a range of different options for the user.
SPDHG#
Stochastic Primal Dual Hybrid Gradient (SPDHG) is a stochastic version of PDHG and deals with optimisation problems of the form:
\[\min_{x} f(Kx) + g(x) = \min_{x} \sum f_i(K_i x) + g(x)\]
where \(f_i\) and the regulariser \(g\) need only be proper, convex and lower semicontinuous ( i.e. do not need to be differentiable). Each iteration considers just one index of the sum, potentially reducing computational cost. For more examples see our [user notebooks]( vaisral/CILDemos).
 class cil.optimisation.algorithms.SPDHG(f=None, g=None, operator=None, tau=None, sigma=None, initial=None, prob=None, gamma=1.0, **kwargs)[source]#
Stochastic Primal Dual Hybrid Gradient
Problem:
\[\min_{x} f(Kx) + g(x) = \min_{x} \sum f_i(K_i x) + g(x)\] Parameters:
f (BlockFunction) – Each must be a convex function with a “simple” proximal method of its conjugate
g (Function) – A convex function with a “simple” proximal
operator (BlockOperator) – BlockOperator must contain Linear Operators
tau (positive float, optional, default=None) – Step size parameter for Primal problem
sigma (list of positive float, optional, default=None) – List of Step size parameters for Dual problem
initial (DataContainer, optional, default=None) – Initial point for the SPDHG algorithm
prob (list of floats, optional, default=None) – List of probabilities. If None each subset will have probability = 1/number of subsets
gamma (float) – parameter controlling the tradeoff between the primal and dual step sizes
**kwargs
norms (list of floats) – precalculated list of norms of the operators
Example
Example of usage: See vaisral/CILDemos
Note
Convergence is guaranteed provided that [2, eq. (12)]:
\[\]sigma[i]^{1/2} * K[i] * tau^{1/2} ^2 < p_i for all i
Note
 Notation for primal and dual stepsizes are reversed with comparison
to PDHG.py
Note
 this code implements serial sampling only, as presented in [2]
(to be extended to more general case of [1] as future work)
References
[1]”Stochastic primaldual hybrid gradient algorithm with arbitrary sampling and imaging applications”, Chambolle, Antonin, Matthias J. Ehrhardt, Peter Richtárik, and CarolaBibiane Schonlieb, SIAM Journal on Optimization 28, no. 4 (2018): 27832808.
[2]”Faster PET reconstruction with nonsmooth priors by randomization and preconditioning”, Matthias J Ehrhardt, Pawel Markiewicz and CarolaBibiane Schönlieb, Physics in Medicine & Biology, Volume 64, Number 22, 2019.
 set_up(f, g, operator, tau=None, sigma=None, initial=None, prob=None, gamma=1.0, norms=None)[source]#
setup of the algorithm :param f: Each must be a convex function with a “simple” proximal method of its conjugate :type f: BlockFunction :param g: A convex function with a “simple” proximal :type g: Function :param operator: BlockOperator must contain Linear Operators :type operator: BlockOperator :param tau: Step size parameter for Primal problem :type tau: positive float, optional, default=None :param sigma: List of Step size parameters for Dual problem :type sigma: list of positive float, optional, default=None :param initial: Initial point for the SPDHG algorithm :type initial: DataContainer, optional, default=None :param prob: List of probabilities. If None each subset will have probability = 1/number of subsets :type prob: list of floats, optional, default=None :param gamma: parameter controlling the tradeoff between the primal and dual step sizes :type gamma: float :param **kwargs: :param norms: precalculated list of norms of the operators :type norms: list of floats
 property objective#
alias of loss
 get_last_loss(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_last_objective(return_all=False)#
Returns the last stored value of the loss function
if update_objective_interval is 1 it is the value of the objective at the current iteration. If update_objective_interval > 1 it is the last stored value.
 get_output()#
Returns the current solution.
 is_provably_convergent()#
Check if the algorithm is convergent based on the provable convergence criterion.
 property iterations#
returns the iterations at which the objective has been evaluated
 property loss#
returns the list of the values of the objective during the iteration
The length of this list may be shorter than the number of iterations run when the update_objective_interval > 1
 property max_iteration#
gets the maximum number of iterations
 run(iterations=None, callbacks: List[Callback]  None = None, verbose=1, **kwargs)#
run upto
iterations
with callbacks/logging. Parameters:
iterations – number of iterations to run. If not set the algorithm will run until
should_stop()
is reachedverbose – 0=quiet, 1=info, 2=debug
callbacks – list of callables which are passed the current Algorithm object each iteration. Defaults to
[ProgressCallback(verbose)]
.
 should_stop()#
default stopping criterion: number of iterations
The user can change this in concrete implementation of iterative algorithms.
Approximate gradient methods#
Alternatively, consider that, in addition to the functions \(f_i\) and the regulariser \(g\) being proper, convex and lower semicontinuous, the \(f_i\) are differentiable. In this case we consider stochastic methods that replace a gradient calculation in a deterministic algorithm with a, potentially cheaper to calculate, approximate gradient. For example, when \(g(x)=0\), the standard Gradient Descent algorithm utilises iterations of the form
\[x_{k+1}=x_k\alpha \nabla f(x_k) =x_k\alpha \sum_{i=0}^{n1}\nabla f_i(x_k).\]
Replacing, \(\nabla f(x_k)=\sum_{i=0}^{n1}\nabla f_i(x_k)\) with \(n \nabla f_i(x_k)\), for an index \(i\) which changes each iteration, leads to the well known stochastic gradient descent algorithm.
In addition, if \(g(x)\neq 0\) and has a calculable proximal ( need not be differentiable) one can consider ISTA iterations:
\[x_{k+1}=prox_{\alpha g}(x_k\alpha \nabla f(x_k) )=prox_{\alpha g}(x_k\alpha \sum_{i=0}^{n1}\nabla f_i(x_k))\]
and again replacing \(\nabla f(x_k)=\sum_{i=0}^{n1}\nabla f_i(x_k)\) with an approximate gradient.
In a similar way, plugging approximate gradient calculations into deterministic algorithms can lead to a range of stochastic algorithms. In the following table, the left hand column has the approximate gradient function subclass, Approximate Gradient base class the header row has one of CIL’s deterministic optimisation algorithm and the body of the table has the resulting stochastic algorithm.
GD 
ISTA 
FISTA 

SGFunction 
SGD 
ProxSGD 
AccProxSGD 
SAGFunction* 
SAG 
ProxSAG 
AccProxSAG 
SAGAFunction* 
SAGA 
ProxSAGA 
AccProxSAGA 
SVRGFunction* 
SVRG 
ProxSVRG 
AccProxSVRG 
LSVRGFunction* 
LSVRG 
ProxLSVRG 
AccProxLSVRG 
*In development
The stochastic gradient functions can be found listed under functions in the documentation.
Stochastic Gradient Descent Example#
The below is an example of Stochastic Gradient Descent built of the SGFunction and Gradient Descent algorithm:
from cil.optimisation.utilities import Sampler
from cil.optimisation.algorithms import GD
from cil.optimisation.functions import LeastSquares, SGFunction
from cil.utilities import dataexample
from cil.plugins.astra.operators import ProjectionOperator
# get the data
data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get()
data.reorder('astra')
data = data.get_slice(vertical='centre')
# create the geometries
ag = data.geometry
ig = ag.get_ImageGeometry()
# partition the data and build the projectors
n_subsets = 10
partitioned_data = data.partition(n_subsets, 'sequential')
A_partitioned = ProjectionOperator(ig, partitioned_data.geometry, device = "cpu")
# create the list of functions for the stochastic sum
list_of_functions = [LeastSquares(Ai, b=bi) for Ai,bi in zip(A_partitioned, partitioned_data)]
#define the sampler and the stochastic gradient function
sampler = Sampler.staggered(len(list_of_functions), stride=2)
f = SGFunction(list_of_functions, sampler=sampler)
#set up and run the gradient descent algorithm
alg = GD(initial=ig.allocate(0), objective_function=f, step_size=1/f.L)
alg.run(300)
Operators#
The two most important methods are direct
and adjoint
methods that describe the result of applying the operator, and its
adjoint respectively, onto a compatible DataContainer
input.
The output is another DataContainer
object or subclass
hereof. An important special case is to represent the tomographic
forward and backprojection operations.
Operator base classes#
All operators extend the Operator
class. A special class is the LinearOperator
which represents an operator for which the adjoint
operation is defined.
A ScaledOperator
represents the multiplication of any operator with a scalar.
 class cil.optimisation.operators.Operator(domain_geometry, **kwargs)[source]#
Operator that maps from a space X > Y
 Parameters:
domain_geometry (ImageGeometry or AcquisitionGeometry) – domain of the operator
range_geometry (ImageGeometry or AcquisitionGeometry, optional, default None) – range of the operator
 direct(x, out=None)[source]#
Calls the operator
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the Operator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the Operator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
 norm(**kwargs)[source]#
Returns the norm of the Operator. On first call the norm will be calculated using the operator’s calculate_norm method. Subsequent calls will return the cached norm.
 Returns:
norm
 Return type:
positive:float
 set_norm(norm=None)[source]#
Sets the norm of the operator to a custom value.
 Parameters:
norm (float, optional) – Positive real valued number or None
Note
The passed values are cached so that when self.norm() is called, the saved value will be returned and not calculated via the power method. If None is passed, the cache is cleared prompting the function to call the power method to calculate the norm the next time self.norm() is called.
 class cil.optimisation.operators.LinearOperator(domain_geometry, **kwargs)[source]#
Linear operator that maps from a space X <> Y
 Parameters:
domain_geometry (ImageGeometry or AcquisitionGeometry) – domain of the operator
range_geometry (ImageGeometry or AcquisitionGeometry, optional, default None) – range of the operator
 adjoint(x, out=None)[source]#
Returns the adjoint/inverse operation evaluated at the point \(x\)
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the Operator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the Operator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
Note
Only available to linear operators
 static PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e05, return_all=False, method='auto')[source]#
Power method or Power iteration algorithm
The Power method computes the largest (dominant) eigenvalue of a matrix in magnitude, e.g., absolute value in the real case and modulus in the complex case.
 Parameters:
operator (LinearOperator)
max_iteration (positive:int, default=10) – Number of iterations for the Power method algorithm.
initial (DataContainer, default = None) – Starting point for the Power method.
tolerance (positive:float, default = 1e5) – Stopping criterion for the Power method. Check if two consecutive eigenvalue evaluations are below the tolerance.
return_all (boolean, default = False) – Toggles the verbosity of the return
method (string one of “auto”, “composed_with_adjoint” and “direct_only”, default = “auto”) – The default auto lets the code choose the method, this can be specified with “direct_only” or “composed_with_adjoint”
 Returns:
dominant eigenvalue (positive:float)
number of iterations (positive:int) – Number of iterations run. Only returned if return_all is True.
eigenvector (DataContainer) – Corresponding eigenvector of the dominant eigenvalue. Only returned if return_all is True.
list of eigenvalues (
list
) – List of eigenvalues. Only returned if return_all is True.convergence (boolean) – Check on wether the difference between the last two iterations is less than tolerance. Only returned if return_all is True.
Note
The power method contains two different algorithms chosen by the method flag.
In the case method=”direct_only”, for operator, \(A\), the power method computes the iterations \(x_{k+1} = A (x_k/\x_{k}\)\) initialised with a random vector \(x_0\) and returning the largest (dominant) eigenvalue in magnitude given by \(\x_k\\).
In the case method=”composed_with_adjoint”, the algorithm computes the largest (dominant) eigenvalue of \(A^{T}A\) returning the square root of this value, i.e. the iterations: \(x_{k+1} = A^TA (x_k/\x_{k}\)\) and returning \(\sqrt{\x_k\}\).
The default flag is method=”auto”, the algorithm checks to see if the operator.domain_geometry() == operator.range_geometry() and if so uses the method “direct_only” and if not the method “composed_with_adjoint”.
Examples
>>> M = np.array([[1.,0],[1.,2.]]) >>> Mop = MatrixOperator(M) >>> Mop_norm = Mop.PowerMethod(Mop) >>> Mop_norm 2.0000654846240296
PowerMethod is called when we compute the norm of a matrix or a LinearOperator.
>>> Mop_norm = Mop.norm() 2.0005647295658866
 calculate_norm()[source]#
Returns the norm of the LinearOperator calculated by the PowerMethod with default values.
 static dot_test(operator, domain_init=None, range_init=None, tolerance=1e06, **kwargs)[source]#
Does a dot linearity test on the operator Evaluates if the following equivalence holds .. math:
Ax\times y = y \times A^Tx
 Parameters:
operator – operator to test the dot_test
range_init – optional initialisation container in the operator range
domain_init – optional initialisation container in the operator domain
seed (int, default = 1) – Seed random generator
tolerance (float, default 1e6) – Check if the following expression is below the tolerance
math:: (..) – Ax\times y  y \times A^Tx/(Axy + 1e12) < tolerance
 Return type:
boolean, True if the test is passed.
 class cil.optimisation.operators.ScaledOperator(operator, scalar, **kwargs)[source]#
A class to represent the scalar multiplication of an Operator with a scalar. It holds an operator and a scalar. Basically it returns the multiplication of the result of direct and adjoint of the operator with the scalar. For the rest it behaves like the operator it holds.
Parameters
a scalar multiplier
Example
The scaled operator behaves like the following:
sop = ScaledOperator(operator, scalar) sop.direct(x) = scalar * operator.direct(x) sop.adjoint(x) = scalar * operator.adjoint(x) sop.norm() = operator.norm() sop.range_geometry() = operator.range_geometry() sop.domain_geometry() = operator.domain_geometry()
 class cil.optimisation.operators.CompositionOperator(*operators, **kwargs)[source]#
Composes one or more operators. For example, CompositionOperator(left, right).direct(x) is equivalent to left.direct(right.direct(x))
 Parameters:
args (`Operator`s) – Operators to be composed. As in mathematical notation, the operators will be applied right to left
 direct(x, out=None)[source]#
Calls the composition operator at the point \(x\).
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the CompositionOperator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the CompositionOperator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
 adjoint(x, out=None)[source]#
Calls the adjoint of the composition operator at the point \(x\).
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the range of the CompositionOperator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the adjoint of the CompositionOperator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
 class cil.optimisation.operators.DiagonalOperator(diagonal, domain_geometry=None)[source]#
Performs an elementwise multiplication, i.e., Hadamard Product of a
DataContainer
x andDataContainer
diagonal, d .\[(D\circ x) = \sum_{i,j}^{M,N} D_{i,j} x_{i, j}\]In matrixvector interpretation, if D is a \(M\times N\) dense matrix and is flattened, we have a \(M*N \times M*N\) vector. A sparse diagonal matrix, i.e.,
DigaonalOperator
can be created if we add the vector above to the main diagonal. If theDataContainer
x is also flattened, we have a \(M*N\) vector. Now, matrixvector multiplcation is allowed and results to a \((M*N,1)\) vector. After reshaping we recover a \(M\times N\)DataContainer
. Parameters:
diagonal (DataContainer) – DataContainer with the same dimensions as the data to be operated on.
domain_geometry (ImageGeometry) – Specifies the geometry of the operator domain. If ‘None’ will use the diagonal geometry directly. default=None .
 class cil.optimisation.operators.ChannelwiseOperator(op, channels, dimension='prepend')[source]#
ChannelwiseOperator: takes in a singlechannel operator op and the number of channels to be used, and creates a new multichannel ChannelwiseOperator, which will apply the operator op independently on each channel for the number of channels specified.
ChannelwiseOperator supports simple operators as input but not BlockOperators. Typically if such behaviour is desired, it can be achieved by creating instead a BlockOperator of ChannelwiseOperators.
 param op:
Singlechannel operator
 param channels:
Number of channels
 param dimension:
‘prepend’ (default) or ‘append’ channel dimension onto existing dimensions
 class cil.optimisation.operators.SumOperator(operator1, operator2)[source]#
Sums two operators. For example, SumOperator(left, right).direct(x) is equivalent to left.direct(x)+right.direct(x)
 Parameters:
operator1 (Operator) – The first Operator in the sum
operator2 (Operator) – The second Operator in the sum
Note
Both operators must have the same domain and range.
 direct(x, out=None)[source]#
Calls the sum operator
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the SumOperator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the SumOperator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
 adjoint(x, out=None)[source]#
Calls the adjoint of the sum operator, evaluated at the point \(x\).
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the range of the SumOperator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the adjoint of the SumOperator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
Trivial operators#
Trivial operators are the following.
 class cil.optimisation.operators.IdentityOperator(domain_geometry, range_geometry=None)[source]#
IdentityOperator: Id: X > Y, Id(x) = xin Y
X : gm_domain Y : gm_range ( Default: Y = X )
 class cil.optimisation.operators.ZeroOperator(domain_geometry, range_geometry=None)[source]#
ZeroOperator: O: X > Y, maps any element of \(x\in X\) into the zero element \(\in Y, O(x) = O_{Y}\)
 Parameters:
gm_domain – domain of the operator
gm_range – range of the operator, default: same as domain
Note:
\[O^{*}: Y^{*} > X^{*} \text{(Adjoint)} < O(x), y > = < x, O^{*}(y) >\]
 class cil.optimisation.operators.MatrixOperator(A)[source]#
Matrix wrapped into a LinearOperator
 Param:
a numpy matrix
 direct(x, out=None)[source]#
Calls the operator
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the Operator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the Operator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
 adjoint(x, out=None)[source]#
Returns the adjoint/inverse operation evaluated at the point \(x\)
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the Operator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the Operator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
Note
Only available to linear operators
 class cil.optimisation.operators.MaskOperator(mask, domain_geometry=None)[source]#
 Parameters:
mask (DataContainer) – Boolean array with the same dimensions as the data to be operated on.
domain_geometry (ImageGeometry) – Specifies the geometry of the operator domain. If ‘None’ will use the mask geometry size and spacing as float32. default = None .
 class cil.optimisation.operators.ProjectionMap(domain_geometry, index, range_geometry=None)[source]#
Projection Map or Canonical Projection (https://en.wikipedia.org/wiki/Projection_(mathematics))
Takes an element \(x = (x_{0},\dots,x_{i},\dots,x_{n})\) from a Cartesian product space \(X_{1}\times\cdots\times X_{n}\rightarrow X_{i}\) and projects it to element \(x_{i}\) specified by the index \(i\).
\[\pi_{i}: X_{1}\times\cdots\times X_{n}\rightarrow X_{i}\]\[\pi_{i}(x_{0},\dots,x_{i},\dots,x_{n}) = x_{i}\]The adjoint operation, is defined as
\[\pi_{i}^{*}(x_{i}) = (0, \cdots, x_{i}, \cdots, 0)\] Parameters:
domain_geometry (BlockGeometry) – The domain of the ProjectionMap. A BlockGeometry is expected.
index (int) – Index to project to the corresponding ImageGeometry
 direct(x, out=None)[source]#
Returns the ith (index) element of the Block data container, \(x\)
 Parameters:
x (BlockDataContainer)
out (DataContainer, default None) – If out is not None the output of the ProjectionMap will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer
 adjoint(x, out=None)[source]#
Returns a BlockDataContainer of zeros with the ith (index) filled with the DataContainer, \(x\)
 Parameters:
x (DataContainer)
out (BlockDataContainer, default None) – If out is not None the output of the adjoint of the ProjectionMap will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
BlockDataContainer
GradientOperator#
 class cil.optimisation.operators.GradientOperator(domain_geometry, method='forward', bnd_cond='Neumann', **kwargs)[source]#
Gradient Operator: Computes firstorder forward/backward differences on 2D, 3D, 4D ImageData under Neumann/Periodic boundary conditions
 Parameters:
domain_geometry (ImageGeometry) – Set up the domain of the function
method (str, default 'forward') – Accepts: ‘forward’, ‘backward’, ‘centered’, note C++ optimised routine only works with ‘forward’
bnd_cond (str, default, 'Neumann') – Set the boundary conditions to use ‘Neumann’ or ‘Periodic’
**kwargs –
 correlation: str, default ‘Space’
’Space’ will compute the gradient on only the spatial dimensions, ‘SpaceChannels’ will include the channel dimension direction
 backend: str, default ‘c’
’c’ or ‘numpy’, defaults to ‘c’ if correlation is ‘SpaceChannels’ or channels = 1
 num_threads: int
If backend is ‘c’ specify the number of threads to use. Default is number of cpus/2
 split: boolean
If ‘True’, and backend ‘c’ will return a BlockDataContainer with grouped spatial domains. i.e. [Channel, [Z, Y, X]], otherwise [Channel, Z, Y, X]
 Returns:
a BlockDataContainer containing images of the derivatives order given by dimension_labels i.e. [‘horizontal_y’,’horizontal_x’] will return [d(‘horizontal_y’), d(‘horizontal_x’)]
 Return type:
Example
2D example
\begin{eqnarray} \nabla : X \rightarrow Y\\ u \in X, \nabla(u) &=& [\partial_{y} u, \partial_{x} u]\\ u^{*} \in Y, \nabla^{*}(u^{*}) &=& \partial_{y} v1 + \partial_{x} v2 \end{eqnarray} direct(x, out=None)[source]#
Computes the firstorder forward differences
 Parameters:
x (ImageData)
out (BlockDataContainer, optional) – preallocated output memory to store result
 Returns:
result data if out not specified
 Return type:
 adjoint(x, out=None)[source]#
Computes the firstorder backward differences
 Parameters:
x (BlockDataContainer) – Gradient images for each dimension in ImageGeometry domain
out (ImageData, optional) – preallocated output memory to store result
 Returns:
result data if out not specified
 Return type:
 calculate_norm()[source]#
Returns the analytical norm of the GradientOperator.
\[\begin{split}(\partial_{z}, \partial_{y}, \partial_{x}) &= \sqrt{\\partial_{z}\^{2} + \\partial_{y}\^{2} + \\partial_{x}\^{2} } \\ &= \sqrt{ \frac{4}{h_{z}^{2}} + \frac{4}{h_{y}^{2}} + \frac{4}{h_{x}^{2}}}\end{split}\]Where the voxel sizes in each dimension are equal to 1 this simplifies to:
2D geometries \(norm = \sqrt{8}\)
3D geometries \(norm = \sqrt{12}\)
 class cil.optimisation.operators.FiniteDifferenceOperator(domain_geometry, range_geometry=None, direction=None, method='forward', bnd_cond='Neumann')[source]#
Computes forward/backward/centered finite differences of a DataContainer under Neumann/Periodic boundary conditions
 Parameters:
domain_geometry – Domain geometry for the FiniteDifferenceOperator
direction (string label from domain geometry or integer number) – Direction to evaluate finite differences
method ('forward', 'backward', 'centered') – Method for finite differences
bnd_cond – ‘Neumann’, ‘Periodic’
 direct(x, out=None)[source]#
Calls the operator
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the Operator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the Operator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
 adjoint(x, out=None)[source]#
Returns the adjoint/inverse operation evaluated at the point \(x\)
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the Operator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the Operator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
Note
Only available to linear operators
 class cil.optimisation.operators.SparseFiniteDifferenceOperator(domain_geometry, range_geometry=None, direction=0, bnd_cond='Neumann')[source]#
Create Sparse Matrices for the Finite Difference Operator
 direct(x)[source]#
Calls the operator
 Parameters:
x (DataContainer or BlockDataContainer) – Element in the domain of the Operator
out (DataContainer or BlockDataContainer, default None) – If out is not None the output of the Operator will be filled in out, otherwise a new object is instantiated and returned.
 Return type:
DataContainer or BlockDataContainer containing the result.
 class cil.optimisation.operators.SymmetrisedGradientOperator(domain_geometry, bnd_cond='Neumann', **kwargs)[source]#
The symmetrised gradient is the operator, \(E\), defined by \(E: V \rightarrow W\) where V is BlockGeometry and W is the range of the Symmetrised Gradient and
\[\begin{split}E(v) = 0.5 ( \nabla v + (\nabla v)^{T} ) \\\end{split}\]In 2 dimensions, let \(v(x,y)=(v_1(x,y),v_2(x,y))\) which gives
\[\begin{split}\nabla v =\left( \begin{matrix} \partial_{x} v_1 & \partial_x v_2\\ \partial_{y}v_1 & \partial_y v_2 \end{matrix}\right)\end{split}\]and thus
\[\begin{split}E(v) = 0.5 ( \nabla v + (\nabla v)^{T} ) =\left( \begin{matrix} \partial_{x} v_1 & 0.5 (\partial_{y} v_1 + \partial_{x} v_2) \\ 0.5 (\partial_{x} v_1 + \partial_{y} v_2) & \partial_{y} v_2 \end{matrix}\right)\end{split}\] Parameters:
domain_geometry (BlockGeometry with shape (2,1) or (3,1)) – Set up the domain of the function.
bnd_cond (str, optional, default
Neumann
) – Boundary condition eitherNeumann
orPeriodic
correlation (str, optional, default
Channel
) – Correlation eitherSpaceChannel
orChannel
WaveletOperator#
We utilise PyWavelets (https://pywavelets.readthedocs.io/en/latest/index.html) to build wavelet operators in CIL:
 class cil.optimisation.operators.WaveletOperator(domain_geometry, range_geometry=None, level=None, wname='haar', axes=None, **kwargs)[source]#
Computes forward or inverse (adjoint) discrete wavelet transform (DWT) of the input
 Parameters:
domain_geometry (cil geometry) – Domain geometry for the WaveletOperator
range_geometry (cil geometry, optional) – Output geometry for the WaveletOperator. Default = domain_geometry with the right coefficient array size deduced from pywavelets
level (int, optional, default= log_2(min(shape(axes)))) – integer for decomposition level. Default = log_2(min(shape(axes))), i.e. the maximum number of accurate downsamplings possible
wname (string, optional, default='haar') – label for wavelet used.
axes (list of ints, optional, default=`None`) – Defines the dimensions to decompose along. Note that channel is the first dimension: for example, spatial DWT is given by axes=range(1,3) and channelwise DWT is axes=range(1) Default = None, meaning all dimensions are transformed. Same as axes = range(ndim)
**kwargs#
 correlation: str, default ‘All’. Note: Only applied if axes = None!
‘All’ will compute the wavelet decomposition on every possible dimension. ‘Space’ will compute the wavelet decomposition on only the spatial dimensions. If there are multiple channels, each channel is decomposed independently. ‘Channels’ will compute the wavelet decomposition on only the channels, independently for every spatial point.
 bnd_cond: str, default ‘symmetric’. More commonly known as the padding or extension method used in discrete convolutions. All options supported by PyWavelets are valid.
Most common examples are ‘symmetric’ (padding by mirroring edge values), ‘zero’ (padding with zeros), ‘periodic’ (wrapping values around as in circular convolution). Some padding methods can have unexpected effect on the wavelet coefficients at the edges. See https://pywavelets.readthedocs.io/en/latest/ref/signalextensionmodes.html for more details and all options.
true_adjoint: bool, default True. For biorthogonal wavelets the true mathematical adjoint should no longer produce perfect reconstructions, setting true_adjoint`as `False the reconstruction (using the adjoint) should be (almost) the original input.
Note
The default decomposition level is the theoretical maximum: log_2(min(input.shape)). However, this is not always recommended and pywavelets should give a warning if the coarsest scales are too small to be meaningful.
Note
We currently do not support wavelets that are not orthogonal or biorthogonal.
 direct(x, out=None)[source]#
Returns the value of the WaveletOperator applied to \(x\)
 Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the WaveletOperator applied to \(x\) or None if out
 adjoint(Wx, out=None)[source]#
Returns the value of the adjoint of the WaveletOperator applied to \(x\)
 Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the adjoint of the WaveletOperator applied to \(x\) or None if out
Functions#
A Function
represents a mathematical function of one or more inputs
and is intended to accept DataContainers
as input as well as any
additional parameters.
Fixed parameters can be passed in during the creation of the function object.
The methods of the function reflect the properties of it, for example, if the function
represented is differentiable the function should contain a method gradient
which should return the gradient of the function evaluated at an input point.
If the function is not differentiable but allows a simple proximal operator,
the method proximal
should return the proximal operator evaluated at an
input point. The function value is evaluated by calling the function itself,
e.g. f(x)
for a Function f
and input point x
.
Base classes#
 class cil.optimisation.functions.Function(L=None)[source]#
Abstract class representing a function
 Parameters:
L (number, positive, default None) – Lipschitz constant of the gradient of the function F(x), when it is differentiable.
Note
The Lipschitz of the gradient of the function is a positive real number, such that \(\f'(x)  f'(y)\ \leq L \xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 gradient(x, out=None)[source]#
Returns the value of the gradient of function \(F\) evaluated at \(x\), if it is differentiable
\[F'(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the function at x.
 proximal(x, tau, out=None)[source]#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 convex_conjugate(x)[source]#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal_conjugate(x, tau, out=None)[source]#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 centered_at(center)[source]#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 class cil.optimisation.functions.SumFunction(*functions)[source]#
SumFunction represents the sum of \(n\geq2\) functions
\[(F_{1} + F_{2} + ... + F_{n})(\cdot) = F_{1}(\cdot) + F_{2}(\cdot) + ... + F_{n}(\cdot)\] Parameters:
*functions (Functions) – Functions to set up a
SumFunction
 Raises:
ValueError – If the number of function is strictly less than 2.
Examples
\[F(x) = \x\^{2} + \frac{1}{2}\x  1\^{2}\]>>> from cil.optimisation.functions import L2NormSquared >>> from cil.framework import ImageGeometry >>> f1 = L2NormSquared() >>> f2 = 0.5 * L2NormSquared(b = ig.allocate(1)) >>> F = SumFunction(f1, f2)
\[F(x) = \sum_{i=1}^{50} \x  i\^{2}\]>>> F = SumFunction(*[L2NormSquared(b=i) for i in range(50)])
 property L#
Returns the Lipschitz constant for the SumFunction
\[L = \sum_{i} L_{i}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 property Lmax#
Returns the maximum Lipschitz constant for the SumFunction
\[L = \max_{i}\{L_{i}\}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 gradient(x, out=None)[source]#
Returns the value of the sum of the gradient of functions evaluated at \(x\), if all of them are differentiable.
\[(F'_{1} + F'_{2} + ... + F'_{n})(x) = F'_{1}(x) + F'_{2}(x) + ... + F'_{n}(x)\] Parameters:
x (DataContainer) – Point to evaluate the gradient at.
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the sum of the gradients evaluated at point \(x\).
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 class cil.optimisation.functions.ScaledFunction(function, scalar)[source]#
ScaledFunction represents the scalar multiplication with a Function.
Let a function F then and a scalar \(\alpha\).
If \(G(x) = \alpha F(x)\) then:
\(G(x) = \alpha F(x)\) ( __call__ method )
\(G'(x) = \alpha F'(x)\) ( gradient method )
\(G^{*}(x^{*}) = \alpha F^{*}(\frac{x^{*}}{\alpha})\) ( convex_conjugate method )
\(\text{prox}_{\tau G}(x) = \text{prox}_{(\tau\alpha) F}(x)\) ( proximal method )
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 convex_conjugate(x)[source]#
Returns the convex conjugate of the scaled function.
\[G^{*}(x^{*}) = \alpha F^{*}(\frac{x^{*}}{\alpha})\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the scaled function.
 gradient(x, out=None)[source]#
Returns the gradient of the scaled function evaluated at \(x\).
\[G'(x) = \alpha F'(x)\] Parameters:
x (DataContainer) – Point to evaluate the gradient at.
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the scaled function evaluated at \(x\).
 proximal(x, tau, out=None)[source]#
Returns the proximal operator of the scaled function, evaluated at \(x\).
\[\text{prox}_{\tau G}(x) = \text{prox}_{(\tau\alpha) F}(x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the scaled function evaluated at \(x\) with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)[source]#
This returns the proximal conjugate operator for the function at \(x\), \(\tau\)
 Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal conjugate operator for the function evaluated at \(x\) and \(\tau\).
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 class cil.optimisation.functions.SumScalarFunction(function, constant)[source]#
SumScalarFunction represents the sum a function with a scalar.
\[(F + scalar)(x) = F(x) + scalar\]Although SumFunction has no general expressions for
convex_conjugate
proximal
proximal_conjugate
if the second argument is a ConstantFunction then we can derive the above analytically.
 convex_conjugate(x)[source]#
Returns the convex conjugate of a \((F+scalar)\), evaluated at \(x\).
\[(F+scalar)^{*}(x^{*}) = F^{*}(x^{*})  scalar\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate evaluated at \(x\).
 proximal(x, tau, out=None)[source]#
Returns the proximal operator of \(F+scalar\)
\[ext{prox}_{ au (F+scalar)}(x) = ext{prox}_{ au F}\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the evaluation of the proximal operator evaluated at \(x\) and :math:` au`.
 property L#
Returns the Lipschitz constant for the SumFunction
\[L = \sum_{i} L_{i}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 property Lmax#
Returns the maximum Lipschitz constant for the SumFunction
\[L = \max_{i}\{L_{i}\}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 gradient(x, out=None)#
Returns the value of the sum of the gradient of functions evaluated at \(x\), if all of them are differentiable.
\[(F'_{1} + F'_{2} + ... + F'_{n})(x) = F'_{1}(x) + F'_{2}(x) + ... + F'_{n}(x)\] Parameters:
x (DataContainer) – Point to evaluate the gradient at.
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the sum of the gradients evaluated at point \(x\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 class cil.optimisation.functions.TranslateFunction(function, center)[source]#
TranslateFunction represents the translation of function F with respect to the center b.
Let a function F and consider \(G(x) = F(x  center)\).
Function F is centered at 0, whereas G is centered at point b.
If \(G(x) = F(x  b)\) then:
\(G(x) = F(x  b)\) ( __call__ method )
\(G'(x) = F'(x  b)\) ( gradient method )
\(G^{*}(x^{*}) = F^{*}(x^{*}) + <x^{*}, b >\) ( convex_conjugate method )
\(\text{prox}_{\tau G}(x) = \text{prox}_{\tau F}(x  b) + b\) ( proximal method )
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 gradient(x, out=None)[source]#
Returns the gradient of the translated function.
\[G'(x) = F'(x  b)\] Parameters:
x (DataContainer) – Point to evaluate the gradient at.
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the gradient of the translated function evaluated at \(x\).
 proximal(x, tau, out=None)[source]#
Returns the proximal operator of the translated function.
\[\text{prox}_{\tau G}(x) = \text{prox}_{\tau F}(xb) + b\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the translated function at \(x\) and \(\tau\).
 convex_conjugate(x)[source]#
Returns the convex conjugate of the translated function.
\[G^{*}(x^{*}) = F^{*}(x^{*}) + <x^{*}, b >\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the translated function at \(x\).
Simple functions#
 class cil.optimisation.functions.ConstantFunction(constant=0)[source]#
ConstantFunction: \(F(x) = constant, constant\in\mathbb{R}\)
 gradient(x, out=None)[source]#
Returns the value of the gradient of the function, \(F'(x)=0\) :param x: Point to evaluate the gradient at. :type x: DataContainer :param out: :type out: return DataContainer, if None a new DataContainer is returned, default None.
 Return type:
A DataContainer of zeros, the same size as \(x\).
 convex_conjugate(x)[source]#
The convex conjugate of constant function \(F(x) = c\in\mathbb{R}\) is
\[\begin{split}F(x^{*}) = \begin{cases} c, & if x^{*} = 0\\ \infty, & \mbox{otherwise} \end{cases}\end{split}\]However, \(x^{*} = 0\) only in the limit of iterations, so in fact this can be infinity. We do not want to have inf values in the convex conjugate, so we have to penalise this value accordingly. The following penalisation is useful in the PDHG algorithm, when we compute primal & dual objectives for convergence purposes.
\[F^{*}(x^{*}) = \sum \max\{x^{*}, 0\}\] Parameters:
x (DataContainer)
 Return type:
The maximum of x and 0, summed over the entries of x.
 proximal(x, tau, out=None)[source]#
Returns the proximal operator of the constant function, which is the same element, i.e.,
\[\text{prox}_{\tau F}(x) = x\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, equal to \(x\).
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 class cil.optimisation.functions.ZeroFunction[source]#
ZeroFunction represents the zero function, \(F(x) = 0\)
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
The convex conjugate of constant function \(F(x) = c\in\mathbb{R}\) is
\[\begin{split}F(x^{*}) = \begin{cases} c, & if x^{*} = 0\\ \infty, & \mbox{otherwise} \end{cases}\end{split}\]However, \(x^{*} = 0\) only in the limit of iterations, so in fact this can be infinity. We do not want to have inf values in the convex conjugate, so we have to penalise this value accordingly. The following penalisation is useful in the PDHG algorithm, when we compute primal & dual objectives for convergence purposes.
\[F^{*}(x^{*}) = \sum \max\{x^{*}, 0\}\] Parameters:
x (DataContainer)
 Return type:
The maximum of x and 0, summed over the entries of x.
 gradient(x, out=None)#
Returns the value of the gradient of the function, \(F'(x)=0\) :param x: Point to evaluate the gradient at. :type x: DataContainer :param out: :type out: return DataContainer, if None a new DataContainer is returned, default None.
 Return type:
A DataContainer of zeros, the same size as \(x\).
 proximal(x, tau, out=None)#
Returns the proximal operator of the constant function, which is the same element, i.e.,
\[\text{prox}_{\tau F}(x) = x\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, equal to \(x\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 class cil.optimisation.functions.Rosenbrock(alpha, beta)[source]#
Rosenbrock function
\[\]F(x,y) = (alpha  x)^2 + beta(yx^2)^2
The function has a global minimum at .. math:: (x,y)=(alpha, alpha^2)
 gradient(x, out=None)[source]#
Gradient of the Rosenbrock function
\[\]nabla f(x,y) = left[ 2*((xalpha)  2beta x(yx^2)) ; 2beta (y  x^2) right]
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Composition of operator and a function#
This class allows the user to write a function which does the following:
where \(A\) is an operator. For instance the least squares function l2norm_ Norm2Sq
can
be expressed as
 class cil.optimisation.functions.OperatorCompositionFunction(function, operator)[source]#
Composition of a function with an operator as : \((F \otimes A)(x) = F(Ax)\)
 parameter function:
Function
F parameter operator:
Operator
A
For general operator, we have no explicit formulas for convex_conjugate, proximal and proximal_conjugate
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Indicator box#
 class cil.optimisation.functions.IndicatorBox(lower=None, upper=None, accelerated=True)[source]#
Indicator function for box constraint
\[\begin{split}f(x) = \mathbb{I}_{[a, b]} = \begin{cases} 0, \text{ if } x \in [a, b] \\ \infty, \text{otherwise} \end{cases}\end{split}\] Parameters:
lower (float, DataContainer or numpy array, default None) – Lower bound. If set to None, it is equivalent to
np.inf
.upper (float, DataContainer or numpy array, default None) – Upper bound. If set to None, it is equivalent to
np.inf
.accelerated (bool, default True) – Specifies whether to use the accelerated version or not, using numba or numpy backends respectively.
If
lower
orupper
are passed aDataContainer
(or derived class such asImageData
orAcquisitionData
) or anumpy array
, the bounds can be set to different values for each element.In order to save computing time it is possible to suppress the evaluation of the function. This is achieved by setting
suppress_evaluation
toTrue
.IndicatorBox
evaluated on any input will then return 0.If
accelerated
is set toTrue
(default), the Numba backend is used. Otherwise, the Numpy backend is used. An optional parameter to set the number of threads used by Numba can be set withset_num_threads
. Setting the number of threads whenaccelerate
is set toFalse
will not have any effect. The default number of threads is defined in thecil.utilities.multiprocessing
module, and it is equivalent to half of the CPU cores available.Example:#
In order to save computing time it is possible to suppress the evaluation of the function.
ib = IndicatorBox(lower=0, upper=1) ib.set_suppress_evaluation(True) ib.evaluate(x) # returns 0
Example:#
Set the number of threads used in accelerated mode.
num_threads = 4 ib = IndicatorBox(lower=0, upper=1) ib.set_num_threads(num_threads)
 set_suppress_evaluation(value)[source]#
Suppresses the evaluation of the function
 Parameters:
value (bool) – If True, the function evaluation on any input will return 0, without calculation.
 proximal(x, tau, out=None)[source]#
Proximal operator of IndicatorBox at x
\[prox_{\tau * f}(x)\] Parameters:
x (DataContainer) – Input to the proximal operator
tau (float) – Step size. Notice it is ignored in IndicatorBox
out (DataContainer, optional) – Output of the proximal operator. If not provided, a new DataContainer is created.
Note
tau
is ignored but it is in the signature of the generic Function class
 gradient(x, out=None)[source]#
IndicatorBox is not differentiable, so calling gradient will raise a
ValueError
 property num_threads#
Get the optional number of threads parameter to use for the accelerated version.
Defaults to the value set in the CIL multiprocessing module.
 set_num_threads(value)[source]#
Set the optional number of threads parameter to use for the accelerated version.
This is discarded if
accelerated=False
.
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
KullbackLeibler#
 class cil.optimisation.functions.KullbackLeibler(b, eta=None, mask=None, backend='numba')[source]#
Kullback Leibler
\[\begin{split}F(u, v) = \begin{cases} u \log(\frac{u}{v})  u + v & \mbox{ if } u > 0, v > 0\\ v & \mbox{ if } u = 0, v \ge 0 \\ \infty, & \mbox{otherwise} \end{cases}\end{split}\]where the \(0\log0 := 0\) convention is used.
 Parameters:
b (DataContainer, nonnegative) – Translates the function at point b.
eta (DataContainer, default = 0) – Background noise
mask (DataContainer, default = None) – Mask for the data b
backend ({'numba','numpy'}, optional) – Backend for the KullbackLeibler methods. Numba is the default backend.
Note
The KullbackLeibler function is used in practice as a fidelity term in minimisation problems where the acquired data follow Poisson distribution. If we denote the acquired data with
b
then, we write\[\underset{i}{\sum} F(b_{i}, (v + \eta)_{i})\]where, \(\eta\) is an additional noise.
In the case of Positron Emission Tomography reconstruction \(\eta\) represents scatter and random events contribution during the PET acquisition. Hence, in that case the KullbackLeibler fidelity measures the distance between \(\mathcal{A}v + \eta\) and acquisition data \(b\), where \(\mathcal{A}\) is the projection operator. This is related to PoissonLogLikelihoodWithLinearModelForMean , definition that is used in PET reconstruction in the SIRF software.
Note
The default implementation uses the buildin function kl_div from scipy. The methods of the
KullbackLeibler
are accelerated provided that numba library is installed.Examples
>>> from cil.optimisation.functions import KullbackLeibler >>> from cil.framework import ImageGeometry >>> ig = ImageGeometry(3,4) >>> data = ig.allocate('random') >>> F = KullbackLeibler(b = data)
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 gradient(x, out=None)#
Returns the value of the gradient of function \(F\) evaluated at \(x\), if it is differentiable
\[F'(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the function at x.
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
L1 Norm#
 class cil.optimisation.functions.L1Norm(b=None, weight=None)[source]#
L1Norm function
Consider the following cases:
 \[F(x) = x_{1}\]
 \[F(x) = x  b_{1}\]
In the weighted case, \(w\) is an array of nonnegative weights.
 \[F(x) = x_{L^1(w)}\]
 \[F(x) = x  b_{L^1(w)}\]
with \(x_{L^1(w)} =  x w_1 = \sum_{i=1}^{n} x_i w_i\).
 Parameters:
weight (DataContainer, numpy ndarray, default None) – Array of nonnegative weights. If
None
returns the L1 Norm.b (DataContainer, default None) – Translation of the function.
 convex_conjugate(x)[source]#
Returns the value of the convex conjugate of the L1 Norm function at x.
This is the indicator of the unit \(L^{\infty}\) norm:
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{\infty}\leq1\}}(x^{*})\]
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{\infty}\leq1\}}(x^{*}) + \langle x^{*},b\rangle\]
\[\begin{split}\mathbb{I}_{\{\\cdot\_{\infty}\leq1\}}(x^{*}) = \begin{cases} 0, \mbox{if } \x^{*}\_{\infty}\leq1\\ \infty, \mbox{otherwise} \end{cases}\end{split}\]In the weighted case the convex conjugate is the indicator of the unit \(L^{\infty}(w^{1})\) norm.
See: https://math.stackexchange.com/questions/1533217/convexconjugateofl1normfunctionwithweight
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{L^\infty(w^{1})}\leq 1\}}(x^{*})\]
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{L^\infty(w^{1})}\leq 1\}}(x^{*}) + \langle x^{*},b\rangle\]
with \(\x\_{L^\infty(w^{1})} = \max_{i} \frac{x_i}{w_i}\) and possible cases of 0/0 are defined to be 1..
 Parameters:
x (DataContainer) – where to evaluate the convex conjugate of the L1 Norm function.
 Returns:
the value of the convex conjugate of the WeightedL1Norm function at x
 Return type:
DataContainer.
 proximal(x, tau, out=None)[source]#
Returns the value of the proximal operator of the L1 Norm function at x with scaling parameter tau.
Consider the following cases:
 \[\mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}_\tau(x)\]
 \[\mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}_\tau(x  b) + b\]
where,
\[\mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}(x) = sgn(x) * \max\{ x  \tau, 0 \}\]The weighted case follows from Example 6.23 in Chapter 6 of “FirstOrder Methods in Optimization” by Amir Beck, SIAM 2017 https://archive.siam.org/books/mo25/mo25_ch6.pdf
 \[\mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}_{\tau*w}(x)\]
 \[\mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}_{\tau*w}(x  b) + b\]
 Parameters:
x (DataContainer)
tau (float, real, ndarray, DataContainer)
out (DataContainer, default None) – If not None, the result will be stored in this object.
 Returns:
The value of the proximal operator of the L1 norm function at x
 Return type:
DataContainer.
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 gradient(x, out=None)#
Returns the value of the gradient of function \(F\) evaluated at \(x\), if it is differentiable
\[F'(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the function at x.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
L2 Norm Squared#
 class cil.optimisation.functions.L2NormSquared(**kwargs)[source]#
L2NormSquared function: \(F(x) = \ x\^{2}_{2} = \underset{i}{\sum}x_{i}^{2}\)
Following cases are considered:
\(F(x) = \x\^{2}_{2}\)
\(F(x) = \x  b\^{2}_{2}\)
 Parameters:
b (DataContainer, optional) – Translation of the function
Note
For case b) we can use
F = L2NormSquared().centered_at(b)
, see TranslateFunction.Example
>>> F = L2NormSquared() >>> F = L2NormSquared(b=b) >>> F = L2NormSquared().centered_at(b)
 gradient(x, out=None)[source]#
Returns the value of the gradient of the L2NormSquared function at x.
Following cases are considered:
\(F'(x) = 2x\)
\(F'(x) = 2(xb)\)
 convex_conjugate(x)[source]#
Returns the value of the convex conjugate of the L2NormSquared function at x.
Consider the following cases:
 \[F^{*}(x^{*}) = \frac{1}{4}\x^{*}\^{2}_{2}\]
 \[F^{*}(x^{*}) = \frac{1}{4}\x^{*}\^{2}_{2} + \langle x^{*}, b\rangle\]
 proximal(x, tau, out=None)[source]#
Returns the value of the proximal operator of the L2NormSquared function at x.
Consider the following cases:
 \[\text{prox}_{\tau F}(x) = \frac{x}{1+2\tau}\]
 \[\text{prox}_{\tau F}(x) = \frac{xb}{1+2\tau} + b\]
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 class cil.optimisation.functions.WeightedL2NormSquared(**kwargs)[source]#
WeightedL2NormSquared function: \(F(x) = \x\_{W,2}^2 = \Sigma_iw_ix_i^2 = \langle x, Wx\rangle = x^TWx\) where \(W=\text{diag}(weight)\) if weight is a DataContainer or \(W=\text{weight} I\) if weight is a scalar.
 Parameters:
**kwargs
weight (a scalar or a DataContainer with the same shape as the intended domain of this WeightedL2NormSquared function)
b (a DataContainer with the same shape as the intended domain of this WeightedL2NormSquared function) – A shift so that the function becomes \(F(x) = \ xb\_{W,2}^2 = \Sigma_iw_i(x_ib_i)^2 = \langle xb, W(xb) \rangle = (xb)^TW(xb)\)
 gradient(x, out=None)[source]#
Returns the value of \(F'(x) = 2Wx\) or, if b is defined, \(F'(x) = 2W(xb)\) where \(W=\text{diag}(weight)\) if weight is a DataContainer or \(\text{weight}I\) if weight is a scalar.
 convex_conjugate(x)[source]#
Returns the value of the convex conjugate of the WeightedL2NormSquared function at x.
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 proximal(x, tau, out=None)[source]#
Returns the value of the proximal operator of the WeightedL2NormSquared function at x.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Least Squares#
 class cil.optimisation.functions.LeastSquares(A, b, c=1.0, weight=None)[source]#
(Weighted) Least Squares function
\[F(x) = c\Axb\_2^2\]or if weighted
\[F(x) = c\Axb\_{2,W}^{2}\]where \(W=\text{diag}(weight)\).
 Parameters:
A (LinearOperator)
b (Data, DataContainer)
c (Scaling Constant, float, default 1.0)
weight (DataContainer with all positive elements of size of the range of operator A, default None)
Note
L is the Lipshitz Constant of the gradient of \(F\) which is \(2 c A_2^2 = 2 c \sigma_1(A)^2\), or \(2 c W A_2^2 = 2cW \sigma_1(A)^2\), where \(\sigma_1(A)\) is the largest singular value of \(A\) and \(W=\text{diag}(weight)\).
 gradient(x, out=None)[source]#
Returns the value of the gradient of \(F(x)\):
\[F'(x) = 2cA^T(Axb)\]or
\[F'(x) = 2cA^T(W(Axb))\]where \(W=\text{diag}(weight)\).
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
L1 Sparsity#
 class cil.optimisation.functions.L1Sparsity(Q, b=None, weight=None)[source]#
L1Sparsity function
Calculates the following cases, depending on if the optional parameter weight or data b is passed. For weight=None:
 \[F(x) = Qx_{1}\]
 \[F(x) = Qx  b_{1}\]
In the weighted case, weight = \(w\) is an array of nonnegative weights.
 \[F(x) = Qx_{L^1(w)}\]
 \[F(x) = Qx  b_{L^1(w)}\]
with \(x_{L^1(w)} =  x \cdot w_1 = \sum_{i=1}^{n} x_i w_i\).
In all cases \(Q\) is an orthogonal operator.
 Parameters:
Q (orthogonal Operator) – Note that for the correct calculation of the proximal the provided operator must be orthogonal
b (Data, DataContainer, default is None)
weight (array, optional, default=None) – nonnegative weight array matching the size of the range of operator \(Q\).
 convex_conjugate(x)[source]#
Returns the value of the convex conjugate of the L1Sparsity function at x. Here, we need to use the convex conjugate of L1Sparsity, which is the Indicator of the unit \(\ell^{\infty}\) norm on the range of the (bijective) operator Q.
Consider the nonweighted case:
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{\infty}\leq1\}}(Qx^{*})\]
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{\infty}\leq1\}}(Qx^{*}) + \langle Qx^{*},b\rangle\]
\[\begin{split}\mathbb{I}_{\{\\cdot\_{\infty}\leq1\}}(x^{*}) = \begin{cases} 0, \mbox{if } \x^{*}\_{\infty}\leq1\\ \infty, \mbox{otherwise} \end{cases}\end{split}\]In the weighted case the convex conjugate is the indicator of the unit \(L^{\infty}( w^{1} )\) norm.
See: https://math.stackexchange.com/questions/1533217/convexconjugateofl1normfunctionwithweight
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{L^\infty(w^{1})}\leq 1\}}(Qx^{*})\]
 \[F^{*}(x^{*}) = \mathbb{I}_{\{\\cdot\_{L^\infty(w^{1})}\leq 1\}}(Qx^{*}) + \langle Qx^{*},b\rangle\]
with \(\x\_{L^\infty(w^{1})} = \max_{i} \frac{x_i}{w_i}\) and possible cases of 0 / 0 are defined to be 1.
 proximal(x, tau, out=None)[source]#
Returns the value of the proximal operator of the L1 Norm function at x with scaling parameter tau.
Consider the following cases:
 \[\mathrm{prox}_{\tau F}(x) = Q^T \mathrm{ShinkOperator}_{\tau}(Qx)\]
 \[\mathrm{prox}_{\tau F}(x) = Q^T \left( \mathrm{ShinkOperator}_\tau(Qx b) + b \right)\]
where,
\[\mathrm{prox}_{\tau  \cdot }(x) = \mathrm{ShinkOperator}(x) = sgn(x) * \max\{ x  \tau, 0 \}\]The weighted case follows from Example 6.23 in Chapter 6 of “FirstOrder Methods in Optimization” by Amir Beck, SIAM 2017 https://archive.siam.org/books/mo25/mo25_ch6.pdf
 \[\mathrm{prox}_{\tau F}(x) = Q^T \mathrm{ShinkOperator}_{\tau*w}(Qx)\]
 \[\mathrm{prox}_{\tau F}(x) = Q^T \left( \mathrm{ShinkOperator}_{\tau*w}(Qxb) + b \right)\]
 Parameters:
x (DataContainer)
tau (float, ndarray, DataContainer)
out (DataContainer, default None) – If not None, the result will be stored in this object.
 Returns:
The value of the proximal operator of the L1 norm function at x
 Return type:
DataContainer.
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 gradient(x, out=None)#
Returns the value of the gradient of function \(F\) evaluated at \(x\), if it is differentiable
\[F'(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the function at x.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Mixed L21 norm#
 class cil.optimisation.functions.MixedL21Norm(**kwargs)[source]#
MixedL21Norm function: \(F(x) = x_{2,1} = \sum x_{2} = \sum \sqrt{ (x^{1})^{2} + (x^{2})^{2} + \dots}\)
where x is a BlockDataContainer, i.e., \(x=(x^{1}, x^{2}, \dots)\)
 convex_conjugate(x)[source]#
Returns the value of the convex conjugate of the MixedL21Norm function at x.
This is the Indicator function of \(\mathbb{I}_{\{\\cdot\_{2,\infty}\leq1\}}(x^{*})\),
i.e.,
\[\begin{split}\mathbb{I}_{\{\\cdot\_{2, \infty}\leq1\}}(x^{*}) = \begin{cases} 0, \mbox{if } \x\_{2, \infty}\leq1\\ \infty, \mbox{otherwise} \end{cases}\end{split}\]where,
\[\x\_{2,\infty} = \max\{ \x\_{2} \} = \max\{ \sqrt{ (x^{1})^{2} + (x^{2})^{2} + \dots}\}\]
 proximal(x, tau, out=None)[source]#
Returns the value of the proximal operator of the MixedL21Norm function at x.
\[\mathrm{prox}_{\tau F}(x) = \frac{x}{\x\_{2}}\max\{ \x\_{2}  \tau, 0 \}\]where the convention 0 · (0/0) = 0 is used.
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 gradient(x, out=None)#
Returns the value of the gradient of function \(F\) evaluated at \(x\), if it is differentiable
\[F'(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the function at x.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Smooth Mixed L21 norm#
 class cil.optimisation.functions.SmoothMixedL21Norm(epsilon)[source]#
SmoothMixedL21Norm function: \(F(x) = x_{2,1} = \sum x_{2} = \sum \sqrt{ (x^{1})^{2} + (x^{2})^{2} + \epsilon^2 + \dots}\)
where x is a BlockDataContainer, i.e., \(x=(x^{1}, x^{2}, \dots)\)
Conjugate, proximal and proximal conjugate methods no closedform solution
 gradient(x, out=None)[source]#
Returns the value of the gradient of the SmoothMixedL21Norm function at x.
frac{x}{x}
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Mixed L11 norm#
 class cil.optimisation.functions.MixedL11Norm(**kwargs)[source]#
MixedL11Norm function
\[F(x) = x_{1,1} = \sum x_{1} + x_{2} + \cdots + x_{n}\]Note
MixedL11Norm is a separable function, therefore it can also be defined using the
BlockFunction
.See also
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 gradient(x, out=None)#
Returns the value of the gradient of function \(F\) evaluated at \(x\), if it is differentiable
\[F'(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the function at x.
 proximal(x, tau, out=None)[source]#
Returns the value of the proximal operator of the MixedL11Norm function at x.
\[\mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}(x)\]where,
\[\mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}(x) := sgn(x) * \max\{ x  \tau, 0 \}\]
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Total variation#
 class cil.optimisation.functions.TotalVariation(max_iteration=10, tolerance=None, correlation='Space', backend='c', lower=None, upper=None, isotropic=True, split=False, strong_convexity_constant=0, warm_start=True)[source]#
Total variation Function
\[\mathrm{TV}(u) := \\nabla u\_{2,1} = \sum \\nabla u\_{2},\, (\mbox{isotropic})\]\[\mathrm{TV}(u) := \\nabla u\_{1,1} = \sum \\nabla u\_{1}\, (\mbox{anisotropic})\]Notes
The
TotalVariation
(TV)Function
acts as a composite function, i.e., the composition of theMixedL21Norm
function and theGradientOperator
operator,\[f(u) = \u\_{2,1}, \Rightarrow (f\circ\nabla)(u) = f(\nabla x) = \mathrm{TV}(u)\]In that case, the proximal operator of TV does not have an exact solution and we use an iterative algorithm to solve:
\[\mathrm{prox}_{\tau \mathrm{TV}}(b) := \underset{u}{\mathrm{argmin}} \frac{1}{2\tau}\u  b\^{2} + \mathrm{TV}(u)\]The algorithm used for the proximal operator of TV is the Fast Gradient Projection algorithm (or FISTA) applied to the _dual problem_ of the above problem, see [1], [2], [12].
See also “Multicontrast MRI Reconstruction with StructureGuided Total Variation”, Ehrhardt, Betcke, 2016.
 Parameters:
max_iteration (
int
, default = 5) – Maximum number of iterations for the FGP algorithm to solve to solve the dual problem of the Total Variation Denoising problem (ROF). If warm_start=False, this should be around 100, or larger, with a set tolerance.tolerance (
float
, default = None) –Stopping criterion for the FGP algorithm used to to solve the dual problem of the Total Variation Denoising problem (ROF). If the difference between iterates in the FGP algorithm is less than the tolerance the iterations end before the max_iteration number.
\[\x^{k+1}  x^{k}\_{2} < \mathrm{tolerance}\]correlation (
str
, default = Space) – Correlation between Space and/or SpaceChannels for theGradientOperator
.backend (
str
, default = c) – Backend to compute theGradientOperator
lower (
'float
, default = None) – A constraint is enforced using theIndicatorBox
function, e.g.,IndicatorBox(lower, upper)
.upper (
'float
, default = None) – A constraint is enforced using theIndicatorBox
function, e.g.,IndicatorBox(lower, upper)
.isotropic (
boolean
, default = True) –Use either isotropic or anisotropic definition of TV.
\[x_{2} = \sqrt{x_{1}^{2} + x_{2}^{2}},\, (\mbox{isotropic})\]\[x_{1} = x_{1} + x_{2}\, (\mbox{anisotropic})\]split (
boolean
, default = False) – Splits the Gradient into spatial gradient and spectral or temporal gradient for multichannel data.strong_convexity_constant (
float
, default = 0) –A strongly convex term weighted by the
strong_convexity_constant
(\(\gamma\)) parameter is added to the Total variation. Now theTotalVariation
function is \(\gamma\)  strongly convex and the proximal operator is\[\underset{u}{\mathrm{argmin}} \frac{1}{2\tau}\u  b\^{2} + \mathrm{TV}(u) + \frac{\gamma}{2}\u\^{2} \Leftrightarrow\]\[\underset{u}{\mathrm{argmin}} \frac{1}{2\frac{\tau}{1+\gamma\tau}}\u  \frac{b}{1+\gamma\tau}\^{2} + \mathrm{TV}(u)\]warm_start (
boolean
, default = True) – If set to true, the FGP algorithm used to solve the dual problem of the Total Variation Denoising problem (ROF) is initiated by the final value from the previous iteration and not at zero. This allows the max_iteration value to be reduced to 510 iterations.
Note
With warm_start set to the default, True, the TV function will keep in memory the range of the gradient of the image to be denoised, i.e. N times the dimensionality of the image. This increases the memory requirements. However, during the evaluation of proximal the memory requirements will be unchanged as the same amount of memory will need to be allocated and deallocated.
Note
In the case where the Total variation becomes a \(\gamma\)  strongly convex function, i.e.,
\[\mathrm{TV}(u) + \frac{\gamma}{2}\u\^{2}\]\(\gamma\) should be relatively small, so as the second term above will not act as an additional regulariser. For more information, see [11], [3].
Examples
\[\underset{u}{\mathrm{argmin}} \frac{1}{2}\u  b\^{2} + \alpha\\nabla u\_{2,1}\]>>> alpha = 2.0 >>> TV = TotalVariation() >>> sol = TV.proximal(b, tau = alpha)
Examples
\[\underset{u}{\mathrm{argmin}} \frac{1}{2}\u  b\^{2} + \alpha\\nabla u\_{1,1} + \mathbb{I}_{C}(u)\]where \(C = \{1.0\leq u\leq 2.0\}\).
>>> alpha = 2.0 >>> TV = TotalVariation(isotropic=False, lower=1.0, upper=2.0) >>> sol = TV.proximal(b, tau = alpha)
Examples
\[\underset{u}{\mathrm{argmin}} \frac{1}{2}\u  b\^{2} + (\alpha\\nabla u\_{2,1} + \frac{\gamma}{2}\u\^{2})\]>>> alpha = 2.0 >>> gamma = 1e3 >>> TV = alpha * TotalVariation(isotropic=False, strong_convexity_constant=gamma) >>> sol = TV.proximal(b, tau = 1.0)
 proximal(x, tau, out=None)[source]#
Returns the proximal operator of the TotalVariation function at
x
.
 convex_conjugate(x)[source]#
Returns the value of convex conjugate of the TotalVariation function at
x
.
 property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 gradient(x, out=None)#
Returns the value of the gradient of function \(F\) evaluated at \(x\), if it is differentiable
\[F'(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the gradient of the function at x.
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 property gradient_operator#
GradientOperator is created if it is not instantiated yet. The domain of the _gradient, is created in the __call__ and proximal methods.
Approximate Gradient base class#
 class cil.optimisation.functions.ApproximateGradientSumFunction(functions, sampler=None)[source]#
ApproximateGradientSumFunction represents the following sum
\[\sum_{i=0}^{n1} f_{i} = (f_{0} + f_{2} + ... + f_{n1})\]where there are \(n\) functions. This function class has two ways of calling gradient
full_gradient calculates the gradient of the sum \(\sum_{i=0}^{n1} \nabla f_{i}\)
gradient calls an approximate_gradient function which may be less computationally expensive to calculate than the full gradient
This class is an abstract class.
 Parameters:
functions (list of functions) – A list of functions: \([f_{0}, f_{2}, ..., f_{n1}]\). Each function is assumed to be smooth with an implemented
gradient()
method. All functions must have the same domain. The number of functions (equivalently the length of the list) must be strictly greater than 1.sampler (An instance of a CIL Sampler class (
sampler()
) or of another class which has a__next__
function implemented to output integers in \({0,...,n1}\).) – This sampler is called each timegradient
is called and sets the internalfunction_num
passed to theapproximate_gradient
function. Default isSampler.random_with_replacement(len(functions))
.
Note
We ensure that the approximate gradient is of a similar order of magnitude to the full gradient calculation. For example, in the
SGFunction
we approximate the full gradient by \(n\nabla f_i\) for an index \(i\) given by the sampler. The multiplication by \(n\) is a choice to more easily allow comparisons between stochastic and nonstochastic methods and between stochastic methods with varying numbers of subsets.Note
Each time
gradient
is called the class keeps track of which functions have been used to calculate the gradient. This may be useful for debugging or plotting after using this function in an iterative algorithm.The property
data_passes_indices
is a list of lists holding the indices of the functions that are processed in each call of gradient. This list is updated each time gradient is called by appending a list of the indices of the functions used to calculate the gradient.The property
data_passes
is a list of floats that holds the amount of data that has been processed up until each call of gradient. This list is updated each time gradient is called by appending the proportion of the data used when calculating the approximate gradient since the class was initialised (a full gradient calculation would be 1 full data pass). Warning: if your functions do not contain an equal amount of data, for example your data was not partitioned into equal batches, then you must first use the `set_data_partition_weights” function for this to be accurate.
Note
The
gradient()
returns the approximate gradient depending on an index provided by thesampler
method.Example
This class is an abstract base class, so we give an example using the SGFunction child class.
Consider the objective is to minimise:
\[\sum_{i=0}^{n1} f_{i}(x) = \sum_{i=0}^{n1}\A_{i} x  b_{i}\^{2}\]>>> list_of_functions = [LeastSquares(Ai, b=bi)] for Ai,bi in zip(A_subsets, b_subsets)) >>> f = ApproximateGradientSumFunction(list_of_functions)
>>> list_of_functions = [LeastSquares(Ai, b=bi)] for Ai,bi in zip(A_subsets, b_subsets)) >>> sampler = Sampler.sequential(len(list_of_functions)) >>> f = SGFunction(list_of_functions, sampler=sampler) >>> f.full_gradient(x) This will return :math:`\sum_{i=0}^{n1} \nabla f_{i}(x)` >>> f.gradient(x) As per the approximate gradient implementation in the SGFunction this will return :math:`\nabla f_{0}`. The choice of the `0` index is because we chose a `sequential` sampler and this is the first time we called `gradient`. >>> f.gradient(x) This will return :math:`\nabla f_{1}` because we chose a `sequential` sampler and this is the second time we called `gradient`.
 full_gradient(x, out=None)[source]#
Returns the value of the full gradient of the sum of functions at \(x\).
\[\nabla_x(f_{0} + f_{1} + ... + f_{n1})(x) = \nabla_xf_{0}(x) + \nabla_xf_{1}(x) + ... + \nabla_xf_{n1}(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Returns:
The value of the gradient of the sum function at x or nothing if out
 Return type:
 abstract approximate_gradient(x, function_num, out=None)[source]#
Returns the approximate gradient at a given point
x
given a function_number in {0,…,len(functions)1}. Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
function_num (int) – Between 0 and the number of functions in the list
 Returns:
the value of the approximate gradient of the sum function at
x
given a function_number in {0,…,len(functions)1} Return type:
 gradient(x, out=None)[source]#
Selects a random function using the sampler and then calls the approximate gradient at
x
 Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Returns:
the value of the approximate gradient of the sum function at
x
 Return type:
 set_data_partition_weights(weights)[source]#
Setter for the partition weights used to calculate the data passes
 Parameters:
weights (list of positive floats that sum to one.) – The proportion of the data held in each function. Equivalent to the proportions that you partitioned your data into.
 property data_passes_indices#
The property
data_passes_indices
is a list of lists holding the indices of the functions that are processed in each call of gradient. This list is updated each time gradient is called by appending a list of the indices of the functions used to calculate the gradient.
 property data_passes#
if your functions do not contain an equal amount of data, for example your data was not partitioned into equal batches, then you must first use the `set_data_partition_weights” function for this to be accurate.
 Type:
The property
data_passes
is a list of floats that holds the amount of data that has been processed up until each call of gradient. This list is updated each time gradient is called by appending the proportion of the data used when calculating the approximate gradient since the class was initialised (a full gradient calculation would be 1 full data pass). Warning
 property L#
Returns the Lipschitz constant for the SumFunction
\[L = \sum_{i} L_{i}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 property Lmax#
Returns the maximum Lipschitz constant for the SumFunction
\[L = \max_{i}\{L_{i}\}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
Stochastic Gradient function#
 class cil.optimisation.functions.SGFunction(functions, sampler=None)[source]#
Stochastic gradient function, a child class of
ApproximateGradientSumFunction
, which defines from a list of functions, \({f_0,...,f_{n1}}\) a SumFunction, \(f_0+...+f_{n1}\) where each time the gradient is called, thesampler
provides an index, \(i \in {0,...,n1}\) and thegradient
method returns the approximate gradient \(n \nabla_x f_i(x)\). This can be used with thecil.optimisation.algorithms
algorithmGD
to give a stochastic gradient descent algorithm. Parameters:
functions (list of functions) – A list of functions:
[f_{0}, f_{1}, ..., f_{n1}]
. Each function is assumed to be smooth with an implementedgradient()
method. Although CIL does not define a domain of aFunction
, all functions are supposed to have the same domain. The number of functions must be strictly greater than 1.sampler (An instance of a CIL Sampler class (
sampler()
) or of another class which has a__next__
function implemented to output integers in {0,…,n1}.) – This sampler is called each time gradient is called and sets the internalfunction_num
passed to theapproximate_gradient
function. Default isSampler.random_with_replacement(len(functions))
.
 approximate_gradient(x, function_num, out=None)[source]#
Returns the gradient of the function at index function_num at
x
. Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
function_num (int) – Between 0 and the number of functions in the list
 Returns:
the value of the approximate gradient of the sum function at
x
given a function_number in {0,…,len(functions)1} Return type:
 property L#
Returns the Lipschitz constant for the SumFunction
\[L = \sum_{i} L_{i}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 property Lmax#
Returns the maximum Lipschitz constant for the SumFunction
\[L = \max_{i}\{L_{i}\}\]where \(L_{i}\) is the Lipschitz constant of the smooth function \(F_{i}\).
 centered_at(center)#
 Returns a translated function, namely if we have a function \(F(x)\) the center is at the origin.
TranslateFunction is \(F(x  b)\) and the center is at point b.
 Parameters:
center (DataContainer) – The point to center the function at.
 Return type:
The translated function.
 convex_conjugate(x)#
Evaluation of the function F* at x, where F* is the convex conjugate of function F,
\[F^{*}(x^{*}) = \underset{x}{\sup} \langle x^{*}, x \rangle  F(x)\] Parameters:
x (DataContainer)
 Return type:
The value of the convex conjugate of the function at x.
 property data_passes#
if your functions do not contain an equal amount of data, for example your data was not partitioned into equal batches, then you must first use the `set_data_partition_weights” function for this to be accurate.
 Type:
The property
data_passes
is a list of floats that holds the amount of data that has been processed up until each call of gradient. This list is updated each time gradient is called by appending the proportion of the data used when calculating the approximate gradient since the class was initialised (a full gradient calculation would be 1 full data pass). Warning
 property data_passes_indices#
The property
data_passes_indices
is a list of lists holding the indices of the functions that are processed in each call of gradient. This list is updated each time gradient is called by appending a list of the indices of the functions used to calculate the gradient.
 full_gradient(x, out=None)#
Returns the value of the full gradient of the sum of functions at \(x\).
\[\nabla_x(f_{0} + f_{1} + ... + f_{n1})(x) = \nabla_xf_{0}(x) + \nabla_xf_{1}(x) + ... + \nabla_xf_{n1}(x)\] Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Returns:
The value of the gradient of the sum function at x or nothing if out
 Return type:
 gradient(x, out=None)#
Selects a random function using the sampler and then calls the approximate gradient at
x
 Parameters:
x (DataContainer)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Returns:
the value of the approximate gradient of the sum function at
x
 Return type:
 proximal(x, tau, out=None)#
Returns the proximal operator of function \(\tau F\) evaluated at x
\[\text{prox}_{\tau F}(x) = \underset{z}{\text{argmin}} \frac{1}{2}\z  x\^{2} + \tau F(z)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the proximal operator of the function at x with scalar \(\tau\).
 proximal_conjugate(x, tau, out=None)#
Returns the proximal operator of the convex conjugate of function \(\tau F\) evaluated at \(x^{*}\)
\[\text{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\text{argmin}} \frac{1}{2}\z^{*}  x^{*}\^{2} + \tau F^{*}(z^{*})\]Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate \(F^{*}\)
\[\text{prox}_{\tau F^{*}}(x) = x  \tau\text{prox}_{\tau^{1} F}(\tau^{1}x)\] Parameters:
x (DataContainer)
tau (scalar)
out (return DataContainer, if None a new DataContainer is returned, default None.)
 Return type:
DataContainer, the value of the proximal operator of the convex conjugate at point \(x\) for scalar \(\tau\) or None if out.
 set_data_partition_weights(weights)#
Setter for the partition weights used to calculate the data passes
 Parameters:
weights (list of positive floats that sum to one.) – The proportion of the data held in each function. Equivalent to the proportions that you partitioned your data into.
Utilities#
Contains utilities for the CIL optimisation framework.
Samplers#
Here, we define samplers that select from a list of indices {0, 1, …, N1} either randomly or by some deterministic pattern.
The cil.optimisation.utilities.sampler
class defines a function next()
which gives the next sample. It also has utility to get_samples
to access which samples have or will be drawn.
For ease of use we provide the following static methods in cil.optimisation.utilities.sampler that allow you to configure your sampler object rather than initialising the classes directly:
 static Sampler.from_function(num_indices, function, prob_weights=None)[source]#
Instantiate a sampler that wraps a function for index selection.
 num_indices: int
The sampler will select from a range of indices 0 to num_indices.
 functioncallable
A deterministic function that takes an integer as an argument, representing the iteration number, and returns an integer between 0 and num_indices. The function signature should be function(iteration_number: int) > int
 prob_weights: list of floats of length num_indices that sum to 1. Default is [1 / num_indices] * num_indices
Consider that the sampler is incremented a large number of times this argument holds the expected number of times each index would be outputted, normalised to 1.
 Sampler
An instance of the Sampler class which samples from a function.
This example creates a sampler that always outputs 2. The probability weights are passed to the sampler as they are not uniform.
>>> num_indices=3 >>> >>> def my_sampling_function(iteration_number): >>> return 2 >>> >>> sampler = Sampler.from_function(num_indices=num_indices, function=my_sampling_function, prob_weights=[0, 0, 1]) >>> print(list(sampler.get_samples(12))) [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
This example creates a sampler that outputs sequential indices, starting from 1. The probability weights are not passed to the sampler as they are uniform.
>>> num_indices=10 >>> >>> def my_sampling_function(iteration_number): >>> return (iteration_number+1)%10 >>> >>> sampler = Sampler.from_function(num_indices=num_indices, function=my_sampling_function) >>> print(list(sampler.get_samples(25))) [1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5]
This example creates a sampler that samples in order from a custom list. The num_indices is 6, although note that the index 5 is never output by the sampler. The number of indices must be at least one greater than any of the elements in the custom_list. The probability weights are passed to the sampler as they are not uniform.
>>> custom_list = [0,0,0,0,0,0,3,2,1,4] >>> num_indices = 6 >>> >>> def my_sampling_function(iteration_number, custom_list=custom_list]): >>> return(custom_list[iteration_number%len(custom_list)]) >>> >>> sampler = Sampler.from_function(num_indices=num_indices, function=my_sampling_function, prob_weights=[0.6, 0.1, 0.1, 0.1, 0.1, 0.0]) >>> print(list(sampler.get_samples(25))) [0, 0, 0, 0, 0, 0, 3, 2, 1, 4, 0, 0, 0, 0, 0, 0, 3, 2, 1, 4, 0, 0, 0, 0, 0] >>> print(sampler) Sampler that wraps a function that takes an iteration number and selects from a list of indices {0, 1, …, N1}, where N is the number of indices. Type : from_function Current iteration number : 0 number of indices : 6 Probability weights : [0.6, 0.1, 0.1, 0.1, 0.1, 0.0]
 static Sampler.sequential(num_indices)[source]#
Instantiates a sampler that outputs sequential indices.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
 Returns:
An instance of the Sampler class that will generate indices sequentially.
 Return type:
Example
>>> sampler = Sampler.sequential(10) >>> print(sampler.get_samples(5)) >>> print(next(sampler)) 0 [0 1 2 3 4] >>> print(sampler.next()) 1
 static Sampler.staggered(num_indices, stride)[source]#
Instantiates a sampler which outputs in a staggered order.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
stride (int) – The stride between returned indices. The stride should be less than the num_indices.
 Returns:
An instance of the Sampler class that will generate indices in a staggered pattern.
 Return type:
Example
>>> sampler = Sampler.staggered(num_indices=21, stride=4) >>> print(next(sampler)) 0 >>> print(sampler.next()) 4 >>> print(sampler.get_samples(5)) [ 0 4 8 12 16]
Example
>>> sampler = Sampler.staggered(num_indices=17, stride=8) >>> print(next(sampler)) 0 >>> print(sampler.next()) 8 >>> print(sampler.get_samples(10)) [ 0 8 16 1 9 2 10 3 11 4]
 static Sampler.herman_meyer(num_indices)[source]#
Instantiates a sampler which outputs in a Herman Meyer order.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices. For HermanMeyer sampling this number should not be prime.
Reference

in (The sampling method was introduced)
I (Singh)
MIC (et al. Deep Image Prior PET Reconstruction using a SIRFBased Objective  IEEE)
https (NSS & RTSD 2022.)
in
GT (Herman)
doi (Meyer LB. Algebraic reconstruction techniques can be made computationally efficient. IEEE Trans Med Imaging.)
 Returns:
An instance of the Sampler class which outputs in a Herman Meyer order.
 Return type:
Example
>>> sampler=Sampler.herman_meyer(12) >>> print(sampler.get_samples(16)) [ 0 6 3 9 1 7 4 10 2 8 5 11 0 6 3 9]
 static Sampler.random_with_replacement(num_indices, prob=None, seed=None)[source]#
Instantiates a sampler which outputs an index between 0  num_indices with a given probability.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices
prob (list of floats, optional) – The probability for each index to be selected by the ‘next’ operation. If not provided, the indices will be sampled uniformly. The list should have a length equal to num_indices, and the values should sum to 1
seed (int, optional) – Used to initialise the random number generator where repeatability is required.
 Returns:
An instance of the RandomSampler class that will generate indices randomly with replacement
 Return type:
RandomSampler
Example
>>> sampler = Sampler.random_with_replacement(5) >>> print(sampler.get_samples(10)) [3 4 0 0 2 3 3 2 2 1] >>> print(next(sampler)) 3 >>> print(sampler.next()) 4
>>> sampler = Sampler.random_with_replacement(num_indices=4, prob=[0.7,0.1,0.1,0.1]) >>> print(sampler.get_samples(10)) [0 1 3 0 0 3 0 0 0 0]
 static Sampler.random_without_replacement(num_indices, seed=None)[source]#
Instantiates a sampler which outputs an index between 0  num_indices. Once sampled the index will not be sampled again until all indices have been returned.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
seed (int, optional) – Used to initialise the random number generator where repeatability is required.
 Returns:
An instance of the RandomSampler class that will generate indices randomly without replacement
 Return type:
RandomSampler
Example
>>> sampler=Sampler.randomWithoutReplacement(num_indices=7, seed=1) >>> print(sampler.get_samples(16)) [6 2 1 0 4 3 5 1 0 4 2 5 6 3 3 2]
They will all instantiate a Sampler defined in the following class:
 class cil.optimisation.utilities.Sampler(num_indices, function, sampling_type=None, prob_weights=None)[source]#
Initialises a sampler that returns and then increments indices from a sequence defined by a function.
Static methods to easily configure several samplers are provided, such as sequential, staggered, HermanMayer, random with and without replacement.
Custom deterministic samplers can be created by using the from_function static method or by subclassing this sampler class.
 Parameters:
function (Callable[[int], int]) – A function that takes an integer iteration number and returns an integer between 0 and num_indices.
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
sampling_type (str, optional, default = None) – The sampling type used. This is recorded for reference and printed when print is called.
prob_weights (list of floats of length num_indices that sum to 1. Default is [1 / num_indices] * num_indices) – Consider that the sampler is incremented a large number of times this argument holds the expected number of times each index would be outputted, normalised to 1.
 Returns:
An instance of the Sampler class representing the desired configuration.
 Return type:
Example
>>> sampler = Sampler.random_with_replacement(5) >>> print(sampler.get_samples(20)) [3 4 0 0 2 3 3 2 2 1 1 4 4 3 0 2 4 4 2 4] >>> print(next(sampler)) 3 >>> print(sampler.next()) 4
>>> sampler = Sampler.staggered(num_indices=21, stride=4) >>> print(next(sampler)) 0 >>> print(sampler.next()) 4 >>> print(sampler.get_samples(5)) [ 0 4 8 12 16]
Example
>>> sampler = Sampler.sequential(10) >>> print(sampler.get_samples(5)) >>> print(next(sampler)) 0 [0 1 2 3 4] >>> print(sampler.next()) 1
Example
>>> sampler = Sampler.herman_meyer(12) >>> print(sampler.get_samples(16)) [ 0 6 3 9 1 7 4 10 2 8 5 11 0 6 3 9]
Example
This example creates a sampler that outputs sequential indices, starting from 1.
>>> num_indices=10 >>> >>> def my_sampling_function(iteration_number): >>> return (iteration_number+1)%10 >>> >>> sampler = Sampler.from_function(num_indices=num_indices, function=my_sampling_function) >>> print(list(sampler.get_samples(25))) [1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5]
Note
The optimal choice of sampler depends on the data and the number of calls to the sampler. Note that a low number of calls to a random sampler won’t give an even distribution. For a small number of samples (e.g. <5*num_indices) the user may wish to consider another sampling method e.g. random without replacement, which, when calling num_indices samples is guaranteed to draw each index exactly once.
 next()[source]#
Returns a sample from the list of indices `{0, 1, …, N1}, where N is the number of indices and increments the sampler.
 get_samples(num_samples)[source]#
Generates a list of the first num_samples output by the sampler. Calling this does not increment the sampler index or affect the behaviour of the sampler .
 Parameters:
num_samples (int) – The number of samples to return.
 Returns:
The first `num_samples” output by the sampler.
 Return type:
List
 static sequential(num_indices)[source]#
Instantiates a sampler that outputs sequential indices.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
 Returns:
An instance of the Sampler class that will generate indices sequentially.
 Return type:
Example
>>> sampler = Sampler.sequential(10) >>> print(sampler.get_samples(5)) >>> print(next(sampler)) 0 [0 1 2 3 4] >>> print(sampler.next()) 1
 static staggered(num_indices, stride)[source]#
Instantiates a sampler which outputs in a staggered order.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
stride (int) – The stride between returned indices. The stride should be less than the num_indices.
 Returns:
An instance of the Sampler class that will generate indices in a staggered pattern.
 Return type:
Example
>>> sampler = Sampler.staggered(num_indices=21, stride=4) >>> print(next(sampler)) 0 >>> print(sampler.next()) 4 >>> print(sampler.get_samples(5)) [ 0 4 8 12 16]
Example
>>> sampler = Sampler.staggered(num_indices=17, stride=8) >>> print(next(sampler)) 0 >>> print(sampler.next()) 8 >>> print(sampler.get_samples(10)) [ 0 8 16 1 9 2 10 3 11 4]
 static random_with_replacement(num_indices, prob=None, seed=None)[source]#
Instantiates a sampler which outputs an index between 0  num_indices with a given probability.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices
prob (list of floats, optional) – The probability for each index to be selected by the ‘next’ operation. If not provided, the indices will be sampled uniformly. The list should have a length equal to num_indices, and the values should sum to 1
seed (int, optional) – Used to initialise the random number generator where repeatability is required.
 Returns:
An instance of the RandomSampler class that will generate indices randomly with replacement
 Return type:
RandomSampler
Example
>>> sampler = Sampler.random_with_replacement(5) >>> print(sampler.get_samples(10)) [3 4 0 0 2 3 3 2 2 1] >>> print(next(sampler)) 3 >>> print(sampler.next()) 4
>>> sampler = Sampler.random_with_replacement(num_indices=4, prob=[0.7,0.1,0.1,0.1]) >>> print(sampler.get_samples(10)) [0 1 3 0 0 3 0 0 0 0]
 static random_without_replacement(num_indices, seed=None)[source]#
Instantiates a sampler which outputs an index between 0  num_indices. Once sampled the index will not be sampled again until all indices have been returned.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
seed (int, optional) – Used to initialise the random number generator where repeatability is required.
 Returns:
An instance of the RandomSampler class that will generate indices randomly without replacement
 Return type:
RandomSampler
Example
>>> sampler=Sampler.randomWithoutReplacement(num_indices=7, seed=1) >>> print(sampler.get_samples(16)) [6 2 1 0 4 3 5 1 0 4 2 5 6 3 3 2]
 static from_function(num_indices, function, prob_weights=None)[source]#
Instantiate a sampler that wraps a function for index selection.
 num_indices: int
The sampler will select from a range of indices 0 to num_indices.
 functioncallable
A deterministic function that takes an integer as an argument, representing the iteration number, and returns an integer between 0 and num_indices. The function signature should be function(iteration_number: int) > int
 prob_weights: list of floats of length num_indices that sum to 1. Default is [1 / num_indices] * num_indices
Consider that the sampler is incremented a large number of times this argument holds the expected number of times each index would be outputted, normalised to 1.
 Sampler
An instance of the Sampler class which samples from a function.
This example creates a sampler that always outputs 2. The probability weights are passed to the sampler as they are not uniform.
>>> num_indices=3 >>> >>> def my_sampling_function(iteration_number): >>> return 2 >>> >>> sampler = Sampler.from_function(num_indices=num_indices, function=my_sampling_function, prob_weights=[0, 0, 1]) >>> print(list(sampler.get_samples(12))) [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
This example creates a sampler that outputs sequential indices, starting from 1. The probability weights are not passed to the sampler as they are uniform.
>>> num_indices=10 >>> >>> def my_sampling_function(iteration_number): >>> return (iteration_number+1)%10 >>> >>> sampler = Sampler.from_function(num_indices=num_indices, function=my_sampling_function) >>> print(list(sampler.get_samples(25))) [1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5]
This example creates a sampler that samples in order from a custom list. The num_indices is 6, although note that the index 5 is never output by the sampler. The number of indices must be at least one greater than any of the elements in the custom_list. The probability weights are passed to the sampler as they are not uniform.
>>> custom_list = [0,0,0,0,0,0,3,2,1,4] >>> num_indices = 6 >>> >>> def my_sampling_function(iteration_number, custom_list=custom_list]): >>> return(custom_list[iteration_number%len(custom_list)]) >>> >>> sampler = Sampler.from_function(num_indices=num_indices, function=my_sampling_function, prob_weights=[0.6, 0.1, 0.1, 0.1, 0.1, 0.0]) >>> print(list(sampler.get_samples(25))) [0, 0, 0, 0, 0, 0, 3, 2, 1, 4, 0, 0, 0, 0, 0, 0, 3, 2, 1, 4, 0, 0, 0, 0, 0] >>> print(sampler) Sampler that wraps a function that takes an iteration number and selects from a list of indices {0, 1, …, N1}, where N is the number of indices. Type : from_function Current iteration number : 0 number of indices : 6 Probability weights : [0.6, 0.1, 0.1, 0.1, 0.1, 0.0]
 static herman_meyer(num_indices)[source]#
Instantiates a sampler which outputs in a Herman Meyer order.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices. For HermanMeyer sampling this number should not be prime.
Reference

in (The sampling method was introduced)
I (Singh)
MIC (et al. Deep Image Prior PET Reconstruction using a SIRFBased Objective  IEEE)
https (NSS & RTSD 2022.)
in
GT (Herman)
doi (Meyer LB. Algebraic reconstruction techniques can be made computationally efficient. IEEE Trans Med Imaging.)
 Returns:
An instance of the Sampler class which outputs in a Herman Meyer order.
 Return type:
Example
>>> sampler=Sampler.herman_meyer(12) >>> print(sampler.get_samples(16)) [ 0 6 3 9 1 7 4 10 2 8 5 11 0 6 3 9]
In addition, we provide a random sampling class which is a child class of cil.optimisation.utilities.sampler and provides options for sampling with and without replacement:
 class cil.optimisation.utilities.SamplerRandom(num_indices, seed=None, replace=True, prob=None, sampling_type='random_with_replacement')[source]#
The user is recommended to not instantiate this class directly but instead use one of the static methods in the parent Sampler class that will return instances of different samplers.
This class produces Samplers that output random samples with and without replacement from the set {0, 1, …, N1} where N=num_indices.
Custom random samplers can be created by subclassing this sampler class.
 Parameters:
num_indices (int) – The sampler will select from a range of indices 0 to num_indices.
sampling_type (str, optional, default = 'random_with_replacement") – The sampling type used. This is recorded for reference and printed when print is called.
prob_weights (list of floats of length num_indices that sum to 1. Default is [1 / num_indices] * num_indices) – Consider that the sampler is incremented a large number of times this argument holds the expected number of times each index would be outputted, normalised to 1.
replace (bool, default is True) – If True, sample with replace, otherwise sample without replacement
seed (int, optional) – Used to initialise the random number generator where repeatability is required.
 Returns:
An instance of the Sampler class representing the desired configuration.
 Return type:
Example
>>> sampler = Sampler.random_with_replacement(5) >>> print(sampler.get_samples(20)) [3 4 0 0 2 3 3 2 2 1 1 4 4 3 0 2 4 4 2 4]
Example
>>> sampler=Sampler.randomWithoutReplacement(num_indices=7, seed=1) >>> print(sampler.get_samples(16)) [6 2 1 0 4 3 5 1 0 4 2 5 6 3 3 2]
 get_samples(num_samples)[source]#
Generates a list of the first num_samples output by the sampler. Calling this does not increment the sampler index or affect the behaviour of the sampler .
 Parameters:
num_samples (int) – The number of samples to return.
 Returns:
The first num_samples produced by the sampler
 Return type:
list
Callbacks#
A list of Callback
s to be executed each iteration can be passed to `Algorithms`_ run
method.
from cil.optimisation.utilities.callbacks import LogfileCallback
...
algorithm.run(..., callbacks=[LogfileCallback("log.txt")])
 class cil.optimisation.utilities.callbacks.Callback(verbose=1)[source]#
Base Callback to inherit from for use in
Algorithm.run(callbacks: list[Callback])
. Parameters:
verbose (int, choice of 0,1,2, default 1) – 0=quiet, 1=info, 2=debug.
Builtin callbacks include:
 class cil.optimisation.utilities.callbacks.ProgressCallback(verbose=1, tqdm_class=<class 'tqdm.asyncio.tqdm_asyncio'>, **tqdm_kwargs)[source]#
tqdm
based progress bar. Parameters:
tqdm_class (default
tqdm.auto.tqdm
)**tqdm_kwargs – Passed to
tqdm_class
.
 class cil.optimisation.utilities.callbacks.TextProgressCallback(verbose=1, *, tqdm_class=<class 'cil.optimisation.utilities.callbacks._TqdmText'>, **tqdm_kwargs)[source]#
ProgressCallback
but printed on separate lines to screen. Parameters:
miniters (int, default
Algorithm.update_objective_interval
) – Number of algorithm iterations between screen prints.
 class cil.optimisation.utilities.callbacks.LogfileCallback(log_file, mode='a', **kwargs)[source]#
TextProgressCallback
but to a file instead of screen. Parameters:
log_file (FileDescriptorOrPath) – Passed to
open()
.mode (str) – Passed to
open()
.
Users can also write custom callbacks.
Below is an example of a custom callback implementing early stopping.
In each iteration of the TestAlgo
, the objective \(x\) is reduced by \(5\). The EarlyStopping
callback terminates the algorithm when \(x \le 15\). The algorithm thus terminates after \(3\) iterations.
from cil.optimisation.algorithms import Algorithm
from cil.optimisation.utilities import callbacks
class TestAlgo(Algorithm):
def __init__(self, *args, **kwargs):
self.x = 0
super().__init__(*args, **kwargs)
self.configured = True
def update(self):
self.x = 5
def update_objective(self):
self.loss.append(2 ** self.x)
class EarlyStopping(callbacks.Callback):
def __call__(self, algorithm: Algorithm):
if algorithm.x <= 15: # arbitrary stopping criterion
raise StopIteration
algo = TestAlgo()
algo.run(20, callbacks=[callbacks.ProgressCallback(), EarlyStopping()])
Output:
15%███  3/20 [00:00<00:00, 11770.73it/s, objective=3.05e5]
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 stepsize.
Currently in CIL we have a base class:
 class cil.optimisation.utilities.StepSizeMethods.StepSizeRule[source]#
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.
We also have a number of example classes:
 class cil.optimisation.utilities.StepSizeMethods.ConstantStepSize(step_size)[source]#
Stepsize rule that always returns a constant stepsize.
 Parameters:
step_size (float) – The stepsize to be returned with each call.
 class cil.optimisation.utilities.StepSizeMethods.ArmijoStepSizeRule(alpha=1000000.0, beta=0.5, max_iterations=None)[source]#
Applies the Armijo rule to calculate the step size (step_size).
The Armijo rule runs a while loop to find the appropriate step_size by starting from a very large number (alpha). The step_size is found by reducing the step size (by a factor beta) in an iterative way until a certain criterion is met. To avoid infinite loops, we add a maximum number of times (max_iterations) the while loop is run.
 Parameters:
alpha (float, optional, default=1e6) – The starting point for the step size iterations
beta (float between 0 and 1, optional, default=0.5) – The amount the step_size is reduced if the criterion is not met
max_iterations (integer, optional, default is numpy.ceil (2 * numpy.log10(alpha) / numpy.log10(2))) – The maximum number of iterations to find a suitable step size
Reference

Optimization (Algorithm 3.1 (Numerical) – https://projecteuclid.org/download/pdf_1/euclid.pjm/1102995080
Nocedal (//www.math.uci.edu/~qnie/Publications/NumericalOptimization.pdf)) – https://projecteuclid.org/download/pdf_1/euclid.pjm/1102995080
(https (Wright)) – https://projecteuclid.org/download/pdf_1/euclid.pjm/1102995080
Preconditioners#
A preconditioner is a class which acts on an algorithm and can be passed to cil.optimisation.algorithm.GD, cil.optimisation.algorithm.ISTA or cil.optimisation.algorithm.FISTA and it’s method apply is called after the calculation of the gradient before the gradient descent step is taken. It modifies and returns a passed gradient.
Currently in CIL we have a base class:
 class cil.optimisation.utilities.preconditioner.Preconditioner[source]#
Abstract base class for Preconditioner objects. The apply method of this class takes an initialised CIL function as an argument and modifies a provided gradient.
 abstract apply(algorithm, gradient, out)[source]#
Abstract method to apply the preconditioner.
 Parameters:
algorithm (Algorithm) – The algorithm object.
gradient (DataContainer) – The calculated gradient to modify
out (DataContainer,) – Container to fill with the modified gradient.
Note
In CIL algorithms, the preconditioners are used inplace. Make sure this method is safe to use in place.
 Returns:
The modified gradient
 Return type:
We also have a number of already provided preconditioners
 class cil.optimisation.utilities.preconditioner.Sensitivity(operator)[source]#
Sensitivity preconditioner class.
In each call to the preconditioner the gradient is multiplied by \(1/(A^T \mathbf{1})\) where \(A\) is an operator, \(\mathbf{1}\) is an object in the range of the operator filled with ones.
 Parameters:
operator (CIL Operator) – The operator used for sensitivity computation.
 compute_preconditioner_matrix()[source]#
Compute the sensitivity. \(A^T \mathbf{1}\) where \(A\) is the operator and \(\mathbf{1}\) is an object in the range of the operator filled with ones. Then perform safe division by the sensitivity to store the preconditioner array \(1/(A^T \mathbf{1})\)
 apply(algorithm, gradient, out=None)[source]#
Update the preconditioner.
 Parameters:
algorithm (object) – The algorithm object.
gradient (DataContainer) – The calculated gradient to modify
out (DataContainer,) – Container to fill with the modified gradient
 Returns:
The modified gradient
 Return type:
 class cil.optimisation.utilities.preconditioner.AdaptiveSensitivity(operator, delta=1e06, max_iterations=100, reference=None)[source]#
Adaptive Sensitivity preconditioner class.
In each call to the preconditioner the gradient is multiplied by \((x+\delta) /(A^T \mathbf{1})\) where \(A\) is an operator, \(\mathbf{1}\) is an object in the range of the operator filled with ones. The point \(x\) is the current iteration, or a reference image, and \(\delta\) is a small positive float.
 Parameters:
operator (CIL object) – The operator used for sensitivity computation.
delta (float, optional) – The delta value for the preconditioner.
reference (DataContainer e.g. ImageData, default is None) – Reference data, an object in the domain of the operator. Recommended to be a best guess reconstruction. If reference data is passed the preconditioner is always fixed.
max_iterations (int, default = 100) – The maximum number of iterations before the preconditoner is frozen and nolonger updates. Note that if reference data is passed the preconditioner is always frozen and iterations is set to 1.
Note
A reference for the freezing of the preconditioner can be found: Twyman R., Arridge S., Kereta Z., Jin B., Brusaferri L., Ahn S., Stearns CW., Hutton B.F., Burger I.A., Kotasidis F., Thielemans K.. An Investigation of Stochastic Variance Reduction Algorithms for Relative Difference Penalized 3D PET Image Reconstruction. IEEE Trans Med Imaging. 2023 Jan;42(1):2941. doi: 10.1109/TMI.2022.3203237. Epub 2022 Dec 29. PMID: 36044488.
 apply(algorithm, gradient, out=None)[source]#
Update the preconditioner.
 Parameters:
algorithm (object) – The algorithm object.
gradient (DataContainer) – The calculated gradient to modify
out (DataContainer,) – Container to fill with the modified gradient
 Returns:
The modified gradient
 Return type:
Block Framework#
To be able to express more advanced optimisation problems we developed the Block Framework, which provides a generic strategy to treat variational problems in the following form:
The block framework consists of:
The block framework allows writing more advanced optimisation problems. Consider the typical Tikhonov regularisation:
where,
\(A\) is the projection operator
\(b\) is the acquired data
\(u\) is the unknown image to be solved for
\(\alpha\) is the regularisation parameter
\(L\) is a regularisation operator
The first term measures the fidelity of the solution to the data. The second term measures the fidelity to the prior knowledge we have imposed on the system, operator \(L\).
This can be rewritten equivalently in the block matrix form:
With the definitions:
\(\tilde{A} = \binom{A}{\alpha L}\)
\(\tilde{b} = \binom{b}{0}\)
this can now be recognised as a least squares problem which can be solved by any algorithm in the cil.optimisation
which can solve least squares problem, e.g. CGLS.
To be able to express our optimisation problems in the matrix form above, we developed the socalled,
Block Framework comprising 4 main actors: BlockGeometry
, BlockDataContainer
,
BlockFunction
and BlockOperator
.
BlockDataContainer#
BlockDataContainer holds DataContainer as column vector. It is possible to do basic algebra between BlockDataContainer s and with numbers, list or numpy arrays.
 class cil.framework.BlockDataContainer(*args, **kwargs)[source]#
Class to hold DataContainers as column vector
Provides basic algebra between BlockDataContainer’s, DataContainer’s and subclasses and Numbers
algebra between `BlockDataContainer`s will be elementwise, only if the shape of the 2 `BlockDataContainer`s is the same, otherwise it will fail
algebra between BlockDataContainer`s and `list or numpy array will work as long as the number of rows and element of the arrays match, independently on the fact that the BlockDataContainer could be nested
algebra between BlockDataContainer and one DataContainer is possible. It will require all the DataContainers in the block to be compatible with the DataContainer we want to operate with.
algebra between BlockDataContainer and a Number is possible and it will be done with each element of the BlockDataContainer even if nested
A = [ [B,C] , D] A * 3 = [ 3 * [B,C] , 3* D] = [ [ 3*B, 3*C] , 3*D ]
 add(other, *args, **kwargs)[source]#
Algebra: add method of BlockDataContainer with number/DataContainer or BlockDataContainer
 Param:
other (number, DataContainer or subclasses or BlockDataContainer
 Param:
out (optional): provides a placehold for the resul.
 subtract(other, *args, **kwargs)[source]#
Algebra: subtract method of BlockDataContainer with number/DataContainer or BlockDataContainer
 Param:
other (number, DataContainer or subclasses or BlockDataContainer
 Param:
out (optional): provides a placeholder for the result.
 multiply(other, *args, **kwargs)[source]#
Algebra: multiply method of BlockDataContainer with number/DataContainer or BlockDataContainer
 Param:
other (number, DataContainer or subclasses or BlockDataContainer)
 Param:
out (optional): provides a placeholder for the result.
 divide(other, *args, **kwargs)[source]#
Algebra: divide method of BlockDataContainer with number/DataContainer or BlockDataContainer
 Param:
other (number, DataContainer or subclasses or BlockDataContainer)
 Param:
out (optional): provides a placeholder for the result.
 power(other, *args, **kwargs)[source]#
Algebra: power method of BlockDataContainer with number/DataContainer or BlockDataContainer
 Param:
other (number, DataContainer or subclasses or BlockDataContainer
 Param:
out (optional): provides a placeholder for the result.
 maximum(other, *args, **kwargs)[source]#
Algebra: power method of BlockDataContainer with number/DataContainer or BlockDataContainer
 Param:
other (number, DataContainer or subclasses or BlockDataContainer)
 Param:
out (optional): provides a placeholder for the result.
 minimum(other, *args, **kwargs)[source]#
Algebra: power method of BlockDataContainer with number/DataContainer or BlockDataContainer
 Param:
other (number, DataContainer or subclasses or BlockDataContainer)
 Param:
out (optional): provides a placeholder for the result.
 sapyb(a, y, b, out, num_threads=2)[source]#
performs axpby elementwise on the BlockDataContainer containers
Does the operation .. math:: a*x+b*y and stores the result in out, where x is self
 Parameters:
a – scalar
b – scalar
y – compatible (Block)DataContainer
out – (Block)DataContainer to store the result
Example:#
>>> a = 2 >>> b = 3 >>> ig = ImageGeometry(10,11) >>> x = ig.allocate(1) >>> y = ig.allocate(2) >>> bdc1 = BlockDataContainer(2*x, y) >>> bdc2 = BlockDataContainer(x, 2*y) >>> out = bdc1.sapyb(a,bdc2,b)
 axpby(a, b, y, out, dtype=<class 'numpy.float32'>, num_threads=2)[source]#
Deprecated method. Alias of sapyb
 binary_operations(operation, other, *args, **kwargs)[source]#
Algebra: generic method of algebric operation with BlockDataContainer with number/DataContainer or BlockDataContainer
Provides commutativity with DataContainer and subclasses, i.e. this class’s reverse algebraic methods take precedence w.r.t. direct algebraic methods of DataContainer and subclasses.
This method is not to be used directly
 unary_operations(operation, *args, **kwargs)[source]#
Unary operation on BlockDataContainer:
generic method of unary operation with BlockDataContainer: abs, sign, sqrt and conjugate
This method is not to be used directly
 __radd__(other)[source]#
Reverse addition
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
 __rsub__(other)[source]#
Reverse subtraction
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
 __rmul__(other)[source]#
Reverse multiplication
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
 __rdiv__(other)[source]#
Reverse division
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
 __rtruediv__(other)[source]#
Reverse truedivision
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
 __rpow__(other)[source]#
Reverse power
to make sure that this method is called rather than the __mul__ of a numpy array the class constant __array_priority__ must be set > 0 https://docs.scipy.org/doc/numpy1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
 __weakref__#
list of weak references to the object
Block Function#
BlockFunction acts on BlockDataContainer as a separable sum function:
\[ \begin{align}\begin{aligned}f = [f_1,...,f_n] \newline\\f([x_1,...,x_n]) = f_1(x_1) + .... + f_n(x_n)\end{aligned}\end{align} \]
 class cil.optimisation.functions.BlockFunction(*functions)[source]#
BlockFunction represents a separable sum function \(F\) defined as
\[F:X_{1}\times X_{2}\cdots\times X_{m} \rightarrow (\infty, \infty]\]where \(F\) is the separable sum of functions \((f_{i})_{i=1}^{m}\),
\[F(x_{1}, x_{2}, \cdots, x_{m}) = \overset{m}{\underset{i=1}{\sum}}f_{i}(x_{i}), \mbox{ with } f_{i}: X_{i} \rightarrow (\infty, \infty].\]A nice property (due to it’s separability structure) is that the proximal operator can be decomposed along the proximal operators of each function \(f_{i}\).
\[\mathrm{prox}_{\tau F}(x) = ( \mathrm{prox}_{\tau f_{i}}(x_{i}) )_{i=1}^{m}\]In addition, if \(\tau := (\tau_{1},\dots,\tau_{m})\), then
\[\mathrm{prox}_{\tau F}(x) = ( \mathrm{prox}_{\tau_{i} f_{i}}(x_{i}) )_{i=1}^{m}\] property L#
Lipschitz of the gradient of function f.
L is positive real number, such that \(\f'(x)  f'(y)\ \leq L\xy\\), assuming \(f: IG \rightarrow \mathbb{R}\)
 __call__(x)[source]#
Returns the value of the BlockFunction \(F\)
\[F(x) = \overset{m}{\underset{i=1}{\sum}}f_{i}(x_{i}), \mbox{ where } x = (x_{1}, x_{2}, \cdots, x_{m}), \quad i = 1,2,\dots,m\]Parameter:
x : BlockDataContainer and must have as many rows as self.length
returns ..math:: sum(f_i(x_i))
 convex_conjugate(x)[source]#
Returns the value of the convex conjugate of the BlockFunction at \(x^{*}\).
\[F^{*}(x^{*}) = \overset{m}{\underset{i=1}{\sum}}f_{i}^{*}(x^{*}_{i})\]Parameter:
x : BlockDataContainer and must have as many rows as self.length
 proximal(x, tau, out=None)[source]#
Proximal operator of the BlockFunction at x:
\[\mathrm{prox}_{\tau F}(x) = (\mathrm{prox}_{\tau f_{i}}(x_{i}))_{i=1}^{m}\]Parameter:
x : BlockDataContainer and must have as many rows as self.length
 gradient(x, out=None)[source]#
Returns the value of the gradient of the BlockFunction function at x.
\[F'(x) = [f_{1}'(x_{1}), ... , f_{m}'(x_{m})]\]Parameter:
x : BlockDataContainer and must have as many rows as self.length
 proximal_conjugate(x, tau, out=None)[source]#
Proximal operator of the convex conjugate of BlockFunction at x:
\[\mathrm{prox}_{\tau F^{*}}(x) = (\mathrm{prox}_{\tau f^{*}_{i}}(x^{*}_{i}))_{i=1}^{m}\]Parameter:
x : BlockDataContainer and must have as many rows as self.length
 __rmul__(other)[source]#
Define multiplication with a scalar
 Parameters:
other – number
Returns a new BlockFunction containing the product of the scalar with all the functions in the block
Block Operator#
BlockOperator represent a block matrix with operators
Column: Share the same domains \(X_{1}, X_{2}\)
Rows: Share the same ranges \(Y_{1}, Y_{2}, Y_{3}\)
\(A_{1}, A_{3}, A_{5}\): share the same domain \(X_{1}\) and \(A_{2}, A_{4}, A_{6}\): share the same domain \(X_{2}\)
For instance with these ingredients one may write the following objective function,
where \(g\) represent the measured values, \(u\) the solution \(\nabla\) is the gradient operator, \( ~~ _{2,1}\) is a norm for the output of the gradient operator and \( xg ^2_2\) is least squares fidelity function as
Then we have rewritten the problem as
Which in Python would be like
op1 = GradientOperator(ig, correlation=GradientOperator.CORRELATION_SPACE)
op2 = IdentityOperator(ig, ag)
# Create BlockOperator
K = BlockOperator(op1, op2, shape=(2,1) )
# Create functions
F = BlockFunction(alpha * MixedL21Norm(), 0.5 * L2NormSquared(b=noisy_data))
 class cil.optimisation.operators.BlockOperator(*args, **kwargs)[source]#
A Block matrix containing Operators
 Parameters:
*args (Operator) – Operators in the block.
**kwargs (dict) – shape (
tuple
, optional): If shape is passed the Operators in vararg are considered input in a rowbyrow fashion.
Note
The Block Framework is a generic strategy to treat variational problems in the following form:
\[\min Regulariser + Fidelity\]BlockOperators have a generic shape M x N, and when applied on an Nx1 BlockDataContainer, will yield and Mx1 BlockDataContainer.
Note
BlockDatacontainer are only allowed to have the shape of N x 1, with N rows and 1 column.
User may specify the shape of the block, by default is a row vector
Operators in a Block are required to have the same domain columnwise and the same range rowwise.
Examples
BlockOperator(op0,op1) results in a row block
BlockOperator(op0,op1,shape=(1,2)) results in a column block
 get_item(row, col)[source]#
Returns the Operator at specified row and col :param row: The row index required. :type row: int :param col: The column index required. :type col: int
 norm()[source]#
Returns the Euclidean norm of the norms of the individual operators in the BlockOperators
 get_norms_as_list()[source]#
Returns a list of the individual norms of the Operators in the BlockOperator
 set_norms(norms)[source]#
Uses the set_norm() function in Operator to set the norms of the operators in the BlockOperator from a list of custom values.
 Parameters:
norms (list) – A list of positive real values the same length as the number of operators in the BlockOperator.
 direct(x, out=None)[source]#
Direct operation for the BlockOperator
Note
BlockOperators work on BlockDataContainers, but they will also work on DataContainers and inherited classes by simple wrapping the input in a BlockDataContainer of shape (1,1)
 adjoint(x, out=None)[source]#
Adjoint operation for the BlockOperator
Note
BlockOperator may contain both LinearOperator and Operator This method exists in BlockOperator as it is not known what type of Operator it will contain.
BlockOperators work on BlockDataContainers, but they will also work on DataContainers and inherited classes by simple wrapping the input in a BlockDataContainer of shape (1,1)
Raises: ValueError if the contained Operators are not linear
 get_output_shape(xshape, adjoint=False)[source]#
Returns the shape of the output BlockDataContainer :param xshape: :type xshape: BlockDataContainer :param adjoint: :type adjoint: bool
Examples
A(N,M) direct u(M,1) > N,1
A(N,M)^T adjoint u(N,1) > M,1
 __rmul__(scalar)[source]#
Defines the left multiplication with a scalar. Returns a block operator with Scaled Operators inside.
 Parameters:
scalar (number or iterable containing numbers)
 property T#
Returns the transposed of self.
Recall the input list is shaped in a rowbyrow fashion
References#
Amir Beck and Marc Teboulle. A fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009. URL: https://doi.org/10.1137/080716542, arXiv:https://doi.org/10.1137/080716542, doi:10.1137/080716542.
Amir Beck and Marc Teboulle. Fast gradientbased algorithms for constrained total variation image denoising and deblurring problems. IEEE Transactions on Image Processing, 18(11):2419–2434, 2009. doi:10.1109/TIP.2009.2028250.
Antonin Chambolle and Thomas Pock. A firstorder primaldual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120–145, May 2011. URL: https://doi.org/10.1007/s1085101002511, doi:10.1007/s1085101002511.
Patrick L. Combettes and Valérie R. Wajs. Signal recovery by proximal forwardbackward splitting. Multiscale Modeling & Simulation, 4(4):1168–1200, 2005. URL: https://doi.org/10.1137/050626090, arXiv:https://doi.org/10.1137/050626090, doi:10.1137/050626090.
Ernie Esser, Xiaoqun Zhang, and Tony F. Chan. A general framework for a class of first order primaldual algorithms for convex optimization in imaging science. SIAM Journal on Imaging Sciences, 3(4):1015–1046, 2010. URL: https://doi.org/10.1137/09076934X, arXiv:https://doi.org/10.1137/09076934X, doi:10.1137/09076934X.
Osman Güler. New proximal point algorithms for convex minimization. SIAM Journal on Optimization, 2(4):649–664, 1992.
J. S. Jörgensen, E. Ametova, G. Burca, G. Fardell, E. Papoutsellis, E. Pasca, K. Thielemans, M. Turner, R. Warr, W. R. B. Lionheart, and P. J. Withers. Core imaging library  part i: a versatile python framework for tomographic imaging. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379(2204):20200192, 2021. URL: https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2020.0192, arXiv:https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2020.0192, doi:10.1098/rsta.2020.0192.
Avinash C. Kak and Malcolm Slaney. Principles of Computerized Tomographic Imaging. Society for Industrial and Applied Mathematics, January 2001. URL: https://doi.org/10.1137/1.9780898719277, doi:10.1137/1.9780898719277.
Yurii Nesterov. Introductory lectures on convex optimization: A basic course. Volume 87. Springer Science & Business Media, 2003.
Evangelos Papoutsellis, Evelina Ametova, Claire Delplancke, Gemma Fardell, Jakob S. Jörgensen, Edoardo Pasca, Martin Turner, Ryan Warr, William R. B. Lionheart, and Philip J. Withers. Core imaging library  part ii: multichannel reconstruction for dynamic and spectral tomography. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379(2204):20200193, 2021. URL: https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2020.0193, arXiv:https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2020.0193, doi:10.1098/rsta.2020.0193.
Julian Rasch and Antonin Chambolle. Inexact firstorder primal–dual algorithms. Computational Optimization and Applications, 76(2):381–430, Jun 2020. URL: https://doi.org/10.1007/s1058902000186y, doi:10.1007/s1058902000186y.
Mingqiang Zhu, Stephen J. Wright, and Tony F. Chan. Dualitybased algorithms for totalvariationregularized image restoration. Computational Optimization and Applications, 47(3):377–400, Nov 2010. URL: https://doi.org/10.1007/s1058900892252, doi:10.1007/s1058900892252.