∎ 11institutetext: Birgit Komander22institutetext: Institute of Analysis and
Algebra, TU Braunschweig, 38092 Braunschweig, Germany, 22email: b.komander@tu-
braunschweig.de 33institutetext: D. A. Lorenz44institutetext: Institute of
Analysis and Algebra, TU Braunschweig, 38092 Braunschweig, Germany, 44email:
<EMAIL_ADDRESS>55institutetext: Lena Vestweber 66institutetext:
Institut Computational Mathematics, AG Numerik, TU Braunschweig, 38092
Braunschweig, Germany, 66email<EMAIL_ADDRESS>
# Denoising of image gradients and total generalized variation denoising
††thanks: This material was based upon work partially supported by the
National Science Foundation under Grant DMS-1127914 to the Statistical and
Applied Mathematical Sciences Institute. Any opinions, findings, and
conclusions or recommendations expressed in this material are those of the
author(s) and do not necessarily reflect the views of the National Science
Foundation.

B. Komander  D. A. Lorenz  L. Vestweber

###### Abstract

We revisit total variation denoising and study an augmented model where we
assume that an estimate of the image gradient is available. We show that this
increases the image reconstruction quality and derive that the resulting model
resembles the total generalized variation denoising method, thus providing a
new motivation for this model. Further, we propose to use a constraint
denoising model and develop a variational denoising model that is basically
parameter free, i.e. all model parameters are estimated directly from the
noisy image.

Moreover, we use Chambolle-Pock’s primal dual method as well as the Douglas-
Rachford method for the new models. For the latter one has to solve large
discretizations of partial differential equations. We propose to do this in an
inexact manner using the preconditioned conjugate gradients method and derive
preconditioners for this. Numerical experiments show that the resulting method
has good denoising properties and also that preconditioning does increase
convergence speed significantly. Finally we analyze the duality gap of
different formulations of the TGV denoising problem and derive a simple
stopping criterion.

###### Keywords:

image denoising, gradient estimate, total generalized variation, Douglas-
Rachford method, preconditioning

##  1 Introduction

In this work we revisit variational denoising of images with total variation
penalties, dating back to the classical Rudin-Osher-Fatemi total variation
denoising method rudin1992rofmodel . We start by augmenting the model with an
estimate of the image gradient and analyze, how this helps for image
denoising. This is related to the method of first estimating image normals and
then using this estimate for a better image denoising, an approach proposed by
Lysaker et al. Lysaker2004 . As we will see, a combined approach, which tries
to estimate the gradient of the denoised image and the denoised image itself
simultaneously, is very close to the successful total generalized variation
denoising from Bredies et al. bredies2010total . A brief introduction to this
idea was already proposed in komander2017denoising . Further, we will propose
different (in some sense equivalent) versions of the total general variation
denoising method (one of these, CTGV, already introduced in
komander2017denoising ) which have several advantages over the classical one:
First, we are going to work with constraints in contrast to penalties, which,
in some cases, allows for a simple, clean and effective parameter choice.
Second, different formulations of these problems lead to different dual
problems and hence, different algorithms and some of these turn out to be a
little simpler regarding duality gaps and stopping. Moreover, the different
models show slightly different numerical performance. Finally, we will make
use of the Douglas-Rachford method to solve these minimization problems. This
involves the solution of large discretizations of linear partial differential
equations and we will develop simple and effective preconditioners for these
equations. In contrast to bredies2015preconditioned  where the authors use
classical linear splitting methods for the inexact solution of the linear
equations, we propose to use a few iterations of the preconditioned conjugate
gradient method.

The paper is organized as follows. Section 2 motivates denoising of gradients
as a mean to improve total variation denoising and derives several new
variational methods. Then, section 3 investigates the corresponding duality
gaps of the problems and section 4 deals with the numerical treatment and
especially with the Douglas-Rachford method and efficient preconditioners for
the respective linear subproblems. Section 5 reports numerical experiments and
section 6 draws some conclusions.

##  2 Total variation denoising with estimates of the gradient

Since its introduction in 1992, the Rudin-Osher-Fatemi model rudin1992rofmodel
, also known as total variation denoising, has found numerous applications.
One way to put this model is that the total variation of an image is used as a
regularizer for an image denoising optimization problem, in general
$\min_{x}F(x)+G(Kx)$, with $u_{0}$ as the input image, defined on a domain
$\Omega$ as

$\displaystyle\begin{split}&\min_{u}\lambda\int_{\Omega}|\operatorname{\nabla}u|\,\mathrm{d}x+\frac{1}{2}\int_{\Omega}|u(x)-u_{0}(x)|^{2}\,\mathrm{d}x\\\
=&\min_{u}\lambda{\|{|\operatorname{\nabla}u|}\|}_{1}+\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}.\end{split}$
(1)

One problem in the resulting denoised images is the occurring staircasing
effect, i.e. the creation of flat areas separated by jumps. One way to
overcome this staircasing, proposed by Lysaker et al. Lysaker2004 , is an
image denoising technique in two steps. There, in a first step, a total
variation filter was used to smooth the normal vectors of the level sets of a
given noisy image and then, as a second step, a surface was fitted to the
resulting normal vectors. The method was formulated in a dynamic way, i.e. by
solving a certain partial differential equation to steady state. A similar
approach has been taken in komander2014  for a problem of deflectometric
surface measurements where the measurement device does not only produce
approximate point coordinates but also approximate surface normals. It turned
out that the incorporation of the surface normals results in an effective, but
fairly complicated and non-linear problem. Switching from surface normals to
image gradients, however, turns the problem into a “more linear” one and leads
to an equally effective method, see komander2014 .

In this section we follow the idea of introducing additional information, i.e.
gradient information, into the ROF-model (1) in order to prevent or reduce the
staircasing effect.

###  2.1 Denoising with prior knowledge on the gradient

Consider the image model $u_{0}=u^{\dagger}+\eta$, where $u_{0}$ is the given
noisy image, $u^{\dagger}$ is the ground truth, i.e. the noise-free image, and
$\eta$ is the additional Gaussian white noise. In the situation of images,
there are methods to obtain a reasonable estimate of the amount of noise, i.e.
an estimate on ${\|{u^{\dagger}-u_{0}}\|}_{2}={\|{\eta}\|}_{2}$ is available.
One can use, for example the techniques from liu2012noise ; liu2013single  to
estimate the noise level of Gaussian white noise quite accurately from a
single image. Using this information, it seems that

$\displaystyle{\|{u-u_{0}}\|}_{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}={\|{\eta}\|}_{2}\eqqcolon\delta_{1}$

is a sensible condition for the denoised image, since one should not look for
an image $u$ further away from $u^{\dagger}$ than $u_{0}$. This motivates to
consider a variant of the ROF model (1) where the discrepancy
${\|{u-u_{0}}\|}_{2}$ is not a penalty in the objective, but taken into
account as a constraint. This leads to a reformulation of the total variation
problem as

$\displaystyle\min_{u}{\|{|\operatorname{\nabla}u|}\|}_{1}\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}={\|{\eta}\|}_{2}.$

By estimating $\eta$ as Gaussian noise from the given image $u_{0}$ (e.g.,
using the method from liu2012noise ; liu2013single ) one obtains a parameter
free denoising method.

Next, assume that we have some additional information on the original image
$u^{\dagger}$ available, namely some estimate $v$ of its gradient. This could
be taken into account as

$\displaystyle\min_{u}{\|{|\operatorname{\nabla}u-v|}\|}_{1}\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq{\|{\eta}\|}_{2}.$ (2)

It turns out, that this information can be quite powerful. The next simple
lemma shows that if we would know the gradient of $u^{\dagger}$ and the noise
level exactly, our model would recover $u^{\dagger}$ perfectly, even for
arbitrary large noise (and also independent of the type of noise).

###### Lemma 1

Assume that $u^{\dagger}$ and $u_{0}$ fulfill
$\int_{\Omega}u^{\dagger}=\int_{\Omega}u_{0}$ and let
$v=\operatorname{\nabla}u^{\dagger}$ and
$\delta_{1}={\|{u^{\dagger}-u_{0}}\|}_{2}$. Then it holds that

$u^{\dagger}=\operatorname*{argmin}_{u}{\|{|\operatorname{\nabla}u-v|}\|}_{1}\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq\delta_{1},$ (3)

i.e. $u^{\dagger}$ is the unique solution of the denoising problem.

###### Proof

The set of minimizers is

$\displaystyle\operatorname*{argmin}_{u}{\|{|\operatorname{\nabla}u-\operatorname{\nabla}u^{\dagger}|}\|}_{1}\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}.$

Clearly, $u^{\dagger}$ is within this set, since the optimal value is $0$ and
$u^{\dagger}$ is feasible, because the constraint is trivially fulfilled.

To show that $u^{\dagger}$ is indeed the unique solution, consider any other
$u$ that also produces an objective value of zero. This implies
$\operatorname{\nabla}u=\operatorname{\nabla}u^{\dagger}$, i.e.
$u=u^{\dagger}+c$ for some constant $c$. Thus, $u$ fulfills the constraint
${\|{u-u_{0}}\|}_{2}^{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}^{2}$ if
${\|{u^{\dagger}-u_{0}+c}\|}_{2}^{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}^{2}$. We
expand the left hand side and get, writing $|\Omega|$ for the measure of
$\Omega$,

${\|{u^{\dagger}-u_{0}}\|}_{2}^{2}+2c\int_{\Omega}(u^{\dagger}-u_{0})+c^{2}|\Omega|\leq{\|{u^{\dagger}-u_{0}}\|}_{2}^{2}.$

Since the middle integral vanishes by assumption, we see that $c=0$.
∎

The next lemma shows, that $v=\operatorname{\nabla}u^{\dagger}$ is also
necessary for perfect reconstruction.

###### Lemma 2

If $\delta_{1}={\|{u^{\dagger}-u_{0}}\|}_{2}$,
$\int_{\Omega}u^{\dagger}=\int_{\Omega}u_{0}$ and $u^{\dagger}$ solves (3),
then $v=\operatorname{\nabla}u^{\dagger}$.

###### Proof

Let
$u^{\dagger}\in\operatorname*{argmin}_{u}{\|{|\operatorname{\nabla}u-v|}\|}_{1}$
s.t. ${\|{u-u_{0}}\|}_{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}$. We denote by
$\operatorname{\mathcal{I}}_{C}(x)$ the indicator function of a set $C$, and
set $K=\nabla$,
$F(u)=\operatorname{\mathcal{I}}_{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}}(u)$,
$G(q)={\|{|q-v|}\|}_{1}$ (i.e. the Fenchel conjugate of $G$ is $G^{*}(\phi)\\\
=\int_{\Omega}v\phi\,\mathrm{d}x+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{1}\leq
1}(\phi)$). The characterization of optimality by the Fenchel-Rockafellar
optimality system (ekeland1976convex, , Remark 4.2) shows that
$(u^{*},\phi^{*})$ is a primal-dual optimal pair if and only if

$\begin{cases}0\in K^{*}\phi^{*}+\partial F(u^{*}),\\\ 0\in-Ku^{*}+\partial
G^{*}(\phi^{*}),\end{cases}$

which amounts to the inclusions

$\displaystyle\begin{cases}0\in K^{*}\phi^{*}+\partial F(u^{*}),\\\ 0\in-
Ku^{*}+\partial G^{*}(\phi^{*}),\end{cases}$ $\displaystyle\Leftrightarrow$
$\displaystyle\begin{cases}0\in
K^{*}\phi^{*}+\partial\operatorname{\mathcal{I}}_{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}}(u^{*}),\\\
0\in-Ku^{*}+v+\partial\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{1}\leq
1}(\phi^{*}).\end{cases}$

That means, that $(u^{\dagger},\phi^{*})$ is optimal if and only if

$\displaystyle\begin{cases}0\in-\operatorname{div}\phi^{*}+\partial\operatorname{\mathcal{I}}_{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq{\|{u^{\dagger}-u_{0}}\|}_{2}}(u^{\dagger}),\\\
0\in-\operatorname{\nabla}u^{\dagger}+v+\partial\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{1}\leq
1}(\phi^{*}).\end{cases}$

Since $u^{\dagger}$ is on the boundary of the domain of the indicator function
in the first inclusion, the subgradient there is the normal cone, which
implies

$\displaystyle\begin{cases}\exists\,t\geq
0\,:\,0=\operatorname{div}\phi^{*}+t(u^{\dagger}-u_{0}),\\\
\begin{cases}\operatorname{\nabla}u^{\dagger}(x)=v(x),\\\ \qquad\text{if}\
|\phi^{*}(x)|<1,\\\ \exists\,s(x)\geq
0\,:\,\operatorname{\nabla}u^{\dagger}(x)-v(x)=s(x)\phi^{*}(x),\\\
\qquad\text{if}\ |\phi^{*}(x)|=1.\end{cases}\end{cases}$

By a result of Bourgain and Brezis (bourgain2003divYf, , Proposition 1) there
is an $L^{\infty}$ solution $\phi$ of
$-\operatorname{div}\phi=u^{\dagger}-u_{0}$, i.e. there exists $L>0$ such that
$|\phi|\leq L$ a.e. and hence for $\tilde{\phi}=\phi/(L+1)$ it holds that
$-\operatorname{div}\tilde{\phi}=\tfrac{1}{L+1}(u^{\dagger}-u_{0})$ and
$|\tilde{\phi}|<1$ a.e. Hence, $v=\operatorname{\nabla}u^{\dagger}$ a.e. ∎

###### Remark 1

If the condition $\int_{\Omega}u^{\dagger}=\int_{\Omega}u_{0}$ in Lemma 1 does
not hold, but $\int_{\Omega}(u^{\dagger}-u_{0})=\varepsilon$, then the proof
of that lemma still shows that all solutions of (3) are of the form
$u^{\dagger}+c$ with
$|c-\tfrac{\varepsilon}{|\Omega|}|\leq\tfrac{\varepsilon}{|\Omega|}$.

If $v\neq\operatorname{\nabla}u^{\dagger}$, then any solution $\tilde{u}$ of
(3) will usually be different from $u^{\dagger}$, although, it will fulfill
the trivial estimate

$\displaystyle{\|{\tilde{u}-u^{\dagger}}\|}_{2}$
$\displaystyle\leq{\|{\tilde{u}-u_{0}}\|}_{2}+{\|{u_{0}-u^{\dagger}}\|}_{2}$
$\displaystyle\leq 2{\|{u_{0}-u^{\dagger}}\|}_{2}=\delta_{1}.$

However, the following lemma shows that $\tilde{u}\to u^{\dagger}$ for
$v\to\operatorname{\nabla}u^{\dagger}$ (for constant noise level
$\delta_{1}$):

###### Lemma 3

Assume $\Omega$ is convex, $u^{\dagger}$ and $u_{0}$ fulfill
$\int_{\Omega}u^{\dagger}=\int_{\Omega}u_{0}$, let $v$ fulfill
${\|{|v-\operatorname{\nabla}u^{\dagger}|}\|}_{1}\leq\varepsilon$ and assume
that there exists a solution $\tilde{u}$ of (3) with
$\delta_{1}={\|{u^{\dagger}-u_{0}}\|}_{2}$ for which the constraint is active
(i.e. ${\|{\tilde{u}-u_{0}}\|}_{2}=\delta_{1}$).

Then there exists another solution $\bar{u}$ of (3) with
$\delta_{1}={\|{u^{\dagger}-u_{0}}\|}_{2}$ that fulfills
$\int_{\Omega}\bar{u}=\int_{\Omega}u^{\dagger}$ and moreover, it holds that

${\|{\bar{u}-u^{\dagger}}\|}_{2}\leq\operatorname{diam}(\Omega)\varepsilon$

where $\operatorname{diam}(\Omega)$ denotes the diameter of $\Omega$.

###### Proof

To obtain $\bar{u}$ we consider $\bar{u}=\tilde{u}+c$ for a suitable constant
$c$. The equality $\int_{\Omega}\bar{u}=\int_{\Omega}u^{\dagger}$ is achieved
for $c=\int_{\Omega}(u^{\dagger}-\tilde{u})/|\Omega|$. Since
$\operatorname{\nabla}\bar{u}=\operatorname{\nabla}\tilde{u}$ holds, $\bar{u}$
is optimal for (3) as soon as it is feasible. To check feasibility we
calculate

$\displaystyle{\|{\bar{u}-u_{0}}\|}_{2}^{2}$
$\displaystyle={\|{\tilde{u}-u_{0}+c}\|}_{2}^{2}$
$\displaystyle={\|{\tilde{u}-u_{0}}\|}_{2}^{2}+2c\int_{\Omega}(\tilde{u}-u_{0})+c^{2}|\Omega|.$

Since

$\displaystyle\int_{\Omega}(\tilde{u}-u_{0})$
$\displaystyle=\int_{\Omega}(\tilde{u}-u^{\dagger})+\int_{\Omega}(u^{\dagger}-u_{0})$
$\displaystyle=-c|\Omega|+0,$

we get that

$\displaystyle{\|{\bar{u}-u_{0}}\|}_{2}^{2}$
$\displaystyle={\|{\tilde{u}-u_{0}}\|}_{2}^{2}-c^{2}|\Omega|$

which shows feasibility of $\bar{u}$.

Now we use the Poincaré-Wirtinger inequality in $L^{1}$ for which the optimal
constant is known from acosta2004optimal  to be
$\operatorname{diam}(\Omega)/2$, i.e. it holds that

$\displaystyle{\|{\bar{u}-u^{\dagger}}\|}_{1}$
$\displaystyle\leq\tfrac{\operatorname{diam}(\Omega)}{2}{\|{|\operatorname{\nabla}(\bar{u}-u^{\dagger})|}\|}_{1}$
$\displaystyle\leq\tfrac{\operatorname{diam}(\Omega)}{2}({\|{|\operatorname{\nabla}\bar{u}-v|}\|}_{1}+{\|{|v-\operatorname{\nabla}u^{\dagger}|}\|}_{1}).$

By optimality of $\bar{u}$ and feasibility of $u^{\dagger}$ we get
${\|{|\operatorname{\nabla}\bar{u}-v|}\|}_{1}\\\
\leq{\|{|\operatorname{\nabla}u^{\dagger}-v|}\|}_{1}$ and hence

${\|{\bar{u}-u^{\dagger}}\|}_{1}\leq\operatorname{diam}(\Omega)\varepsilon.\hfill\qed$

###  2.2 Denoising of image gradients

The previous lemma shows that any approximation $v$ of the true gradient
$\operatorname{\nabla}u^{\dagger}$ is helpful for total variation denoising
according to (3). In order to determine such a $v$, one way is to denoise the
gradient of the input image, specifically, by a variational method with a
smoothness penalty for the gradient and some discrepancy term. Naturally, a
norm of the derivate of the gradients can be used. A first candidate could be
the Jacobian of the gradient, i.e.

$\displaystyle
J(\operatorname{\nabla}u)=\begin{pmatrix}\partial_{1}(\partial_{1}u)&\partial_{2}(\partial_{1}u)\\\
\partial_{1}(\partial_{2}u)&\partial_{2}(\partial_{2}u)\end{pmatrix}$

which amounts to the Hessian of $u$. Thus, the matrix is symmetric as soon as
$u$ is twice continuously differentiable. However, notice, that the Jacobian
of an arbitrary vector field is not necessarily symmetric and hence using
${\|{J(v)}\|}$ as smoothness penalty seems unnatural. Instead, we could use
the symmetrized Jacobian,

$\displaystyle\operatorname{\mathcal{E}}(v)=\begin{pmatrix}\partial_{1}v_{1}&\tfrac{1}{2}(\partial_{1}v_{2}+\partial_{2}v_{1})\\\
\tfrac{1}{2}(\partial_{1}v_{2}+\partial_{2}v_{1})&\partial_{2}v_{2}\end{pmatrix},$

where $v_{1}$ and $v_{2}$ are the components of $v$. Note that for twice
differentiable $u$ we have

$\displaystyle\operatorname{\mathcal{E}}(\operatorname{\nabla}u)=J(\operatorname{\nabla}u)=\operatorname{Hess}(u),$

i.e. in both cases we obtain the Hessian of $u$. Imitating the
$\operatorname{TV}$-seminorm (and also followoing the idea of total
generalized variation), we take
$F(v)={\|{|\operatorname{\mathcal{E}}v|}\|}_{1}$.

Similar to the constraint in (3) the denoised gradient should not differ more
from the true gradient than $\operatorname{\nabla}u_{0}$, thus, we consider
the minimization problem with a constraint

$\displaystyle\min_{v}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}\text{\quad
s.t.\quad}{\|{|\operatorname{\nabla}u_{0}-v|}\|}_{1}\leq\delta_{2}$

The parameter $\delta_{2}$ can be chosen as follows: If we set
$\delta_{2}\coloneqq c{\|{|\operatorname{\nabla}u_{0}|}\|}_{1}$, then $c=1$
would allow the trivial minimizer $v=0$ and any $c<1$ will enforce some
structure of $\operatorname{\nabla}u_{0}$ onto the minimizer and smaller $c$
leads to less denoising.

Putting the pieces together, we arrive at a two-stage denoising method:

  1. 1.

Choose $0<c<1$ and calculate a denoised gradient by solving

$\displaystyle\begin{split}\hat{v}\in&\operatorname*{argmin}_{v}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}\\\
&\text{\quad s.t.\quad}{\|{|\operatorname{\nabla}u_{0}-v|}\|}_{1}\leq
c{\|{|\operatorname{\nabla}u_{0}|}\|}_{1}\eqqcolon\delta_{2}.\end{split}$ (4)

  2. 2.

Denoise $u_{0}$ by solving

$\displaystyle\begin{split}\hat{u}\in&\operatorname*{argmin}_{u}{\|{|\operatorname{\nabla}u-\hat{v}|}\|}_{1}\\\
&\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq{\|{\eta}\|}_{2}\eqqcolon\delta_{1}.\end{split}$
(5)

Instead of using the constrained formulation of the first problem, we can also
use a penalized formulation. Thus, the gradient denoising problem writes as:

  1. 1.

Choose $\alpha>0$ and calculate a denoised gradient by solving

$\displaystyle\hat{v}$
$\displaystyle=\operatorname*{argmin}_{v}{\|{|\operatorname{\nabla}u_{0}-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}$
(6)

  2. 2.

Denoise $u_{0}$ by solving (5).

Since we use a denoised gradient prior to apply total variation denoising, we
term the method (4) and (5) Denoised Gradient Total Variation (DGTV). Due to
the similarity with total generalized variation, we call the second method
based on (6) and (5) Denoised Gradient Total Generalized Variation (DGTGV).

###  2.3 Constrained and Morozov total generalized variation

Both two-stage denoising methods DGTV and DGTGV for the gradient resemble
previously known methods: The latter is related to total generalized variation
(TGV) bredies2010total  while the former to constrained total generalized
variation (CTGV) komander2017denoising . The TGV of second order, defined in
bredies2010total  has been shown to be equal to

$\displaystyle\operatorname{TGV}^{2}_{(\alpha_{0},\alpha_{1})}(u)=\min_{v}\,\alpha_{1}{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha_{0}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}$
(7)

in bredies2013properties ; bredies2010total  while CTGV from
komander2017denoising  is defined as

$\displaystyle\operatorname{CTGV}_{\delta}(u)=\min_{v}\Big{\\{}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}\text{\quad
s.t.\quad}{\|{|\operatorname{\nabla}u-v|}\|}_{1}\leq\delta\Big{\\}}.$ (8)

Considering $\operatorname{\nabla}u$ to be the given data in these problems,
one could say, following the notion from lorenz2013neccond , that TGV is a
Tikhonov-type estimation of $\operatorname{\nabla}u$ while CTGV is a Morozov-
type estimation of $\operatorname{\nabla}u$.

Now, combining the two steps of DGTGV into one optimization problem, where in
each step the image as well as the gradient is updated simultaneously, we get

$\displaystyle\begin{split}&\min_{u}\Big{\\{}\operatorname{TGV}^{2}_{(\alpha,1)}(u)\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq\delta_{1}\Big{\\}}\\\
=&\min_{u,v}\Big{\\{}{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}\\\
&\qquad\qquad\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq\delta_{1}\Big{\\}}.\end{split}$ (9)

This formulation is a Morozov type formulation of the TGV problem

$\min_{u}\tfrac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}+\operatorname{TGV}^{2}_{(\alpha_{0},\alpha_{1})}(u)$
(10)

and thus, in the following referred to as MTGV.

Another approach to form a combined optimization problem out of DGTGV is to
preserve both constraints, i.e. taking (4) and (5) to obtain

$\displaystyle\begin{split}\min_{u,v}\,{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}\text{\quad
s.t.\quad}&{\|{u-u_{0}}\|}_{2}\leq\delta_{1},\\\
&{\|{|\operatorname{\nabla}u-v|}\|}_{1}\leq\delta_{2}.\end{split}$ (11)

In both problem formulations $\delta_{1}$ again is the noise level. Using (8)
we see that problem (11) becomes

$\displaystyle\min_{u}\,\operatorname{CTGV}_{\delta_{2}}(u)\text{\quad
s.t.\quad}{\|{u-u_{0}}\|}_{2}\leq\delta_{1}.$ (12)

Obviously the TGV, MTGV and the CTGV denoising problems are implicitly
equivalent in the sense that, knowing the solution to one of the problems
allows to calculate the respective parameters of one of the other problems
such that the solution stays one (cf. (lorenz2013neccond, , Theorem 2.3) for a
general result and (komander2017denoising, , Lemma 1, Lemma 2) for a result in
the case of CTGV and TGV).

###  2.4 Parameter choice

A few words on parameter choice for all methods are in order. Frequently, a
constraint ${\|{u-u_{0}}\|}_{2}\leq\delta_{1}$ appears in the problem and the
parameter $\delta_{1}$ has a large influence on the denoising result. A
natural choice is to adapt the parameter to the noise level, i.e.

$\displaystyle\delta_{1}={\|{\eta}\|}_{2}={\|{u_{0}-u^{{\dagger}}}\|}_{2}.$

For a given discrete image, this number can be estimated as follows: Under the
assumption that the noise is additive Gaussian white noise, the methods from
liu2012noise ; liu2013single  allow to estimate the standard deviation
$\sigma$ of the noise $\eta$. Then, it is well known that the 2-norm of $\eta$
is estimated ny ${\|{\eta}\|}_{2}\approx\sigma\sqrt{k}$, where $k$ is the
number of pixels.

The second parameter in DGTV is
$\delta_{2}=c{\|{|\operatorname{\nabla}u_{0}|}\|}_{1}$, the constraint
parameter in (4). Setting $c=1$ leads to $v=0$ as a feasible and also optimal
solution, hence, the gradient we insert into the second optimization problem
is zero and thus the second problem becomes pure total variation denoising
without any additional information. Therefore, $c\in(0,1)$ is a reasonable
choice. Experiments showed (cf. section 5) that a lot of gradient denoising,
i.e. smoothening of $\operatorname{\nabla}u_{0}$, leads to good
reconstructions of the image in the second step. Thus, we set $c\approx 0.99$.

The method DGTGV, the penalized variant of the two-stage method, includes the
parameter $\alpha$, the penalization parameter in (6). Again experiments
showed that $\alpha=1$ leads to a good image reconstruction (cf. section 5),
independent of the noise level, the image size or the image type.

Numerical experiments on the performance of the two-stage methods with regard
to image quality and computational speed can be found in section 5.

Also the problems CTGV-, TGV-, and MTGV-denoising come with two parameters
each that have to be chosen. For CTGV- and MTGV-denoising the choice for
$\delta_{1}$ is the same as above. For the parameter $\alpha$ in MTGV-
denoising (9) we have experimental experience that hints that $\alpha\approx
2$ is a good universal parameter (cf. section 5). This inline with the usual
recommendation that $\alpha_{1}=2\alpha_{0}$ is a good choice for TGV
denoising from (10) cf. KnollBrediesTGV2011 . For the remaining parameter
$\delta_{2}$ for CTGV there is following heuristic from komander2017denoising
: We denote with $u^{\operatorname{TV}}$ the $\operatorname{TV}$ denoised
image with $\delta_{1}$ chosen according to the noise level estimate from
Section 2.4 and set
$\delta_{2}=c{\|{|\operatorname{\nabla}u^{\operatorname{TV}}|}\|}_{1}$ with
$0<c<1$. However, CTGV will not be included in the experiment in section 5 and
the main reason is, that the parameter choice here does lead to inferior
results compared to MTGV.

##  3 Duality gaps

Besides the choice of the problem parameters, controlling the denoising, one
has to choose a reasonable stopping criterion. Since all considered problems
are convex, the duality gap is a natural candidate. Fenchel-Rockafellar
optimality ekeland1976convex  states that, under appropriate conditions, the
primal problem $\min F(x)+G(Kx)$ has the dual problem
$\max_{y}-F^{*}(-K^{*}y)-G^{*}(y)$, that the duality gap, defined as

$\operatorname{gap}(x,y)=F(x)+G(Kx)+F^{*}(-K^{*}y)+G^{*}(y)$

is an upper bound on the difference from the current objective value to the
optimal one, i.e. $\operatorname{gap}(x,y)\geq F(x)+G(Kx)-F(x^{*})-G(Kx^{*})$
where $x^{*}$ is a solution of the primal problem, and that the duality gap
vanishes exactly at primal-dual optimal pairs.

For the MTGV denoising problem (9) we have the gap function:

$\begin{split}\operatorname{gap}&{}_{\operatorname{MTGV}}(u,v,p,q)=\\\
&{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\operatorname{\mathcal{I}}_{\\{{\|{\cdot-
u_{0}}\|}_{2}\leq\delta_{1}\\}}(u)\\\
&+\delta_{1}{\|{\operatorname{\nabla}^{*}p}\|}_{2}-\langle
p,\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{0\\}}(p-\operatorname{\mathcal{E}}^{*}q)\\\
&+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha\\}}(p)+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq
1\\}}(q).\end{split}$ (13)

As also noted in bredies2015preconditioned  (for the TGV denoising problem
(10)) this gap is not helpful: The problem is, that the indicator function
$\operatorname{\mathcal{I}}_{\\{0\\}}(p-\operatorname{\mathcal{E}}^{*}q)$ is
usually not finite as $p=\operatorname{\mathcal{E}}^{*}q$ is usually not
fulfilled. To circumvent this problem, one could simply replace $p$ by
$\operatorname{\mathcal{E}}^{*}q$ in the gap function. If we do this, we
obtain a gap that only depends on $q$:

$\begin{split}\operatorname{gap}&{}_{\operatorname{MTGV}}(u,v,q)=\\\
&{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\operatorname{\mathcal{I}}_{\\{{\|{\cdot-
u_{0}}\|}_{2}\leq\delta_{1}\\}}(u)\\\
&+\delta_{1}{\|{\operatorname{\nabla}^{*}\operatorname{\mathcal{E}}^{*}q}\|}_{2}-\langle\operatorname{\mathcal{E}}^{*}q,\operatorname{\nabla}u_{0}\rangle\\\
&+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha\\}}(\operatorname{\mathcal{E}}^{*}q)+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq
1\\}}(q).\end{split}$ (14)

We still have the problem, that
$\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha\\}}(\operatorname{\mathcal{E}}^{*}q)$
might be infinite. So replacing $q$ and $p$ in $\eqref{eq:gapmtgv}$ by

$\displaystyle\tilde{q}$
$\displaystyle\coloneqq\frac{q}{\max\left(1,\frac{{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}}{\alpha}\right)}$
$\displaystyle\tilde{p}$
$\displaystyle\coloneqq\frac{\operatorname{\mathcal{E}}^{*}q}{\max\left(1,\frac{{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}}{\alpha}\right)}=\operatorname{\mathcal{E}}^{*}\tilde{q},$

we obtain a finite gap function which is a valid stopping criterion for the
problem:

###### Lemma 4

Let
$\Phi(u,v)\coloneqq{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\operatorname{\mathcal{I}}_{\\{{\|{\cdot-
u_{0}}\|}_{2}\leq\delta_{1}\\}}(u)$ be the primal functional of the MTGV
denoising problem and $(u^{*},v^{*})$ be a solution of $\min_{u,v}\Phi(u,v)$.
Then the error of the primal energy can be estimated by

$\begin{split}\Phi&(u,v)-\Phi(u^{*},v^{*})\\\
\leq&{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\operatorname{\mathcal{I}}_{\\{{\|{\cdot-
u_{0}}\|}_{2}\leq\delta_{1}\\}}(u)\\\
&+\delta_{1}{\|{\operatorname{\nabla}^{*}(\operatorname{\mathcal{E}}^{*}\tilde{q})}\|}_{2}-\langle\operatorname{\mathcal{E}}^{*}\tilde{q},\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq
1\\}}(\tilde{q}),\end{split}$

where
$\tilde{q}\coloneqq\frac{q}{\max\left(1,\frac{{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}}{\alpha}\right)}$
and $q$ is any feasible dual variable, i.e. ${\|{|q|}\|}_{\infty}\leq 1$.
Moreover, if $(u^{n},v^{n},q^{n})$ converge to a primal-dual solution, the
upper bound converges to zero.

###### Proof

The estimate is clear since any dual value is smaller or equal than any primal
value. So we can change the dual variables in the gap as we like.

Since the gap function is continuous on its domain, the sequence
$(u^{n},v^{n},q^{n})$ converges to a primal-dual solution and has the same
limit as $(u^{n},v^{n},\tilde{q}^{n})$, the second claim follows.∎

The remaining indicator functions that are left in the modified gap are
guaranteed to be finite by several algorithms (namely ones that use
projections onto the respective constraints). As we will see later, this
holds, for example, for classical methods like Douglas-Rachford and Chambolle-
Pock.

Similar problems with an infinite duality gap appear in other problems, too,
e.g. in the CTGV denoising (12). A closer look at the MTGV denoising problem
(9) reveals that the above construction of the modified gap can indeed be
avoided by a different choice of variables and that this holds for a broad
class of problems. First, we illustrate this for the MTGV problem. The MTGV
problem does not change if we introduce a new variable
$w=\operatorname{\nabla}u-v$ and replace $v$ in the formulation. We obtain:

$\begin{split}\min_{u,w}\,{\|{|w|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(\operatorname{\nabla}u-w)|}\|}_{1}\text{\quad
s.t.\quad}&{\|{u-u_{0}}\|}_{2}\leq\delta_{1}.\end{split}$ (15)

For this problem, the gap function is

$\begin{split}\operatorname{gap}&{}_{\operatorname{MTGV}}(u,w,q)=\\\
&{\|{|w|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(\operatorname{\nabla}u-w)|}\|}_{1}+\operatorname{\mathcal{I}}_{\\{{\|{\cdot-
u_{0}}\|}_{2}\leq\delta_{1}\\}}(u)\\\
&+\delta_{1}{\|{\operatorname{\nabla}^{*}(\operatorname{\mathcal{E}}^{*}(q))}\|}_{2}-\langle\operatorname{\mathcal{E}}^{*}q,\operatorname{\nabla}u_{0}\rangle\\\
&+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha\\}}(\operatorname{\mathcal{E}}^{*}q)+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq
1\\}}(q),\end{split}$ (16)

which is exactly the same as $\eqref{eq:gapmtgvmod}$. As above the indicator
function
$\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{0}\\}}(\operatorname{\mathcal{E}}^{*}q)$
might be infinite and therefore $q$ should be replaced by $\tilde{q}$ from
above. So replacing $p$ by $\operatorname{\mathcal{E}}^{*}q$ in
$\eqref{eq:gapmtgv}$ results in the gap function of
$\eqref{eq:mtgvvarchange}$. Both gap functions can be used as valid stopping
criteria for both problem formulations.

This method for the stopping criterion does apply to general problems of the
form

$\displaystyle\begin{split}\min_{u,v}F(u)+G(Av)+H(Bu-v),\end{split}$ (17)

where $A$ and $B$ are linear (and standard regularity conditions, implying
Fenchel-Rockafellar duality is fulfilled). The dual problem is

$\displaystyle\max_{p,q}-F^{*}(-B^{*}p)-H^{*}(p)-G^{*}(q)+\operatorname{\mathcal{I}}_{\\{0\\}}(p-A^{*}q).$

If we replace $p$ by $A^{*}q$, we obtain

$\displaystyle\max_{q}-F^{*}(-B^{*}A^{*}q)-H^{*}(A^{*}q)-G^{*}(q),$

which is the dual problem of $\eqref{eq:minprob}$ with variable change
$w=Bu-v$:

$\displaystyle\min_{u,w}F(u)+G(A(Bu-w))+H(w).$

The TGV denoising problem for example is very similar to MTGV. The gap is:

$\displaystyle\operatorname{gap}$ ${}_{\operatorname{TGV}}(u,v,p,q)=$
$\displaystyle\alpha_{1}{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha_{0}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}$
$\displaystyle+\frac{1}{2}{\|{\operatorname{\nabla}^{*}p}\|}_{2}^{2}-\langle
p,\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{0\\}}(p-\operatorname{\mathcal{E}}^{*}q)$
$\displaystyle+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{0}\\}}(p)+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{1}\\}}(q).$

With $\tilde{p}$ and $\tilde{q}$ as before one gets a simple gap for TGV:

$\displaystyle\operatorname{gap}$ ${}_{\operatorname{TGV}}(u,v,\tilde{q})=$
$\displaystyle\alpha_{1}{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha_{0}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}$
$\displaystyle+\frac{1}{2}{\|{\operatorname{\nabla}^{*}(\operatorname{\mathcal{E}}^{*}\tilde{q})}\|}_{2}^{2}-\langle\operatorname{\mathcal{E}}^{*}\tilde{q},\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{1}\\}}(\tilde{q}).$

###### Corollary 1

Let
$\Phi(u,v)\coloneqq{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}$
be the primal functional of the TGV denoising problem and $(u^{*},v^{*})$ be a
solution of $\min_{u,v}\Phi(u,v)$. Then the error of the primal energy can be
estimated by

$\begin{split}\Phi&(u,v)-\Phi(u^{*},v^{*})\\\
&\leq\alpha_{1}{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha_{0}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}\\\
&+\frac{1}{2}{\|{\operatorname{\nabla}^{*}(\operatorname{\mathcal{E}}^{*}\tilde{q})}\|}_{2}^{2}-\langle\operatorname{\mathcal{E}}^{*}\tilde{q},\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{1}\\}}(\tilde{q}),\end{split}$

where
$\tilde{q}\coloneqq\frac{q}{\max\left(1,\frac{{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}}{\alpha}\right)}$.
Moreover, if $(u^{n},v^{n},q^{n})$ converge to a primal-dual solution, the
upper bound converges to zero.

##  4 Numerics

In this section we describe methods to solve the convex optimization problems
related to the various denoising methods from the previous sections. We will
work with standard discretizations of the images and the derivative operators,
but state them for the sake of completeness in Appendix A.1. In this section
we focus on the optimization methods.

###  4.1 Douglas-Rachford’s method

The Douglas-Rachford algorithm (see douglasrachford1956  and
lions1979splitting ) is a splitting algorithm to solve monotone inclusions of
the form $0\in\mathcal{A}(x)+\mathcal{B}(x)$, which requires only the
resolvents $\mathcal{R}_{t\mathcal{A}}=(I+t\mathcal{A})^{-1}$ and
$\mathcal{R}_{t\mathcal{B}}=(I+t\mathcal{B})^{-1}$ but not the resolvent of
the sum $\mathcal{A}+\mathcal{B}$.

One way to write down the Douglas-Rachford iteration is as a fixed point
iteration

$\displaystyle z^{k}$ $\displaystyle=F(z^{k-1}),\quad\text{ with}$
$\displaystyle F(z)$
$\displaystyle=z+\mathcal{R}_{t\mathcal{A}}(2\mathcal{R}_{t\mathcal{B}}(z)-z)-\mathcal{R}_{t\mathcal{B}}(z)$

for some step-size $t>0$. It is possible to employ relaxation for the Douglas-
Rachford iteration as

$z^{k}=z^{k-1}+\rho(F(z^{k-1})-z^{k-1}),$ (18)

where $1<\rho<2$ is overrelaxation and $0<\rho<1$ is underrelaxation. For
$\rho$ in the whole range from $0$ to $2$ and any $t>0$ it holds that the
iteration converges to some fixed point $z$ of $F$ such that
$\mathcal{R}_{t\mathcal{B}}(z)$ is a zero of $\mathcal{A}+\mathcal{B}$, see,
e.g. eckstein1992douglas .

The Douglas-Rachford method can be used to solve the saddle point problem

$\displaystyle\begin{split}\min_{x}\max_{y}\ \langle
y,Kx\rangle+F(x)-G^{*}(y).\end{split}$ (19)

To that end, in oconnor2014primal  the authors propose to use the splitting

$\displaystyle\begin{split}\mathcal{B}(x,y)&=\left[\begin{tabular}[]{c}$\partial
F(x)$\\\ $\partial G^{*}(y)$\end{tabular}\right],\\\
\mathcal{A}(x,y)&=\left[\begin{tabular}[]{l r}$0$&$K^{*}$\\\
$-K$&$0$\end{tabular}\right]\left[\begin{tabular}[]{c}$x$\\\
$y$\end{tabular}\right].\end{split}$ (20)

If $\mathcal{A}$ is a matrix the resolvent $\mathcal{R}_{t\mathcal{A}}$ simply
is the matrix inverse $(I+t\mathcal{A})^{-1}$. The resolvent of $\mathcal{B}$
is given by the proximal mappings of $F$ and $G^{*}$, namely

$\displaystyle\mathcal{R}_{t\mathcal{B}}(x,y)=\left[\begin{array}[]{c}(I+t\partial
F)^{-1}(x)\\\ (I+t\partial
G^{*})^{-1}(y)\end{array}\right]=\left[\begin{array}[]{c}\operatorname{prox}_{tF}(x)\\\
\operatorname{prox}_{tG^{*}}(y)\end{array}\right].$

For further flexibility, it is proposed in oconnor2014primal  to rescale the
problem by replacing $G$ with $\tilde{G}(y)\coloneqq G(\frac{y}{\beta})$ and
$K$ with $\tilde{K}\coloneqq\beta K$. Then one can replace the dual variable
by $\tilde{y}^{k}\coloneqq\beta y^{k}$ and finally use
$\operatorname{prox}_{t\tilde{G}^{*}}(y)=\frac{1}{\beta}\operatorname{prox}_{t\beta^{2}G^{*}}(\beta
y)$ to obtain a second stepsize $s\coloneqq\beta^{2}t$, which can be chosen
independently from the stepsize $t$. In total, the resolvents then become

$\displaystyle\mathcal{R}_{t\mathcal{B}}(x,y)=\left[\begin{array}[]{c}\operatorname{prox}_{tF}(x)\\\
\operatorname{prox}_{sG^{*}}(y)\end{array}\right],$
$\displaystyle\mathcal{R}_{t\mathcal{A}}=\left[\begin{tabular}[]{l
r}$I$&$tK^{*}$\\\ $-sK$&$I$\end{tabular}\right]^{-1}.$

In each step of Douglas-Rachford we need the inverse of a fairly large block
matrix. However, as also noted in oconnor2014primal  this can be done
efficiently with the help of the Schur complement. If $K$ is a matrix,
$K^{*}=K^{T}$ and

$\displaystyle\left[\begin{array}[]{cc}I&tK^{T}\\\
-sK&I\end{array}\right]^{-1}$ $\displaystyle=\left[\begin{array}[]{cc}0&0\\\
0&I\end{array}\right]+\left[\begin{array}[]{c}I\\\
sK\end{array}\right](I+stK^{T}K)^{-1}\left[\begin{array}[]{c}I\\\
-tK\end{array}\right]^{T}.$

To use this, we need to solve equations with the positive definite coefficient
matrix $I+\lambda K^{T}K$ efficiently.

In the following we interpret the linear operators as matrices, without
changing notation, i.e. the matrix $\operatorname{\mathcal{E}}$ also stands
for the matrix realizing the linear operation $\operatorname{\mathcal{E}}$.
With the help of the vectorization operation $\operatorname{vec}$, it holds
that
$\operatorname{\mathcal{E}}\cdot\operatorname{vec}(u)=\operatorname{vec}(\operatorname{\mathcal{E}}u)$,
where on the left hand side $\operatorname{\mathcal{E}}$ is a matrix and on
the right hand side $\operatorname{\mathcal{E}}$ is the linear operator from
A.1.

As discussed is section 3, we have different possibilities to choose the
primal variables (and thus, also for the dual variables): Namely, we could use
the primal variables $u$ and $v$, as, e.g. in the MTGV problem (9), or the
primal variables $u$ and $w$, as in the corresponding problem (15). This
choice does not only influence the dual problem and the duality gap, but also
the involved linear map $K^{T}K$. The formulation with $u$ and $v$ as primal
variables leads to

$\displaystyle K^{T}K=\begin{pmatrix}-\operatorname{div}&0\\\
-I&\mathcal{E}^{T}\end{pmatrix}\begin{pmatrix}\nabla&-I\\\
0&\mathcal{E}\end{pmatrix}=\begin{pmatrix}-\Delta&\operatorname{div}\\\
-\nabla&I+\mathcal{E}^{T}\mathcal{E}\end{pmatrix}$ (21)

while the variable change to $u$ and $w$ gives

$\displaystyle K^{T}K$
$\displaystyle=\begin{pmatrix}-\operatorname{div}(\mathcal{E}^{T})\\\
-\mathcal{E}^{T}\end{pmatrix}\begin{pmatrix}\mathcal{E}(\nabla)&-\mathcal{E}\end{pmatrix}=\begin{pmatrix}H^{T}\\\
-\mathcal{E}^{T}\end{pmatrix}\begin{pmatrix}H&-\mathcal{E}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}H^{T}H&-H^{T}\mathcal{E}\\\
-\mathcal{E}^{T}H&\mathcal{E}^{T}\mathcal{E}\end{pmatrix},$ (22)

where $H\coloneqq\mathcal{E}(\nabla)$ is the Hessian matrix, cf. Appendix
A.2.. The same is true for CTGV from (11) and the usual TGV problem (10). For
problem (6) we have

$\displaystyle K^{T}K$
$\displaystyle=\operatorname{\mathcal{E}}^{T}\operatorname{\mathcal{E}}.$ (23)

###  4.2 Inexact Douglas-Rachford

For the solution of the linear equation with coefficient matrix $I+\lambda
K^{T}K$ we propose to use the preconditioned conjugate gradient method (PCG).
For preconditioning we do several approximations: First we replace complicated
discrete difference operators by simpler one and then we use a block-diagonal
preconditioner and use the incomplete Cholesky decomposition in each block. As
we will see, with these preconditioners we only need one or two iterations of
PCG to obtain good convergence results of the Douglas-Rachford iteration. If
we only use one iteration of PCG and denote $A=(I+stK^{T}K)$, the linear step
is approximated by

$\displaystyle A^{-1}b\approx x_{k}-aM(Ax_{k}+b),$

where the PCG stepsize is given by

$\displaystyle a=\frac{\langle b-Ax_{k},M(b-Ax_{k})\rangle}{\langle
AM(b-Ax_{k}),M(b-Ax_{k})\rangle}$

and $M$ is the preconditioner for $A$, i.e. $M^{-1}\approx A$.

In the following the coefficient matrices and preconditioners are given in
detail. With the matrices $D_{1}$ and $D_{2}$ representing the derivatives in
the first and second direction, we have

$\displaystyle\operatorname{\nabla}=\begin{pmatrix}D_{1}\\\
D_{2}\end{pmatrix},\quad\mathcal{E}=\begin{pmatrix}D_{1}&0\\\
\frac{1}{2}D_{2}&\frac{1}{2}D_{1}\\\ \frac{1}{2}D_{2}&\frac{1}{2}D_{1}\\\
0&D_{2}\end{pmatrix},$ $\displaystyle J=\begin{pmatrix}D_{1}&0\\\ D_{2}&0\\\
0&D_{1}\\\ 0&D_{2}\end{pmatrix},\quad H=\begin{pmatrix}D_{1}^{2}\\\
D_{1}D_{2}\\\ D_{1}D_{2}\\\ D_{2}^{2}\end{pmatrix},$
$\displaystyle-\Delta=\nabla^{T}\nabla=D_{1}^{T}D_{1}+D_{2}^{T}D_{2},$
$\displaystyle\mathcal{E}^{T}\operatorname{\mathcal{E}}=\begin{pmatrix}D_{1}^{T}D_{1}+\frac{1}{2}D_{2}^{T}D_{2}&\frac{1}{2}D_{2}^{T}D_{1}\\\
\frac{1}{2}D_{1}^{T}D_{2}&\frac{1}{2}D_{1}^{T}D_{1}+D_{2}^{T}D_{2}\end{pmatrix},$
$\displaystyle J^{T}J=\begin{pmatrix}D_{1}^{T}D_{1}+D_{2}^{T}D_{2}&0\\\
0&D_{1}^{T}D_{1}+D_{2}^{T}D_{2}\end{pmatrix},$ $\displaystyle
H^{T}H=(D_{1}^{2})^{T}D_{1}^{2}+2(D_{1}D_{2})^{T}(D_{1}D_{2})+(D_{2}^{2})^{T}D_{2}^{2},$
$\displaystyle-\operatorname{\mathcal{E}}^{T}H=\begin{pmatrix}-D_{1}^{T}D_{1}^{2}-D_{2}^{T}(D_{1}D_{2})\\\
-D_{1}^{T}(D_{1}D_{2})-D_{2}^{T}D_{2}^{2}\end{pmatrix}.$

The operator $J^{T}J$ is in fact the negative (discrete) Laplace operator
applied component-wise and also called vector Laplacian. Note that the
operator $\operatorname{\mathcal{E}}^{T}\operatorname{\mathcal{E}}$ decomposes
as

$\operatorname{\mathcal{E}}^{T}\operatorname{\mathcal{E}}=\tfrac{1}{2}J^{T}J+\tfrac{1}{2}\begin{bmatrix}D_{1}^{T}\\\
D_{2}\end{bmatrix}\begin{bmatrix}D_{1}D_{2}\end{bmatrix}.$

Note that the boundary conditions are implicitely contained in the
discretization operators and adjoints (e.g. we use Neumann boundary conditions
for the gradient and thus, equations $-\Delta u=f$ always contain Neumann
boundary conditions).

Our linear operators are dicretizations of continuous differential operators.
Hence, the resolvent steps correspond to solutions of certain differential
equations. In this context it has been shown to be beneficial, to motivate
preconditioners for the linear systems by their continuous counterparts, see
mardal2011preconditioning . Block diagonal preconditioners are a natural
choice for these operators. The conjugate gradient method is usually used for
linear systems of the form $\mathcal{A}x=f$, where $f$ is an element of a
finite dimensional space $X$. As discussed in mardal2011preconditioning  it
can also be used in the case where $X$ is infinite dimensional and
$\mathcal{A}:X\to X$ is a symmetric and positive-definite isomorphism.
Consider for example the negative Laplace operator
$\mathcal{A}:X=H_{0}^{1}(\Omega)\to X^{*}$ defined by

$\displaystyle\langle\mathcal{A}u,v\rangle=\int_{\Omega}\nabla u\cdot\nabla
v\;dx,\quad u,v\in X.$

The standard weak formulation of the Dirichlet problem is now
$\mathcal{A}x=f$, where $f\in X^{*}$. Since $X\neq X^{*}$ the linear operator
$\mathcal{A}$ maps $x$ out of the space and the conjugate gradient method is
not well defined. To overcome this problem, we introduce a preconditioner
$\mathcal{M}$, which is a symmetric and positive definite isomorphism mapping
$X^{*}$ to $X$. The preconditioned system
$\mathcal{M}\mathcal{A}x=\mathcal{M}f$ can then be solved by the conjugate
gradient method. We consider now the corresponding continuous linear operator
from (21):

$\displaystyle\mathcal{A}_{u,v}=I+stK^{*}K=I+st\begin{pmatrix}-\Delta&-\nabla^{*}\\\
-\nabla&I+\operatorname{\mathcal{E}}^{*}\operatorname{\mathcal{E}}\end{pmatrix}.$

The operator $\mathcal{A}_{u,v}$ is an isomorphism mapping
$X=H^{1}(\Omega)\times(H^{1}(\Omega))^{2}$ into its dual
$X^{*}=H^{1}(\Omega)^{\prime}\times(H^{1}(\Omega)^{\prime})^{2}$. The
canonical choice of a preconditioner, in the sense of
mardal2011preconditioning , is therefore given as the blockdiagonal operator

$\displaystyle\mathcal{M}_{u,v}=\begin{pmatrix}(I-st\Delta)^{-1}&0\\\
0&\left(I+st\left(I+J^{*}J\right)\right)^{-1}\end{pmatrix}$

(the inverses of the respective operators also exist in the continuous
setting, see (alt2016linfuncana, , Theorem 6.6)). To see that
$\mathcal{A}_{u,v}$ has a bounded inverse, we can check coercivity, i.e. for
some $k>0$ we have

$\displaystyle\left\langle\mathcal{A}_{u,v}(u,v),(u,v)\right\rangle$
$\displaystyle=\langle u,u\rangle_{L^{2}(\Omega)}+st\langle\nabla u,\nabla
u\rangle_{(L^{2}(\Omega))^{2}}$ $\displaystyle\quad-2st\langle\nabla
u,v\rangle_{(L^{2}(\Omega))^{2}}$ $\displaystyle\quad+(1+st)\langle
v,v\rangle_{(L^{2}(\Omega))^{2}}+st\langle\operatorname{\mathcal{E}}v,\operatorname{\mathcal{E}}v\rangle_{(L^{2}(\Omega))^{4}}$
$\displaystyle=\|u\|^{2}_{L^{2}(\Omega)}+st\|\nabla
u\|^{2}_{(L^{2}(\Omega))^{2}}-2st\langle\nabla
u,v\rangle_{(L^{2}(\Omega))^{2}}$
$\displaystyle\quad+(1+st)\|v\|^{2}_{(L^{2}(\Omega))^{2}}+st\|\operatorname{\mathcal{E}}v\|^{2}_{(L^{2}(\Omega))^{4}}$
$\displaystyle\geq\|u\|^{2}_{L^{2}(\Omega)}+st\|\nabla
u\|^{2}_{(L^{2}(\Omega))^{2}}$ $\displaystyle\quad-2st\|\nabla
u\|_{(L^{2}(\Omega))^{2}}\|v\|_{(L^{2}(\Omega))^{2}}$
$\displaystyle\quad+(1+st)\|v\|^{2}_{(L^{2}(\Omega))^{2}}+st\|\operatorname{\mathcal{E}}v\|^{2}_{(L^{2}(\Omega))^{4}}$
$\displaystyle\geq\|u\|^{2}_{L^{2}(\Omega)}+st\|\nabla
u\|^{2}_{(L^{2}(\Omega))^{2}}$ $\displaystyle\quad-2st\left(\frac{\|\nabla
u\|^{2}_{(L^{2}(\Omega))^{2}}}{2\varepsilon}+\frac{\varepsilon\|v\|^{2}_{(L^{2}(\Omega))^{2}}}{2}\right)$
$\displaystyle\quad+(1+st)\|v\|^{2}_{(L^{2}(\Omega))^{2}}+st\|\operatorname{\mathcal{E}}v\|^{2}_{(L^{2}(\Omega))^{4}}$
$\displaystyle=\|u\|^{2}_{L^{2}(\Omega)}+st(1-\frac{1}{\varepsilon})\|\nabla
u\|^{2}_{(L^{2}(\Omega))^{2}}$
$\displaystyle\quad+(1+st(1-\varepsilon))\|v\|^{2}_{(L^{2}(\Omega))^{2}}+st\|\operatorname{\mathcal{E}}v\|^{2}_{(L^{2}(\Omega))^{4}}$
$\displaystyle\geq\tilde{k}\left(\|u\|^{2}_{L^{2}(\Omega)}+\|\nabla
u\|^{2}_{(L^{2}(\Omega))^{2}}\right)$
$\displaystyle\quad+\tilde{k}\left(\|v\|^{2}_{(L^{2}(\Omega))^{2}}+\|\operatorname{\mathcal{E}}v\|^{2}_{(L^{2}(\Omega))^{4}}\right)$
$\displaystyle\geq
k\left(\|u\|_{H^{1}(\Omega)}^{2}+\|v\|^{2}_{(H^{1}(\Omega))^{2}}\right).$

The second to last estimate holds for some $\tilde{k}>0$ since for all $st>0$
we can find $\varepsilon>0$, such that
$st\left(1-\frac{1}{\varepsilon}\right)>0$ and $(1+st(1-\varepsilon))>0$. We
can chose $\tilde{k}$ as the minimum of both of them. The last estimate
follows for some $k>0$ as a consequence of Korn’s inequality ciarlet2013linear
.

The corresponding continuous linear operator from (23) is

$\displaystyle\mathcal{A}_{v}=I+st\operatorname{\mathcal{E}}^{*}\operatorname{\mathcal{E}}.$

The operator $\mathcal{A}_{v}$ is an isomorphism mapping
$X=(H^{1}(\Omega))^{2}$ into its dual $X^{*}=(H^{1}(\Omega)^{\prime})^{2}$.
The canonical choice of a preconditioner is therefore given as the
blockdiagonal operator

$\displaystyle\mathcal{M}_{v}=(I+stJ^{*}J)^{-1}.$

The corresponding continuous linear operator from (22) is

$\displaystyle\mathcal{A}_{u,w}=I+stK^{*}K=I+st\begin{pmatrix}H^{*}H&-H^{*}\operatorname{\mathcal{E}}\\\
-\operatorname{\mathcal{E}}^{*}H&\operatorname{\mathcal{E}}^{*}\operatorname{\mathcal{E}}\end{pmatrix},$

where $H=\operatorname{\mathcal{E}}(\nabla)$ is the Hessian matrix. In our
experiments we chose the discrete version of the blockdiagonal operator

$\displaystyle\mathcal{M}_{u,w}=\begin{pmatrix}(I+stH^{*}H)^{-1}&0\\\
0&\left(I+stJ^{*}J\right)^{-1}\end{pmatrix}$

for preconditioning. It gives good numerical results, but is not as nicely
accompanied by the theory as the previous operator. To obtain fast algorithms
the inverses in the preconditioners can be well approximated by the incomplete
Cholesky decomposition.111We used the MATLAB function ichol with default
parameters in the experiments. Varying the parameters did not lead to
significantly different results. For the denoising problems with the primal
variables $u$ und $v$ we use the preconditioner

$\displaystyle C_{u,v}=\begin{pmatrix}\operatorname{ichol}(I-st\Delta)&0\\\
0&\operatorname{ichol}\left(I+st\left(I+J^{T}J\right)\right)\end{pmatrix}$
(24)

and for the denoising problems with the primal variables $u$ and $w$ we use
the preconditioner

$\displaystyle C_{u,w}=\begin{pmatrix}\operatorname{ichol}(I+H^{T}H)&0\\\
0&\operatorname{ichol}\left(I+stJ^{T}J\right)\end{pmatrix}.$ (25)

For problem (6) we use the preconditioner

$\displaystyle C_{v}=\operatorname{ichol}(I+stJ^{T}J).$ (26)

Here the preconditioners are given in the form $M^{-1}=C^{T}C$.

##  5 Experiments

In the previous sections we introduced several different models for
variational image denoising and also described two algorithms, each applicable
to each model. This leaves us with a large number of parameters that have to
be chosen. In this section, we try to illustrate the effects of these
parameters and to give some guidance on how these parameters should be chosen.
To that end, we divide the set of parameters into two groups:

Problem parameters:

These are the parameters of the model itself. For example in the case of TGV
denoising (10), the two problem parameters are the two regularization
parameters $\alpha_{0}$ and $\alpha_{1}$ while for MTGV (9) we have the
parameters $\alpha$ and $\delta_{1}$ and the DGTGV method consisting of (6)
and (5) also has the two parameters $\alpha$ and $\delta_{1}$.

Algorithmic parameters:

These are parameters, that influence only the algorithm, but not the
theoretical minimizers. These can be for example: One or more step-sizes, the
relaxation parameter, the stopping criterion (e.g. a tolerance for the duality
gap), the parameters of the CG iteration, or the preconditioner.

The problem parameters influence the quality of the denoising, while the
algoritmic parameters influence the performance (or speed) of the method.
Moreover, there is a trade-off between speed and quality: if the algorithm is
stopped too early, the minimizer may not be approximated well. Note, however,
that sometimes early stopping may increase reconstruction quality (which often
indicates that the model can be improved), but we shall not deal with this
question here but rather focus on the analysis of the problem and algorithmic
parameters separately.

###  5.1 Optimal constants in DGTV, DGTGV and MTGV

The two-stage denoising method DGTV from section 2.2 use the parameter
$\delta_{2}$ which we motivated to be
$\delta_{2}=c{\|{|\operatorname{\nabla}u^{0}|}\|}_{2}$ in section 2.4. We
calculated the optimal constants $c$ for various images with different noise-
levels, cf. figure 1. As expected, all $c\geq 1$ lead to the same result (in
this case $v=0$ is a feasible solution and optimal and therefor, the two-stage
method becomes pure TV-denoising). If we choose $c<1$, we transfer a bit of
structure of the input image into the gradient as an additional information
for the image denoising step. Hence, $c\approx 0.99$ seems like a sensible
choice.

Figure 1: Optimal values for various images and different noise-level for the
DGTV method.

A similar experiment for the parameter $\alpha$ in the DGTGV method from
section 2.2 revealed that all these optimal values are close to $1$, cf.
figure 2. For $\alpha>1$, the change in the norm distance is minimal, hence,
the denoised image will be similar to the denoised image with $\alpha=1$.
Smaller values $\alpha<1$ lead to worse reconstructions. Therefore, we use a
default value of $\alpha=1$. For the MTGV method from section 2.3 we report
the results of the optimization of the parameter $\alpha$ in figure 2 and we
see that values around $\alpha=2$ seem optimal (while the variance is larger
than for DGTGV).

(a) Optimal $\alpha$ for various images with different noise-levels in the
DGTGV method.

(b) Optimal $\alpha$ for various images with different noise-levels in the
MTGV method. Figure 2: Optimal $\alpha$ values according to relative norm
distance between $u$ and the original image $u^{\dagger}$ for various images
and different noise-levels.

###  5.2 Quality and runtime of MTGV and DGTGV

Table 1 collects results on the reconstruction quality (measured in PSNR) and
the runtime of the two-stage method DGTGV and the MTGV method. In table 1(a)
we compared the two-stage method DGTGV and the MTGV denoising method with the
best possible $\alpha$ (i.e. we calculated an optimal value of this parameter
for each image, noise-level and method) and the default $\delta_{1}$ as
estimated from the images. In table 1(b) we compared both methods with the
default $\alpha$ values, i.e. $\alpha=1$ for DGTGV and $\alpha=2$ for MTGV
(again with $\delta_{1}$ estimated from the image).

In both comparisons the PSNR values differ only slightly from each other
showing that the default values always lead to results close to the optimal
ones. Moreover, the two-stage method DGTGV consistently gives slightly lower
PSNR values than MTGV. The difference in PSNR is so small, that the resulting
images are very similar to each other (cf. figures 3, 4, 5 for some examples).

Image(noise factor)
---
affine(0.05)
affine(0.1)
affine(0.25)
eye(0.05)
eye(0.1)
eye(0.25)
cameraman(0.05)
cameraman(0.1)
cameraman(0.25)
moonsurface(0.05)
moonsurface(0.1)
moonsurface(0.25)
barbara(0.05)
barbara(0.1)
barbara(0.25)

DGTGV | MTGV
---|---
best $\alpha$
37.26 | 38.30
32.58 | 33.97
27.74 | 28.41
31.34 | 31.64
28.94 | 29.37
26.26 | 26.98
30.57 | 30.79
27.14 | 27.32
23.20 | 23.34
30.69 | 30.87
28.46 | 28.73
26.08 | 26.38
28.35 | 28.66
24.91 | 25.14
22.34 | 22.48
(a) PSNR values for best possible $\alpha$ according to each method.

DGTGV | MTGV
---|---
default $\alpha$
37.07 | 38.25
32.66 | 33.65
26.80 | 27.87
31.32 | 31.49
28.95 | 29.28
26.09 | 26.89
30.53 | 30.78
27.09 | 27.32
23.15 | 23.26
30.68 | 30.77
28.46 | 28.65
26.04 | 26.36
28.35 | 28.49
24.91 | 25.06
22.34 | 22.46
(b) PSNR values for default $\alpha$ values according to each method.

DGTGV | DGTGV | DGTGV | MTGV
---|---|---|---
(grad.) | (img.) | (all) |
0.12 | 0.06 | 0.17 | 0.20
0.15 | 0.05 | 0.19 | 0.25
0.27 | 0.04 | 0.31 | 0.39
0.14 | 0.16 | 0.30 | 0.91
0.23 | 0.18 | 0.41 | 0.93
0.67 | 0.16 | 0.83 | 1.53
0.30 | 0.15 | 0.46 | 1.28
0.38 | 0.16 | 0.54 | 1.38
0.73 | 0.16 | 0.89 | 1.81
0.14 | 0.16 | 0.30 | 0.84
0.28 | 0.17 | 0.45 | 0.93
0.65 | 0.15 | 0.80 | 1.52
1.19 | 0.58 | 1.77 | 10.98
1.03 | 0.60 | 1.63 | 7.46
1.94 | 0.75 | 2.69 | 7.33
(c) Time in seconds (Chambolle-Pock).

DGTGV | MTGV
---|---
0.18 | 0.28
0.20 | 0.30
0.44 | 0.34
0.46 | 1.21
0.83 | 1.08
1.66 | 1.26
0.80 | 1.79
1.03 | 1.55
1.51 | 1.43
0.52 | 1.17
0.72 | 1.08
1.46 | 1.26
3.40 | 9.67
4.66 | 7.53
8.03 | 6.04
(d) Time in seconds (Douglas-Rachford).

Table 1: 1(a), 1(b): PSNR values of DGTGV and MTGV methods with best possible
$\alpha$ value and default values for each method; 1(c): Time in seconds for
DGTGV and MTGV methods implemented with Chambolle-Pock, for the DGTGV both
steps also separately; 1(d): Time in seconds for both methods implemented with
Douglas Rachford.

In tables 1(c) and 1(d) we compare different run-times: First, we used
Chabolle-Pock’s primal dual method for both steps of the DGTGV method
separately, i.e. the gradient denoising as a pre-step, after that the actual
image denoising and MTGV in table 1(c). It can be seen that the two-stage
method DGTGV is faster (usually by a factor of 2 or 3). In all cases we used a
tolerance value of $10^{-3}$ for the relative primal-dual gap as a stopping
criterion. Table 1(d) also shows runtimes for inexact preconditioned Douglas-
Rachford method. Here we set the tolerance for the primal-dual gap to
$10^{-2}$ since, this already gave better or comparable PSRN values (a
phenomenon which we can not explain theoretically). However, even with this
larger tolerance, the DR-method is only faster for large noise level and the
MTGV method.

(a) $u^{\dagger}$

(b) $u_{0}$

(c) DGTGV

(d) MTGV

Figure 3: Comparison of DGTGV and MTGV method of affine image with default
$\alpha$ value; noise factor $0.1$.

(a) $u^{\dagger}$

(b) $u_{0}$

(c) DGTGV

(d) MTGV

Figure 4: Comparison of DGTGV and MTGV method of affine image with best
possible $\alpha$ value; noise factor $0.25$.

(a) $u^{\dagger}$

(b) $u_{0}$

(c) DGTGV

(d) MTGV

Figure 5: Comparison of DGTGV and MTGV method with best possible $\alpha$
value; noise factor $0.1$.

###  5.3 Comparison of different formulations

The TGV minimization problem (10), the MTGV minimization problem (9) and the
MTGV minimization problem with variable change $\eqref{eq:mtgvvarchange}$ have
the same solution parameter $u$ if we choose
$\delta_{1}={\|{u_{0}-u^{\operatorname{TGV}}}\|}_{2}$ and $\alpha=2$, where
$u^{\operatorname{TGV}}$ is the solution of (10). For the first two problems,
we used the algorithms of Douglas-Rachford and Chambolle-Pock. The latter one
is an $\mathcal{O}(1/k)$ primal-dual algorithm with constant step sizes (see
chambolle2011first ), namely

$\displaystyle\begin{split}y^{n+1}&=(I+\sigma\partial G^{*})^{-1}(y^{n}+\sigma
K(2x^{n}-x^{n-1}))\\\ x^{n+1}&=(I+\tau\partial F)^{-1}(x^{n}-\tau
K^{*}y^{n+1})\end{split}$

For MTGV with variable change we only show the results for the Douglas-
Rachford algorithm, since the Chambolle-Pock algorithm was much slower for
this problem. In figure 6 we compare iteration number and time needed to
obtain the desired accuracy of the image. Since the duality gap is not
suitable to compare different minimization problems, the tolerance is given by
$\frac{{\|{u-u^{\operatorname{TGV}}}\|}_{2}}{{\|{u^{\operatorname{TGV}}}\|}_{2}}$.
The reference value $u^{\operatorname{TGV}}$ is obtained by solving (10) with
1000000 iterations of Chambolle-Pock ($\tau=0.004$,
$\sigma=\frac{1}{\tau\|K\|^{2}}$, $\|K\|^{2}=12$). The algorithms are tested
with the image eye (256x256) corrupted by Gaussian white noise of mean $0$ and
variance $0.1$. The stepsizes of the algorithms are chosen by trial and error:

  * •

CP TGV: $\tau=0.008$, $\sigma=\frac{1}{\tau\|K\|^{2}}$, $\|K\|^{2}=12$,

  * •

CP MTGV: $\tau=0.004$, $\sigma=\frac{1}{\tau\|K\|^{2}}$, $\|K\|^{2}=12$,

  * •

DR TGV: $s=60$, $t=0.28$,

  * •

DR MTGV: $s=60$, $t=0.1$,

  * •

DR MTGV (var. change): $s=60$, $t=0.04$.

The optimal stepsize depends on the accuracy needed. The Douglas-Rachford
algorithm is used in an inexact manner. The linear operator is approximated by
two iterations of the preconditioned conjugate gradients method, where the
preconditioners are given by $\eqref{preconditioner_uv}$ and
$\eqref{preconditioner_uw}$. In figure 6 we can see that the algorithms for
TGV and MTGV are competitive. The variable change leads to much slower
algorithms.

Figure 6: Runtime TGV vs. MTGV.

###  5.4 Inexactness for the Douglas-Rachford method

The Douglas-Rachford iteration in general allows inexact evaluation of the
operators as long as the error stays summable, see e.g. combettes2004averaged
. We made experiments with MTGV, using a few iterates of the conjugate
gradient method with and without preconditioning (see figure 7).The algorithms
are tested with the image eye (256x256) corrupted by Gaussian white noise of
mean $0$ and variance $0.1$. The stepsizes of the Douglas-Rachford iteration
are chosen by trial and error: $s=120$, $t=0.1$. The Douglas-Rachford
iteration is very slow for one or two iterations of CG without
preconditioning. It does not even converge for three iterations of CG without
preconditioning. Using the preconditioners as proposed in section 4.2 we
obtain very good convergence of Douglas-Rachford for any number of iterations
of PCG. In figure 7 we can see that preconditioning is crucial for the method
to converge. After 300 iterations all inexact algorithms have caught up with
the exact one regarding the iteration count, while only one or two iterations
perform best regarding the computational time.

(a) Convergence of MTGV without preconditioning according to number of
iterations.

(b) Convergence of MTGV with preconditioning according to number of
iterations.

(c) Convergence of MTGV without preconditioning according to time in seconds.

(d) Convergence of MTGV with preconditioning according to time in seconds.

Figure 7: A comparisson of the MTGV method with and without preconditioning
according to iterations and time.

##  6 Conclusion

We investigated variants of variational denoising methods using total-
variation penalties and an estimate of the image gradient. First, this
provides a natural and, at least to us, new interpretation of the successful
TGV method. The reformulation with a constraint for the discrepancy term
${\|{u-u^{0}}\|}_{2}$ (which also works for all other norms) together with our
empirical observations allows for variational denoising methods that are
basically parameter free and we mainly investigated the methods DGTGV and
MTGV. Since the gradient of $u^{0}$ is even more noisy that $u^{0}$ itself,
our experiments in Section 5 show, that it is still useful to use a denoised
version of $\nabla u^{0}$ as estimate for $\nabla u^{{\dagger}}$. Indeed, we
obtained that the two-stage method DGTGV sacrifices only little denoising
performance for a substantial gain in speed. Put differently, solving the two
denoising problems for the gradient and the image is significantly easier,
than solving the combined MTGV problem and the denoising result is still good.
This qualifies the two-stage DGTGV as an alternative to MTGV (and hence, TGV)
as a denoising method.

Another part of the investigation involved the Douglas-Rachford method for
these problems. Here we could derive natural preconditioners for the linear
sub-problems and show that they greatly improve the overall speed of the
method, especially in the inexact case. However, the overall runtime was only
better than the simpler primal-dual method by Chambolle and Pock in case of
large noise.

##  Appendix A Appendix

###  A.1 Discretization

We equip the space of $M\times N$ images with the inner product and induced
norm

$\displaystyle\langle
u,v\rangle=\sum_{i=1}^{M}\sum_{j=1}^{N}u_{i,j}v_{i,j},\quad\|u\|_{2}=\left(\sum_{i=1}^{M}\sum_{j=1}^{N}u_{i,j}^{2}\right)^{\frac{1}{2}}.$

For $u\in\mathbf{R}^{M\times N}$ we define the discrete partial forward
derivative (with constant boundary extension $u_{M+1,j}=u_{M,j}$ and
$u_{i,N+1}=u_{i,N}$ as

$\displaystyle(\partial_{1}u)_{i,j}$ $\displaystyle=u_{i+1,j}-u_{i,j},$
$\displaystyle(\partial_{2}u)_{i,j}$ $\displaystyle=u_{i,j+1}-u_{i,j}.$

The discrete gradient $\nabla:\mathbf{R}^{M\times N}\to\mathbf{R}^{M\times
N\times 2}$ is defined by

$\displaystyle(\nabla u)_{i,j,k}=(\partial_{k}u)_{i,j}.$

The symmetrized gradient $\mathcal{E}$ maps from $\mathbf{R}^{M\times N\times
2}$ to $\mathbf{R}^{M\times N\times 4}$. For simplification of notation, the 4
blocks are written in one plane:

$\displaystyle\mathcal{E}v$ $\displaystyle=\frac{1}{2}\left(\nabla v+(\nabla
v)^{T}\right)$
$\displaystyle=\begin{pmatrix}\partial_{1}v_{1}&\frac{1}{2}(\partial_{1}v_{2}+\partial_{2}v_{1})\\\
\frac{1}{2}(\partial_{1}v_{2}+\partial_{2}v_{1})&\partial_{2}v_{2}\end{pmatrix}.$

The norm $\||\cdot|\|_{1}$ in the space $\mathbf{R}^{M\times N\times K}$
reflects that for $v\in\mathbf{R}^{M\times N\times K}$ we consider $v_{i,j}$
as a vector in $\mathbf{R}^{K}$ on which we use the Euclidean norm:

$\displaystyle\||v|\|_{1}\coloneqq\sum_{i=1}^{M}\sum_{j=1}^{N}|v_{i,j}|\quad\text{with
}|v_{i,j}|\coloneqq\left(\sum_{k=1}^{K}v_{i,j,k}^{2}\right)^{\frac{1}{2}}.$

The discrete divergence is the negative adjoint of $\nabla$, i.e. the unique
linear mapping $\operatorname{div}:\mathbf{R}^{M\times N\times
2}\to\mathbf{R}^{M\times N}$, which satisfies

$\displaystyle\langle\nabla u,v\rangle=-\langle
u,\operatorname{div}v\rangle,\quad\text{for all }u,v.$

The adjoint of the symmetriced gradient is the unique linear mapping
$\operatorname{\mathcal{E}}^{*}:\mathbf{R}^{M\times N\times
4}\to\mathbf{R}^{M\times N\times 2}$, which satisfies

$\displaystyle\langle\operatorname{\mathcal{E}}v,p\rangle=\langle
v,\operatorname{\mathcal{E}}^{*}p\rangle,\quad\text{for all }v,p.$

###  A.2 Prox operators and duality gaps for considered problems

In order to calculate experiments with the methods proposed in the previous
sections, in this section we will state all primal and dual functionals
according to a general optimization problem $\min_{x}F(x)+G(Kx)$ along with a
study of the corresponding primal-dual gaps and possibilities to ensure
feasibility of the iterates throughout the program by reformulation of the
problems by introducing a substitution variable. We also give the proximal
operators needed for the Chambolle-Pock and Douglas-Rachford algorithm.

####  A.2.1 DGTV

In section 2 we formulated a two- staged denoising method in two ways. First,
as a constrained version (4) and (5). In this formulation we get the primal
functionals for the first problem (4)

$\displaystyle\begin{split}F(v)&=\operatorname{\mathcal{I}}_{{\|{|\operatorname{\nabla}u_{0}-\,\cdot\,|}\|}_{1}\leq\delta_{2}}(v),\\\
G(\psi)&={\|{|\psi|}\|}_{1}\end{split}$ (27)

with operator $K=\operatorname{\mathcal{E}}$. The dual problems, in general
written as

$\displaystyle\max_{y}-F^{*}(-K^{*}y)-G^{*}(y)$

are

$\displaystyle\begin{split}F^{*}(t)&=\delta_{2}{\|{|t|}\|}_{\infty}+\langle
t,\operatorname{\nabla}u_{0}\rangle,\\\
G^{*}(q)&=\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(q)\end{split}$ (28)

with operator $K^{*}=\operatorname{\mathcal{E}}^{*}$. The primal-dual gap
writes as

$\displaystyle\begin{split}\operatorname{gap}&{}_{\operatorname{DGTV}}^{(1)}(v,q)=\operatorname{\mathcal{I}}_{{\|{|\operatorname{\nabla}u_{0}-\,\cdot\,|}\|}_{1}\leq
1}(v)+{\|{|\operatorname{\mathcal{E}}v|}\|}_{1}\\\
&+\delta_{2}{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}-\langle\operatorname{\mathcal{E}}^{*}q,\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(q).\end{split}$ (29)

The proximal operators are

$\displaystyle\begin{split}\operatorname{prox}_{\tau
F}(v)&=\operatorname{proj}_{{\|{|\operatorname{\nabla}u_{0}-\,\cdot\,|}\|}_{1}\leq\delta_{2}}(v),\\\
\operatorname{prox}_{\sigma
G^{*}}(q)&=\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(q).\end{split}$ (30)

For the second problem within DGTV, we have the denoising problem of the image
with respect to the denoised gradient $\widehat{v}$ as output of the previous
problem, cf. problem (5). There, the primal functionals with
$K=\operatorname{\nabla}$ are

$\displaystyle\begin{split}F(u)&=\operatorname{\mathcal{I}}_{{\|{u_{0}-\,\cdot\,}\|}_{2}\leq\delta_{1}}(u)\\\
G(\phi)&={\|{|\phi-\widehat{v}|}\|}_{1}\end{split}$ (31)

and the corresponding dual functionals write as

$\displaystyle\begin{split}F^{*}(s)&=\delta_{1}{\|{s}\|}_{2}+\langle
u_{0},s\rangle\\\
G^{*}(p)&=\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(p)+\langle\widehat{v},p\rangle.\end{split}$ (32)

Hence, the primal-dual gap for this problem is

$\displaystyle\begin{split}\operatorname{gap}_{\operatorname{DGTGV}}^{(2)}&(u,p)=\operatorname{\mathcal{I}}_{{\|{u_{0}-\,\cdot\,}\|}_{2}\leq\delta_{1}}(u)+{\|{|\operatorname{\nabla}u-\widehat{v}|}\|}_{1}\\\
&-\langle
u_{0},\operatorname{\nabla}^{*}p\rangle+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(p)+\langle\widehat{v},p\rangle.\end{split}$ (33)

The proximal operators are given by

$\displaystyle\begin{split}\operatorname{prox}_{\tau
F}(u)&=\operatorname{proj}_{{\|{u_{0}-\,\cdot\,}\|}_{2}\leq\delta_{1}}(u)\\\
\operatorname{prox}_{\sigma
G^{*}}(p)&=\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(p-\sigma\widehat{v}).\end{split}$ (34)

####  A.2.2 DGTGV

We reformulate problem (6) by using a substitution
$w=\operatorname{\nabla}u_{0}-v$, also considered in section 3, where we
calculated another duality gap. Hence, we get the primal functionals for the
gradient denoising problem with operator $K=\operatorname{\mathcal{E}}$ as

$\displaystyle\begin{split}F(w)&={\|{|w|}\|}_{1}\\\
G(\psi)&=\alpha{\|{|\operatorname{\mathcal{E}}\operatorname{\nabla}u_{0}-w|}\|}_{1}.\end{split}$
(35)

Therefore, the dual functionals write as

$\displaystyle\begin{split}F^{*}(t)&=\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(t)\\\
G^{*}(q)&=\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha}(q)+\langle\operatorname{\mathcal{E}}\operatorname{\nabla}u_{0},q\rangle.\end{split}$
(36)

With that the primal-dual gap is

$\displaystyle\begin{split}\operatorname{gap}^{(1)}&{}_{\operatorname{DGTGV}}(w,q)={\|{|w|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}\operatorname{\nabla}u_{0}-w|}\|}_{1}\\\
&+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(\operatorname{\mathcal{E}}^{*}q)+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha}(q)+\langle\operatorname{\mathcal{E}}\operatorname{\nabla}u_{0},q\rangle.\end{split}$
(37)

The proximal operators are, with Moreau’s identity for the first one,

$\displaystyle\begin{split}\operatorname{prox}_{\tau
F}(w)&=w-\tau\operatorname{prox}_{\tau^{-1}F^{*}}(\tau^{-1}w)\\\
&=w-\tau\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(\tau^{-1}w)\\\ \operatorname{prox}_{\sigma
G^{*}}(q)&=\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha}(q-\tau\operatorname{\mathcal{E}}\operatorname{\nabla}u_{0}).\end{split}$
(38)

For the second problem, we already derived all functionals, gaps and proximal
operators in subsection A.2.1, equations (31)–(34).

####  A.2.3 CTGV

The Morozov type constrained total generalized variation denoising problem was
formulated in section 2.3 (cf. (12)). For this formulation the primal
functionals are

$\displaystyle\begin{split}F(u,w)&=\operatorname{\mathcal{I}}_{{\|{u-u_{0}}\|}_{2}\leq\delta_{1}}(u)+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{1}\leq\delta_{2}}(w),\\\
G(\phi)&={\|{|\phi|}\|}_{1}\end{split}$ (39)

with block operator

$\displaystyle
K=\operatorname{\mathcal{E}}\begin{pmatrix}\operatorname{\nabla}&-\operatorname{Id}\end{pmatrix}.$
(40)

The corresponding dual functionals are

$\displaystyle\begin{split}F^{*}(s,t)&=\delta_{1}{\|{s}\|}_{2}+\langle
s,u_{0}\rangle+\delta_{2}{\|{|t|}\|}_{\infty},\\\
G^{*}(q)&=\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(q)\end{split}$ (41)

with dual block operator

$\displaystyle K^{*}=\begin{pmatrix}\operatorname{\nabla}^{*}\\\
-\operatorname{Id}\end{pmatrix}\operatorname{\mathcal{E}}^{*}.$ (42)

The proximal operators are accordingly given by

$\displaystyle\begin{split}&\operatorname{prox}_{\tau
F}(u,w)=\begin{pmatrix}\operatorname{proj}_{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq\delta_{1}}(u)\\\
\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{1}\leq\delta_{2}}(w)\end{pmatrix},\\\
&\operatorname{prox}_{\sigma
G^{*}}(q)=\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(q)\end{split}$ (43)

and the primal-dual gap writes as

$\displaystyle\begin{split}\operatorname{gap}&{}_{\operatorname{CTGV}}(u,w,q)=\\\
&{\|{|\operatorname{\mathcal{E}}(\operatorname{\nabla}u-w)|}\|}_{1}+\operatorname{\mathcal{I}}_{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq\delta_{1}}(u)+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{1}\leq\delta_{2}}(w)\\\
&+\delta_{1}{\|{\operatorname{\nabla}^{*}\operatorname{\mathcal{E}}^{*}q}\|}_{2}-\langle\operatorname{\nabla}^{*}\operatorname{\mathcal{E}}^{*}q,u_{0}\rangle\\\
&+\delta_{2}{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(q).\end{split}$ (44)

####  A.2.4 MTGV

In section 2.3 we defined the MTGV optimization problem (9) as a sort of mixed
version between the TGV (10) and the CTGV (12) problems. Thus, the primal
functionals are

$\displaystyle\begin{split}F(u,v)&=\operatorname{\mathcal{I}}_{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq\delta_{1}}(u),\\\
G(\phi,\psi)&={\|{|\phi|}\|}_{1}+\alpha{\|{|\psi|}\|}_{1}\end{split}$ (45)

with the block operator

$\displaystyle K=\begin{pmatrix}\operatorname{\nabla}&-\operatorname{Id}\\\
0&\operatorname{\mathcal{E}}\end{pmatrix}.$ (46)

Accordingly, the dual functionals are

$\displaystyle\begin{split}F^{*}(s,t)&=\delta_{1}{\|{s}\|}_{2}+\langle
s,u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{0\\}}(t),\\\
G^{*}(p,q)&=\delta_{2}{\|{|p|}\|}_{\infty}+\operatorname{\mathcal{I}}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(q)\end{split}$ (47)

and the dual operator

$\displaystyle K=\begin{pmatrix}\operatorname{\nabla}^{*}&0\\\
-\operatorname{Id}&\operatorname{\mathcal{E}}^{*}\end{pmatrix}.$ (48)

The proximal operators are given by

$\displaystyle\operatorname{prox}_{sF}(u,v)=(\operatorname{proj}_{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq\delta_{1}}(u),v)$
$\displaystyle\operatorname{prox}_{tG^{*}}(p,q)=(\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1}(p),\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha}(q))$

and the gap function is given by

$\displaystyle\operatorname{gap}$ ${}_{\operatorname{MTGV}}(u,v,p,q)=$
$\displaystyle{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\operatorname{\mathcal{I}}_{\\{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq\delta_{1}\\}}(u)$
$\displaystyle+\delta_{1}{\|{\operatorname{\nabla}^{*}p}\|}_{2}-\langle
p,\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{0\\}}(p-\operatorname{\mathcal{E}}^{*}q)$
$\displaystyle+\operatorname{\mathcal{I}}_{\\{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha\\}}(p)+\operatorname{\mathcal{I}}_{\\{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1\\}}(q).$

To circumvent the feasibility problem, as introduced in 3 one can use the
modified gap function

$\displaystyle\operatorname{gap}$ ${}_{\operatorname{MTGV}}(u,v,\tilde{q})=$
$\displaystyle\leq{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}+\operatorname{\mathcal{I}}_{\\{{\|{\,\cdot\,-u_{0}}\|}_{2}\leq\delta_{1}\\}}(u)$
$\displaystyle+\delta_{1}{\|{\operatorname{\nabla}^{*}(\operatorname{\mathcal{E}}^{*}\tilde{q})}\|}_{2}-\langle\operatorname{\mathcal{E}}^{*}\tilde{q},\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{{\|{|\,\cdot\,|}\|}_{\infty}\leq
1\\}}(\tilde{q}),$

where
$\tilde{q}\coloneqq\frac{q}{\max\left(1,\frac{{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}}{\alpha}\right)}$.

####  A.2.5 TGV:

For TGV (10) we have the primal functionals

$\displaystyle F(u,v)=\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}$ $\displaystyle
G(\phi,\psi)=\alpha_{1}{\|{|\phi|}\|}_{1}+\alpha_{0}{\|{|\psi|}\|}_{1},$

with the same block operator (46). The corresponding dual functionals are

$\displaystyle F^{*}(s,t)=\frac{1}{2}{\|{s}\|}_{2}^{2}+\langle
s,u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{0\\}}(t),$ $\displaystyle
G^{*}(p,q)=\operatorname{\mathcal{I}}_{\\{{\|{\cdot}\|}_{\infty}\leq\alpha_{1}\\}}(p)+\operatorname{\mathcal{I}}_{\\{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha_{0}\\}}(q)$

and the dual operator is (48). The proximal operators are given by

$\displaystyle\operatorname{prox}_{tF}(u,v)=\left(\frac{u+tu_{0}}{1+t},v\right)$
$\displaystyle\operatorname{prox}_{tG^{*}}(p,q)=(\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha_{1}}(p),\operatorname{proj}_{{\|{|\,\cdot\,|}\|}_{\infty}\leq\alpha_{0}}(q))$

and the gap function is given by

$\displaystyle\operatorname{gap}$
${}_{\operatorname{TGV}}(u,v,p,q)=\alpha_{1}{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha_{0}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}$
$\displaystyle+\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}+\frac{1}{2}{\|{\operatorname{\nabla}^{*}p}\|}_{2}^{2}-\langle
p,\operatorname{\nabla}u_{0}\rangle+\operatorname{\mathcal{I}}_{\\{0\\}}(p-\operatorname{\mathcal{E}}^{*}q)$
$\displaystyle+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{0}\\}}(p)+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{1}\\}}(q).$

To circumvent the feasability problem, as introduced in 3 one can use the
modified gap function

$\displaystyle\operatorname{gap}$
${}_{\operatorname{TGV}}(u,v,\tilde{q})=\alpha_{1}{\|{|\operatorname{\nabla}u-v|}\|}_{1}+\alpha_{0}{\|{|\operatorname{\mathcal{E}}(v)|}\|}_{1}$
$\displaystyle+\frac{1}{2}{\|{u-u_{0}}\|}_{2}^{2}+\frac{1}{2}{\|{\operatorname{\nabla}^{*}(\operatorname{\mathcal{E}}^{*}\tilde{q})}\|}_{2}^{2}-\langle\operatorname{\mathcal{E}}^{*}\tilde{q},\operatorname{\nabla}u_{0}\rangle$
$\displaystyle+\operatorname{\mathcal{I}}_{\\{{\|{|\cdot|}\|}_{\infty}\leq\alpha_{1}\\}}(\tilde{q}).$

where
$\tilde{q}\coloneqq\frac{q}{\max\left(1,\frac{{\|{|\operatorname{\mathcal{E}}^{*}q|}\|}_{\infty}}{\alpha}\right)}$.

###  A.3 Projections

The projections used for the algorithms are

$\displaystyle\operatorname{proj}_{{\|{\cdot-u_{0}}\|}_{2}\leq\delta}(u)$
$\displaystyle=u_{0}+\max\left(\frac{\delta}{{\|{u-u_{0}}\|}_{2}},1\right)\cdot(u-u_{0}),$
$\displaystyle\operatorname{proj}_{{\|{|\cdot|}\|}_{\infty}\leq\delta}(v)$
$\displaystyle=\frac{v}{\max(1,\frac{|v|}{\delta})},$

with $|v|=\left(\sum_{k=1}^{d}v_{k}^{2}\right)^{\frac{1}{2}}$ in the pointwise
sense.

The idea how to project onto an mixed norm ball, i.e.
${\|{|\,\cdot\,|}\|}_{1}$, can be found with in Songsiri2011 . There the
author developed an algorithm to project onto an $l_{1}$-norm ball and after
that onto a sum of $l_{2}$-norm balls. The projection itself cannot be stated
in a simple closed form.

## References

  * (1) Gabriel Acosta and Ricardo G Durán.  An optimal Poincaré inequality in $L^{1}$ for convex domains.  Proceedings of the American Mathematical Society, 132(1):195–202, 2004.
  * (2) Hans Wilhelm Alt.  Linear functional analysis.  Universitext. Springer-Verlag London, Ltd., London, 2016.  An application-oriented introduction, Translated from the German edition by Robert Nürnberg.
  * (3) Jean Bourgain and Haïm Brezis.  On the equation $\mathrm{div}{Y}=f$ and application to control of phases.  Journal of the American Mathematical Society, 16(2):393–426, 2003\.
  * (4) Kristian Bredies, Karl Kunisch, and Thomas Pock.  Total generalized variation.  SIAM Journal on Imaging Sciences, 3(3):492–526, 2010.
  * (5) Kristian Bredies, Karl Kunisch, and Tuomo Valkonen.  Properties of $L^{1}-TGV^{2}$: The one-dimensional case.  Journal of Mathematical Analysis and Applications, 398(1):438–454, 2013.
  * (6) Kristian Bredies and Hong Peng Sun.  Preconditioned Douglas–Rachford algorithms for Tv- and TGV-regularized variational imaging problems.  Journal of Mathematical Imaging and Vision, 52(3):317–344, 2015\.
  * (7) 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, 2011\.
  * (8) Philippe G Ciarlet.  Linear and nonlinear functional analysis with applications, volume 130.  SIAM, 2013.
  * (9) Patrick L. Combettes.  Solving monotone inclusions via compositions of nonexpansive averaged operators.  Optimization, 53(5-6):475–504, 2004.
  * (10) Jim Douglas, Jr. and H. H. Rachford, Jr.  On the numerical solution of heat conduction problems in two and three space variables.  Trans. Amer. Math. Soc., 82:421–439, 1956.
  * (11) Jonathan Eckstein and Dimitri P Bertsekas.  On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators.  Mathematical Programming, 55(1-3):293–318, 1992.
  * (12) Ivar Ekeland and Roger Temam.  Convex analysis and variational problems.  SIAM, 1976.
  * (13) Florian Knoll, Kristian Bredies, Thomas Pock, and Rudolf Stollberger.  Second order total generalized variation (TGV) for MRI.  Magnetic Resonance in Medicine, 65(2):480–491, 2011.
  * (14) Birgit Komander, Dirk Lorenz, Marc Fischer, Marcus Petz, and Rainer Tutsch.  Data fusion of surface normals and point coordinates for deflectometric measurements.  J. Sens. Sens. Syst., pages 281–290, 2014.
  * (15) Birgit Komander and Dirk A. Lorenz.  Denoising of image gradients and constrained total generalized variation.  In François Lauze, Yiqiu Dong, and Anders Bjorholm Dahl, editors, Scale Space and Variational Methods in Computer Vision: 6th International Conference, SSVM 2017, Kolding, Denmark, June 4-8, 2017, Proceedings, pages 435–446, Cham, 2017. Springer International Publishing.
  * (16) P.-L. Lions and B. Mercier.  Splitting algorithms for the sum of two nonlinear operators.  SIAM J. Numer. Anal., 16(6):964–979, 1979.
  * (17) Xinhao Liu, Masayuki Tanaka, and Masatoshi Okutomi.  Noise level estimation using weak textured patches of a single noisy image.  In Image Processing (ICIP), 2012 19th IEEE International Conference on, pages 665–668. IEEE, 2012.
  * (18) Xinhao Liu, Masayuki Tanaka, and Masatoshi Okutomi.  Single-image noise level estimation for blind denoising.  IEEE transactions on image processing, 22(12):5226–5237, 2013.
  * (19) Dirk A. Lorenz and Nadja Worliczek.  Necessary conditions for variational regularization schemes.  Inverse Problems, 29:075016pp, 2013.
  * (20) M. Lysaker, S. Osher, and Xue-Cheng Tai.  Noise removal using smoothed normals and surface fitting.  IEEE Trans. Img. Proc., 13(10):1345–1357, October 2004.
  * (21) Kent-Andre Mardal and Ragnar Winther.  Preconditioning discretizations of systems of partial differential equations.  Numer. Linear Algebra Appl., 18(1):1–40, 2011.
  * (22) Daniel O’Connor and Lieven Vandenberghe.  Primal-dual decomposition by operator splitting and applications to image deblurring.  SIAM Journal on Imaging Sciences, 7(3):1724–1754, 2014.
  * (23) Leonid I. Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: Nonlinear Phenomena, 60(1):259 – 268, 1992.
  * (24) Jitkomut Songsiri.  Projection onto an $l_{1}$-norm ball with application to identification of sparse autoregressive models.  Asean Symposium on Automatic Control, 2011.  Vietnam.