Optimisation framework#

This package allows rapid prototyping of optimisation-based reconstruction problems, i.e. defining and solving different optimization problems to enforce different properties on the reconstructed image.

Firstly, it provides an object-oriented framework for defining mathematical operators and functions as well a collection of useful example operators and functions. Both smooth and non-smooth functions can be used.

Further, it provides a number of high-level 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 algorithm
update is the actual iteration updating the solution
update_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 and update_objective methods

A courtesy method run is available to run n iterations. The method accepts a callbacks 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. The run method will stop when the stopping criterion is met or StopIteration is raised.

set_up(*args, **kwargs)[source]#: Set up the algorithm

update()[source]#: A single iteration of the algorithm

should_stop()[source]#

default stopping criterion: number of iterations

The user can change this in concrete implementation of iterative algorithms.

get_output()[source]#: Returns the current solution.

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.

update_objective()[source]#: calculates the objective with the current solution

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: Optional[List[cil.optimisation.utilities.callbacks.Callback]] = 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 reached
verbose – 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, alpha=1000000.0, beta=0.5, rtol=1e-05, atol=1e-08, **kwargs)[source]#

Gradient Descent algorithm

set_up(initial, objective_function, step_size)[source]#

initialisation of the algorithm

Parameters

initial – initial guess
objective_function – objective function to be minimised
step_size – step size

update()[source]#: Single iteration

update_objective()[source]#: calculates the objective with the current solution

armijo_rule()[source]#

Applies the Armijo rule to calculate the step size (step_size)

https://projecteuclid.org/download/pdf_1/euclid.pjm/1102995080

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 dividing alpha by 2 in an iterative way until a certain criterion is met. To avoid infinite loops, we add a maximum number of times the while loop is run.

This rule would allow to reach a minimum step_size of 10^-alpha.

if alpha = numpy.power(10,gamma) delta = 3 step_size = numpy.power(10, -delta) with armijo rule we can get to step_size from initial alpha by repeating the while loop k times where alpha / 2^k = step_size 10^gamma / 2^k = 10^-delta 2^k = 10^(gamma+delta) k = gamma+delta / log10(2) approx 3.3 * (gamma+delta)

if we would take by default delta = gamma kmax = numpy.ceil ( 2 * gamma / numpy.log10(2) ) kmax = numpy.ceil (2 * numpy.log10(alpha) / numpy.log10(2))

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: Optional[List[cil.optimisation.utilities.callbacks.Callback]] = 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 reached
verbose – 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()[source]#

default stopping criterion: number of iterations

The user can change this in concrete implementation of iterative algorithms.

CGLS#

class cil.optimisation.algorithms.CGLS(initial=None, operator=None, data=None, tolerance=1e-06, **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

Reference:: https://web.stanford.edu/group/SOL/software/cgls/

set_up(initial, operator, data, tolerance=1e-06)[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

update()[source]#: single iteration

update_objective()[source]#: calculates the objective with the current solution

should_stop()[source]#

default stopping criterion: number of iterations

The user can change this in concrete implementation of iterative algorithms.

flag()[source]#: returns whether the tolerance has been reached

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: Optional[List[cil.optimisation.utilities.callbacks.Callback]] = 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 reached
verbose – 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 [7].

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 constraint
upper (float, default = None) – Upper bound constraint
constraint (Function, default = None) – A function with proximal method, e.g., IndicatorBox function and IndicatorBox.proximal(), or TotalVariation function and TotalVariation.proximal().
kwargs – Keyword arguments used from the base class Algorithm.

Note

If constraint is not passed, lower and upper are used to create an IndicatorBox and apply its proximal.

If constraint is passed, proximal method is required to be implemented.

Note

The preconditioning arrays (weights) M and D 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) ) ) )\]

update_objective()[source]#: Returns the objective

\[\frac{1}{2}\|A x - b\|^{2}\]

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: Optional[List[cil.optimisation.utilities.callbacks.Callback]] = 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 reached
verbose – 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, **kwargs)[source]#

Iterative Shrinkage-Thresholding Algorithm, see [1], [2].

Iterative Shrinkage-Thresholding 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 [].

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, default = None) – Step size for the gradient step of ISTA. The default step_size is \(\frac{1}{L}\) or 1 if f=None.
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()

FISTA#

class cil.optimisation.algorithms.FISTA(initial, f, g, step_size=None, **kwargs)[source]#

Fast Iterative Shrinkage-Thresholding Algorithm, see [1], [2].

Fast Iterative Shrinkage-Thresholding 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_{k-1}}(x_{k} - x_{k-1}) \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, default = None) – Step size for the gradient step of FISTA. The default step_size is \(\frac{1}{L}\) or 1 if f=None.
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 ([8] 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 ([5] 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()

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], [4].

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_iterationint, optional, default=0
Maximum number of iterations of the PDHG.

update_objective_intervalint, optional, default=1
Evaluates objectives, e.g., primal/dual/primal-dual gap every update_objective_interval.

check_convergenceboolean, default=True
Checks scalar sigma and tau values satisfy convergence criterion

Example

In our CIL-Demos 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 [6], [9].

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 first-order primal-dual algorithm for convex optimization problems with known saddle-point structure with applications in imaging.

The general problem considered in the PDHG algorithm is the generic saddle-point 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 finite-dimensional 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 over-relaxation 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 is None and tau is None:
\[\sigma = \frac{1}{\|K\|}, \tau = \frac{1}{\|K\|}\]
- If tau is None:
\[\tau = \frac{1}{\sigma\|K\|^{2}}\]
- If sigma is None:
\[\sigma = \frac{1}{\tau\|K\|^{2}}\]
To monitor the convergence of the algorithm, we compute the primal/dual objectives and the primal-dual gap in update_objective().

The primal objective is

\[f(Kx) + g(x)\]

and the dual objective is

\[- g^{*}(-K^{*}y) - f^{*}(y)\]

The primal-dual gap (or duality gap) is

\[f(Kx) + g(x) + g^{*}(-K^{*}y) + f^{*}(y)\]

and measures how close is the primal-dual pair (x,y) to the primal-dual solution. It is always non-negative 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 step-sizes \(\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 1-strongly 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 step-sizes \(\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) –

get_output()[source]#: Returns the current solution.

update()[source]#: Performs a single iteration of the PDHG algorithm

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 step-sizes for the PDHG algorithm. The step sizes can be either scalar or array-objects.

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 primal-dual 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: Optional[List[cil.optimisation.utilities.callbacks.Callback]] = 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 reached
verbose – 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 x-update.

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/s10107-018-1321-1.pdf

x^{k} = prox_{ au f } (x^{k-1} - tau/sigma A^{T}(Ax^{k-1} - z^{k-1} + u^{k-1} )

z^{k} = prox_{sigma g} (Ax^{k} + u^{k-1})

u^{k} = u^{k-1} + Ax^{k} - z^{k}

set_up(f, g, operator, tau=None, sigma=1.0, initial=None)[source]#: Set up the algorithm

update()[source]#: A single iteration of the algorithm

update_objective()[source]#: calculates the objective with the current solution

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: Optional[List[cil.optimisation.utilities.callbacks.Callback]] = 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 reached
verbose – 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:

\[\min_{x} f(x)+g(x) = \min_{x} \sum_{i=0}^{n-1} f_{i}(x) + g(x) = \min_{x} (f_{0}(x) + f_{1}(x) + ... + f_{n-1}(x))+g(x)\]

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 semi-continuous ( 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]( vais-ral/CIL-Demos).

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 trade-off between the primal and dual step sizes
**kwargs –
norms (list of floats) – precalculated list of norms of the operators

Example

Example of usage: See vais-ral/CIL-Demos

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 step-sizes 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 primal-dual hybrid gradient algorithm with arbitrary sampling and imaging applications”, Chambolle, Antonin, Matthias J. Ehrhardt, Peter Richtárik, and Carola-Bibiane Schonlieb, SIAM Journal on Optimization 28, no. 4 (2018): 2783-2808.

[2]”Faster PET reconstruction with non-smooth priors by randomization and preconditioning”, Matthias J Ehrhardt, Pawel Markiewicz and Carola-Bibiane 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]#: set-up 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 trade-off 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

update()[source]#: A single iteration of the algorithm

update_objective()[source]#: calculates the objective with the current solution

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: Optional[List[cil.optimisation.utilities.callbacks.Callback]] = 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 reached
verbose – 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 semi-continuous, 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}^{n-1}\nabla f_i(x_k).\]

Replacing, \(\nabla f(x_k)=\sum_{i=0}^{n-1}\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}^{n-1}\nabla f_i(x_k))\]

and again replacing \(\nabla f(x_k)=\sum_{i=0}^{n-1}\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	Prox-SGD	Acc-Prox-SGD
SAGFunction*	SAG	Prox-SAG	Acc-Prox-SAG
SAGAFunction*	SAGA	Prox-SAGA	Acc-Prox-SAGA
SVRGFunction*	SVRG	Prox-SVRG	Acc-Prox-SVRG
LSVRGFunction*	LSVRG	Prox-LSVRG	Acc-Prox-LSVRG

*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

is_linear()[source]#: Returns if the operator is linear :rtype: Bool

is_orthogonal()[source]#: Returns if the operator is orthogonal :rtype: Bool

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, or None if out is passed.

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.

calculate_norm()[source]#: Returns the norm of the Operator. Note that this gives a NotImplementedError if the SumOperator is not linear. :returns: Scalar :rtype: the norm of the Operator

range_geometry()[source]#: Returns the range of the Operator: Y space

domain_geometry()[source]#: Returns the domain of the Operator: X space

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

is_linear()[source]#: Returns if the operator is linear

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, or None if out is passed.

Note

Only available to linear operators

static PowerMethod(operator, max_iteration=10, initial=None, tolerance=1e-05, 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 = 1e-5) – 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=1e-06, **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 1e-6) – Check if the following expression is below the tolerance
math:: (..) – |Ax\times y - y \times A^Tx|/(|A||x||y| + 1e-12) < 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

operator: a Operator or LinearOperator scalar: Number

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

direct(x, out=None)[source]#: direct method

adjoint(x, out=None)[source]#: adjoint method

norm(**kwargs)[source]#: norm of the operator

is_linear()[source]#: returns a boolean indicating whether the operator is linear

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, or None if out is passed.

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, or None if out is passed.

is_linear()[source]#: Returns if the operator is linear :rtype: Bool

class cil.optimisation.operators.DiagonalOperator(diagonal, domain_geometry=None)[source]#

Performs an element-wise multiplication, i.e., Hadamard Product of a DataContainer x and DataContainer diagonal, d .

\[(D\circ x) = \sum_{i,j}^{M,N} D_{i,j} x_{i, j}\]

In matrix-vector 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 the DataContainer x is also flattened, we have a \(M*N\) vector. Now, matrix-vector 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 .

direct(x, out=None)[source]#: Returns \(D\circ x\)

adjoint(x, out=None)[source]#: Returns \(D^*\circ x\)

calculate_norm(**kwargs)[source]#: Returns the operator norm of DiagonalOperator which is the \(\infty\) norm of diagonal

\[\|D\|_{\infty} = \max_{i}\{|D_{i}|\}\]

class cil.optimisation.operators.ChannelwiseOperator(op, channels, dimension='prepend')[source]#

ChannelwiseOperator: takes in a single-channel operator op and the number of channels to be used, and creates a new multi-channel 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

Single-channel operator

param channels

Number of channels

param dimension

‘prepend’ (default) or ‘append’ channel dimension onto existing dimensions

direct(x, out=None)[source]#: Returns D(x)

adjoint(x, out=None)[source]#: Returns D^{*}(y)

calculate_norm(**kwargs)[source]#: Evaluates operator norm of DiagonalOperator

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, or None if out is passed.

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, or None if out is passed.

is_linear()[source]#: Returns if the operator is linear :rtype: Bool

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 )

direct(x, out=None)[source]#: Returns Id(x)

adjoint(x, out=None)[source]#: Returns Id(x)

calculate_norm(**kwargs)[source]#: Evaluates operator norm of IdentityOperator

is_orthogonal()[source]#: Returns if the operator is orthogonal :rtype: Bool

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:

\[ \begin{align}\begin{aligned}O^{*}: Y^{*} -> X^{*} \text{(Adjoint)}\\< O(x), y > = < x, O^{*}(y) >\end{aligned}\end{align} \]

direct(x, out=None)[source]#: Returns O(x)

adjoint(x, out=None)[source]#: Returns O^{*}(y)

calculate_norm(**kwargs)[source]#: Evaluates operator norm of ZeroOperator

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, or None if out is passed.

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, or None if out is passed.

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 first-order 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

BlockDataContainer

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 first-order forward differences

Parameters

x (ImageData) –
out (BlockDataContainer, optional) – pre-allocated output memory to store result

Returns

result data if out not specified

Return type

BlockDataContainer

adjoint(x, out=None)[source]#

Computes the first-order backward differences

Parameters

x (BlockDataContainer) – Gradient images for each dimension in ImageGeometry domain
out (ImageData, optional) – pre-allocated output memory to store result

Returns

result data if out not specified

Return type

ImageData

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, or None if out is passed.

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, or None if out is passed.

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, or None if out is passed.

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 either Neumann or Periodic
correlation (str, optional, default Channel) – Correlation either SpaceChannel or Channel

direct(x, out=None)[source]#

Returns \(E(v) = 0.5 * ( \nabla v + (\nabla v)^{T} )\)

x: BlockDataContainer out: BlockDataContainer, default None

If out is not None the output of direct will be filled in out, otherwise a new object is instantiated and returned.

adjoint(x, out=None)[source]#

Returns the adjoint of the symmetrised gradient operator

x: BlockDataContainer out: BlockDataContainer, default None

If out is not None the output of adjoint will be filled in out, otherwise a new object is instantiated and returned.

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)

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/signal-extension-modes.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 bi-orthogonal.

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

calculate_norm()[source]#

Returns the norm of WaveletOperator, which is equal to 1.0 if the wavelet is orthogonal

Returns: norm
Return type: float

is_orthogonal()[source]#

Returns if the operator is orthogonal

Return type: Bool

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 \|x-y\|\), 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 or None if out

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\) or None if out.

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\|x-y\|\), 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\) or None if out.

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\) or None if out.

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\|x-y\|\), 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\) or None if out.

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\) or None if out.

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\) or None if out.

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` or None if out.

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\) or None if out.

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\|x-y\|\), 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\) or None if out.

proximal(x, tau, out=None)[source]#

Returns the proximal operator of the translated function.

\[\text{prox}_{\tau G}(x) = \text{prox}_{\tau F}(x-b) + 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\) or None if out.

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\) or None if out

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\) or None if out.

property L#

Lipschitz of the gradient of function f.

L is positive real number, such that \(\|f'(x) - f'(y)\| \leq L\|x-y\|\), 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\|x-y\|\), 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}\]

\[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)#

Return type: A DataContainer of zeros, the same size as \(x\) or None if out

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\) or None if out.

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(y-x^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*((x-alpha) - 2beta x(y-x^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\|x-y\|\), 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\) or None if out.

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:

\[F ( x ) = G ( Ax )\]

where \(A\) is an operator. For instance the least squares function l2norm_ Norm2Sq can be expressed as

\[F(x) = || Ax - b ||^2_2\]

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\|x-y\|\), assuming \(f: IG \rightarrow \mathbb{R}\)

gradient(x, out=None)[source]#

Return the gradient of F(Ax),

..math :: (F(Ax))’ = A^{T}F’(Ax)

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\) or None if out.

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 or upper are passed a DataContainer (or derived class such as ImageData or AcquisitionData) or a numpy 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 to True. IndicatorBox evaluated on any input will then return 0.

If accelerated is set to True (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 with set_num_threads. Setting the number of threads when accelerate is set to False will not have any effect. The default number of threads is defined in the cil.utilities.multiprocessing module, and it is equivalent to half of the CPU cores available.

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

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\|x-y\|\), 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, non-negative) – 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 Kullback-Leibler 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 build-in 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\|x-y\|\), 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 or None if out

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\) or None if out.

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 non-negative 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 non-negative 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/convex-conjugate-of-l1-norm-function-with-weight

\[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 “First-Order 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\|x-y\|\), 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 or None if out

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(x-b)\)

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{x-b}{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\|x-y\|\), 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) = \| x-b\|_{W,2}^2 = \Sigma_iw_i(x_i-b_i)^2 = \langle x-b, W(x-b) \rangle = (x-b)^TW(x-b)\)

gradient(x, out=None)[source]#: Returns the value of \(F'(x) = 2Wx\) or, if b is defined, \(F'(x) = 2W(x-b)\) 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\|x-y\|\), 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\|Ax-b\|_2^2\]

or if weighted

\[F(x) = c\|Ax-b\|_{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 = 2c||W|| \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(Ax-b)\]

\[F'(x) = 2cA^T(W(Ax-b))\]

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\|x-y\|\), 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\) or None if out.

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 non-negative 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) – non-negative 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 non-weighted 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/convex-conjugate-of-l1-norm-function-with-weight

\[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 “First-Order 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}(Qx-b) + 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\|x-y\|\), 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 or None if out

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\|x-y\|\), 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 or None if out

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 closed-form 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\|x-y\|\), 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\) or None if out.

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.

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 the MixedL21Norm function and the GradientOperator 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], [11].

See also “Multicontrast MRI Reconstruction with Structure-Guided 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 the GradientOperator.
backend (str, default = c) – Backend to compute the GradientOperator
lower ('float, default = None) – A constraint is enforced using the IndicatorBox function, e.g., IndicatorBox(lower, upper).
upper ('float, default = None) – A constraint is enforced using the IndicatorBox 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 the TotalVariation 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 5-10 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 [10], [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 = 1e-3
>>> 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\|x-y\|\), assuming \(f: IG \rightarrow \mathbb{R}\)

calculate_Lipschitz()[source]#: Default value for the Lipschitz constant.

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 or None if out

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}^{n-1} f_{i} = (f_{0} + f_{2} + ... + f_{n-1})\]

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}^{n-1} \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_{n-1}]\). 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,...,n-1}\).) – This sampler is called each time gradient is called and sets the internal function_num passed to the approximate_gradient function. Default is Sampler.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 non-stochastic 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 the sampler 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}^{n-1} f_{i}(x) = \sum_{i=0}^{n-1}\|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}^{n-1} \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_{n-1})(x) = \nabla_xf_{0}(x) + \nabla_xf_{1}(x) + ... + \nabla_xf_{n-1}(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

DataContainer

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

DataContainer

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

DataContainer

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\) or None if out.

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_{n-1}}\) a SumFunction, \(f_0+...+f_{n-1}\) where each time the gradient is called, the sampler provides an index, \(i \in {0,...,n-1}\) and the gradient method returns the approximate gradient \(n \nabla_x f_i(x)\). This can be used with the cil.optimisation.algorithms algorithm GD to give a stochastic gradient descent algorithm.

Parameters

functions (list of functions) – A list of functions: [f_{0}, f_{1}, ..., f_{n-1}]. Each function is assumed to be smooth with an implemented gradient() method. Although CIL does not define a domain of a Function, 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,…,n-1}.) – This sampler is called each time gradient is called and sets the internal function_num passed to the approximate_gradient function. Default is Sampler.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

DataContainer

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#

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_{n-1})(x) = \nabla_xf_{0}(x) + \nabla_xf_{1}(x) + ... + \nabla_xf_{n-1}(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

DataContainer

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

DataContainer

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\) or None if out.

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, …, N-1} 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, …, N-1}, 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: Sampler

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

Sampler

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 Herman-Meyer 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 SIRF-Based 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

Sampler

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, Herman-Mayer, 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

Sampler

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, …, N-1}, 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: Sampler

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

Sampler

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

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]

>>> 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, …, N-1}, 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 Herman-Meyer 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 SIRF-Based 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

Sampler

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, …, N-1} 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

Sampler

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.

Built-in 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.05e-5]

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:

\[\min \text{Regulariser} + \text{Fidelity}\]

The block framework consists of:

BlockDataContainer
BlockFunction
BlockOperator

The block framework allows writing more advanced optimisation problems. Consider the typical Tikhonov regularisation:

\[\underset{u}{\mathrm{argmin}}\begin{Vmatrix}A u - b \end{Vmatrix}^2_2 + \alpha^2\|Lu\|^2_2\]

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 re-written equivalently in the block matrix form:

\[\underset{u}{\mathrm{argmin}}\begin{Vmatrix}\binom{A}{\alpha L} u - \binom{b}{0}\end{Vmatrix}^2_2\]

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.

\[\underset{u}{\mathrm{argmin}}\begin{Vmatrix}\tilde{A} u - \tilde{b}\end{Vmatrix}^2_2\]

To be able to express our optimisation problems in the matrix form above, we developed the so-called, 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.

\[ \begin{align}\begin{aligned}x = [x_{1}, x_{2} ]\in (X_{1}\times X_{2})\\y = [y_{1}, y_{2}, y_{3} ]\in(Y_{1}\times Y_{2} \times Y_{3})\end{aligned}\end{align} \]

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 element-wise, 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 ]

__iter__()[source]#: BlockDataContainer is Iterable

next()[source]#: python2 backwards compatibility

is_compatible(other)[source]#: basic check if the size of the 2 objects fit

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 element-wise 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

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

copy()[source]#: alias of clone

__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/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

__rsub__(other)[source]#

Reverse subtraction

__rmul__(other)[source]#

Reverse multiplication

__rdiv__(other)[source]#

Reverse division

__rtruediv__(other)[source]#

Reverse truedivision

__rpow__(other)[source]#

Reverse power

__iadd__(other)[source]#: Inline addition

__isub__(other)[source]#: Inline subtraction

__imul__(other)[source]#: Inline multiplication

__idiv__(other)[source]#: Inline division

__itruediv__(other)[source]#: Inline truedivision

__neg__()[source]#: Return - self

__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} \]

\[ \begin{align}\begin{aligned}\begin{split}Y = \begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \end{bmatrix}, \quad F = [ f_{1}, f_{2}, f_{3} ]\end{split}\\F(Y) : = f_{1}(y_{1}) + f_{2}(y_{2}) + f_{3}(y_{3})\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}\]

__init__(*functions)[source]#

property L#

Lipschitz of the gradient of function f.

L is positive real number, such that \(\|f'(x) - f'(y)\| \leq L\|x-y\|\), 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

\[\begin{split}K = \begin{bmatrix} A_{1} & A_{2} \\ A_{3} & A_{4} \\ A_{5} & A_{6} \end{bmatrix}_{(3,2)} * \quad \underbrace{\begin{bmatrix} x_{1} \\ x_{2} \end{bmatrix}_{(2,1)}}_{\textbf{x}} = \begin{bmatrix} A_{1}x_{1} + A_{2}x_{2}\\ A_{3}x_{1} + A_{4}x_{2}\\ A_{5}x_{1} + A_{6}x_{2}\\ \end{bmatrix}_{(3,1)} = \begin{bmatrix} y_{1}\\ y_{2}\\ y_{3} \end{bmatrix}_{(3,1)} = \textbf{y}\end{split}\]

Column: Share the same domains \(X_{1}, X_{2}\)

Rows: Share the same ranges \(Y_{1}, Y_{2}, Y_{3}\)

\[K : (X_{1}\times X_{2}) \rightarrow (Y_{1}\times Y_{2} \times 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}\)

\[\begin{split}A_{1}: X_{1} \rightarrow Y_{1} \\ A_{3}: X_{1} \rightarrow Y_{2} \\ A_{5}: X_{1} \rightarrow Y_{3} \\ A_{2}: X_{2} \rightarrow Y_{1} \\ A_{4}: X_{2} \rightarrow Y_{2} \\ A_{6}: X_{2} \rightarrow Y_{3}\end{split}\]

For instance with these ingredients one may write the following objective function,

\[\alpha ||\nabla u||_{2,1} + ||u - g||_2^2\]

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 \(|| x-g ||^2_2\) is least squares fidelity function as

\[ \begin{align}\begin{aligned}\begin{split}K = \begin{bmatrix} \nabla \\ \mathbb{1} \end{bmatrix}\end{split}\\F(x) = \Big[ \alpha \lVert ~x~ \rVert_{2,1} ~~ , ~~ || x - g||_2^2 \Big]\\w = [ u ]\end{aligned}\end{align} \]

Then we have rewritten the problem as

\[F(Kw) = \alpha \left\lVert \nabla u \right\rVert_{2,1} + ||u-g||^2_2\]

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 row-by-row 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 column-wise and the same range row-wise.

Examples

BlockOperator(op0,op1) results in a row block

BlockOperator(op0,op1,shape=(1,2)) results in a column block

__init__(*args, **kwargs)[source]#

column_wise_compatible()[source]#: Operators in a Block should have the same domain per column

row_wise_compatible()[source]#: Operators in a Block should have the same range per row

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

is_linear()[source]#: Returns whether all the elements of the BlockOperator are 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 row-by-row fashion

domain_geometry()[source]#

Returns the domain of the BlockOperator

If the shape of the BlockOperator is (N,1) the domain is a ImageGeometry or AcquisitionGeometry. Otherwise it is a BlockGeometry.

range_geometry()[source]#: Returns the range of the BlockOperator

__getitem__(index)[source]#: Returns the index-th operator in the block irrespectively of it’s shape

get_as_list()[source]#: Returns the list of operators

Return Home

References#

1(1,2,3): Amir Beck and Marc Teboulle. A fast iterative shrinkage-thresholding 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.
2(1,2,3): Amir Beck and Marc Teboulle. Fast gradient-based 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.
3(1,2): Antonin Chambolle and Thomas Pock. A first-order primal-dual 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/s10851-010-0251-1, doi:10.1007/s10851-010-0251-1.
4: Ernie Esser, Xiaoqun Zhang, and Tony F. Chan. A general framework for a class of first order primal-dual 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.
5: Osman Güler. New proximal point algorithms for convex minimization. SIAM Journal on Optimization, 2(4):649–664, 1992.
6: 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.
7: 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.
8: Yurii Nesterov. Introductory lectures on convex optimization: A basic course. Volume 87. Springer Science & Business Media, 2003.
9: 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.
10: Julian Rasch and Antonin Chambolle. Inexact first-order primal–dual algorithms. Computational Optimization and Applications, 76(2):381–430, Jun 2020. URL: https://doi.org/10.1007/s10589-020-00186-y, doi:10.1007/s10589-020-00186-y.
11: Mingqiang Zhu, Stephen J. Wright, and Tony F. Chan. Duality-based algorithms for total-variation-regularized image restoration. Computational Optimization and Applications, 47(3):377–400, Nov 2010. URL: https://doi.org/10.1007/s10589-008-9225-2, doi:10.1007/s10589-008-9225-2.

Table of Contents

Optimisation framework#

Algorithms (Deterministic)#

Base class#

GD#

CGLS#

SIRT#

ISTA#

FISTA#

PDHG#

LADMM#

Algorithms (Stochastic)#

SPDHG#

Approximate gradient methods#

Stochastic Gradient Descent Example#

Operators#

Operator base classes#

Trivial operators#

GradientOperator#

WaveletOperator#

Functions#

Base classes#

Simple functions#

Composition of operator and a function#

Indicator box#

KullbackLeibler#

L1 Norm#

L2 Norm Squared#

Least Squares#

L1 Sparsity#

Mixed L21 norm#

Smooth Mixed L21 norm#

Mixed L11 norm#

Total variation#

Approximate Gradient base class#

Stochastic Gradient function#

Utilities#

Samplers#

Callbacks#

Block Framework#

BlockDataContainer#

Block Function#

Block Operator#

References#