sporco.admm.tvl2

Classes for ADMM algorithms for Total Variation (TV) optimisation with an \(\ell_2\) data fidelity term

Classes

TVL2Denoise(*args, **kwargs)	ADMM algorithm for \(\ell_2\)-TV denoising problem [41], [27], [8].
TVL2Deconv(*args, **kwargs)	ADMM algorithm for \(\ell_2\)-TV deconvolution problem.

Class Descriptions

class sporco.admm.tvl2.TVL2Denoise(*args, **kwargs)

Bases: ADMM

ADMM algorithm for \(\ell_2\)-TV denoising problem [41], [27], [8].

Solve the optimisation problem

\[\mathrm{argmin}_\mathbf{x} \; (1/2) \| W_{\mathrm{df}}(\mathbf{x} - \mathbf{s}) \|_2^2 + \lambda \left\| W_{\mathrm{tv}} \sqrt{(G_r \mathbf{x})^2 + (G_c \mathbf{x})^2} \right\|_1\]

via the ADMM problem

\[\begin{split}\mathrm{argmin}_{\mathbf{x},\mathbf{y}_r,\mathbf{y}_c} \; (1/2) \| W_{\mathrm{df}}(\mathbf{x} - \mathbf{s}) \|_2^2 + \lambda \left\| W_{\mathrm{tv}} \sqrt{(\mathbf{y}_r)^2 + (\mathbf{y}_c)^2} \right\|_1 \;\text{such that}\; \left( \begin{array}{c} G_r \\ G_c \end{array} \right) \mathbf{x} - \left( \begin{array}{c} \mathbf{y}_r \\ \mathbf{y}_c \end{array} \right) = \left( \begin{array}{c} \mathbf{0} \\ \mathbf{0} \end{array} \right) \;\;,\end{split}\]

where \(G_r\) and \(G_c\) are gradient operators along array rows and columns respectively, and \(W_{\mathrm{df}}\) and \(W_{\mathrm{tv}}\) are diagonal weighting matrices.

While these equations describe the default behaviour of regularisation in two dimensions, this class supports an arbitrary number of dimensions. For example, for 3D TV regularisation in a 3D array, the object should be initialised with parameter axes set to (0, 1, 2).

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

ObjFun : Objective function value

DFid : Value of data fidelity term \((1/2) \| W_{\mathrm{df}} (\mathbf{x} - \mathbf{s}) \|_2^2\)

RegTV : Value of regularisation term \(\| W_{\mathrm{tv}} \sqrt{(G_r \mathbf{x})^2 + (G_c \mathbf{x})^2} \|_1\)

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

GSIter : Number of Gauss-Seidel iterations

GSRelRes : Relative residual of Gauss-Seidel solution

Time : Cumulative run time

Call graph

Parameters:

Sarray_like

Signal vector or matrix

lmbdafloat

Regularisation parameter

optTVL2Denoise.Options object

Algorithm options

axestuple, optional (default (0, 1))

Axes on which TV regularisation is to be applied

caxisint or None, optional (default None)

Axis on which channels of a multi-channel image are stacked. If None, TV regularisation is applied indepdendently to each channel, otherwise Vector TV [8] regularisation is applied jointly to all channels.

class Options(opt=None)

Bases: Options

TVL2Denoise algorithm options

Options include all of those defined in sporco.admm.admm.ADMM.Options, together with additional options:

gEvalY : Flag indicating whether the \(g\) component of the objective function should be evaluated using variable Y (True) or X (False) as its argument.

MaxGSIter : Maximum Gauss-Seidel iterations.

GSTol : Gauss-Seidel stopping tolerance.

DFidWeight : Data fidelity weight matrix.

TVWeight : TV term weight matrix.

Parameters:

optdict or None, optional (default None)

TVL2Denoise algorithm options

defaults = {'AbsStopTol': 0.0, 'AutoRho': {'AutoScaling': True, 'Enabled': False, 'Period': 1, 'RsdlRatio': 1.2, 'RsdlTarget': None, 'Scaling': 1000.0, 'StdResiduals': False}, 'Callback': None, 'DFidWeight': 1.0, 'DataType': None, 'FastSolve': False, 'GSTol': 0.0, 'IterTimer': 'solve', 'MaxGSIter': 2, 'MaxMainIter': 1000, 'RelStopTol': 0.001, 'RelaxParam': 1.8, 'StatusHeader': True, 'TVWeight': 1.0, 'U0': None, 'Verbose': False, 'Y0': None, 'gEvalY': True, 'rho': None}

Default content and allowed dict keys

itstat_fields_objfn = ('ObjFun', 'DFid', 'RegTV')

Fields in IterationStats associated with the objective function; see eval_objfn

itstat_fields_extra = ('GSIter', 'GSRelRes')

Non-standard fields in IterationStats; see itstat_extra

hdrtxt_objfn = ('Fnc', 'DFid', 'RegTV')

Display column headers associated with the objective function; see eval_objfn

hdrval_objfun = {'DFid': 'DFid', 'Fnc': 'ObjFun', 'RegTV': 'RegTV'}

Dictionary mapping display column headers in hdrtxt_objfn to IterationStats entries

uinit(ushape)

Return initialiser for working variable U.

xstep()

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

ystep()

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

obfn_gvar()

Variable to be evaluated in computing regularisation term, depending on ‘gEvalY’ option value.

eval_objfn()

Compute components of objective function as well as total contribution to objective function. Data fidelity term is \((1/2) \| \mathbf{x} - \mathbf{s} \|_2^2\) and regularisation term is \(\| W_{\mathrm{tv}} \sqrt{(G_r \mathbf{x})^2 + (G_c \mathbf{x})^2}\|_1\).

itstat_extra()

Non-standard entries for the iteration stats record tuple.

cnst_A(X)

Compute \(A \mathbf{x}\) component of ADMM problem constraint. In this case \(A \mathbf{x} = (G_r^T \;\; G_c^T)^T \mathbf{x}\).

cnst_AT(X)

Compute \(A^T \mathbf{x}\) where \(A \mathbf{x}\) is a component of ADMM problem constraint. In this case \(A^T \mathbf{x} = (G_r^T \;\; G_c^T) \mathbf{x}\).

cnst_B(Y)

Compute \(B \mathbf{y}\) component of ADMM problem constraint. In this case \(B \mathbf{y} = -\mathbf{y}\).

cnst_c()

Compute constant component \(\mathbf{c}\) of ADMM problem constraint. In this case \(\mathbf{c} = \mathbf{0}\).

LaplaceCentreWeight()

Centre weighting matrix for TV Laplacian.

GaussSeidelStep(S, X, ATYU, rho, lcw, W2)

Gauss-Seidel step for linear system in TV problem.

solve()

Call graph

class sporco.admm.tvl2.TVL2Deconv(*args, **kwargs)

Bases: ADMM

ADMM algorithm for \(\ell_2\)-TV deconvolution problem.

Solve the optimisation problem

\[\mathrm{argmin}_\mathbf{x} \; (1/2) \| H \mathbf{x} - \mathbf{s} \|_2^2 + \lambda \left\| W_{\mathrm{tv}} \sqrt{(G_r \mathbf{x})^2 + (G_c \mathbf{x})^2} \right\|_1 \;\;,\]

where \(H\) denotes the linear operator corresponding to a convolution, \(G_r\) and \(G_c\) are gradient operators along array rows and columns respectively, and \(W_{\mathrm{df}}\) and \(W_{\mathrm{tv}}\) are diagonal weighting matrices, via the ADMM problem

\[\begin{split}\mathrm{argmin}_{\mathbf{x},\mathbf{y}_r,\mathbf{y}_c} \; (1/2) \| H \mathbf{x} - \mathbf{s} \|_2^2 + \lambda \left\| W_{\mathrm{tv}} \sqrt{(\mathbf{y}_r)^2 + (\mathbf{y}_c)^2} \right\|_1 \;\text{such that}\; \left( \begin{array}{c} G_r \\ G_c \end{array} \right) \mathbf{x} - \left( \begin{array}{c} \mathbf{y}_r \\ \mathbf{y}_c \end{array} \right) = \left( \begin{array}{c} \mathbf{0} \\ \mathbf{0} \end{array} \right) \;\;.\end{split}\]

While these equations describe the default behaviour of regularisation in two dimensions, this class supports an arbitrary number of dimensions. For example, for 3D TV regularisation in a 3D array, the object should be initialised with parameter axes set to (0, 1, 2).

Note that the convolution is implemented in the frequency domain, having the same phase offset as fftconv, which differs from that of scipy.ndimage.convolve with the default origin parameter.

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

ObjFun : Objective function value

DFid : Value of data fidelity term \((1/2) \| H \mathbf{x} - \mathbf{s} \|_2^2\)

RegTV : Value of regularisation term \(\| W_{\mathrm{tv}} \sqrt{(G_r \mathbf{x})^2 + (G_c \mathbf{x})^2} \|_1\)

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

Call graph

Parameters:

Aarray_like

Filter kernel corresponding to operator \(H\) above

Sarray_like

Signal vector or matrix

lmbdafloat

Regularisation parameter

optTVL2Deconv.Options object

Algorithm options

axestuple, optional (default (0, 1))

Axes on which TV regularisation is to be applied

caxisint or None, optional (default None)

Axis on which channels of a multi-channel image are stacked. If None, TV regularisation is applied indepdendently to each channel, otherwise Vector TV [8] regularisation is applied jointly to all channels.

class Options(opt=None)

Bases: Options

TVL2Deconv algorithm options

Options include all of those defined in sporco.admm.admm.ADMM.Options, together with additional options:

gEvalY : Flag indicating whether the \(g\) component of the objective function should be evaluated using variable Y (True) or X (False) as its argument.

LinSolveCheck : If True, compute relative residual of X step solver.

TVWeight : TV term weight matrix.

Parameters:

optdict or None, optional (default None)

TVL2Deconv algorithm options

defaults = {'AbsStopTol': 0.0, 'AutoRho': {'AutoScaling': True, 'Enabled': True, 'Period': 1, 'RsdlRatio': 1.2, 'RsdlTarget': None, 'Scaling': 1000.0, 'StdResiduals': False}, 'Callback': None, 'DataType': None, 'FastSolve': False, 'IterTimer': 'solve', 'LinSolveCheck': False, 'MaxMainIter': 1000, 'RelStopTol': 0.001, 'RelaxParam': 1.8, 'StatusHeader': True, 'TVWeight': 1.0, 'U0': None, 'Verbose': False, 'Y0': None, 'gEvalY': True, 'rho': None}

Default content and allowed dict keys

itstat_fields_objfn = ('ObjFun', 'DFid', 'RegTV')

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 = ('Fnc', 'DFid', 'RegTV')

Display column headers associated with the objective function; see eval_objfn

hdrval_objfun = {'DFid': 'DFid', 'Fnc': 'ObjFun', 'RegTV': 'RegTV'}

Dictionary mapping display column headers in hdrtxt_objfn to IterationStats entries

uinit(ushape)

Return initialiser for working variable U.

xstep()

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

ystep()

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

obfn_gvar()

Variable to be evaluated in computing regularisation term, depending on ‘gEvalY’ option value.

eval_objfn()

Compute components of objective function as well as total contribution to objective function. Data fidelity term is \((1/2) \| H \mathbf{x} - \mathbf{s} \|_2^2\) and regularisation term is \(\| W_{\mathrm{tv}} \sqrt{(G_r \mathbf{x})^2 + (G_c \mathbf{x})^2}\|_1\).

itstat_extra()

Non-standard entries for the iteration stats record tuple.

cnst_A(X, Xf=None)

Compute \(A \mathbf{x}\) component of ADMM problem constraint. In this case \(A \mathbf{x} = (G_r^T \;\; G_c^T)^T \mathbf{x}\).

solve()

Call graph

cnst_AT(X)

Compute \(A^T \mathbf{x}\) where \(A \mathbf{x}\) is a component of ADMM problem constraint. In this case \(A^T \mathbf{x} = (G_r^T \;\; G_c^T) \mathbf{x}\).

cnst_B(Y)

Compute \(B \mathbf{y}\) component of ADMM problem constraint. In this case \(B \mathbf{y} = -\mathbf{y}\).

cnst_c()

Compute constant component \(\mathbf{c}\) of ADMM problem constraint. In this case \(\mathbf{c} = \mathbf{0}\).