\documentclass[11pt]{ansarticle}

%\usepackage[pagebackref=true]{hyperref} %zeigt auf welcher Seite die Referenz verwendet wird

% definitions used by included articles, reproduced here for 
% educational benefit, and to minimize alterations needed to be made
% in developing this sample file.
\usepackage{a4wide,amssymb,euscript,stmaryrd,epsfig}
\usepackage{paralist}
\usepackage{epsfig}
\usepackage{amsfonts}
\usepackage{latexsym}
\usepackage{xspace}
\usepackage{graphicx}
\usepackage{subcaption}
\usepackage{booktabs}
\usepackage{pgfplots}
\pgfplotsset{compat=1.9}

\usepackage{url}
\usepackage{lscape}

\usepackage{dune}

\newcommand{\pe}{\psi}
\def\d{\delta} 
\def\ds{\displaystyle} 
\def\e{{\epsilon}} 
\def\eb{\bar{\eta}}  
\def\enorm#1{\|#1\|_2} 
\def\Fp{F^\prime}  
\def\fishpack{{FISHPACK}} 
\def\fortran{{FORTRAN}} 
\def\gmres{{GMRES}} 
\def\gmresm{{\rm GMRES($m$)}} 
\def\Kc{{\cal K}} 
\def\norm#1{\|#1\|} 
\def\wb{{\bar w}} 
\def\zb{{\bar z}} 
\def\Prj{\mathbb{P}}
\def\nsm{\code{non-smooth-minization}}
\def\U{\mathcal{U}}
\def\P{\mathcal{P}}
\def\cimg{\textsc{CImg} }
\def\div{\operatorname{div}}


%\usepackage{fullpage}
% to activate colors, uncomment the following two lines
%\def\nnew{\color{red}}
%\def\mnew{\color{black}}
% to de-activate colors, uncomment the following two lines
\def\nnew{}
\def\mnew{}


%EXTRA MACROS FOR NONSIAM
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{alg}{Algorithm}
\newtheorem{assumption}{Assumption}
\newtheorem{property}{Property}
\newtheorem{proof}{Proof}

%MACROS FOR SIAM
%\newtheorem{remark}[theorem]{Remark}
%\newtheorem{alg}{Algorithm}
%\newtheorem{example}[theorem]{Example}
%\newtheorem{examples}[theorem]{Examples}



% some definitions of bold math italics to make typing easier.
% They are used in the corollary.

\def\bfE{\mbox{\boldmath$E$}}
\def\bfG{\mbox{\boldmath$G$}}


\newcommand{\supp}{\operatorname{supp}}
\newcommand{\sign}{\operatorname{sign}}
%\renewcommand{\R}{\mathbb R}
\newcommand{\I}{\mathbb I}
\newcommand{\Rb}{\overline{\R}}
\newcommand{\C}{\mathbb C}
\newcommand{\Cb}{\overline{\C}}
\newcommand{\Z}{\mathbb Z}
\newcommand{\N}{\mathbb N}
\newcommand{\D}{\mathbb D}

\newcommand{\F}{{\mathcal F}}
\newcommand{\J}{{\mathcal J}}
\newcommand{\V}{{\mathcal V}}
\def\H{{\mathcal H}}
\newcommand{\Var}{\operatorname{Var}}
\newcommand{\EAV}{\operatorname{EAV}}
\newcommand{\E}{\mathbb{E}}
\newcommand{\PSNR}{\operatorname{PSNR}}
\newcommand{\mcT}{{\mathcal T}}
\newcommand{\diam}{\operatorname{diam}}
\newcommand{\mcL}{{\mathcal L}}
\newcommand{\mcS}{{\mathcal S}}
\newcommand{\mcI}{{\mathcal I_h}}

\newcommand{\argmin}{\arg\min}
\def\NL{{\operatorname{NL}}}
\def\Tr{{\operatorname{Tr}}}

\title{Using \dune[ACFem] for Non-smooth Minimization of Bounded Variation Functions}
%\title{Properties of the $L^1$-$L^2$-$TV$ Model and Automated Parameter Selection for Removing Mixed Gaussian-Impulse Noise}
\author{Martin Alk\"amper and Andreas Langer\thanks{ Institute for Applied Analysis and Numerical Simulation, University of Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany, {\tt martin.alkaemper@mathematik.uni-stuttgart.de, andreas.langer@mathematik.uni-stuttgart.de} 
} 
}%\and Michael Hinterm\"uller\thanks{Department of Mathematics, Humboldt-University of Berlin, Unter den Linden 6, 10099 Berlin, Germany Email: {\tt hint@math.hu-berlin.de}}} 

\begin{document}

\maketitle
\pagestyle{myheadings}
\thispagestyle{plain}

\begin{abstract}
The utility of \dune[ACFem] is demonstrated to work well for solving a non-smooth minimization problem over bounded variation functions by implementing a primal-dual algorithm. The implementation is based on the simplification provided by \dune[ACFem]. Moreover, the convergence of the discrete minimizer to the continuous one is shown theoretically. 
\end{abstract}

\section{Introduction}
The \dune -framework \cite{BaBlDeEnKlKuOhSa06} and in particular the discretization module \dune[Fem] \cite{DeKlNoOh10} provides the means to handle discrete functions, operators and solvers on different grids. Still, implementing complicated partial differential equations (PDEs) and their solvers is cumbersome and tedious. The \dune - module \dune[ACFem] aims to simplify the usage of \dune[Fem] by defining expression templates for discrete functions and PDE models. This allows to linearly combine discrete functions and PDE models. Additionally it supports parallel and adaptive finite-element schemes on continuous discrete functions for the predefined and combined models. \cite{Heine14}\\
The flexibility of \dune (and \dune[Fem], \dune[ACFem]) allows e.g. to exchange discrete spaces by a single line of code or to change the gridtype and linear solvers.  \\
We will demonstrate the usefulness of \dune[ACFem] by minimizing a non-smooth functional consisting of a combined $L^1$/$L^2$-data fidelity term and a total variation term. Such an optimization problem has been shown to effectively remove Gaussian and salt-and-pepper noise, see \cite{HinLan2013, Lan2016}. In order to compute an approximate solution we use the primal-dual algorithm proposed in \cite{ChaPoc}.  For the numerical implementation we discretize using finite-element spaces defined over locally refined conforming grids. Motivated by the works \cite{Bar2012, Bar2013, Bar2015, HinRin}, where the considered functional is composed solely of an $L^2$-data term and a total variation term, we refine the grid adaptively using an a priori criterion. Similar as in \cite{Bar2012} we show for the considered minimization problem, that a minimizer over a finite element space converges to a minimizer in the space of functions of bounded variation as the mesh-size goes to 0. \\
In contrast to previous works \cite{Bar2012, Bar2013, Bar2015, HinRin}, we consider an additional non-smooth $L^1$-data term in the objective, which has to be treated carefully. Moreover due to the use of \dune[ALUGrid] \cite{AlDeKlNo16}  and the capabilities of \dune[ACFem] the resulting algorithm is intrinsically parallelized by domain decomposition.
 
%using an adaptation of the adaptive Steepest-Descent algorithm for non-smooth minimization introduced by Bartels in \cite{Bartels}. Non-smooth minimization of bounded variation functions is a problem that arises in imaging, e.g. from denoising images. (\cite{something})  To the functional of \cite{Bartels} we add an $L_1$-regularization term to additionally cope with salt-and-pepper noise. This has been shown to be effective \cite{blah}.\\

\section{Problem Formulation}
We consider the following problem
\begin{equation}
\min_{v\in BV(\Omega)\cap L^2(\Omega)} J_{\alpha_1,\alpha_2}(v):=\alpha_1 \| v - g\|_{L^1(\Omega)} + \alpha_2 \| v - g \|_{L^2(\Omega)}^2 + |Dv|(\Omega)
\label{eq:minProb}
\end{equation}
where $\Omega \subset \R^d$, $d\in \N$, is an open bounded set with Lipschitz boundary, $g \in L^2(\Omega)$ is a given datum, $\alpha_i \geq 0$ for $i=1,2$ with $\alpha_1 + \alpha_2 >0$  and $BV(\Omega)\subset L^1(\Omega)$ denotes the space of functions with bounded variation. that is, $v\in BV(\Omega)$ if and only if
\begin{equation}\label{eq:tv}
|Dv|(\Omega) := \sup\left\{ \int_{\Omega} v \div \vec\phi dx, \ \vec\phi\in C_0^{\infty}(\Omega,\R^d), \ \| \ |\vec\phi|_2\|_{C^0(\Omega)}\leq 1\right\}
\end{equation}
is finite, see \cite{AmbFusPal, Giu}. The space $BV(\Omega)$ endowed with the norm $\|v\|_{BV(\Omega)} = \|v\|_{L^1(\Omega)} + |Dv|(\Omega)$ is a Banach space \cite{Giu}. For $\alpha_2>0$ the minimization problem \eqref{eq:minProb} admits a unique solution owing to the strict convexity of the quadratic term \cite{Lan2016}. Note, that if $\alpha_1=0$ in \eqref{eq:minProb}, then we obtain the functional used in \cite{Bar2012, Bar2013, Bar2015, HinRin}. 


\subsection{Discretization}
Let $(\mcT_h)_{h>0}$ be a sequence of regular triangulations of $\Omega$ with (maximal) diameter $h = \max_{T\in \mcT_h} \diam(T)$. We define the following finite element spaces
\begin{equation*}
\begin{split}
&\mcL^0(\mcT_h)=\{q_h \in L^1(\Omega) : q_h\mid_T \mbox{ is constant for each } T\in \mcT_h\} \\
&\mcS^1(\mcT_h)=\{v_h \in C(\overline{\Omega}) : v_h\mid_T \mbox{ is affine for each } T\in \mcT_h\}.
\end{split}
\end{equation*}
The nodal interpolant $\mcI v \in \mcS^1(\mcT_h)$ of a function $v\in W^{2,p}$, with $\frac{d}{2}< p \leq \infty$ or $p=1$ if $d=2$, satisfies
$$
\| v - \mcI v \|_{L^p(\Omega)} + h \| \nabla (v - \mcI v) \|_{L^p(\Omega)} \leq c_\mcI h^2 \|D^2 v\|_{L^{p}(\Omega)},
$$ 
where $c_\mcI>0$; cf. \cite{BreSco}.

We recall, that the space $BV(\Omega)$ is continuously embedded in $L^p(\Omega)$ for $1\leq p\leq \frac{d}{d-1}$, i.e., there is a constant $c_{BV}>0$ such that $\|v\|_{L^p(\Omega)} \leq c_{BV}\|v\|_{BV} = c_{BV} \left( \|v\|_{L^1(\Omega)} + |Dv|(\Omega) \right)$ for any $v\in BV(\Omega)$. For $1\leq p < \frac{d}{d-1}$ this embedding is compact; cf. \cite{AmbFusPal}. Smooth functions are dense in $BV(\Omega) \cap L^p(\Omega)$, $1\leq p <\infty$, with respect to strict convergence in the sense that for $v\in BV(\Omega) \cap L^p(\Omega)$ and $\delta >0$ there exists $\varepsilon_0 > 0$ and functions $(v_\varepsilon)_{\varepsilon>0} \subset C^{\infty} \cap BV(\Omega) \cap L^p(\Omega)$ such that  for all $\varepsilon\leq \varepsilon_0$ we have
\begin{eqnarray}
&&\|\nabla v_\varepsilon\|_{L^1(\Omega)} \leq |Dv|(\Omega) + c_0\delta ,\label{Eq1}\\
&&\|v-v_\varepsilon\|_{L^p(\Omega)} \leq c_1 \delta ,\label{Eq2}\\
&&\|D^2 v_\varepsilon\|_{L^p(\Omega)} \leq c \varepsilon^{-2} \|v\|_{L^p(\Omega)}. \label{Eq3}
%&&\|D^2 v_\varepsilon\|_{L^1(\Omega)} \leq c \varepsilon^{-1} |Dv|(\Omega)\label{Eq4}.
\end{eqnarray}

Now we are able to show the following convergence result, which follows similar ideas as the proof of \cite[Theorem 3.1]{Bar2012}. 
\begin{theorem}
Let $u\in BV(\Omega)\cap L^2(\Omega)$ be a minimizer of the function $J_{\alpha_1,\alpha_2}$ and $u_h\in \argmin_{v\in \mcS^1(\mcT_h)} J_{\alpha_1,\alpha_2}(v)$. Then we have that $J_{\alpha_1,\alpha_2}(u_h) \to J_{\alpha_1,\alpha_2}(u)$ as $h\to 0$. If additionally $\alpha_2>0$, then $u_h\to u$ in $L^2(\Omega)$ as $h\to0$.
\end{theorem}

\begin{proof}
By the optimality of $u$ we have $J_{\alpha_1,\alpha_2}(u_h) - J_{\alpha_1,\alpha_2}(u) \geq 0$. For $\delta>0$ let $u_\varepsilon \in C^\infty(\Omega) \cap L^2(\Omega)$ as above and $\mcI u_\varepsilon$ its nodal interpolant, i.e., $\mcI u_\varepsilon \in \mcS^1(\mcT_h)$. Then we deduce
\begin{equation*}
\begin{split}
J_{\alpha_1,\alpha_2}(u_h) - J_{\alpha_1,\alpha_2}(u)  &\leq J_{\alpha_1,\alpha_2}(\mcI u_\varepsilon) - J_{\alpha_1,\alpha_2}(u) \\
 &= \| \nabla \mcI u_\varepsilon \|_{L^1(\Omega)} + \alpha_1 \| \mcI u_\varepsilon - g \|_{L^1(\Omega)} + \alpha_2 \| \mcI u_\varepsilon - g \|_{L^2(\Omega)} \\
 & \phantom{\| \nabla \mcI u_\varepsilon \|_{L^1(\Omega)} +}- |Du|(\Omega) - \alpha_1 \| u - g \|_{L^1(\Omega)} - \alpha_2 \| u - g \|_{L^2(\Omega)}.
\end{split}
\end{equation*}
Using \eqref{Eq1} and $\| \mcI u_\varepsilon - g\|^2_{L^2(\Omega)} - \|u-g\|^2_{L^2(\Omega)} = \int_{\Omega} (\mcI u\varepsilon - u) (\mcI u_\varepsilon + u -2g)$ we obtain
\begin{equation*}
\begin{split}
J_{\alpha_1,\alpha_2}(u_h) - J_{\alpha_1,\alpha_2}(u)  \leq  \| \nabla \mcI &u_\varepsilon \|_{L^1(\Omega)}  - \| \nabla u_\varepsilon\|_{L^1(\Omega)} + c_0 \delta + \alpha_1 \| \mcI u_\varepsilon - g \|_{L^1(\Omega)} \\
& - \alpha_1 \| u - g \|_{L^1(\Omega)} +\alpha_2 \int_{\Omega} (\mcI u\varepsilon - u) (\mcI u_\varepsilon + u -2g).
\end{split}
\end{equation*}
By the triangle-inequality and the Cauchy-Schwarz inequality we get
\begin{equation*}
\begin{split}
J_{\alpha_1,\alpha_2}(u_h) - J_{\alpha_1,\alpha_2}(u)  &\leq  \| \nabla( \mcI u_\varepsilon - u_\varepsilon)\|_{L^1(\Omega)} + c_0 \delta + \alpha_1 \| \mcI u_\varepsilon - g -(u - g)\|_{L^1(\Omega)}\\
&  \phantom{\leq  \| \nabla}+\alpha_2 \| \mcI u\varepsilon - u\|_{L^2(\Omega)} (\|\mcI u_\varepsilon\|_{L^2(\Omega)} + \|u\|_{L^2(\Omega)} + 2\|g\|_{L^2(\Omega)})\\
&=\| \nabla( \mcI u_\varepsilon - u_\varepsilon)\|_{L^1(\Omega)} + c_0 \delta + \alpha_1 \| \mcI u_\varepsilon - u_\varepsilon + u_\varepsilon-u\|_{L^1(\Omega)}\\
& \phantom{\leq  \| \nabla} +\alpha_2 \| \mcI u_\varepsilon -u_\varepsilon + u_\varepsilon- u\|_{L^2(\Omega)} (\|\mcI u_\varepsilon\|_{L^2(\Omega)} + \|u\|_{L^2(\Omega)} + 2\|g\|_{L^2(\Omega)})\\
&\leq\| \nabla( \mcI u_\varepsilon - u_\varepsilon)\|_{L^1(\Omega)} + c_0 \delta + \alpha_1 \| \mcI u_\varepsilon - u_\varepsilon\|_{L^1(\Omega)} + \|u_\varepsilon - u\|_{L^1(\Omega)}\\
& \phantom{\leq  \| \nabla} +\alpha_2 (\| \mcI u_\varepsilon - u_\varepsilon\|_{L^2(\Omega)} + \|  u_\varepsilon - u\|_{L^2(\Omega)}) (\|\mcI u_\varepsilon\|_{L^2(\Omega)} + \|u\|_{L^2(\Omega)} + 2\|g\|_{L^2(\Omega)})\\
\end{split}
\end{equation*}
The bound $\|\mcI u_\varepsilon\|_{L^2(\Omega)} + \|u\|_{L^2(\Omega)} + 2\|g\|_{L^2(\Omega)} \leq \tilde{c}$, which holds provided that $h\leq c\varepsilon$, and the nodal interpolant estimate yield
\begin{equation*}
\begin{split}
J_{\alpha_1,\alpha_2}(u_h) - J_{\alpha_1,\alpha_2}(u)  \leq c_\mcI h \|D^2 u_\varepsilon\|_{L^1(\Omega)} + c_0 \delta + c_\mcI \alpha_1 h^2 \|D^2 u_\varepsilon\|_{L^1(\Omega)} + \alpha_1  \|  u_\varepsilon - u\|_{L^1(\Omega)}\\
+ \tilde{c} \alpha_2 (c_\mcI h^2 \|D^2u_\varepsilon\|_{L^2(\Omega)} +  \|  u_\varepsilon - u\|_{L^2(\Omega)} ).
\end{split}
\end{equation*}
Using \eqref{Eq2}, \eqref{Eq3}, and the bound $ \|u\|_{L^2(\Omega)} \leq \tilde{c}$ we get
\begin{equation*}
\begin{split}
J_{\alpha_1,\alpha_2}(u_h) - J_{\alpha_1,\alpha_2}(u)  \leq C(C_1\frac{h}{\varepsilon^2} + C_2 \delta + C_3 \frac{h^2}{\varepsilon^2} + C_4 \delta + C_5 \frac{h^2}{\varepsilon^2} + C_6 \delta).
\end{split}
\end{equation*}
We choose $\delta=\frac{h}{\varepsilon^2}$ and deduce for $h\to 0$ that $\delta \to 0$ and hence $J_{\alpha_1,\alpha_2}(u_h) \to J_{\alpha_1,\alpha_2}(u)$.



If $\alpha_2>0$, then $J_{\alpha_1,\alpha_2}$ is strictly convex and it follows that
$$
J_{\alpha_1,\alpha_2}(u_h) - J_{\alpha_1,\alpha_2}(u) \geq \alpha_2 \|u_h-u\|_{L^2(\Omega)}^2,
$$
cf. \cite[Lemma 3.8]{Lan2016}. Hence $u_h\to u$ for $h\to 0$. \hfill $\square$
\end{proof}



\subsection{Steepest Descent Algorithm}\label{sec:steepdesc}

By the definition of the total variation \eqref{eq:tv} an equivalent formulation of \eqref{eq:minProb} reads
\begin{equation}
\min_{v\in BV(\Omega)\cap L^2(\Omega)} \max_{\vec p \in C_0^\infty(\Omega,\R^d), |\vec p|_2\leq 1}\alpha_1 \| v - g\|_{L^1(\Omega)} + \alpha_2 \| v - g \|_{L^2(\Omega)}^2 + \int_\Omega v \div \vec p.
\label{eq:minmaxProb}
\end{equation}
We discretize using finite element spaces $\U$ and $\P$ (e.g. $\mcS^1(\mcT_h)$ or $\mcL^0(\mcT_h)$) for the primal variable $u$ and the dual variable $\vec p$ respectively.

Following the ideas of Chambolle and Pock \cite[Algorithm 1]{ChaPoc} we formulate our steepest-descent algorithm by identifying 
$$F^*(\vec p)=\begin{cases} 0 & \text{if } |\vec p(x)|_2\leq 1\\ \infty & \text{else}\end{cases}, \ G(v)=\alpha_1 \| v - g\|_{L^1(\Omega)} + \alpha_2 \| v - g \|_{L^2(\Omega)}^2, \text{ and }\langle K v, \vec p \rangle = \int_\Omega v \div \vec p,
$$
where we assume that each element in $\P$ has a weak derivative. If this is not the case, we require it for $\U$ and define $\langle K v , \vec p \rangle = \int_\Omega \nabla v \vec p dx$.

We assume, that the datum $g \in L^2(\Omega)$ and the cost parameters $\alpha_i$ are given. Then in our algorithm we initialize $u_0 = \bar{u}_0 \in \U$, e.g., $u_0 \equiv 0$ or $u_0 =\Prj g$, where $\Prj : L^2(\Omega) \rightarrow \U$ denotes the $L^2$-Projection onto the discrete space $\U$, and $\vec p_0 \in \P$ to be the zero function. As parameters we need
 the step sizes $\sigma, \tau >0$,
 the coefficient $\beta = \frac{\tau\alpha_2}{1+\tau\alpha_1}$,
 and the overrelaxation parameter $\theta \in [0,1]$.\\
Then our steepest-descent algorithm iterates starting with $k=0$ as follows:
\begin{enumerate}
\item Set $\bar{p}$ according to the descent direction with step size $\sigma$ as
\begin{equation}
\bar{p} = \vec p_k + \sigma \nabla \bar{u}_k.
\label{eq:updateP}
\end{equation}
This is to be understood in the weak sense using test-functions from $\P$.
As the algorithm is derived from formulation $\eqref{eq:minmaxProb}$, we need $|p|_2 \leq 1$, so we have to project $p$ using
\begin{equation}
\vec p_{k+1}  = (I + \sigma \partial F^*)^{-1}(\bar{p}) \quad \Leftrightarrow  \quad \vec p_{k+1} = \frac{\bar{p}}{\max{(|\bar{p}|_2,1)}}.
\label{eq:projectP}
\end{equation}
\item Update $u_k$ using
\begin{equation}
u_{k+1} = (I + \tau \partial G)^{-1}(u_k + \div \vec p_{k+1}) \quad \Leftrightarrow \quad u_{k+1} = \begin{cases}
 z-\beta & \text{if } z-\beta \geq \Prj g \\
   z+\beta       & \text{if } z+\beta \leq \Prj g \\
   \Prj g &\text{else },
\end{cases}
\label{eq:updateU}
\end{equation}
where  
\begin{equation}
z = \frac{1}{1+\tau\alpha_1}\left(u_k + \tau \div \vec p_{k+1} + \tau \alpha_2 \Prj g \right) \in \U
\label{eq:calcZ}
\end{equation}
in the weak sense with test-space $\U$.
In equation \eqref{eq:updateU} and \eqref{eq:calcZ} instead of using the projection $\Prj g$ we could use the real datum $g$, which is equivalent to the above formulation, as the update is understood in the weak sense.
\item As proposed in \cite{ChaPoc} we overrelaxate to speed up the algorithm:
\begin{equation}
\bar{u}_{k+1} = u_{k+1} + \theta  (u_{k+1} - u_k)
\label{eq:updateUBar}
\end{equation}
\item If the stopping criterion
\begin{equation}
\frac{\|u_{k+1} - u_k\|_{L^2}}{\|u_k\|_{L^2}} < TOL
\label{eq:brkCrit}
\end{equation}
holds, we terminate the algorithm, otherwise we set $k \rightarrow k+1$ and repeat.
\end{enumerate}
 
Convergence of this algorithm to the discrete minimum of problem \eqref{eq:minProb} is guaranteed, if $\theta =1 $ and $\sigma \tau B^2<1 $, by using the standard convergence analysis from \cite{ChaPoc}, where $B:= \|\nabla\|_{L^2}$  is the operator norm of the gradient operator. This norm is bounded for discrete functions and of order $O(1/h)$, which implies that using finer meshes will always result in more steps of the steepest-descent algorithm. For such algorithms this is a typical and well-known behavior. Even using an adaptive mesh is not beneficial to the stepsize, as $B$ depends on the minimum cell-width. 


\section{Implementation Details}

The algorithm is implemented in the \dune -project \nsm. This project is compatible with the 2.4-release and available on the website \url{http://www.ians.uni-stuttgart.de/nmh/downloads}. It requires the modules \dune[ACFem] and all modules that are required by it, available at \url{gitlab.dune-project.org}. The installation works the standard \dune -way. (see Appendix A)

\subsection{Image Read-In and Noise Simulation} To read-in images we use The \cimg  Library -  C++ Template Image Processing Toolkit \cite{CIMG}. The single header file library \cimg  is an open source project to  provide easy handling and processing of images. It is capable of reading-in standard image formats  (e.g. \code{.png, .tif, .jpg}) and of adding certain types of noise.  Each pixel is accessible by its coordinates utilizing the  data type \code{image}, whose constructor gets the filename containing the image location. Moreover, \code{image} provides a method to add noise to the data with a given noise level for different noise-types. 
In particular, Gaussian noise with noise level (variance) $\gamma$ and mean $0$ and salt-and-pepper noise $S$ with noise level $s\in [0,1]$, which corrupts an original image $\hat{g}$ by
\[
S\hat{g}(x) = \begin{cases}  \min{\hat{g}}  & \text{with probability }s/2 \\ \max{\hat{g}} & \text{with probability }s/2  \\ \hat{g}(x) & \text{with probability }1-s\end{cases},
\]
are available.\\
We define two adapter classes to convert the data type \code{image} and one adapter class to convert a algebraic expression into a discrete function within the header file \code{src/datafunction.hh}.
\begin{itemize}
\item The \code{LocalDataAdapter} constructs an $\mcL^0(\mcT_h)$ function that approximates an inserted algebraic expression.
\item The \code{LocalImageAdapter} converts an instance of class \code{image} into an $\mcL^0(\mcT_h)$ function using the \code{GridFunctionAdapter} provided by \dune[ACFem]. 
\item The \code{LocalNoiseAdapter} enhances an instance of class \code{image} by noise. More precisely it adds two independent types of noise, where every noise implemented by \cimg can be applied. Noise type and noise level are controlled by a set of parameters defined in the parameter file. For instance, setting \code{nsm.noiseType1} to 0 $0$ leads to Gaussian noise and setting it to $2$ corresponds to salt-and-pepper noise. Setting any noise level parameter to 0 disables the respective noise. For other noise types and more information consult the documentation of \cimg. After applying the noise the corrupted instance of the class \code{image} is then converted into an $\mcL^0(\mcT_h)$ as in the \code{LocalImageAdapter}.
\end{itemize}
The artificial addition of noise (\code{LocalNoiseAdapter}) is only needed for demonstration purposes as in real-world application the image is already corrupted by noise due to certain physical processes.\\
The most important method of the adapter classes is the \code{init()} method, that sets the value on each mesh cell. For example, using the \code{LocalDataAdapter} the following code results in an approximation of a circle of radius $0.3$ and center $(0.5,0.5)$ with value $1$ inside and $0$ outside, see Figure \ref{subfig:circle}.
\lstinputlisting{imageread.code}
The adapters convert the image into discrete functions on the unit square $[0,1] \times [0,1]$,  as we can scale any rectangular image accordingly. The image resolution and aspect ratio relate to the level of initial refinement and to the number of cells in each direction. So $n$ uniform refinements with $l\times k$ initial cells lead to an image resolution of size $l2^n \times k2^n$ square pixels, where every square is composed of two triangles and $l/k$ corresponds to the aspect ratio. This is done in the gridfile \code{data/unitsquare2d.dgf} which describes the unit square using the \dune \code{DGF - Interval} \cite{BaBlDeEnKlKoOhSa08} nomenclature. 
\subsection{Implementation of the steepest-descent algorithm} The steepest-descent algorithm from section \ref{sec:steepdesc} is implemented in the file \code{src/nsm.cc}, where some of its subroutines, which depend on the continuity of the discrete spaces $\U$ and $\P$, are outsourced into the file \code{src/phc.hh}. The template class \code{ProjectionHelperClass} uses partial template specialization to correctly define the methods \code{calculateZ()} and \code{entitywiseProjection()} . The main part of the steepest-descent algorithm reads:\\
\lstinputlisting{nsm.code}

 The method \code{phc.entitywiseProjection()} realizes equations \eqref{eq:updateP} and \eqref{eq:projectP}.  Solving eq. \eqref{eq:updateP} is relatively easy, as all steps are provided by \dune[ACFem], e.g. for continuous $\U$ and $\P$:
\lstinputlisting{updateP.code}
 \code{Models} (e.g. \code{massModel}, \code{gradientModel}, ...) of \dune[ACFem] are objects that define the integral kernels to be multiplied with  test-functions and thus describe the equation to be solved. \code{Models} can be linearly combined over the field $\R$ and discrete functions can be added (or subtracted). \code{FemSchemes} (e.g. \code{EllipticFemScheme}, ...) are objects that represent the complete finite-element method.  A \code{FemScheme} always solves \code{Model} $=0$. In the above case, for given $p$ and $\bar{u}$ the equation
\[
0 = \langle \bar{p}, \phi \rangle - \sigma \langle \bar{u} , -\div \phi \rangle - \langle p, \phi \rangle \quad \forall \phi \in \P
\]
with dirichlet-zero boundary conditions is solved with respect to $\bar{p}$. The solution is written into the variable $p$.\\

Projecting $p$ to be feasible (eq. \eqref{eq:projectP})  is not in the features of \dune[ACFem]. So we have to use the features of \dune[Fem] directly. We choose to do the projection entity-wise, i.e. we iterate over all entities of the grid. On each entity we iterate over the degrees of freedom and if $|p|_2 > 1$ holds, we restrict the corresponding value to length $1$ in the same direction. If the polynomial order of $\P$ is less or equal than $1$, this implies the necessary condition $\| |p|_2\|_{L^\infty} \leq 1$.\\

 The calculation of $z$ (eq. \eqref{eq:calcZ}) translates very nicely into code. We use the \code{L2Projection} of \dune[ACFem] and its ability to linearly combine discrete functions. If $\P$ is continuous, the weak form of  eq. \eqref{eq:calcZ} is
\[
\langle z_k , \psi \rangle = \frac{1}{1+\tau\alpha_1}\left\langle u_k +\tau \div p_{k+1} + \tau \alpha_2 g, \psi \right\rangle \quad \forall \psi \in \U
\]
and translates into code in the following way:
\lstinputlisting{l2projection.code}
\dune[ACFem] currently does not implement DG-methods. Hence, if $\P$ is not continuous, we have to use the weak divergence to shift the derivative to the test-space $\U$, which means that $\U$ has to be continuous.  If neither $\U$ nor $\P$ is continuous, the program fails.  As above we use the basic \code{Models} of \dune[ACFem], in particular the \code{massModel} and the \code{weakDivergenceModel} to implement eq. \eqref{eq:calcZ} for discontinuous $\P$ as follows:
\lstinputlisting{calculateZ.code}
Note that in contrast to the implementation of method \code{entitywiseProjection()} we construct the \code{massModel} from a discrete space instead of a discrete function. This is explicitly allowed by \dune[ACFem] as the object simply has to provide the data type for the test-function space.\\

\subsection{Adaptive Refinement} Adaptive refinement is initiated by a positive parameter \code{nsm.localRefine}. This parameter denotes the number of additional local refinements to be done in the initialization phase in addition to the uniform refinements. If it is set to $\leq 0$, no local refinement will be performed. The grid is locally refined at discontinuities of the piecewise constant datum $g \in \mcL^0(\mcT_h)$ in the following way: Given an Entity $E$ and its neighbour $N$, if
\begin{equation}
\label{eq:adaptCrit}
|g_{|N} - g_{|E} | > TOL
\end{equation}
holds, then $E$ is marked for refinement, where the value of $TOL$ is given by the parameter \code{nsm.adaptTolerance}. We do not refine during the iterations of the steepest-descent algorithm, but keep the mesh static. So the refinement increases the resolution of discontinuities and hence drastically improves the projection $\Prj g$. This is implemented in the method \code{adaptGrid()} in the file \code{nsm.cc}. We use the grid manager \dune[ALUGrid] \cite{AlDeKlNo16}, which is among other things capable of handling the needed conforming, adaptive triangular grids.



\section{Numerical Experiments}\label{sec:experiments}

To reproduce the data of the experiments of this section, consult Appendix A and follow the described steps.


\subsection{Comparison of Discrete Spaces}\label{subsec:exp1}
We investigate for different discrete spaces $\U$ and $\P$ the behaviour of our algorithm with respect to the non-smooth minimization problem \eqref{eq:minProb} considering the following set-ups
\[
\quad \U = \mcS^1(\mcT_h) =\P,
\quad \U = \mcL^0(\mcT_h), \P = \mcS^1(\mcT_h),
\quad \U = \mcS^1(\mcT_h), \P = \mcL^0(\mcT_h).
\]  For this set of experiments we choose a similar example as in \cite{Bar2015}, i.e., the observed data is 1 on a disk of radius $0.3$ and 0 elsewhere. Hence, we choose the discrete datum $g$ to be given as a piecewise constant approximation of the function
\[
  f(x) = \begin{cases}
  1, & \text{if }|x-(0.5,0.5)| \leq 0.3 \\
  0, & \text{else}.
  \end{cases}
\]    The cost parameters are set to $\alpha_1 = 10$ and $\alpha_2 = 20$, the stopping tolerance to $10^{-5}$ and the adaptation tolerance to $0.1$. We do  $5$ uniform refinements and $3$ additional local refinements. As we use the unit square as our domain the operator norm of the gradient is $B=2^8$. Consequently the step sizes are automatically set to $\tau=\sigma=2^{-8}$. Then, in order to guarantee convergence, the overrelaxation parameter $\theta$ is chosen to be 1.

%%%%%%%%%%%
%% figure Pg
%%%%%%%%%%
\begin{figure}[H]
\begin{subfigure}{.5\textwidth}
\includegraphics[width=\linewidth]{./udiscont0pcont1/prjgudiscont0.png}
\caption{$g = \Prj g$, as $\U = \mcL^0(\mcT_h)$ }
\label{subfig:circle}
\end{subfigure}
\begin{subfigure}{.5\textwidth}
\includegraphics[width=\linewidth]{./ucont1pcont1/prjgucont1.png}
\caption{$\Prj g$ for $\U = \mcS^1(\mcT_h)$}
\end{subfigure}
\caption{The Projection of the piecewise constant datum $g$ onto the space $\U$.}
\label{fig:ex1PrjG}
\end{figure}
By imposing all equations weakly the steepest-descent algorithm de facto uses the $L^2$-projection $\Prj g \in \U$ instead of the datum $g \in \mcL^0(\mcT_h)$.  Figure \ref{fig:ex1PrjG} shows that in the case of $\U = \mcS^1(\mcT_h)$  by projecting a discontinuous function onto a continuous space we introduce an additional error in contrast to the case $\U=\mcL^0(\mcT_h)$. So intuitively it would be a good idea to choose $\U=\mcL^0(\mcT_h)$ to execute the algorithm, but we observe that the other setting performs better. This is due to the fact that in general, approximating functions of bounded variation is not possible with piecewise constant functions on a triangulation, since the length of discontinuities may not diminish over refinement \cite{Bar2012}. The over- and undershoots of the continuous projected data do not pose a problem, because they are regularized over the course of the algorithm, as we can see in Figure \ref{fig:ex1u}.


%%%%%%%%%%%
%% figure u
%%%%%%%%%%
\begin{figure}[H]
\begin{center}
\begin{subfigure}{.3\textwidth}
\includegraphics[width=\linewidth]{./ucont1pcont1/ucont1pcont1.png}
\caption{$\U = \mcS^1(\mcT_h) = \P$}
\end{subfigure}
\begin{subfigure}{.3\textwidth}
\includegraphics[width=\linewidth]{./ucont1pdiscont0/ucont1pdiscont0.png}
\caption{$\U = \mcS^1(\mcT_h)$, $\P = \mcL^0$}
\end{subfigure}
\begin{subfigure}{.3\textwidth}
\includegraphics[width=\linewidth]{./udiscont0pcont1/udiscont0pcont1.png}
\caption{$\U =\mcL^0(\mcT_h)$, $\P= \mcS^1(\mcT_h)$}
\label{subfig:ul0}
\end{subfigure}
\caption{The discrete minima $u^*$ of the functional $J_{10,20}$}
\label{fig:ex1u}
\end{center}
\end{figure}

%%%%%%%%%%%
%% figure Energy
%%%%%%%%%%

\begin{figure}[H]\begin{center}
\begin{subfigure}{0.3\textwidth}
\begin{tikzpicture}[scale=.5]
\begin{axis}[xlabel={No. Iteration}, ymin=2, ymax=6,xmin=0,xmax=11000]
\addplot table[x index=0, y index=4, col sep=comma, mark={} ] {./ucont1pcont1/graph.txt};
\end{axis}
\end{tikzpicture}
\caption{$\U = \mcS^1(\mcT_h) =\P$,\\ \# Iterations: 9907,\\ $J_{\alpha_1,\alpha_2}(u^*) = 2.41776$}
\end{subfigure}
\begin{subfigure}{0.3\textwidth}
\begin{tikzpicture}[scale=.5]
\begin{axis}[ xlabel={No. Iteration}, ymin=2, ymax=6,xmin=0,xmax=11000]
\addplot table[x index=0, y index=4, col sep=comma, mark={} ] {./ucont1pdiscont0/graph.txt};
\end{axis}
\end{tikzpicture}
\caption{$\U = \mcS^1(\mcT_h)$, $\P = \mcL^0(\mcT_h)$,\\ \# Iterations: 10234,\\ $J_{\alpha_1,\alpha_2}(u^*) = 2.06217$}
\end{subfigure}
\begin{subfigure}{0.3\textwidth}
\begin{tikzpicture}[scale=.5]
\begin{axis}[xlabel={No. Iteration}, ymin=2, ymax=6,xmin=0,xmax=11000]
\addplot table[x index=0, y index=4, col sep=comma, mark={} ] {./udiscont0pcont1/graph.txt};
\end{axis}
\end{tikzpicture}
\caption{$\U = \mcL^0(\mcT_h)$, $\P = \mcS^1(\mcT_h)$,\\ \# Iterations: 8021,\\ $J_{\alpha_1,\alpha_2}(u^*) = 3.49582$}
\end{subfigure}
\caption{The value of the functional $J_{\alpha_1,\alpha_2}$ over the course of the algorithm.}
\label{fig:energy}\end{center}
\end{figure}


Figure \ref{fig:energy} depicts the performance of the energy $J_{\alpha_1,\alpha_2}$ of the algorithm during the iterations.  We observe that the algorithm does not converge monotonically and also rather slowly, which is typical for steepest descent type methods. The combination $\U = \mcS^1(\mcT_h)$, $\P = \mcL^0(\mcT_h)$ yields the best result, while for $\U = \mcL^0(\mcT_h)$, $\P = \mcS^1(\mcT_h)$ the energy of the solution stays far from the minimum as already expected and also visible in Figure \ref{subfig:ul0}, where one may still observe over- and undershoots.

\subsection{Image Denoising}
Here we demonstrate the denoising capability of the algorithm. Motivated by the above experiments we choose the spaces to be  $\U = \mcS^1(\mcT_h)$, $\P = \mcL^0(\mcT_h)$ again. The considered image is first corrupted by Gaussian white noise with variance $0.1$ and then salt-and-pepper noise with $s=0.1$ is added. The obtained image is shown in Figure \ref{fig:cameramana}. In order to reconstruct the image we choose $\alpha_1 =250$, $\alpha_2=150$ in \eqref{eq:minProb}, and perform our steepest-descent algorithm with the tolerance set to $10^{-4}$, $\theta=1$ and $\sigma=\tau=2^{-8}$. To resolve the $256 \times 256$ pixel-sized picture we apply $8$ uniform refinements and no additional local refinement. In Figure \ref{fig:cameramanb} we depict the output of our algorithm. The result is reasonably smoothed due to the total variation regularization, while discontinuities are still preserved. 


\begin{figure}[hpbt]\begin{center}
\begin{subfigure}{0.3\textwidth}
\includegraphics[width=\linewidth]{./ex2/noised.png}
\caption{Noisy image}\label{fig:cameramana}
\end{subfigure}
\begin{subfigure}{0.3\textwidth}
\includegraphics[width=\linewidth]{./ex2/denoised.png}
\caption{Result of the algorithm}\label{fig:cameramanb}
\end{subfigure}
\caption{A corrupted image before and after the algorithm}
\label{fig:cameraman}\end{center}
\end{figure}


%The result looks reasonable, as we did not try to choose optimal parameters and the original image is strongly corrupted.

\subsection{Strong Scaling} 
We do a strong scaling experiment on the same datum as in Section \ref{subsec:exp1}. The computation is executed on a 32 core shared memory system.  The grid manager \dune[ALUGrid] balances the computational load on the initial grid by partitioning it onto the different processors. So we cannot expect the algorithm to scale if the initial mesh is to coarse. E.g. if there are only two initial elements, only two processors can get partitions that are non-empty. Consequently we discretize the unit square by $8 \times 8$ elements, see \code{data/finecube_2d.dgf}. This results in a different operator norm on the initial grid, so we have to set the corresponding parameter \code{nsm.constantL} to $0.125$ for the stepsizes $\tau, \sigma$ to be computed correctly. As the grid is already fine, we do not need as many uniform refinements to reach a good resolution, so we do $3$ uniform and another $3$ local refinements. The tolerance is set to $10^{-4}$ and $\alpha_1 = \alpha_2 = 150$ to reach a short runtime.


\begin{figure}[H]\begin{center}
\begin{tikzpicture}[scale=0.8]
\begin{loglogaxis}[legend entries={simulation, optimal}, xlabel={No. Procs}, ylabel={Runtime in s}]
\addplot table[x index=0, y index=3, col sep=comma ] {./parallelTest/parallelTest.time.graph};
 \addplot[mark=none, red] coordinates {(1,10) (32,0.3125)};
\end{loglogaxis}
\end{tikzpicture}
\caption{Strong scaling of the algorithm.}
\label{fig:scaling}\end{center}
\end{figure}


Figure \ref{fig:scaling} indicates that the strong scaling at least up to $8$ cores is quite good and afterwards it stagnates, which is probably due to the small grid and not enough workload. We would have to increase the number of initial cells to improve the scaling further. The important part about this result is, that the parallelism is done inside \dune[ACFem] and we spent almost no effort to parallelize the code. The only line of code needed initializes \code{MPI}.

\section{Summary and Outlook}

We have shown that the discrete minimizer converges to the continuous one.
The steepest-descent algorithm used here is implemented very conveniently within \dune[ACFem]. The flexibility of \dune[Fem] allows us to easily exchange discrete spaces and easily set parameters. Also the algorithm is now intrinsically parallel by domain decomposition with a decent strong scaling. Future research includes implementing a semi-smooth newton method in \dune[ACFem] to increase convergence speed even for highly adaptive grids.

\section*{Acknowledgements}
The author Martin Alk\"amper would like to thank the German Research Foundation (DFG) for financial support of the project within the Cluster of Excellence in Simulation Technology (EXC 310/2) at the University of Stuttgart.

\section*{Appendix A - Installation Instructions}

The program is designed to work in a Linux-environment. It may as well work on a unix machine.
To install the program, first download the \dune modules, \dune[common], \dune[geometry], \dune[grid], \dune[istl], \dune[localfunctions], \dune[Fem], \dune[ALUGrid] and \dune[ACFem] from \url{gitlab.dune-project.org} and \nsm from \url{http://www.ians.uni-stuttgart.de/nmh/downloads}. For the \dune- modules checkout the branch \code{releases/2.4}. Put all the projects in a directory as direct subfolders. Then run \\ \code{./dune=common/bin/dunecontrol --opts=config.opts --module=non-smooth-minization all}\\ An example config file \code{example.opts} is provided in the main directory of the \dune project \nsm. \\
Now \dune is installed and the executable \code{nsm} from the subdirectory \code{src} should have been built within the \code{cmake} build directory. The scripts \code{example1.sh} to \code{example3.sh} produce the results presented in this paper. They require a \code{python} interpreter and \code{example3.sh} requires \code{mpi}. The \code{cmake} build directory is assumed to be named \code{build-cmake}. The scripts are simple so changes should be easy. The scripts basically consist of changing directories, compiling the executable \code{nsm}, executing it with the right set of parameters and parsing the output into suitable directories. There are two types of output. Console output from \code{std::out} is parsed into a file \code{output.graph} in a subdirectory of \code{src}, which describes the value of the functional and the value of all its parts with respect to both $g$ and $\Prj g$ and a file \code{output.graph}, where the applied parameters and other output can be looked up. Console output from \code{std::err} is displayed on the console. The other type of output is in the \code{cmake} build directory located within the directory \code{output} in suitable subdirectories. These contain a set of \code{.vtu/.vtk} files (readable with \code{Paraview}), that contain the datum $g$ (called "image"), the projection $\Prj g$ (called "projection of g"), the solution (called "u0") and the adjoint variable (componentwise "p0", "p1"). It will take some time to run the scripts. To get faster results, simply edit the \code{LOCREF} or \code{INITREF} variable of the bash scripts or lower the tolerance parameter \code{nsm.tolerance}.

\section*{Appendix B - Parameters}
There are two types of parameters: Runtime parameters and Compile-time Parameters.
\begin{itemize}
 
\item Runtime parameters can be overloaded on the command line \\ \code{./nsm nsm.parameter:value ...}.\\ The default values are set inside the parameter file located at \code{data/parameter}. This file is copied into the \code{cmake} build directory at compile time. So parameters changed in the build directory will be overwritten when recompiling. Parameters with prefix fem, or istl belong to \dune[Fem] and \dune[ISTL] respectively, and are explained in their documentation. The specific parameters of this algorithm are prefixed nsm. The extensive list of runtime parameters reads:\\
\begin{tabular}{l|r}
nsm.theta & overrelaxation parameter $\theta$ \\
nsm.maxIt & Maximum number of iterations of the algorithm\\
nsm.tolerance & The algorithm breaks if the tolerance is reached \\
nsm.outputStep &  Output every nth step \\
nsm.constantL & Norm of the gradient on the initial mesh \\
nsm.tau & Stepsize $\tau$ \\
nsm.sigma & Stepsize $\sigma$ \\
nsm.lambda\_1 & $L^1$-data term cost coefficient $\alpha_1$ \\
nsm.lambda\_2 & $L^2$-data term cost coefficient $\alpha_2$ \\
nsm.image & Filename of the image to be read-in \\
nsm.initialRefinements & Number of initial uniform refinements \\
nsm.localRefine & Number of initial adaptive refinements\\
nsm.adaptTolerance & Adaptive tolerance\\
nsm.noiseType1 &  The \cimg type of noise to be applied first\\
nsm.noiseLevel1 & The \cimg noise level of the first noise\\
nsm.noiseType2 & The \cimg type of noise to be applied second\\
nsm.noiseLevel2 & The \cimg noise level of the second noise\\
\end{tabular}
\item Compile-time parameters have to be declared when compiling. The location is in the file \code{CMakelists.txt}. As they are \code{cmake}-cache variables, they can be redefined in the usual way. (see e.g. \code{./example1.sh}) They determine what kind of problem to treat and which discrete spaces are to be used.\\
\begin{tabular}{l|r}
NSM\_SET\_STEPSIZE & If true, tau and sigma are set manually\\
NSM\_USE\_IMAGE & If true, image is used instead of geometric expression\\
NSM\_USE\_NOISE & If true, noised imaged is used, requires NSM\_USE\_IMAGE\\
NSM\_U\_POLORDER & The polynomial order of $u$ (default 1)\\
NSM\_P\_POLORDER & The polynomial order of $p$ (default 1)\\
NSM\_U\_DISCONT & If defined, $\U$ is discontinuous.\\
NSM\_P\_DISCONT & If defined,  $\P$ is discontinuous.\\
\end{tabular}
\end{itemize}

\bibliographystyle{plain}
\bibliography{bibliography} 
\end{document}