sporco.admm.ccmodmd

ADMM algorithms for the Convolutional Constrained MOD problem with Mask Decoupling

Functions

ConvCnstrMODMaskDcpl(*args, **kwargs)	A wrapper function that dynamically defines a class derived from one of the implementations of the Convolutional Constrained MOD with Mask Decoupling problems, and returns an object instantiated with the provided.
ConvCnstrMODMaskDcplOptions([opt, method])	A wrapper function that dynamically defines a class derived from the Options class associated with one of the implementations of the Convolutional Constrained MOD with Mask Decoupling problem, and returns an object instantiated with the provided parameters.

Classes

ConvCnstrMODMaskDcplBase(*args, **kwargs)	Base class for ADMM algorithms for Convolutional Constrained MOD
ConvCnstrMODMaskDcpl_IterSM(*args, **kwargs)	ADMM algorithm for Convolutional Constrained MOD with Mask Decoupling [28] with the \(\mathbf{x}\) step solved via iterated application of the Sherman-Morrison equation [53].
ConvCnstrMODMaskDcpl_CG(*args, **kwargs)	ADMM algorithm for Convolutional Constrained MOD with Mask Decoupling [28] with the \(\mathbf{x}\) step solved via Conjugate Gradient (CG) [53].
ConvCnstrMODMaskDcpl_Consensus(*args, **kwargs)	Hybrid ADMM Consensus algorithm for Convolutional Constrained MOD with Mask Decoupling [26].

Function Descriptions

sporco.admm.ccmodmd.ConvCnstrMODMaskDcpl(*args, **kwargs)

A wrapper function that dynamically defines a class derived from one of the implementations of the Convolutional Constrained MOD with Mask Decoupling problems, and returns an object instantiated with the provided. parameters. The wrapper is designed to allow the appropriate object to be created by calling this function using the same syntax as would be used if it were a class. The specific implementation is selected by use of an additional keyword argument ‘method’. Valid values are:

• 'ism' : Use the implementation defined in ConvCnstrMODMaskDcpl_IterSM. This method works well for a small number of training images, but is very slow for larger training sets.

• 'cg' : Use the implementation defined in ConvCnstrMODMaskDcpl_CG. This method is slower than 'ism' for small training sets, but has better run time scaling as the training set grows.

• 'cns' : Use the implementation defined in ConvCnstrMODMaskDcpl_Consensus. This method is the best choice for large training sets.

The default value is 'cns'.

sporco.admm.ccmodmd.ConvCnstrMODMaskDcplOptions(opt=None, method='cns')

A wrapper function that dynamically defines a class derived from the Options class associated with one of the implementations of the Convolutional Constrained MOD with Mask Decoupling problem, and returns an object instantiated with the provided parameters. The wrapper is designed to allow the appropriate object to be created by calling this function using the same syntax as would be used if it were a class. The specific implementation is selected by use of an additional keyword argument ‘method’. Valid values are as specified in the documentation for ConvCnstrMODMaskDcpl.

Class Descriptions

class sporco.admm.ccmodmd.ConvCnstrMODMaskDcplBase(*args, **kwargs)

Bases: ADMMTwoBlockCnstrnt

Base class for ADMM algorithms for Convolutional Constrained MOD with Mask Decoupling [28].

Solve the optimisation problem

\[\mathrm{argmin}_\mathbf{d} \; (1/2) \left\| W \left(\sum_m \mathbf{d}_m * \mathbf{x}_m - \mathbf{s}\right) \right\|_2^2 \quad \text{such that} \quad \mathbf{d}_m \in C \;\; \forall m\]

where \(C\) is the feasible set consisting of filters with unit norm and constrained support, and \(W\) is a mask array, via the ADMM problem

\[\begin{split}\mathrm{argmin}_{\mathbf{d},\mathbf{g}_0,\mathbf{g}_1} \; (1/2) \| W \mathbf{g}_0 \|_2^2 + \iota_C(\mathbf{g}_1) \;\text{such that}\; \left( \begin{array}{c} X \\ I \end{array} \right) \mathbf{d} - \left( \begin{array}{c} \mathbf{g}_0 \\ \mathbf{g}_1 \end{array} \right) = \left( \begin{array}{c} \mathbf{s} \\ \mathbf{0} \end{array} \right) \;\;,\end{split}\]

where \(\iota_C(\cdot)\) is the indicator function of feasible set \(C\), and \(X \mathbf{d} = \sum_m \mathbf{x}_m * \mathbf{d}_m\).

The implementation of this class is substantially complicated by the support of multi-channel signals. In the following, the number of channels in the signal and dictionary are denoted by C and Cd respectively, the number of signals and the number of filters are denoted by K and M respectively, X, Z, and S denote the dictionary, coefficient map, and signal arrays respectively, and Y0 and Y1 denote blocks 0 and 1 of the auxiliary (split) variable of the ADMM problem. We need to consider three different cases:

1. Single channel signal and dictionary (C = Cd = 1)

2. Multi-channel signal, single channel dictionary (C > 1, Cd = 1)

3. Multi-channel signal and dictionary (C = Cd > 1)

The final three (non-spatial) dimensions of the main variables in each of these cases are as in the following table:

Var.	C = Cd = 1	C > 1, Cd = 1	C = Cd > 1
X	1 x 1 x M	1 x 1 x M	Cd x 1 x M
Z	1 x K x M	C x K x M	1 x K x M
S	1 x K x 1	C x K x 1	C x K x 1
Y0	1 x K x 1	C x K x 1	C x K x 1
Y1	1 x 1 x M	1 x 1 x M	C x 1 x M

In order to combine the block components Y0 and Y1 of variable Y into a single array, we need to be able to concatenate the two component arrays on one of the axes, but the shapes Y0 and Y1 are not compatible for concatenation. The solution for cases 1. and 3. is to swap the K and M axes of Y0` before concatenating, as well as after extracting the Y0 component from the concatenated Y variable. In case 2., since the C and K indices have the same behaviour in the dictionary update equation, we combine these axes in __init__, so that the case 2. array shapes become

Var.	C > 1, Cd = 1
X	1 x 1 x M
Z	1 x C K x M
S	1 x C K x 1
Y0	1 x C K x 1
Y1	1 x 1 x M

making it possible to concatenate Y0 and Y1 using the same axis swapping strategy as in the other cases. See block_sep0 and block_cat for additional details.

After termination of the solve method, attribute itstat is a list of tuples representing statistics of each iteration. The fields of the named tuple IterationStats are:

Iter : Iteration number

DFid : Value of data fidelity term \((1/2) \sum_k \| W (\sum_m \mathbf{d}_m * \mathbf{x}_{k,m} - \mathbf{s}_k) \|_2^2\)

Cnstr : Constraint violation measure

PrimalRsdl : Norm of primal residual

DualRsdl : Norm of dual residual

EpsPrimal : Primal residual stopping tolerance \(\epsilon_{\mathrm{pri}}\)

EpsDual : Dual residual stopping tolerance \(\epsilon_{\mathrm{dua}}\)

Rho : Penalty parameter

XSlvRelRes : Relative residual of X step solver

Time : Cumulative run time

Parameters:

Zarray_like

Coefficient map array

Sarray_like

Signal array

Warray_like

Mask array. The array shape must be such that the array is compatible for multiplication with the internal shape of input array S (see cnvrep.CDU_ConvRepIndexing for a discussion of the distinction between external and internal data layouts) after reshaping to the shape determined by cnvrep.mskWshape.

dsztuple

Filter support size(s)

optConvCnstrMODMaskDcplBase.Options object

Algorithm options

dimK0, 1, or None, optional (default None)

Number of dimensions in input signal corresponding to multiple independent signals

dimNint, optional (default 2)

Number of spatial dimensions

class Options(opt=None)

Bases: Options

ConvCnstrMODMaskDcplBase algorithm options

Options include all of those defined in ADMMTwoBlockCnstrnt.Options, together with additional options:

LinSolveCheck : Flag indicating whether to compute relative residual of X step solver.

ZeroMean : Flag indicating whether the solution dictionary \(\{\mathbf{d}_m\}\) should have zero-mean components.

Parameters:

optdict or None, optional (default None)

ConvCnstrMODMaskDcpl algorithm options

defaults = {'AbsStopTol': 0.0, 'AutoRho': {'AutoScaling': False, 'Enabled': False, 'Period': 10, 'RsdlRatio': 10.0, 'RsdlTarget': None, 'Scaling': 2.0, 'StdResiduals': False}, 'AuxVarObj': False, 'Callback': None, 'DataType': None, 'FastSolve': False, 'IterTimer': 'solve', 'LinSolveCheck': False, 'MaxMainIter': 1000, 'RelStopTol': 0.001, 'RelaxParam': 1.8, 'ReturnVar': 'Y1', 'ReturnX': True, 'StatusHeader': True, 'U0': None, 'Verbose': False, 'Y0': None, 'ZeroMean': False, 'fEvalX': True, 'gEvalY': False, 'rho': 1.0}

Default content and allowed dict keys

itstat_fields_objfn = ('DFid', 'Cnstr')

Fields in IterationStats associated with the objective function; see eval_objfn

itstat_fields_extra = ('XSlvRelRes',)

Non-standard fields in IterationStats; see itstat_extra

hdrtxt_objfn = ('DFid', 'Cnstr')

Display column headers associated with the objective function; see eval_objfn

hdrval_objfun = {'Cnstr': 'Cnstr', 'DFid': 'DFid'}

Dictionary mapping display column headers in hdrtxt_objfn to IterationStats entries

uinit(ushape)

Return initialiser for working variable U.

setcoef(Z)

Set coefficient array.

getdict(crop=True)

Get final dictionary. If crop is True, apply cnvrep.bcrop to returned array.

xstep_check(b)

Check the minimisation of the Augmented Lagrangian with respect to \(\mathbf{x}\) by method xstep defined in derived classes. This method should be called at the end of any xstep method.

ystep()

Minimise Augmented Lagrangian with respect to \(\mathbf{y}\).

relax_AX()

Implement relaxation if option RelaxParam != 1.0.

block_sep0(Y)

Separate variable into component corresponding to \(\mathbf{y}_0\) in \(\mathbf{y}\;\;\). The method from parent class ADMMTwoBlockCnstrnt is overridden here to allow swapping of K (multi-image) and M (filter) axes in block 0 so that it can be concatenated on axis M with block 1. This is necessary because block 0 has the dimensions of S while block 1 has the dimensions of D. Handling of multi-channel signals substantially complicate this issue. There are two multi-channel cases: multi-channel dictionary and signal (Cd = C > 1), and single-channel dictionary with multi-channel signal (Cd = 1, C > 1). In the former case, S and D shapes are (N x C x K x 1) and (N x C x 1 x M) respectively. In the latter case, __init__ has already taken care of combining C (multi-channel) and K (multi-image) axes in S, so the S and D shapes are (N x 1 x C K x 1) and (N x 1 x 1 x M) respectively.

block_cat(Y0, Y1)

Concatenate components corresponding to \(\mathbf{y}_0\) and \(\mathbf{y}_1\) to form \(\mathbf{y}\;\;\). The method from parent class ADMMTwoBlockCnstrnt is overridden here to allow swapping of K (multi-image) and M (filter) axes in block 0 so that it can be concatenated on axis M with block 1. This is necessary because block 0 has the dimensions of S while block 1 has the dimensions of D. Handling of multi-channel signals substantially complicate this issue. There are two multi-channel cases: multi-channel dictionary and signal (Cd = C > 1), and single-channel dictionary with multi-channel signal (Cd = 1, C > 1). In the former case, S and D shapes are (N x C x K x 1) and (N x C x 1 x M) respectively. In the latter case, __init__ has already taken care of combining C (multi-channel) and K (multi-image) axes in S, so the S and D shapes are (N x 1 x C K x 1) and (N x 1 x 1 x M) respectively.

cnst_A(X, Xf=None)

Compute \(A \mathbf{x}\) component of ADMM problem constraint.

obfn_g0var()

Variable to be evaluated in computing ADMMTwoBlockCnstrnt.obfn_g0, depending on the AuxVarObj option value.

cnst_A0(X, Xf=None)

Compute \(A_0 \mathbf{x}\) component of ADMM problem constraint.

cnst_A0T(Y0)

Compute \(A_0^T \mathbf{y}_0\) component of \(A^T \mathbf{y}\) (see ADMMTwoBlockCnstrnt.cnst_AT).

cnst_c0()

Compute constant component \(\mathbf{c}_0\) of \(\mathbf{c}\) in the ADMM problem constraint.

eval_objfn()

Compute components of regularisation function as well as total contribution to objective function.

obfn_g0(Y0)

Compute \(g_0(\mathbf{y}_0)\) component of ADMM objective function.

obfn_g1(Y1)

Compute \(g_1(\mathbf{y_1})\) component of ADMM objective function.

itstat_extra()

Non-standard entries for the iteration stats record tuple.

reconstruct(D=None)

Reconstruct representation.

rsdl_s(Yprev, Y)

Compute dual residual vector.

rsdl_sn(U)

Compute dual residual normalisation term.

class sporco.admm.ccmodmd.ConvCnstrMODMaskDcpl_IterSM(*args, **kwargs)

Bases: ConvCnstrMODMaskDcplBase

ADMM algorithm for Convolutional Constrained MOD with Mask Decoupling [28] with the \(\mathbf{x}\) step solved via iterated application of the Sherman-Morrison equation [53].

Multi-channel signals/images are supported [51]. See ConvCnstrMODMaskDcplBase for interface details.

Call graph

class Options(opt=None)

Bases: Options

ConvCnstrMODMaskDcpl_IterSM algorithm options

Options are the same as those defined in ConvCnstrMODMaskDcplBase.Options.

Parameters:

optdict or None, optional (default None)

ConvCnstrMODMaskDcpl_IterSM algorithm options

defaults = {'AbsStopTol': 0.0, 'AutoRho': {'AutoScaling': False, 'Enabled': False, 'Period': 10, 'RsdlRatio': 10.0, 'RsdlTarget': None, 'Scaling': 2.0, 'StdResiduals': False}, 'AuxVarObj': False, 'Callback': None, 'DataType': None, 'FastSolve': False, 'IterTimer': 'solve', 'LinSolveCheck': False, 'MaxMainIter': 1000, 'RelStopTol': 0.001, 'RelaxParam': 1.8, 'ReturnVar': 'Y1', 'ReturnX': True, 'StatusHeader': True, 'U0': None, 'Verbose': False, 'Y0': None, 'ZeroMean': False, 'fEvalX': True, 'gEvalY': False, 'rho': 1.0}

Default content and allowed dict keys

xstep()

Minimise Augmented Lagrangian with respect to \(\mathbf{x}\).

solve()

Call graph

class sporco.admm.ccmodmd.ConvCnstrMODMaskDcpl_CG(*args, **kwargs)

Bases: ConvCnstrMODMaskDcplBase

ADMM algorithm for Convolutional Constrained MOD with Mask Decoupling [28] with the \(\mathbf{x}\) step solved via Conjugate Gradient (CG) [53].

Multi-channel signals/images are supported [51]. See ConvCnstrMODMaskDcplBase for interface details.

Call graph

class Options(opt=None)

Bases: Options

ConvCnstrMODMaskDcpl_CG algorithm options

Options include all of those defined in ConvCnstrMODMaskDcplBase.Options, together with additional options:

CG : CG solver options

MaxIter : Maximum CG iterations.

StopTol : CG stopping tolerance.

Parameters:

optdict or None, optional (default None)

ConvCnstrMODMaskDcpl_CG algorithm options

defaults = {'AbsStopTol': 0.0, 'AutoRho': {'AutoScaling': False, 'Enabled': False, 'Period': 10, 'RsdlRatio': 10.0, 'RsdlTarget': None, 'Scaling': 2.0, 'StdResiduals': False}, 'AuxVarObj': False, 'CG': {'MaxIter': 1000, 'StopTol': 0.001}, 'Callback': None, 'DataType': None, 'FastSolve': False, 'IterTimer': 'solve', 'LinSolveCheck': False, 'MaxMainIter': 1000, 'RelStopTol': 0.001, 'RelaxParam': 1.8, 'ReturnVar': 'Y1', 'ReturnX': True, 'StatusHeader': True, 'U0': None, 'Verbose': False, 'Y0': None, 'ZeroMean': False, 'fEvalX': True, 'gEvalY': False, 'rho': 1.0}

Default content and allowed dict keys

itstat_fields_extra = ('XSlvRelRes', 'XSlvCGIt')

Non-standard fields in IterationStats; see itstat_extra

xstep()

Minimise Augmented Lagrangian with respect to \(\mathbf{x}\).

itstat_extra()

Non-standard entries for the iteration stats record tuple.

solve()

Call graph

class sporco.admm.ccmodmd.ConvCnstrMODMaskDcpl_Consensus(*args, **kwargs)

Bases: ConvCnstrMOD_Consensus

Hybrid ADMM Consensus algorithm for Convolutional Constrained MOD with Mask Decoupling [26].

Solve the optimisation problem

\[\mathrm{argmin}_\mathbf{d} \; (1/2) \left\| W \left(\sum_m \mathbf{d}_m * \mathbf{x}_m - \mathbf{s} \right) \right\|_2^2 \quad \text{such that} \quad \mathbf{d}_m \in C \;\; \forall m\]

where \(C\) is the feasible set consisting of filters with unit norm and constrained support, and \(W\) is a mask array, via a hybrid ADMM Consensus problem.

See the documentation of ConvCnstrMODMaskDcplBase for a detailed discussion of the implementational complications resulting from the support of multi-channel signals.

Call graph

Parameters:

Zarray_like

Coefficient map array

Sarray_like

Signal array

Warray_like

Mask array. The array shape must be such that the array is compatible for multiplication with input array S (see cnvrep.mskWshape for more details).

dsztuple

Filter support size(s)

optConvCnstrMOD_Consensus.Options object

Algorithm options

dimK0, 1, or None, optional (default None)

Number of dimensions in input signal corresponding to multiple independent signals

dimNint, optional (default 2)

Number of spatial dimensions

setcoef(Z)

Set coefficient array.

var_y1()

Get the auxiliary variable that is constrained to be equal to the dictionary. The method is named for compatibility with the method of the same name in ConvCnstrMODMaskDcpl_IterSM and ConvCnstrMODMaskDcpl_CG (it is not variable Y1 in this class).

relax_AX()

The parent class method that this method overrides only implements the relaxation step for the variables of the baseline consensus algorithm. This method calls the overridden method and then implements the relaxation step for the additional variables required for the mask decoupling modification to the baseline algorithm.

xstep()

The xstep of the baseline consensus class from which this class is derived is re-used to implement the xstep of the modified algorithm by replacing self.ZSf, which is constant in the baseline algorithm, with a quantity derived from the additional variables self.Y1 and self.U1. It is also necessary to set the penalty parameter to unity for the duration of the x step.

ystep()

The parent class ystep method is overridden to allow also performing the ystep for the additional variables introduced in the modification to the baseline algorithm.

ustep()

The parent class ystep method is overridden to allow also performing the ystep for the additional variables introduced in the modification to the baseline algorithm.

obfn_dfd()

Compute data fidelity term \((1/2) \| W \left( \sum_m \mathbf{d}_m * \mathbf{x}_m - \mathbf{s} \right) \|_2^2\).

solve()

Call graph

compute_residuals()

Compute residuals and stopping thresholds. The parent class method is overridden to ensure that the residual calculations include the additional variables introduced in the modification to the baseline algorithm.