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 genericlly formulated optimisation problems constructed from operator and function objects. The fundamental components are: + :code:`Operator`: A class specifying a (currently linear) operator + :code:`Function`: A class specifying mathematical functions such as a least squares data fidelity. + :code:`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. 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: .. math:: \min \text{Regulariser} + \text{Fidelity} The block framework consists of: + BlockDataContainer + BlockFunction + BlockOperator `BlockDataContainer`_ holds `DataContainer` as column vector. It is possible to do basic algebra between ``BlockDataContainer`` s and with numbers, list or numpy arrays. `BlockFunction`_ acts on ``BlockDataContainer`` as a separable sum function: .. math:: f = [f_1,...,f_n] \newline f([x_1,...,x_n]) = f_1(x_1) + .... + f_n(x_n) `BlockOperator`_ represent a block matrix with operators .. math:: 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} Column: Share the same domains :math:`X_{1}, X_{2}` Rows: Share the same ranges :math:`Y_{1}, Y_{2}, Y_{3}` .. math:: K : (X_{1}\times X_{2}) \rightarrow (Y_{1}\times Y_{2} \times Y_{3}) :math:`A_{1}, A_{3}, A_{5}`: share the same domain :math:`X_{1}` and :math:`A_{2}, A_{4}, A_{6}`: share the same domain :math:`X_{2}` .. math:: 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} For instance with these ingredients one may write the following objective function, .. math:: \alpha ||\nabla u||_{2,1} + ||u - g||_2^2 where :math:`g` represent the measured values, :math:`u` the solution :math:`\nabla` is the gradient operator, :math:`|| ~~ ||_{2,1}` is a norm for the output of the gradient operator and :math:`|| x-g ||^2_2` is least squares fidelity function as .. math:: K = \begin{bmatrix} \nabla \\ \mathbb{1} \end{bmatrix} F(x) = \Big[ \alpha \lVert ~x~ \rVert_{2,1} ~~ , ~~ || x - g||_2^2 \Big] w = [ u ] Then we have rewritten the problem as .. math:: F(Kw) = \alpha \left\lVert \nabla u \right\rVert_{2,1} + ||u-g||^2_2 Which in Python would be like .. code-block:: python op1 = Gradient(ig, correlation=Gradient.CORRELATION_SPACE) op2 = Identity(ig, ag) # Create BlockOperator K = BlockOperator(op1, op2, shape=(2,1) ) # Create functions F = BlockFunction(alpha * MixedL21Norm(), 0.5 * L2NormSquared(b=noisy_data)) Algorithm ========= 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) and Fast Iterative Shrinkage Thresholding Algorithm (FISTA). An algorithm is designed for a particular generic optimisation problem accepts and number of :code:`Function`s and/or :code:`Operator`s 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 :code:`Algorithm` class. The user is required to implement only 4 methods: set_up, __init__, update and update_objective. + :code:`set_up` and :code:`__init__` are used to configure the algorithm + :code:`update` is the actual iteration updating the solution + :code:`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: .. code-block :: python 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 :code:`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. .. autoclass:: ccpi.optimisation.algorithms.Algorithm :members: :private-members: :special-members: .. autoclass:: ccpi.optimisation.algorithms.GradientDescent :members: .. autoclass:: ccpi.optimisation.algorithms.CGLS :members: .. autoclass:: ccpi.optimisation.algorithms.FISTA :members: :special-members: .. autoclass:: ccpi.optimisation.algorithms.PDHG :members: .. autoclass:: ccpi.optimisation.algorithms.SIRT :members: Operator ======== The two most important methods are :code:`direct` and :code:`adjoint` methods that describe the result of applying the operator, and its adjoint respectively, onto a compatible :code:`DataContainer` input. The output is another :code:`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 :code:`Operator` class. A special class is the :code:`LinearOperator` which represents an operator for which the :code:`adjoint` operation is defined. A :code:`ScaledOperator` represents the multiplication of any operator with a scalar. .. autoclass:: ccpi.optimisation.operators.Operator :members: :special-members: .. autoclass:: ccpi.optimisation.operators.LinearOperator :members: :special-members: .. autoclass:: ccpi.optimisation.operators.ScaledOperator :members: :special-members: Trivial operators ----------------- Trivial operators are the following. .. autoclass:: ccpi.optimisation.operators.Identity :members: :special-members: .. autoclass:: ccpi.optimisation.operators.ZeroOperator :members: :special-members: .. autoclass:: ccpi.optimisation.operators.LinearOperatorMatrix :members: :special-members: Gradient ----------------- In the following the required classes for the implementation of the :code:`Gradient` operator. .. autoclass:: ccpi.optimisation.operators.Gradient :members: :special-members: .. autoclass:: ccpi.optimisation.operators.FiniteDiff :members: :special-members: .. autoclass:: ccpi.optimisation.operators.SparseFiniteDiff :members: :special-members: .. autoclass:: ccpi.optimisation.operators.SymmetrizedGradient :members: :special-members: Shrinkage operator ------------------ .. autoclass:: ccpi.optimisation.operators.ShrinkageOperator :members: :special-members: Function ======== A :code:`Function` represents a mathematical function of one or more inputs and is intended to accept :code:`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 :code:`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 :code:`proximal` should return the proximal operator evaluated at an input point. The function value is evaluated by calling the function itself, e.g. :code:`f(x)` for a :code:`Function f` and input point :code:`x`. Base classes ------------ .. autoclass:: ccpi.optimisation.functions.Function :members: :special-members: .. autoclass:: ccpi.optimisation.functions.ScaledFunction :members: :special-members: .. autoclass:: ccpi.optimisation.functions.SumFunction :members: :special-members: Zero Function ------------- .. autoclass:: ccpi.optimisation.functions.ZeroFunction :members: :special-members: .. autoclass:: ccpi.optimisation.functions.SumFunctionScalar :members: :special-members: Constant Function ----------------- .. autoclass:: ccpi.optimisation.functions.ConstantFunction :members: :special-members: Composition of operator and a function -------------------------------------- This class allows the user to write a function which does the following: .. math:: F ( x ) = G ( Ax ) where :math:`A` is an operator. For instance the least squares function l2norm_ :code:`Norm2Sq` can be expressed as .. math:: F(x) = || Ax - b ||^2_2 .. code::python F1 = Norm2Sq(A, b) # or equivalently F2 = FunctionOperatorComposition(L2NormSquared(b=b), A) .. autoclass:: ccpi.optimisation.functions.FunctionOperatorComposition :members: :special-members: Indicator box ------------- .. autoclass:: ccpi.optimisation.functions.IndicatorBox :members: :special-members: KullbackLeibler --------------- .. autoclass:: ccpi.optimisation.functions.KullbackLeibler :members: :special-members: L1 Norm ------- .. autoclass:: ccpi.optimisation.functions.L1Norm :members: :special-members: Squared L2 norm --------------- .. l2norm: .. autoclass:: ccpi.optimisation.functions.L2NormSquared :members: :special-members: Least Squares ------------- .. autoclass:: ccpi.optimisation.functions.LeastSquares :members: :special-members: Mixed L21 norm -------------- .. autoclass:: ccpi.optimisation.functions.MixedL21Norm :members: :special-members: Smooth Mixed L21 norm --------------------- .. autoclass:: ccpi.optimisation.functions.smoothMixedL21Norm :members: :special-members: Block Framework *************** Block Operator Block Function ============== A Block vector of functions, Size of vector coincides with the rows of :math:`K`: .. math:: Kx = \begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \end{bmatrix}, \quad f = [ f_{1}, f_{2}, f_{3} ] f(Kx) : = f_{1}(y_{1}) + f_{2}(y_{2}) + f_{3}(y_{3}) .. autoclass:: ccpi.optimisation.functions.BlockFunction :members: :special-members: Block DataContainer ============== .. math:: 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}) .. autoclass:: ccpi.framework.BlockDataContainer :members: :special-members: :ref:`Return Home ` .. _BlockDataContainer: framework.html#ccpi.framework.BlockDataContainer .. _BlockFunction: optimisation.html#ccpi.optimisation.functions.BlockFunction .. _BlockOperator: optimisation.html#ccpi.optimisation.operators.BlockOperators