\section{Introduction}

The task of single image Super-Resolution (SR) consists in improving the spatial resolution of an observed Low-Resolution (LR) imaging data so as to obtain a High-Resolution (HR) version which, typically, can be used as a reference for subsequent analysis. Image resolution is limited in many applications due to the optical characteristics and the physical limitations of the acquisition devices. Some standard examples are biomedical and astronomic imaging where, due to light aberration phenomena, close objects (molecules, stars\ldots) on LR images cannot be correctly distinguished/detected, see, e.g.~\cite{Galbraith2011,Willett2004}.
SR techniques are often employed also in image recognition problems. This is the case, for instance, of QR code recognition where images are often captured by scanning tools (e.g. cell-phones) from relatively large distances which may affect the accuracy of the recognition \cite{Kato2011}.
Analogously, in remote sensing applications such as land-cover classification, only LR measurements are often available, which may limit significantly the classification precision \cite{Ling2016,wickham2014}. 


Mathematically, the task can be formulated as an ill-posed inverse problem: for a given vectorised LR image ${\bm{g}} \in \mathbb{R}^{M}$, we look for its HR version ${\bm{u}} \in \mathbb{R}^{N}$ defined on a space of dimension $N=L^{2}M$ with magnification factor $L >1$ which satisfies the following linear degradation model: 
\begin{equation}\label{eq:model}
    {\bm{g}} = \mathbf{S} \mathbf{H} {\bm{u}} + \boldsymbol{\eta}.
\end{equation}
Here, $\mathbf{S} \in \mathbb{R}^{M \times N}$ stands for the down-sampling operator,  $\mathbf{H} \in \mathbb{R}^{N \times N}$ describes blur degradation and $\boldsymbol{\eta}$ denotes the realisation of an Additive White Gaussian Noise (AWGN) r.v. with zero mean and standard deviation $\sigma_{\mathbf{\eta}}>0$.

Due to the ill-posedness of the operator $\mathbf{S}\mathbf{H}$, a standard approach for solving \eqref{eq:model} consists in encoding prior knowledge about the solution ${\bm{u}}$ and on the data statistics via an energy minimisation approach, so that an approximated solution ${\bm{u}}^{*}\in\mathbb{R}^{N}$ is computed by solving
\begin{equation}  \label{eq:reg_mod}
    {\bm{u}}^{*} \in  \argmin_{{\bm{u}} \in \mathbb{R}^{N}} \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \rVert^{2}_{2} + \mu R({\bm{u}}),
\end{equation}
where the quadratic fidelity term models the presence of AWGN while the (possibly non-convex) regularisation term $R: \mathbb{R}^{N}\to \mathbb{R} \cup \left\{+\infty\right\}$ encodes prior information on the target image ${\bm{u}}$, thus ensuring the stability of the inversion process. The two terms are balanced by a regularisation parameter $\mu>0$. We refer the reader to \cite{Chaudhuri2001,Yue16} and to the references therein for a review on variational approaches for SR problems.

In this work, we choose $R$ so as to promote gradient sparsity, which is often desirable in image segmentation applications whenever a simplified, edge-preserving version of the original data ${\bm{g}}$ is required for further analysis.
In recent years, sparse and non-convex gradient-based regularisation approaches have become very popular in the context of image reconstruction due to their better ability of preserving sharp edges even in low-contrast scenarios. A significant contribution has been made by Storath et al. in a series of papers \cite{Storath14,Storath15,Storath17} where sparsity on the image gradient $\mathbf{D} {\bm{u}}\in\mathbb{R}^{2N}$ is promoted by $\ell^0$ regularisation which reads
\begin{equation}
\lVert \mathbf{D} {\bm{u}}\rVert_0:=\# \left\{ (\mathbf{D} {\bm{u}})_i, i=1,\ldots 2N: (\mathbf{D} {\bm{u}})_i\neq 0 \right\}.
\end{equation}
This choice has been thoroughly studied for several imaging problems such as deconvolution, sparse recovery, joint reconstruction and segmentation, see, e.g., \cite{Xu11}. Moreover, it has has been shown to be very useful in many situations where a further classification/labelling step is required. In this work we propose a novel numerical scheme endowed with convergence guarantees which justifies the use of this type of regularisation in the context of SR problems with high blur and noise degradation.


\subsection{Related work}
The vast majority of sparse optimisation approaches for SR problems enforces sparsity either on the signal itself \cite{Candes} or its representation w.r.t. to some basis/overcomplete dictionary \cite{Yang2010}. These methods and their non-convex extensions have been shown to be very powerful in several applications such as image microscopy \cite{Gazagnes} where signal-sparsity can be assumed. However, for non-point-like objects (such as piece-wise constant regions), this type of regularisation is not the appropriate choice. 
Other classical approaches to the SR problem are based on the use of least squares, Fourier series and Tikhonov-type gradient regularisations \cite{Chaudhuri2001}, which favour noise removal at the price of creating smoothing and ringing artefacts which are undesirable in many applications such as object detection, where images with sharper edges are preferable for better classification.
To overcome this drawback, the use of edge-preserving convex regularisations based on the idea of gradient sparsity, such as Total Variation (TV) \cite{Chan2007,MarquinaOsher2008,tao2009alternating,Gao2016}, its fractional  \cite{Yao2020} and  $\ell^0$ extension \cite{Storath14,Storath15,Storath17}, has been proposed. Such methods have shown good performances in many applications, although their convexity (in the case of TV) or their challenging numerical realisation (in the case of non-convex approaches) often limit their practical use and precision. Different approaches for solving the SR problem make use of deep architectures encoding prior information on the desired HR solution from a training set of examples \cite{yang2019,wang2020,sun2020learned}. In particular, in \cite{zhang2017} the authors present a Plug-and-Play (PnP) framework \cite{venkatakrishnan2013} which exploits deep convolutional neural network denoisers embedded in standard optimisation algorithms, such as Alternating Direction of Multipliers (ADMM) or Half-Quadratic Splitting (HQS). Differently from model-based variational approaches, deep learning-based methods do not require an explicit expression of the regularisation term $R$, since this can be learned directly from the data and adapted to the particular application considered. Those methods have currently reached state-of-the-art performances in several image reconstruction problems, although their theoretical foundation and their stability to noise perturbations still limits their practical use in the case of highly-degraded image data. 


\subsection{Contribution}

We consider a variational model for solving problem \eqref{eq:reg_mod} where a quadratic data fidelity is combined with an $\ell^0$-gradient regularisation term both in a coupled (isotropic) and decoupled (anisotropic) form, the latter being better suited for directionally-biased images, such as QR scans.
To solve the model efficiently, we propose to use an ADMM algorithm which decomposes the original problem into substeps cheaply solved by means of direct hard-thresholding and standard iterative Conjugate Gradient (CG) linear solvers. Our variable splitting differs from the one introduced by  Storath et al. in \cite{Storath14,Storath15,Storath17}, where the non-convex substeps are solved by means either of approximate graph-cuts approaches \cite{Boykov2001} or
dynamic programming algorithms. As well as for these different numerical schemes, we prove in this paper fixed-point convergence for the proposed ADMM algorithm.
Up to our knowledge, the same variable splitting has been used only in the case of convex regularisation functions, such as TV, in \cite{tao2009alternating,Gao2016} where convergence to the global minimum is proved.

We test our SR model on real-world applications (QR scanning, cell detection and land-cover labelling) where a simplified HR version of the given LR image ${\bm{g}}$ is required in view of further analysis, showing that the proposed model improves significantly
segmentation and labelling precision. 


\subsection{Organisation of the paper} 

In Section \ref{sec:SRmodel} we provide a review of gradient-sparse variational methods for single-image SR. In Section \ref{sec:ADMM} we present a novel converging ADMM scheme for solving the proposed model along with details on its practical realisation. In Section \ref{sec:modelling} we report some numerical tests on model parameter sensitivity performed on synthetic data. Finally, in Section \ref{sec:experiments} we apply our model to some real-world applications such as QR scanning, cell detection, compressed JPG SR and land-cover classification
We report the convergence proofs of the proposed ADMM schemes in Appendix \ref{sec:convergence} to improve the flow of the manuscript.



\input{section/model}
\input{section/method}
\input{section/results}





\section{Conclusions}

We considered a variational model with $\ell^0$ gradient-sparsity-promoting regularisation combined with a quadratic data fidelity for single-image super-resolution of images corrupted by blur and Gaussian noise. The use of non-convex $\ell^0$ jump-sparse regularisations has been originally proposed in \cite{Storath14} in the context of general 1D inverse problems and subsequently applied in \cite{Storath15,Storath17} to image segmentation and reconstruction problems. To overcome the computational limitations required by the use of ADMM splitting strategies considered in these works, we propose a novel ADMM algorithm allowing for the efficient solution of its subproblems by means of direct hard-thresholding or standard CG solvers. For the proposed scheme we prove fixed-point convergence results assuming specific growth conditions on the sequence of penalty parameters. We validate our model on synthetic data and test it on real-world examples where gradient-sparse super-resolved outputs are required in view of an accurate recognition/classification step (such as QR code recognition \cite{Kato2011}, cell detection and land-cover classification \cite{cihlar2000}).  By numerous comparisons with convex and non-convex variational approaches, and with state-of-the-art deep learning methods \cite{zhang2017,sun2020learned}, we show that the proposed approach significantly improves classification precision, while limiting at the same times smoothing and loss-of-contrast artefacts in comparison with classical convex regularisations.

Further work should address the use of analogous regularisations and algorithms for the joint modelling of SR and segmentation problems via, e.g., Mumford-Shah functionals \cite{Storath15}. Furthermore, the extension of the convergence results to other gradient discretisations and to less restrictive growth conditions for the sequence of penalty parameters is envisaged. 










\appendices


\input{section/Theorem1}





\section*{Acknowledgments}

LC and PC acknowledge the support received by the Academy "Complex Systems" of the JEDI IDEX of the Université Côte d'Azur.
ELP and PC acknowledge the support received by the INDAM-GNCS (Research projects 2020).

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi






\input{main.bbl}
\vspace{-1cm}







\begin{IEEEbiography}
[{\includegraphics[width=1in,height=1in,clip=true]{./bio/pasquale_pic.png}}]{Pasquale Cascarano}
received his master degree in Mathematics in 2018 at University of Bologna. He is currently a Ph.D.\ student in Applied Mathematics at the University of Bologna (IT). His research focuses on variational and deep learning methods for imaging inverse problems.
\end{IEEEbiography}
\vspace{-1.5cm}
\begin{IEEEbiography}
[{\includegraphics[width=1in,height=1.3in,trim=1.5cm 7cm 2.5cm 0cm,clip=true]{./bio/luca_pic.pdf}}]{Luca Calatroni}
completed his Ph.D.\ in Applied Mathematics in 2015 as part of the Cambridge Image Analysis research group (UK). He was then a Marie Skłowdoska-Curie research fellow at the University of Genova (Italy) and \emph{Lecteur Hadamard} FMJH fellow at the \'Ecole Polytechnique (France). From October 2019, he is permanent CNRS researcher at the I3S laboratory in Sophia Antipolis (France) within the Morpheme research group. His research focuses on variational methods and non-smooth optimisation algorithms for imaging and vision.
\end{IEEEbiography}
\vspace{-1.2cm}
\begin{IEEEbiography}
[{\includegraphics[width=0.8in,height=1.3in,clip=true]{./bio/loli_photo.jpg}}]{Elena Loli Piccolomini} is Professor of Numerical Analysis at the University of Bologna. Her research topics are regularisation methods for inverse problems in imaging, with particular focus on medical imaging and tomographic image reconstruction.
\end{IEEEbiography}







\end{document}

\section{Convergence analysis} \label{sec:convergence}

We report here a complete convergence proof of Theorem \ref{th:conv_ADMManiso} and a sketch of the proof of Theorem \ref{th:conv_ADMMiso}, which is based on similar arguments.

\subsection{Proof of Theorem \ref{th:conv_ADMManiso}}

\begin{proof} 
We consider the ADMM sequences $({\bm{u}}^{k}),({\bm{t}}^{k}),({\bm{s}}^{k})$, defined in \eqref{eq:ADMM_aniso_1}-\eqref{eq:ADMM_aniso_3}. We want to show that there exists ${\bm{u}}^*$ such that:
\begin{align}
{\bm{u}}^{k} \rightarrow {\bm{u}}^{*}, \quad {\bm{t}}^{k} \rightarrow\mathbf{D}_{h}{\bm{u}}^{*}\nonumber, \quad {\bm{s}}^{k} \rightarrow  \mathbf{D}_{v}{\bm{u}}^{*}. 
\end{align}
To shorten the proof, we remark that everything proved for the sequences $({\bm{t}}^{k})$,$(\beta_{t}^{k})$, $(\boldsymbol{\lambda}_{t}^{k})$ and $(\mathbf{D}_{h}{\bm{u}}^{k})$ can be deduced for $({\bm{s}}^{k})$, $(\beta_{s}^{k})$, $(\boldsymbol{\lambda}_{s}^{k})$ and $(\mathbf{D}_{v}{\bm{u}}^{k})$ in the same way.

We start defining the following functionals:
\begin{align*}
&G^{h}_{k}({\bm{t}}):= \mu \lVert {\bm{t}} \rVert_{0} + \dfrac{\beta_{t}^{k}}{2} \lVert {\bm{t}} - (\mathbf{D}_{h}{\bm{u}}^{k}+\dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})\rVert_{2}^{2}, \\
& F_{k}({\bm{u}}):= \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \rVert_{2}^{2}  + \dfrac{\beta_{t}^{k}}{2}\lVert \mathbf{D}_{h}{\bm{u}} - ({\bm{t}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})\rVert^{2}_{2} + \nonumber \\  & + \dfrac{\beta_{s}^{k}}{2}\lVert \mathbf{D}_{v}{\bm{u}} - ({\bm{s}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta_{s}^{k}})\rVert^{2}_{2}.
\end{align*}
\paragraph*{Step 1}
There holds: 
\begin{align}
&\lVert {\bm{t}}^{k+1} -\mathbf{D}_{h}{\bm{u}}^{k} - \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{2} \leq \sqrt{\dfrac{2 \mu N}{\beta_{t}^{k}}}. \label{proof:step1_limitato_t}
\end{align}
This inequality can  be trivially shown by the minimality of ${\bm{t}}^{k+1}$ in \eqref{eq:ADMM_aniso_1} which entails   $G^{h}_{k}({\bm{t}}^{k+1}) \leq  G^{h}_{k}(\mathbf{D}_{h}{\bm{u}}^{k} + \frac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})$, therefore we get:
\begin{align*}
& \mu \lVert {\bm{t}}^{k+1}\rVert_{0} + \dfrac{\beta_{t}^{k}}{2} \lVert {\bm{t}}^{k+1} - (\mathbf{D}_{h}{\bm{u}}^{k}+\dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})\rVert_{2}^{2} \notag\\
& \leq \mu \lVert \mathbf{D}_{h}{\bm{u}}^{k} + \frac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{0} \leq \mu N,
\end{align*}
by definition of $\lVert \cdot \rVert_0$,where we recall $N$ is the dimension of the vector ${\bm{u}}^{k}$. By neglecting the first term on the Left Hand Side (LHS)  of the above inequality, we deduce \eqref{proof:step1_limitato_t}.

\paragraph*{Step 2}
From the minimality of ${\bm{u}}^{k+1}$ in \eqref{eq:ADMM_aniso_3} we have: $F_{k}({\bm{u}}^{k+1}) \leq F_{k}({\bm{u}}^{k})$ for every $k$. By definition of $F_{k}$ and applying \eqref{proof:step1_limitato_t} and its analogous related to the sequences $({\bm{s}}^{k})$, $(\beta_{s}^{k})$, $(\boldsymbol{\lambda}_{s}^{k})$ and $(\mathbf{D}_{v}{\bm{u}}^{k})$, we deduce: 
\begin{align}
    & \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k+1} - {\bm{g}} \rVert_{2}^{2}  + \dfrac{\beta_{t}^{k}}{2} \lVert \mathbf{D}_{h} {\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{2}^{2}  \label{proof:riferimento} \\
    & + \dfrac{\beta_{s}^{k}}{2} \lVert \mathbf{D}_{v} {\bm{u}}^{k+1} - {\bm{s}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta_{s}^{k}}  \rVert_{2}^{2} \leq  
   
   
   \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} - {\bm{g}} \rVert_{2}^{2} + 2\mu N. \nonumber
\end{align}
Since the all the terms on the LHS of \eqref{proof:riferimento} are nonnegative, the following inequality holds: 
\begin{align}
    & \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k+1} -{\bm{g}} \rVert_{2}^{2} \leq \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} -{\bm{g}} \rVert_{2}^{2} + 2 \mu N \leq \ldots \\
    & \leq  \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{0} - {\bm{g}} \rVert_{2}^{2} +  2 \mu N k \nonumber
\end{align}
From \eqref{proof:riferimento} and by the sub-additivity property of the square root we can also derive the following inequality:
\begin{align}\label{ineq:series_convergent2}
    &\lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}\rVert_{2}    \leq \sqrt{\frac{1}{{\beta_{t}^{k}}}}\lVert \mathbf{S} \mathbf{H} {\bm{u}}^{0} - {\bm{g}} \rVert_{2} + \sqrt{4\mu N \frac{k}{\beta_{t}^{k}}}
\end{align}


   

\paragraph*{Step 3}
We show that the sequences $\mathbf{D}_{h}{\bm{u}}^{k}$ and $\mathbf{D}_{v}{\bm{u}}^{k}$ are Cauchy sequences, hence they converge. We prove this for $\mathbf{D}_{h}{\bm{u}}^{k}$, the proof for $\mathbf{D}_{v}{\bm{u}}^{k}$ is identical.
\begin{align*}
    & \lVert \mathbf{D}_{h} {\bm{u}}^{k+1} - \mathbf{D}_{h} {\bm{u}}^{k} \rVert_{2} \leq \nonumber \\
    &\leq \lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} +  \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}  \rVert_{2}  +   \lVert \mathbf{D}_{h}{\bm{u}}^{k} - {\bm{t}}^{k+1} +  \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}  \rVert_{2} . 
\end{align*}
By assumption \ref{cond:1_ADMManiso} applied on the RHS of \eqref{ineq:series_convergent2} we deduce:
\begin{align} 
    &\lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}\rVert_{2} \rightarrow 0, \label{prof:bound1
   
\end{align}
which, combined with \eqref{proof:step1_limitato_t} and \eqref{ineq:series_convergent2} entails that $\mathbf{D}_{h}{\bm{u}}^{k}$ is a Cauchy sequence. Hence  it converges to a point ${\bm{t}}^{*}$. Similarly, $\mathbf{D}_{v}{\bm{s}}^{k}$  converges to a point ${\bm{s}}^{*}$.


\paragraph*{Step 4}
We prove now the convergence of the sequences ${\bm{t}}^{k}$ and $\mathbf{D}_{h} {\bm{u}}^{k}$.
By writing \eqref{eq:ADMM_aniso_4} as: 
\begin{align}
 &\frac{\boldsymbol{\lambda}_{t}^{k+1}}{\beta_{t}^{k}} = \mathbf{D}_{h} {\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}, \label{eq:proof_dual1
\end{align}
and  from \eqref{prof:bound1} we deduce that $\frac{\lVert \boldsymbol{\lambda}_{t}^{k+1} \rVert_{2}}{\sqrt{\beta_{t}^{k}}} \rightarrow 0 $. By monotonicity of the $(\beta^k_t)$ we then deduce that  $\frac{\lVert \boldsymbol{\lambda}_{t}^{k} \rVert_{2}}{\sqrt{\beta_{t}^{k}}} \rightarrow 0 $. Hence:
\begin{align*}
& \lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} \rVert_{2} \leq  \dfrac{\lVert \boldsymbol{\lambda}_{t}^{k+1}\rVert_{2} + \lVert \boldsymbol{\lambda}_{t}^{k}\rVert_{2}}  {\sqrt{\beta_{t}^{k}}}, 
\end{align*}
where both quantities on the RHS tend to $0$ as $k\to\infty$. Therefore, by the uniqueness of the limit, ${\bm{t}}^{k} \longrightarrow {\bm{t}}^{*}$ and $\mathbf{D}_{h}{\bm{u}}^{k} \longrightarrow {\bm{t}}^{*}$.  

\paragraph*{Step 5}
We can now prove convergence of the sequence $({\bm{u}}^{k})$ . For simplicity, let us define the quantities $\mathbf{A}:=\mathbf{S}\mathbf{H}$ and  $\mathbf{M}_{k}:=\dfrac{1}{\beta_t^{k}} \mathbf{A}^{T}\mathbf{A} + \mathbf{D}_{h}^{T}\mathbf{D}_{h} + \dfrac{\beta_{s}^{k}}{\beta_{t}^{k}}\mathbf{D}_{v}^{T}\mathbf{D}_{v}$, for every $k$. By \ref{cond:2_ADMManiso}, we observe that the matrix $\mathbf{M}_k$ is invertible for all $k$ and that the optimality condition of \eqref{eq:ADMM_aniso_3} reads:
\begin{equation}
\mathbf{M}_{k}{\bm{u}}^{k}= \mathbf{D}_{h}^{T} ({\bm{t}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta^{k}}) + \dfrac{\beta_{s}^{k}}{\beta_{t}^{k}}\mathbf{D}_{v}^{T}({\bm{s}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta^{k}}) + \dfrac{1}{\beta_{t}^{k}}\mathbf{A}^{T}{\bm{g}}. \nonumber
\end{equation}
Since ${\bm{t}}^{k+1} \to {\bm{t}}^{*}$, ${\bm{s}}^{k+1} \to {\bm{s}}^{*}$,  $\frac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \to 0$, $\frac{\boldsymbol{\lambda}_{s}^{k}}{\beta_{s}^{k}} \to 0$,  and by Assumptions \ref{cond:1_ADMManiso} and \ref{cond:2_ADMManiso},  we have that   $\frac{1}{\beta_{t}^{k}}\mathbf{A}^{T}{\bm{g}} \to \mathbf{0}$ so that the RHS converges pointwise to ${\bm{z}}^{*}=\mathbf{D}_{h}^{T}t^{*} + c \mathbf{D}_{v}^{T}s^{*} $. Additionally, the sequence $\mathbf{M}_{k}^{-1}$ converges pointwise to $\mathbf{M}^{*}$. We thus have that ${\bm{u}}^k=\mathbf{M}_{k}^{-1}\mathbf{M}_{k}{\bm{u}}^{k} \to \mathbf{M}^{*}{\bm{z}}^{*}:={\bm{u}}^{*}$. \\
We now want to show that ${\bm{t}}^*=\mathbf{D}_h {\bm{u}}^*$ and, similarly, that ${\bm{s}}^*=\mathbf{D}_v {\bm{u}}^*$. We show the details only for the former case. By the triangle inequality we get:
\begin{align}
\lVert {\bm{t}}^* -  \mathbf{D}_h {\bm{u}}^*\rVert_2  \leq & \lVert {\bm{t}}^* -  \mathbf{D}_h {\bm{u}}^k\rVert_2 + \lVert \mathbf{D}_h {\bm{u}}^k - \mathbf{D}_h {\bm{u}}^* \rVert_2 \nonumber \\
& \leq \lVert {\bm{t}}^* -  \mathbf{D}_h {\bm{u}}^k\rVert_2 + \lVert \mathbf{D}_h \rVert_2 \lVert {\bm{u}}^k - {\bm{u}}^* \rVert_2, \nonumber
\end{align}
where both terms tend to $\mathbf{0}$ since $\mathbf{D}_h {\bm{u}}^k \to {\bm{t}}^* $ and ${\bm{u}}^k \to {\bm{u}}^*$. 
\end{proof}




\subsection{Proof of Theorem \eqref{th:conv_ADMMiso}}

\begin{proof}
The proof of Theorem \eqref{th:conv_ADMMiso} follows the same steps as the previous one. The only main difference in it is the definiton of $\mathbf{M}_{k}$, which reads in this case:
\begin{align}
    \mathbf{M}_{k}:= \dfrac{1}{\beta^{k}} \mathbf{A}^{T}\mathbf{A} + \mathbf{D}^{T}\mathbf{D} =
    \dfrac{1}{\beta^{k}} \mathbf{A}^{T}\mathbf{A}+
    \mathbf{D}_{h}^{T}\mathbf{D}_{h} + \mathbf{D}_{v}^{T}\mathbf{D}_{v}. \nonumber
\end{align}
By proceeding similarly as above the conclusion holds.
\end{proof}
\section{An efficient ADMM splitting} \label{sec:ADMM}
 
For $p\in\left\{1,2\right\}$, we consider the non-smooth and non-convex SR model \eqref{eq:reg_mod} with the choice \eqref{eq:zero_norm}, that is:
\begin{equation}\label{eq:min_TV0_unconstrained}
\mathbf{u^*} \in \argmin_{{\bm{u}} \in \mathbb{R}^{N}}~ \left\{\Phi({\bm{u}};\mu,p):=\frac{1}{2}\rVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \lVert_{2}^{2} + \mu \lVert \mathbf{D} {\bm{u}} \rVert_{0,p} \right\}.
\end{equation}
Existence of solutions for \eqref{eq:min_TV0_unconstrained} is guaranteed by the following theorem whose proof can be found in \cite[Theorem 1]{Storath14} for a general forward operator $\mathbf{A}$.

\smallskip

\begin{theorem}  \label{th:existence}
The solution set of both the anisotropic ($p=1$) and  isotropic ($p=2$)  problem \eqref{eq:min_TV0_unconstrained} is non-empty.
\end{theorem}

To solve numerically  problem \eqref{eq:min_TV0_unconstrained} we propose an iterative \emph{alternating direction method of multipliers} (ADMM) based on a suitable variable splitting. 
We separate the description for the anisotropic and isotropic case. For both cases, fixed-point convergence of the ADMM iterates upon suitable conditions is proved.

\subsection{ADMM for the anisotropic regularisation}

For $p=1$, we can rewrite the unconstrained minimisation problem \eqref{eq:min_TV0_unconstrained} in the following equivalent constrained form:
\begin{align*
   \argmin_{{\bm{u}}}  & ~\frac{1}{2}\rVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \lVert_{2}^{2}~ + ~\mu (\lVert {\bm{t}} \rVert_{0}  + \lVert {\bm{s}} \rVert_{0} ) \\  \nonumber
    s.t. \quad           & ~ {\bm{t}}:=\mathbf{D}_{h}{\bm{u}},\quad {\bm{s}}:=\mathbf{D}_{v}{\bm{u}}
\end{align*}
where ${\bm{t}},{\bm{s}} \in \mathbb{R}^{N}$ represent the horizontal/vertical gradient components, respectively.

We can then define the augmented Lagrangian function:
\begin{align} \label{eq:lagrangian_anis}
& L_{\beta_t,\beta_s}({\bm{u}};{\bm{t}},{\bm{s}},\boldsymbol{\lambda}_{t},\boldsymbol{\lambda}_{s}):=   \frac{1}{2}\rVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \lVert_{2}^{2}  + \mu \lVert {\bm{t}} \rVert_{0} + \mu\lVert {\bm{s}} \rVert_{0}\nonumber \\
                          & + \langle \boldsymbol{\lambda}_{t}, \mathbf{D}_{h} {\bm{u}} - {\bm{t}} \rangle + \dfrac{\beta_{t}}{2}\lVert \mathbf{D}_{h}{\bm{u}} - {\bm{t}} \rVert_{2}^{2} + \langle \boldsymbol{\lambda}_{s}, \mathbf{D}_{v}{\bm{u}} - {\bm{s}} \rangle  \nonumber\\
                          & + \dfrac{\beta_{s}}{2}\lVert \mathbf{D}_{v}{\bm{u}} - {\bm{s}} \rVert_{2}^{2}
\end{align}
where $\beta_{t}$ and $\beta_{s}$ are two positive penalty parameters and $\boldsymbol{\lambda}_{t}$ and $\boldsymbol{\lambda}_{s}$ are the vectors of Lagrange multipliers related to the auxiliary variables ${\bm{t}}$ and ${\bm{s}}$, respectively. By letting the two parameters $\beta_t, \beta_s$ increase along the iterations (we will provide specific growth conditions in the following Theorem \ref{th:conv_ADMManiso}), we can then minimise \eqref{eq:lagrangian_anis} with respect to ${\bm{t}},{\bm{s}}$ and ${\bm{u}}$ by iterating the following scheme:
 \begin{numcases}{}
  \scalebox{0.9}{${\bm{t}}^{k+1}  \in  \underset{{\bm{t}}}{\argmin}~ \mu \lVert {\bm{t}} \rVert_{0} + \dfrac{\beta^{k}_{t}}{2} \lVert {\bm{t}} - (\mathbf{D}_{h}{\bm{u}}^{k}+\dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta^{k}_{t}})\rVert_{2}^{2}$} &  \label{eq:ADMM_aniso_1} \\
  \scalebox{0.9}{${\bm{s}}^{k+1}  \in  \underset{{\bm{s}}}{\argmin}~ \mu \lVert {\bm{s}} \rVert_{0} + \dfrac{\beta^{k}_{s}}{2} \lVert {\bm{s}} - (\mathbf{D}_{v}{\bm{u}}^{k}+\dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta^{k}_{s}})\rVert_{2}^{2}$} & \label{eq:ADMM_aniso_2} \\
  \scalebox{0.9}{${\bm{u}}^{k+1}  = \underset{{\bm{u}}}{\argmin} ~\dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \rVert_{2}^{2} +$ } \nonumber \\
  \hspace{0cm} \scalebox{0.85}{$+ \dfrac{\beta^{k}_{t}}{2}\lVert \mathbf{D}_{h}{\bm{u}} - ({\bm{t}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta^{k}_{t}})\rVert^{2}_{2} + \dfrac{\beta^{k}_{s}}{2}\lVert \mathbf{D}_{v}{\bm{u}} - ({\bm{s}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta^{k}_{s}})\rVert^{2}_{2}$} & \label{eq:ADMM_aniso_3}\\
  \scalebox{0.9}{$\boldsymbol{\lambda}_{t}^{k+1} = \boldsymbol{\lambda}_{t}^{k} - \beta^{k}_{t}({\bm{t}}^{k+1}-\mathbf{D}_{h}{\bm{u}}^{k+1})$} &\label{eq:ADMM_aniso_4}\\
  \scalebox{0.9}{$\boldsymbol{\lambda}_{s}^{k+1} = \boldsymbol{\lambda}_{s}^{k} - \beta^{k}_{s}({\bm{s}}^{k+1}-\mathbf{D}_{v}{\bm{u}}^{k+1})$,} & \label{eq:ADMM_aniso_5}
  \end{numcases}
where a gradient ascent update of $\boldsymbol{\lambda}_{t}$ and $\boldsymbol{\lambda}_{s}$ is also applied. 

Under suitable growth assumptions, the sequences \eqref{eq:ADMM_aniso_1}, \eqref{eq:ADMM_aniso_2}, \eqref{eq:ADMM_aniso_3} converge to a fixed point (see Appendix \ref{sec:convergence} for the proof). 




\begin{theorem}   \label{th:conv_ADMManiso}
Let the ADMM iterations \eqref{eq:ADMM_aniso_1}-\eqref{eq:ADMM_aniso_5} be defined under the following conditions:
\begin{enumerate}[label=\textbf{A.\arabic*}]
    \item $(\beta_t^k)$,$(\beta_s^k)$ are increasing sequences such that $\sum_{k=1}^{+\infty}\sqrt{\frac{k}{\beta_t^{k}}} < + \infty$, $\sum_{k=1}^{+\infty} \sqrt{\frac{k}{\beta_s^{k}}} < + \infty$ and $  \frac{\beta_s^{k}}{\beta_t^{k}} \to c\neq 0$. \label{cond:1_ADMManiso}
    \item $\mathbf{D}_{h}$ and $\mathbf{D}_{v}$ are full rank.\label{cond:2_ADMManiso}
\end{enumerate}
Then, the sequences $({\bm{t}}^{k}), ({\bm{s}}^{k}), ({\bm{u}}^{k})$ converge, i.e.:
$$
{\bm{t}}^{k} \longrightarrow {\bm{t}}^{*},\ \ 
{\bm{s}}^{k} \longrightarrow {\bm{s}}^{*},\ \
{\bm{u}}^{k} \longrightarrow {\bm{u}}^{*},
$$
with ${\bm{t}}^{*}=\mathbf{D}_{h}{\bm{u}}^{*}$ and ${\bm{s}}^{*}=\mathbf{D}_{v}{\bm{u}}^{*}$.
\end{theorem}
\medskip

We remark that the full rank assumption on the operators $\mathbf{D}_h$ and $\mathbf{D}_v$ is verified, for instance, if Dirichlet boundary conditions are assumed. A sufficient condition which guarantees the required growth of the penalty sequences is $\beta_t^k = \beta_s^k= O(k(1+\epsilon)^{k}),~ 0<\epsilon\ll 1$.


\subsection{ADMM for the isotropic regularisation}

For $p=2$ we can write problem \eqref{eq:min_TV0_unconstrained} in the following equivalent constrained form:
\begin{align}\label{eq:iso_constrained}
    \argmin_{u} &~ \frac{1}{2}\rVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \lVert_{2}^{2}~ +~ \mu  \sum_{i=1}^{N} \big{|} \lVert {\bm{z}}_{i} \rVert \big{|}_{0} \\ \nonumber
    s.t. \quad           & ~ {\bm{z}} := \mathbf{D} {\bm{u}}
\end{align}
where ${\bm{z}}_{i}:= \big{(}(\mathbf{D}_{h}{\bm{u}})_{i},(\mathbf{D}_{v}{\bm{u}})_{i}\big{)}\in \mathbb{R}^2$, for each $i=1,\ldots,N$. The augmented Lagrangian function reads in this case:
\begin{align} \label{eq:lagrangian_iso}
L_{\beta}({\bm{u}};{\bm{z}},\boldsymbol{\lambda}):=  & \frac{1}{2}\rVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \lVert_{2}^{2}  ~+~\mu  \sum_{i=1}^{N} \big{|} \lVert {\bm{z}}_{i} \rVert \big{|}_{0} \nonumber \\
                          & + \langle \boldsymbol{\lambda}, \mathbf{D} {\bm{u}} - {\bm{z}} \rangle + \dfrac{\beta}{2}\lVert \mathbf{D} {\bm{u}} - {\bm{z}} \rVert_{2}^{2}    
\end{align}
where $\beta>0$ is a scalar penalty parameter and $\boldsymbol{\lambda}\in \mathbb{R}^{2\times N}$  is the Lagrange multiplier vector. As above, by letting the penalty parameter increases along the iterations at a certain growth (see the following Theorem \ref{th:conv_ADMMiso}), we seek for minimisers  of \eqref{eq:iso_constrained} by iterating the following scheme: 
\begin{numcases}{}
 \scalebox{0.85}{${\bm{z}}^{k+1} \in \underset{{\bm{z}}}{\argmin} \ \mu\sum_{i=1}^{N} \big{|} \lVert {\bm{z}}_{i} \rVert \big{|}_{0} + \dfrac{\beta^{k}}{2} \lVert {\bm{z}} - (\mathbf{D} {\bm{u}}^{k}+\dfrac{\boldsymbol{\lambda}^{k}}{\beta^{k}})\rVert_{2}^{2}$ }  \label{eq:ADMM_iso_1} &\\
\scalebox{0.85}{${\bm{u}}^{k+1} = \underset{{\bm{u}}}{\argmin} \ \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \rVert_{2}^{2} + \dfrac{\beta^{k}}{2}\lVert \mathbf{D} {\bm{u}} - ({\bm{z}}^{k+1} - \dfrac{\boldsymbol{\lambda}^{k}}{\beta^{k}})\rVert^{2}_{2}$} &  \label{eq:ADMM_iso_2} \\
\scalebox{0.85}{$\boldsymbol{\lambda}^{k+1} = \boldsymbol{\lambda}^{k} - \beta^{k}({\bm{z}}^{k+1}-\mathbf{D} {\bm{u}}^{k+1})$}. &  \label{eq:ADMM_iso_3}
\end{numcases}
For this scheme, a similar result as the one in Theorem \ref{th:conv_ADMManiso} holds (see Appendix \ref{sec:convergence} for a sketch of the proof).

\begin{theorem}   \label{th:conv_ADMMiso}
Let the ADMM iterations \eqref{eq:ADMM_iso_1}-\eqref{eq:ADMM_iso_2} be defined under the following conditions:
\begin{enumerate}[label=\textbf{I.\arabic*}]
    \item $(\beta^k)$ is an increasing sequence  such that $\sum_{k}^{+\infty}\sqrt{\frac{k}{\beta^{k}}} < + \infty$\label{cond:1_ADMMiso}
    \item $\mathbf{D}$ is full rank.\label{cond:2_ADMMiso}
\end{enumerate}
Then, $({\bm{z}}^k)\longrightarrow {\bm{z}}^{*}$, $({\bm{u}}^k)\longrightarrow {\bm{u}}^*$ and ${\bm{z}}^{*}=\mathbf{D}{\bm{u}}^*$.
\end{theorem}






We remark that in order to guarantee the convergence of the sequence $({\bm{u}}^{k})$ to  ${\bm{u}}^*$,  Theorems \ref{th:conv_ADMManiso} and \ref{th:conv_ADMMiso} require full rank on the operators $\mathbf{D}_h$ and $\mathbf{D}_v$. This is not very limiting since Dirichlet boundary conditions can always be imposed through an artificial image padding of the image. Our numerical experiments, however, showed numerical convergence even when periodic boundary conditions are used. A theoretical convergence proof in this case is left for future research. As far as the growth condition on the penalty parameters is concerned, we remark that in \cite{Storath14} a geometric growth was assumed. Unfortunately, this is not enough for our theoretical convergence result to hold, as oscillations may appear if this is violated. We comment more on this in Section \ref{sec:numaspects}. 


\subsection{Efficient solution of the ADMM subproblems} \label{subsec:subproblems}

We report here some practical details on the the efficient solutions of the subproblems \eqref{eq:ADMM_aniso_1}-\eqref{eq:ADMM_aniso_3} and  \eqref{eq:ADMM_iso_1}-\eqref{eq:ADMM_iso_2}. 

\subsubsection*{Solution of $\ell^0$ subproblems} Due to decomposability of the $\ell^{0}$ term, solving problems \eqref{eq:ADMM_aniso_1},\eqref{eq:ADMM_aniso_2} corresponds to solve the $N$ one-dimensional $\ell^2-\ell^0$ problems
\begin{align}\label{eq:1Dstandard_l0}
     \underset{{\bm{t}}_{i} \in \mathbb{R}}{\argmin}~  \delta  \lvert {\bm{t}}_{i} \rvert_{0} + ( {\bm{t}}_{i} - {\bm{f}}_{i} )_{2}^{2}
\end{align}
where $\delta=\frac{2\mu}{\beta^{k}_{t}}$ and ${\bm{f}}_i= (\mathbf{D}_{h}{\bm{u}}^{k}+\frac{\boldsymbol{\lambda}_{t}^{k}}{\beta^{k}_{t}})_{i}$ for \eqref{eq:ADMM_aniso_1}, while $\delta=\frac{2\mu}{\beta^{k}_{s}}$, ${\bm{f}}_i= (\mathbf{D}_{v}{\bm{u}}^{k}+\frac{\boldsymbol{\lambda}_{s}^{k}}{\beta^{k}_{s}})_i$ for \eqref{eq:ADMM_aniso_2}.  As far as the problem \eqref{eq:ADMM_iso_1} is concerned, it similarly reduces to the solution of the $N$ two-dimensional $\ell^{0}$-regularised problems
\begin{equation}\label{eq:2Dstandard_l0}
 \underset{{\bm{z}}_{i} \in \mathbb{R}^{2}}{\argmin} \  \delta  \big{|} \lVert {\bm{z}}_{i} \rVert \big{|}_{0} + \lVert {\bm{z}}_{i} - {\bm{f}}_{i} \rVert_{2}^{2}
\end{equation}
where $\delta=\frac{2\mu}{\beta^{k}}$ and ${\bm{f}}_{i} = (\mathbf{D}_{h} {\bm{u}}^{k}_{i}+\frac{(\boldsymbol{\lambda}^{k})_{1,i}}{\beta^{k}},\mathbf{D}_{v} {\bm{u}}^{k}_{i}+\frac{(\boldsymbol{\lambda}^{k})_{2,i}}{\beta^{k}})$.
Solving \eqref{eq:1Dstandard_l0} and \eqref{eq:2Dstandard_l0} corresponds to compute the proximal map of $|\cdot|_0$ with parameter $\delta$ evaluated in ${\bm{f}}_i$, which is nothing but the 1D \cite{blumensath2009iterative} and 2D \cite{Xu11} hard-thresholding operators, respectively.
\subsubsection*{Solution of the quadratic subproblems} 
The first order optimality conditions of problems \eqref{eq:ADMM_aniso_3} and  \eqref{eq:ADMM_iso_2} lead to the solution of large-size linear systems, whose coefficient matrix is symmetric and positive definite. To solve them efficiently, we make use of Conjugate Gradient (CG) algorithm with a warm-start initialisation at every iteration. We remark that, due to the presence of the downsampling operator $\mathbf{S}$, the use of more efficient solvers based, for instance, on discrete Fourier transforms are here not possible, as the product matrix $\mathbf{S}\mathbf{H}$ does not have a block-circulant structure. However, under suitable assumptions on the down-sampling operator $\mathbf{S}$, the problem admits a closed form solution \cite{zhao2016fast}.







\subsection{Comparisons with previous splittings}\label{sec:comparisons}
The variable splitting and the ADMM iterations considered above are different than the ones considered in \cite{Storath14,Storath15,Storath17} where the choice ${\bm{z}}={\bm{u}}$ in \eqref{eq:iso_constrained} is made. Our choice avoids the presence of the gradient operator in the $\ell^{0}$-based problems \eqref{eq:ADMM_aniso_1}-\eqref{eq:ADMM_aniso_2} and \eqref{eq:ADMM_iso_1}, leading to the faster computation of their solution by direct solvers without requiring the use of approximate solvers based on approximate graph-cut algorithms  \cite{Storath14}. These latter algorithms have  well-known drawbacks such as strong dependence on the initialisation and require an approximate inner solver \cite{Storath15,boykov2004}. As an alternative, in \cite{Storath15}, the isotropic substep is solved by a set of anisotropic problems along the diagonal or knight-move directions, each of which is computed by dynamic programming algorithms with computational cost $O(N^{2})$ compared to $O(N)$ in our approach.

















\begin{comment}

\begin{table}[!htb]
\centering
 \begin{tabular}{l}
 \hline 
 \hspace{20mm} \textbf{Algorithm 1: SRL0-ADMM}\\
 \hline
  \textbf{Input:} LR data ${\bm{g}} \in \mathbb{R}^{M}$, regularization parameter $\mu > 0$, measurement matrices $S \in \mathbb{R}^{M \times N}$, $H \in \mathbb{R}^{N \times N}$,   \\
  {\tt repeat} \\
  \textbf{Output:} $v= x_m -  \frac{ \lambda_k }{r} \sum_{j=1}^r H^TS^T (S H  x_m - g_j ) $\\
  \hspace{5mm} $p=0$ \\
 
  \hspace{5mm}  $W_{\ell, \mu}=D\left (div  \ p - \dfrac{v}{\beta}\right)_{\ell, \mu}$ \\
  \hspace{5mm} $p_{\ell, \mu} = \left (p_{\ell, \mu} + \tau W_{\ell, \mu}\right ) / \left ( 1 + \tau \left\lvert W_{\ell, \mu} \right\rvert \right )$ \\
 
 
  \hspace{5mm}  $x_{m+1} = v - \beta p$ \\
  \hspace{5mm}  $t_{m+1}= 1 + \frac{\sqrt{1 + 4 t^{2}_{m}}}{2}$ \\
  \hspace{5mm}  $y_{m+1}= x_{m} +  \frac{t_{m}-1}{t_{m+1}} \left( x_{m+1} - x_{m} \right)$\\
  \hspace{5mm}  $m=m+1$ \\
  {\tt until} $ |\Phi(x_{m},\lambda_{k})- \Phi(x_{m-1},\lambda_{k}) | \leq \tau_f \Phi(x_{m},\lambda_{k})$ \\
  $u_k=x_m$\\
  \hline
 \end{tabular}
\caption{Forward Backward Algorithm (FBA). The values  of the indices $(\ell,\mu)$ are the  rows and columns of the  HR image.}
\label{tab:FBA}
\end{table}

\end{comment}




\section{Implementation notes}  \label{sec:modelling}



\subsubsection{Operators} For the following synthetic example, we simulate the LR data  from a ground-truth HR image by applying the forward model \eqref{eq:model} where the action of the blur matrix $\mathbf{H}$ is computed by assuming a Gaussian PSF with zero mean and standard deviation $\sigma_H$ which will be specified later on. As $\mathbf{S}$, we consider the discretised 2D Lanczos down-sampling operator \cite{duchon1979} inbuilt in the MATLAB function \texttt{imresize}. Finally, we consider AWGN with zero mean and standard deviation $\sigma_{\bm{\eta}}$ whose values will be made precise in the following.



\subsubsection{Comparisons}  \label{sec:competitors}
We compare our results with the ones obtained by models based on gradient-sparse regularisation such as convex isotropic TV (I-TV) \cite{MarquinaOsher2008}, non-convex capped TV (c-TV) \cite{zhang2009} and anisotropic fractional TV \cite{chen2016computing} which, for consistency, have been implemented within the same ADMM optimisation framework. We further add comparisons with the results obtained by two state-of-the-art Deep Learning-based approaches. The former is the Content Adaptive Resampler (CAR) \cite{sun2020learned} convolutional neural network, which is characterised by  a downsampler-upsampler structure. For that, we use a pre-trained model \footnote{\url{https://github.com/sunwj/CAR}} taking into account only the trained upsampler part. The latter is the Image Restoration Convolutional Neural Network (IRCNN) \cite{zhang2017}, which is a Plug and Play (PnP) method based on HQS optimisation.




\subsubsection{Initialisation, parameters and evaluation metrics}
We initialise ${\bm{u}}^{0}$ in our model as ${\bm{u}}^0 = \mathbf{S}^T {\bm{g}}$. Given the non convexity of problem \eqref{eq:min_TV0_unconstrained}, the choice of a wise initialisation is important. We tested several ones (the aforementioned one, the zero image and the I-TV initialisation) and kept the one providing the best results. 
The variables ${\bm{t}}^0,{\bm{s}}^0,{\bm{z}}^0$ as well as  $\bm{\lambda}^0_{t},\bm{\lambda}^0_{s},\bm{\lambda}^0$ in \eqref{eq:ADMM_aniso_1}-\eqref{eq:ADMM_aniso_5} and in \eqref{eq:ADMM_iso_1}-\eqref{eq:ADMM_iso_3}  were set to $\mathbf{0}$. 
To ensure convergence by Theorems \ref{th:conv_ADMManiso} and \ref{th:conv_ADMMiso}, the penalty sequences are chosen as $(\beta^k)=k(1+\epsilon)^{k}$ with $\epsilon=10^{-4}$. Note that for such small choice of $\epsilon$, $ k(1+\epsilon)^{k}\approx k$, i.e. the growth of $(\beta^k)$ is almost linear. The process is stopped when the relative change between consecutive iterates ${\bm{u}}^k$ is lower than $10^{-3}$.

For simulated data, we evaluate the quality of the SR outputs by means of Peak-Signal-to-Noise-Ratio (PSNR) and  Structure Similarity index (SSIM) as well as the Jaccard index, an evaluation metric in the range $[0,1]$ measuring the ratio between correctly detected points and false detections frequently employed in the context of microscopy imaging. We remark that choosing the right evaluation metric for SR problems is not trivial, see, e.g., \cite{Sage2019} for a review. While PSNR and SSIM are good choices to quantify reconstruction quality, the Jaccard index is more appropriate to assess correct versus false pixel localisation. 







\section{$\ell^0$ gradient-based super-resolution}  \label{sec:SRmodel}

The use of convex gradient-based regularisations for SR problems dates back to \cite{Chan2007,MarquinaOsher2008}, where TV regularisation\footnote{By $\lVert \cdot\rVert$ we denote the standard Euclidean modulus.
 \begin{equation}   \label{eq:TV_reg}
 \lVert \mathbf{D} {\bm{u}} \rVert_{1,p} = \sum_{i=1}^{N} \left( \lVert (\mathbf{D_h} {\bm{u}})_i\rVert^p + \lVert(\mathbf{D_v} {\bm{u}})_i\rVert^p\right)^{1/p},
\end{equation}}
was employed to promote sparsity on the image gradient $\mathbf{D} {\bm{u}} = (\mathbf{D_h} {\bm{u}}, \mathbf{D_v} {\bm{u}})\in \mathbb{R}^{2 \times N} $. 
Note, that for  $p\in\left\{1,2\right\}$ anisotropic/isotropic regularisation is promoted, respectively. We remark that fractional generalisations to exponents $1<p<2$ are also possible \cite{Yao2020}.


Gradient-sparsity can be enforced more severely by means of non-convex $\ell^0$ gradient smoothing, see, e.g., \cite{Xu11} and \cite{Storath14}. Using an analogous notation as in \eqref{eq:TV_reg}, for $p\in\left\{1,2\right\}$ we thus consider the $\ell^0$ gradient regularisation functional defined by:
\begin{align} \label{eq:zero_norm}
R({\bm{u}})& =\lVert \mathbf{D} {\bm{u}} \rVert_{0,p}   \\
& := \sum_{i=1}^{N} \begin{cases}
\big{|}  (\mathbf{D_{h}} {\bm{u}})_{i} \big{|}_0 +  \big{|}  (\mathbf{D_{v}}{\bm{u}})_{i}  \big{|}_0 & \text{ for } p=1,\\
\big{|} \lVert (\mathbf{D_{h}}{\bm{u}})_{i}, (\mathbf{D_{v}}{\bm{u}})_{i} \rVert \big{|}_0  &  \text{ for } p=2, \notag
\end{cases}
\end{align}

where by $|\cdot|_0$ we denote the function:
\begin{align*}
| z | _0:= \begin{cases}
 0 & z = 0  \\
1 & z \neq 0.
 \end{cases}
\end{align*}

The functional \eqref{eq:zero_norm} counts the number of \textit{jumps} of ${\bm{u}}$ in terms of the non-zero values of its gradient magnitude. In particular, in the case $p=1$ the regulariser independently counts the jumps along the two  horizontal and vertical Cartesian directions, whereas for $p=2$ the gradient magnitudes are taken into account jointly. In both cases, the term $\lVert \mathbf{D} {\bm{u}} \rVert_{0,p}$ penalizes low-amplitude structures while preserving edges in the images, thus favouring sharp piece-wise constant reconstructions which are particularly desirable for image segmentation problems. We notice that $0 \leq \lVert \mathbf{D} {\bm{u}} \rVert_{0,p} \leq 2N$ for $p\in\left\{1,2\right\}$.

In the following, we will refer to \eqref{eq:zero_norm} with $p=1$ as the \emph{anisotropic} $\ell^0$-gradient regularisation  (A-TV$^{0}$), while for $p=2$ we will refer to   \emph{isotropic} $\ell^0$-gradient regularisation (I-TV$^{0}$).



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

We report here several experiments performed on synthetic and real data. All the experiments are executed on a PC Intel(R) Core(TM) i5-6200U CPU $@$ 2.30 GHz 2.40GHz with 8.00Gb RAM using Matlab R2018b and Python 3. 
The codes are available at \url{https://github.com/pcascarano}.    

\vspace{-0.3cm}

\subsection{Computational analysis on synthetic data}\label{sec:numaspects}

We first analyse the reconstruction and the convergence properties of the proposed models/algorithms and
comment on their parameter sensitivity.

For this first example, LR data were generated by applying  \eqref{eq:model} to the HR  $428 \times 600$ grayscale image in Figure \ref{fig:butterfly} (a). Gaussian blur with $\sigma_{H}=1$ and down-sampling with factor $L=4$ were applied and AGWN with standard deviation $\sigma_{\eta}=0.01$ was added to get the LR image in Figure \ref{fig:butterfly} (b). 
In Figure \ref{fig:butterfly} (c)-(f) we report the results computed  by the anisotropic (A-$\text{TV}^0$) and isotropic (I-$\text{TV}^0$) $\ell^0$-gradient model for two different values of the regularisation parameter $\mu\in\left\{0.005,0.01\right\}$.
The jump-sparse regularisation flattens out many details in the reconstruction, promoting a  cartoon-like reconstruction which can then be used for subsequent classification and segmentation purposes: the higher the regularisation parameter $\mu$, the more simplified the reconstruction. We further add a close-up of two ROIs: the blue square contains both fine details (filaments, yellow arrows) and corner points (green arrows), the red one textured details.
The directional bias of the A-TV$^{0}$ regularisation along the horizontal and vertical direction is here clearly visible.  
We report in the captions of the Figure \ref{fig:butterfly} (c)-(f) the values $\lVert\mathbf{D}{\bm{u}}^*\rVert_{0,1}$ and $\lVert\mathbf{D}{\bm{u}}^*\rVert_{0,2}$ which corresponds to the number of gradient jumps on the output image. Note that choosing a larger $\mu$, more jump-sparsity is promoted so the number of jumps on ${\bm{u}}^*$ is smaller.

We now validate the algorithmic convergence behaviour w.r.t. to the choice of the penalty sequences $(\beta_t^k), (\beta_s^k), (\beta_k)$.  Namely, in Figure \ref{fig:beta1} (a)  and \ref{fig:beta1} (b) we report the behaviour of the objective functions $\Phi({\bm{u}}^k;\mu,p)$ in \eqref{eq:min_TV0_unconstrained} along the ADMM iterations for different choices of the penalty sequences (left).  For both cases $p=1$ and $p=2$ we choose $\beta^{k}=\beta^k_t=\beta^k_s\equiv10$ for all $k$ (blue line), $\beta^{k}=\beta^k_t=\beta^k_s=k^{0.5}$ (red line) and $\beta^{k}=\beta^k_t=\beta^k_s=k(1 + \epsilon)^{k}$ with $\epsilon=10^{-4}$ (yellow line). On the same plots we further show the decay of the quadratic data term (right). We observe that  when the penalty sequence fulfil the required growth condition then the convergence is nicely monotone. whereas for the other two choices, the decay exhibits oscillations while preserving a globally decreasing trend. Numerically, this suggests that possibly less severe growth conditions may be employed, such as a sufficiently large constant values of the penalty parameters. A further study on this is left for future research. 


To confirm the improved computational performance of our ADMM algorithm w.r.t.~  to the one proposed in \cite{Storath15} and adapted to solve the SR problem \eqref{eq:min_TV0_unconstrained}, we report in Table \ref{tab:4} a comparison table both in terms of number of iterations-to-convergence and  computational times. We stress that the poor performance of the ADMM algorithm in \cite{Storath15} is due here to the large computational cost required to solve the $\ell^0$ gradient steps via inner optimisation routines.  This, combined with the use of CG solvers (required for the SR problem under consideration as no Fourier-based approaches can be used in general) makes the overall cost much higher in comparison to our more explicit splitting.

\begin{table}[t] 
\begin{center}
\caption{\small{Iterations till convergence (\texttt{iter}) and computational time (in seconds) for different methods solving \eqref{eq:min_TV0_unconstrained}.}} 
\scalebox{0.9}{
\begin{tabular}{c|cccc}\hline
\textbf{Method}  & \cite{Storath15} & A-TV$^{0}$ & I-TV$^{0}$  \\ \hline 
\texttt{iter}   & 1905   & 63 & 59 \\
time (s) & 2866.31   & 214.83  & 195.99 \\
\hline
\end{tabular} 
}
\label{tab:4}
\end{center}
\end{table}

\input{tikz_fig/fig_analysis}



\begin{figure}[ht!]
	\centering
	\subfloat[]{{\includegraphics[width=0.24\textwidth]{images/algo_analisi/fig_anis.eps}}}
	\subfloat[]{{\includegraphics[width=0.24\textwidth]{images/algo_analisi/fig_iso.eps}}}
	\caption{\small{Values of the cost function in \eqref{eq:min_TV0_unconstrained} (left $y$-axis) and of the fidelity term (right $y$-axis)  along iterations in the two cases $\Phi ({\bm{u}}^k;\mu,1)$ (a) and $\Phi ({\bm{u}}^k;\mu,2)$ (b). The penalty sequences are chosen as  $\beta^{k}=\beta^k_t=\beta^k_s\equiv10$ (blue), $\beta^{k}=\beta^k_t=\beta^k_s=k^{0.5}$ (red), $\beta^{k}=\beta^k_t=\beta^k_s=k(1 + \epsilon)^{k}$} with $\epsilon=10^{-4}$ (yellow).}
	\label{fig:beta1}
\end{figure}
\vspace{-0.5cm}



\subsection{Real-world applications} \label{sec:Applications}
We now report the results obtained by applying the proposed model to different real-world applications where a SR version of the given LR image is required for further image analysis.

\subsubsection{QR code recognition} \label{sec:QR}

The first application we consider is the problem of QR super-resolution. As described in, e.g., \cite{Kato2011}, images of QR codes are often scanned nowadays by means of portable devices with limited resolution. Furthermore, QR scans are often taken from a distance and in non-optimal optical conditions so that blur and noise  further limit the amount of visible information, thus making SR desirable. 

For our tests, we first generate a binary QR code image of size $250 \times 250$ by using a free QR code generator \footnote{\url{https://www.qrme.co.uk/}}, then we simulate several LR acquisitions for different levels of degradation. We consider three test cases: $\sigma_{\eta}=0.01$ and $\sigma_H=1$ (TEST 1), $\sigma_{\eta}=0.05$ and $\sigma_H=1$ (TEST 2) and $\sigma_{\eta}=0.01$ and $\sigma_H=4$ (TEST 3). We compare the results obtained by our model with the ones obtained by the models in Section \ref{sec:competitors}. For each method, we select the model parameters maximising the Jaccard index. 
Furthermore, to avoid non-binary outputs (required for Jaccard index computations), we post-process the SR results by means of an adaptive Otsu thresholding and re-compute the evaluation metrics on the binarised output, see Table \ref{tab:2}.

In Figure \ref{fig:QRcode1} we report the results obtained by the different methods for the TEST 2 image before (red frame) and after (blue frame) binarisation. We observe that due to the sharp nature of the the TV$^{0}$ regularisations, the results are almost binary so they do not benefit much from the post-processing step in terms of Jaccard index values as the other methods do. In Figure \ref{fig:QRcode2} we report a zoom of the best results obtained before binarisation by all  methods starting from the TEST 3 highly corrupted LR image. 

\input{tikz_fig/fig_QR1}

\input{tikz_fig/fig_QR2}

The quantitative evaluation of the results in terms of PSNR, SSIM and Jaccard index for the three different test cases is reported in Table \ref{tab:2}. Without any binarisation, the A-TV$^{0}$ model outperforms all the others as far as the PSNR, SSIM and Jaccard indices are concerned. The simplified geometry of the QR images considered (i.e. the sole presence of horizontal/vertical edges) makes in fact this kind of data tailored for such geometrically-biased regularisations. Furthermore, the highly non-convex jump-sparsification forces the ouptut to be almost  binary, without the need of any further post-processing binarisation, as it is required by all the other regularisations to achieve comparable (if not better) quality scores. This simple example shows that the image simplification intrinsically favoured by the use of TV$^{0}$ regularisers shall limits the need of post-processing techniques. 

As far as the deep-learning results are concerned,
we remark that the CAR network in this experiment is used in a transfer learning mode, with no noisy nor blurred images observed in the training phase.
For a fairer comparison, we thus consider the IRCNN PnP network which is capable to handle different levels of degradations, although it is shown to fail in the presence of highly-degraded data, see Figure \ref{fig:QRcode2}.


\begin{table}[h] 
\begin{center}
\caption{\small{Quantitative evaluation of SR models performance on QR  for three different TEST images and methods.
By ``-B" we denote results after binarisation. In each column we colour red the best method, blue the second-best.}} 
\scalebox{0.9}{
\begin{tabular}{c|cccccc}\hline
\textbf{LR}                        & \textbf{Method} & \textbf{PSNR} & \textbf{PSNR-B} & \textbf{SSIM}  & \textbf{SSIM-B} & \textbf{Jaccard} \\ \hline 
\multirow{6}{*}{\textbf{TEST 1}}                     & I-TV$^{0}$   & 22.5199   & 29.0809   & 0.9423   & 0.9873  & 0.9980   \\ 
                           & A-TV$^{0}$ & \textcolor{red}{32.5943}   & \textcolor{blue}{35.8478}   & \textcolor{red}{0.9913}   & \textcolor{blue}{0.9989}  & \textcolor{red}{0.9999}   \\ 
                          & I-TV    & 23.3845   & 26.3357   & 0.9489   & 0.9762  & 0.9963   \\ 
                          & c-TV    & 19.4522   & \textcolor{red}{36.7496}   & 0.8849   & 0.9977  & \textcolor{blue}{0.9997}  \\
                         & A-TV$^{1/2}$    & 18.6328   & \textcolor{red}{36.7496}   & 0.8594   & 0.9989  & \textcolor{blue}{0.9997}  \\
                         
                           & CAR   & 20.2460   & 27.8163   &  0.8159  & 0.9801  & 0.9966    \\
                          & IRCNN  & \textcolor{blue}{25.0589}   & 35.3363   & \textcolor{blue}{0.9622}   & \textcolor{red}{0.9992} & 0.9995   \\ \hline
\multirow{6}{*}{\textbf{TEST 2}}                           & I-TV$^{0}$   & 19.3318   & 18.6308   & 0.8766   & 0.9156  & 0.9781   \\ 
                          & A-TV$^{0}$   & \textcolor{red}{22.6887}   & 22.6256   & \textcolor{red}{0.9242}   & 0.9653  & \textcolor{red}{0.9912}   \\ 
                          & I-TV    & 18.1101   & 18.9848   & 0.8012   & 0.9171  & 0.9798   \\ 
                          & c-TV    & 18.7331   & 21.3473   & 0.8211   & 0.9595  & 0.9882   \\
                          & A-TV$^{1/2}$   & 19.2182   & 22.5108   & 0.8664   & 0.9660  & \textcolor{blue}{0.9910}  \\
                         
                          & CAR   & 18.1320   & \textcolor{red}{26.7831}   &  0.7493  & \textcolor{blue}{0.9805}  & 0.9906    \\
                          & IRCNN  & \textcolor{blue}{21.4314}   & \textcolor{blue}{26.3968}   & \textcolor{blue}{0.9057}  & \textcolor{red}{0.9850}  & 0.9902   \\ \hline
                          
\multirow{6}{*}{\textbf{TEST 3}}     & I-TV$^{0}$  & \textcolor{blue}{18.3763}   & 19.7532   & \textcolor{blue}{0.8634}   & 0.9294  & 0.9831  \\ 
                          & A-TV$^{0}$  & \textcolor{red}{19.2908}   & \textcolor{blue}{21.9341}   & \textcolor{red}{0.8861}   & \textcolor{blue}{0.9556}  & \textcolor{blue}{0.9897}  \\ 
                          & I-TV    & 17.9552   & 20.1585   & 0.8222   & 0.9282  & 0.9846   \\ 
                          & c-TV    & 16.9580   & \textcolor{red}{22.4648}   & 0.7915   & \textcolor{red}{0.9605}  & \textcolor{red}{0.9917}   \\
                          &  A-TV$^{1/2}$    & 17.0785  & 20.6874   & 0.7706   & 0.9372  & 0.9863   \\
                         
                          & CAR   & 11.1809   & 11.5412   &  0.4057  & 0.6342  & 0.8887    \\
                          & IRCNN  & 14.2915   & 12.5640   & 0.6342   & 0.6565  & 0.9133   \\ \hline
\end{tabular} 
}
\label{tab:2}
\end{center}
\end{table}


\subsubsection{Land-cover classification} \label{sec:landcover} 
The exploitation of Multi-Spectral Images (MSIs) is fundamental in the field of land-cover mapping and classification \cite{cihlar2000}. MSIs are satellite images whose numerous channels (from 4 to 200) are acquired at a different electromagnetic spectrum bandwidth, such as visible or infrared bands, which quantifies different types of information about the objects in the recorded scene, such as their physical composition and their temperature. Existing  segmentation techniques exploit these properties to label each pixel of the MSI within a class, thus producing a final 2D labelled image. These maps are essential in many sustainability-related applications and monitoring purposes for detecting land-cover changes (e.g. deforestation) over the years at the same geographical location, which cannot be done directly by simply looking at the MSIs (see \cite{malkin2018} and references therein). Among the many existing open-source MSI datasets, we consider here e.g. the National Agriculture Imagery Program (NAIP)  \cite{maxwell2017} dataset and the Hamlin Beach State Park (HBSP)  \cite{kemker2017} dataset. The former is a collection of HR MSIs with 1 meter resolution and three RGB channels. The latter is a database of MSIs with 6 channels, 3 for the RGB and 3 for the infrared bands and is used for validating the performance of deep-learning-based segmentation algorithms aiming to differentiate land objects with analogous characteristics (e.g. a grass from a tree), see \cite{kemker2017}. %

For this problem, we apply SR methods to increase the spatial resolution of the given MSI image so as to produce an output image which could be easily segmented by standard segmentation algorithms. The need of a SR model in this specific application is justified by the physical limitations preventing HR acquisitions, such as the limited spatial resolution in some bandwidths as the infrared band \cite{mandanici2020}. On the other hand, a simplified image where noise and blur artefacts are removed comes very handy for classification purposes.  To compute the land-cover mapping on the output of the SR regularised images we use in the following a standard $k$-Means segmentation and the state-of-the-art U-Net neural network \cite{Ronneberger2015}, specifically developed for segmentation tasks. 




\input{tikz_fig/fig_landcover_1}

In the first experiment we consider a LR test image \footnote{Image identification number: M 4207221 NW 18 1 20120709} from the NAIP dataset (Figure \ref{fig:landcover1} (a)). We first run the $k$-Means algorithm directly on this image, choosing empirically the number of classes to be $k=5$. The classification obtained looks speckled and significant classification errors occur (see \ref{fig:landcover1} (b)). In Figures \ref{fig:landcover1} (c)-(d), we report the classification results obtained by applying  $k$-Means to the I-TV and I-TV$^{0}$ SR reconstructions (with $L=2$).
The segmentation results obtained on these gradient-sparsified images appear much more reliable. We notice, in particular, that some parts of the vegetation are wrongly labelled as water in the I-TV result (red boxes), whereas this is not the case for the I-TV$^{0}$ reconstruction (green boxes), due to its enhanced flattening properties.  

In the second experiment we use the I-TV$^0$ model for SR to pre-process an image from the validation set of the HBSP dataset before giving it as an input to the U-Net \cite{kemker2017}. To do so, we consider a LR MSI acquisition of size $440 \times 350 \times 6$ and apply the SR model (with $L=2$) to each channel. For comparisons, we use the U-Net both on the given LR MSI and on the computed SR reconstruction. We report the results in Figure \ref{fig:landcover2}. Note, that differently from $k$-Means, U-NET does not require the user to specify the number of required classes.  
We observe that the quality of the U-Net segmentation is significantly improved when a pre-processing with SR I-TV$^{0}$ is made. When applied to the given LR image (see Figure \ref{fig:landcover2} (a)), the U-Net is in fact not capable to differentiate the group of trees (blue) from the grass (red). Increasing the resolution and promoting sparsity on the image gradient seems to be of great help for achieving more accurate results.

\begin{figure
	\centering
	\subfloat[]{{\includegraphics[width=0.14\textwidth]{images/landcover/landcover_1.png}}} \quad
	\subfloat[]{{\includegraphics[width=0.14\textwidth]{images/landcover/landcover_2.png}}}
	\caption{\small{Results of MSI segmentation by U-NET. (a) Result on the given LR image (x2). (b) Result on the I-TV$^{0}$ reconstruction.}}
	\label{fig:landcover2}
\end{figure}







\subsubsection{Cell detection}
Standard light microscopes suffer from a limited resolving power which often causes blur artefacts and limits spatial resolution. In such conditions, the good performance of segmentation algorithms allowing for a precise location of isolated cells as well as cell clusters  is very challenging, despite their large use in biomedical applications where a fast segmentation is important for data analysis \cite{kothari2009automated}.
We test our $\ell^0$-gradient SR model on the light-microscope EVICAN data (Figure \ref{fig:cell} (a)) \cite{schwendy2020evican} for which the reference GT image in Figure \ref{fig:cell} (b)  has been found based on star-convexity shape prior \cite{gulshan2010geodesic}.
We apply the I-TV$^0$ model and its competitors on the LR acquisition obtained by \eqref{eq:model} from GT setting $L=4$, $\sigma_{H}=6$ and $\sigma_{\eta}=0.02$. For the different methods, the segmented regions are shown in Figure \ref{fig:cell} (c)-(e), while in Table \ref{tab:3} the PSNR, SSIM and Jaccard index values are reported. The $\ell^0$-gradient sparsity enforced by the I-TV$^0$ method allows for a better detection of the two isolated cells (green boxes) as well as the cell cluster (red boxes). However, when compared to I-TV, such simplification penalises more strongly image reconstruction metrics (PSNR and SSIM).




\begin{table}[h] 
\begin{center}
\caption{\small{Quantitative comparisons on cell image SR between different methods.}}
\scalebox{0.9}{
\begin{tabular}{c|ccc}\hline
\textbf{Method}  & \textbf{PSNR}  & \textbf{SSIM} & \textbf{Jaccard} \\ \hline 
I-TV   & \textcolor{red}{35.6891}   & \textcolor{red}{0.9198}   & 0.6855   \\ 
I-TV$^{1/2}$  & 35.2428                    & 0.9102                    &  0.8753      \\ 
I-TV$^{0}$ & 35.2863   & 0.9135   & \textcolor{red}{0.8778}    \\ 
CAR    & 35.1664   &  0.9044   & 0.8057    \\
\end{tabular} 
}
\label{tab:3}
\end{center}
\end{table}

\input{tikz_fig/fig_cell}




\input{tikz_fig/fig_MM1}

\subsubsection{Compressed JPG images SR}  \label{sec:JPG}
In \cite{Xu11,Ono2017,Storath14}  $\ell^0$-gradient regularisation has been used for JPG compression artefact removal. Here, we consider a scenario where such task is performed along with a resolution improvement. To do so, we consider an RGB LR cartoon-type image of size $170 \times 170$ suffering from JPG compression artefacts and with small, not discernible details, and apply the gradient-sparse SR models. As no ground truth is available for this example, for all models we empirically select the parameters producing the best visual output. 

 In Figure \ref{fig:MM} we report two close-ups of the computed SR reconstructions marked by blue and red boxes. The blue box highlights small details which are poorly discernible in the LR image, while the red box considers a patch of the face with some blunt edges and a small (but meaningful!) face mole (see green box). We see that both A-TV$^{0}$ and I-TV$^{0}$ reconstructions are sharper and more cartoonised than the ones obtained by the other models.  Furthermore, the well-known I-TV and c-TV loss of contrast reconstruction artefact makes small details hardly discernible. Due to the high-level of compression artefacts, we remark that both IRCNN and CAR results are very blurred.














\section{Convergence analysis} \label{sec:convergence}

We report here the convergence proofs of Theorem \ref{th:conv_ADMManiso} and Theorem \ref{th:conv_ADMMiso}. We provide a detailed proof for Theorem \ref{th:conv_ADMManiso} while providing only few details for Theorem \ref{th:conv_ADMMiso} whose proof is based on essentially the same arguments.

\subsection{Proof of Theorem \ref{th:conv_ADMManiso}}
We consider the ADMM sequences $({\bm{u}}^{k}),({\bm{t}}^{k}),({\bm{s}}^{k})$, defined in \eqref{eq:ADMM_aniso_1}-\eqref{eq:ADMM_aniso_3}. We want to show that there exists ${\bm{u}}^*$ such that:
\begin{align}
{\bm{u}}^{k} \rightarrow {\bm{u}}^{*}, \quad {\bm{t}}^{k} \rightarrow\mathbf{D}_{h}{\bm{u}}^{*}\nonumber, \quad {\bm{s}}^{k} \rightarrow  \mathbf{D}_{v}{\bm{u}}^{*}. 
\end{align}
\textcolor{red}{To shorten the proof, we remark that everything proved for the sequences $({\bm{t}}^{k})$,$(\beta_{t}^{k})$, $(\boldsymbol{\lambda}_{t}^{k})$ and $(\mathbf{D}_{h}{\bm{u}}^{k})$ can be deduced for $({\bm{s}}^{k})$, $(\beta_{s}^{k})$, $(\boldsymbol{\lambda}_{s}^{k})$ and $(\mathbf{D}_{v}{\bm{u}}^{k})$, respectively.}

We start defining the following functionals:
\begin{align*}
&G^{h}_{k}({\bm{t}}):= \mu \lVert {\bm{t}} \rVert_{0} + \dfrac{\beta_{t}^{k}}{2} \lVert {\bm{t}} - (\mathbf{D}_{h}{\bm{u}}^{k}+\dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})\rVert_{2}^{2}, \\
& F_{k}({\bm{u}}):= \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}} - {\bm{g}} \rVert_{2}^{2}  + \dfrac{\beta_{t}^{k}}{2}\lVert \mathbf{D}_{h}{\bm{u}} - ({\bm{t}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})\rVert^{2}_{2} + \nonumber \\  & + \dfrac{\beta_{s}^{k}}{2}\lVert \mathbf{D}_{v}{\bm{u}} - ({\bm{s}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta_{s}^{k}})\rVert^{2}_{2}.
\end{align*}


\paragraph*{Step 1}
There holds: 
\begin{align}
&\lVert {\bm{t}}^{k+1} -\mathbf{D}_{h}{\bm{u}}^{k} - \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{2} \leq \sqrt{\dfrac{2 \mu N}{\beta_{t}^{k}}}. \label{proof:step1_limitato_t}
\end{align}
\textcolor{red}{This inequality} can  be trivially shown by the minimality of ${\bm{t}}^{k+1}$ in \eqref{eq:ADMM_aniso_1} which entails   $G^{h}_{k}({\bm{t}}^{k+1}) \leq  G^{h}_{k}(\mathbf{D}_{h}{\bm{u}}^{k} + \frac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})$, therefore we get:
\begin{align*}
& \mu \lVert {\bm{t}}^{k+1}\rVert_{0} + \dfrac{\beta_{t}^{k}}{2} \lVert {\bm{t}}^{k+1} - (\mathbf{D}_{h}{\bm{u}}^{k}+\dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}})\rVert_{2}^{2} \notag\\
& \leq \mu \lVert \mathbf{D}_{h}{\bm{u}}^{k} + \frac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{0} \leq \textcolor{red}{\mu} N,
\end{align*}
by definition of $\lVert \cdot \rVert_0$. We remark that $N$ is the size of $ubold^{k}$. By neglecting the first term on the Left Hand Side (LHS)  of the above inequality \textcolor{red}{and applying the square root} we deduce \eqref{proof:step1_limitato_t}.

\paragraph*{Step 2}
From the minimality of ${\bm{u}}^{k+1}$ in \eqref{eq:ADMM_aniso_3} we have: $F_{k}({\bm{u}}^{k+1}) \leq F_{k}({\bm{u}}^{k})$ for every $k$. By definition of $F_{k}$ and applying \eqref{proof:step1_limitato_t} and \textcolor{red}{its counterpart related to the sequences $({\bm{s}}^{k})$, $(\beta_{s}^{k})$, $(\boldsymbol{\lambda}_{s}^{k})$ and $(\mathbf{D}_{v}{\bm{u}}^{k})$}, we deduce: 
\begin{align}
    & \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k+1} - {\bm{g}} \rVert_{2}^{2}  + \dfrac{\beta_{t}^{k}}{2} \lVert \mathbf{D}_{h} {\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{2}^{2} + \nonumber \\
    & + \dfrac{\beta_{s}^{k}}{2} \lVert \mathbf{D}_{v} {\bm{u}}^{k+1} - {\bm{s}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta_{s}^{k}}  \rVert_{2}^{2} \leq  \label{proof:riferimento} \\
    & \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} - {\bm{g}} \rVert_{2}^{2}  + \dfrac{\beta_{t}^{k}}{2} \lVert \mathbf{D}_{h} {\bm{u}}^{k} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{2}^{2} + \nonumber \\
    & + \dfrac{\beta_{s}^{k}}{2} \lVert \mathbf{D}_{v} {\bm{u}}^{k} - {\bm{s}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta_{s}^{k}}  \rVert_{2}^{2} \leq \dfrac{1}{2}\lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} - {\bm{g}} \rVert_{2}^{2} + 2\mu N \nonumber.
\end{align}


Since the all the terms in the LHS of \eqref{proof:riferimento} are nonnegative, also the following inequality holds: 
\begin{align}
    & \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k+1} -{\bm{g}} \rVert_{2}^{2} \leq \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} -{\bm{g}} \rVert_{2}^{2} + 2 \mu N,  \nonumber
\end{align}
which can be recursively iterated, thus getting the following relation  for every $k\geq 1$:
\begin{align*}
  \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} -{\bm{g}} \rVert_{2}^{2} \leq \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{0} - {\bm{g}} \rVert_{2}^{2} +  2 \mu N k.
\end{align*} 


From \eqref{proof:riferimento} we can also derive the following inequality (due again to the nonnegativity of all the three terms in the  LHS):
\begin{align}
    & \dfrac{\beta_{t}^{k}}{2} \lVert \mathbf{D}_{h} {\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \rVert_{2}^{2} \leq \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} -{\bm{g}} \rVert_{2}^{2} + 2 \mu N,  \nonumber
\end{align}
and using the estimate above, we can deduce the following:

\begin{align*}
    & \dfrac{\beta_{t}^{k}}{2} \lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}\rVert^{2}_{2}  \leq \nonumber \\
    &  \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{k} - {\bm{g}} \rVert_{2}^{2} + 2 \mu N \leq \dfrac{1}{2} \lVert \mathbf{S} \mathbf{H} {\bm{u}}^{0} - {\bm{g}} \rVert_{2}^{2} +  2 \mu N k , 
\end{align*}
which entails:
\begin{align}\label{ineq:series_convergent2}
\textcolor{red}{\lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}\rVert_{2}  \leq \sqrt{\frac{1}{{\beta_{t}^{k}}}}\lVert \mathbf{S} \mathbf{H} {\bm{u}}^{0} - {\bm{g}} \rVert_{2} + \sqrt{4\mu N \frac{k}{\beta_{t}^{k}}}}. 
\end{align}

By assumption \ref{cond:1_ADMManiso} we get  $(\beta_t^k)\nearrow +\infty$ and $\dfrac{k}{\beta^{k}_{t}} \rightarrow 0$, then it follows:
\begin{align} 
    &\lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}\rVert_{2} \rightarrow 0. \label{prof:bound1
   
\end{align}
   

\paragraph*{Step 3}
We show that the sequences $\mathbf{D}_{h}{\bm{u}}^{k}$ ans $\mathbf{D}_{v}{\bm{u}}^{k}$ are Cauchy sequences, hence they converge. We prove this for $\mathbf{D}_{h}{\bm{u}}^{k}$, the proof for $\mathbf{D}_{v}{\bm{u}}^{k}$ is identical.
\begin{align*}
    & \lVert \mathbf{D}_{h} {\bm{u}}^{k+1} - \mathbf{D}_{h} {\bm{u}}^{k} \rVert_{2} \leq \nonumber \\
    &\leq \lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} +  \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}  \rVert_{2}  +   \lVert \mathbf{D}_{h}{\bm{u}}^{k} - {\bm{t}}^{k+1} +  \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}  \rVert_{2} . 
\end{align*}
By assumption \ref{cond:1_ADMManiso} we get $\sum_{k=0}^{+\infty} \sqrt{\frac{1}{\beta^{k}_{t}}} < + \infty$. Therefore, from \eqref{proof:step1_limitato_t} and \eqref{ineq:series_convergent2} $\mathbf{D}_{h}{\bm{u}}^{k}$ is a Cauchy sequence, which thus converges to a point ${\bm{t}}^{*}$. Similarly, $\mathbf{D}_{v}{\bm{s}}^{k}$ which converges to some ${\bm{s}}^{*}$.


\paragraph*{Step 4}
We start rewriting \eqref{eq:ADMM_aniso_4}: 
\begin{align}
 &\frac{\boldsymbol{\lambda}_{t}^{k+1}}{\beta_{t}^{k}} = \mathbf{D}_{h} {\bm{u}}^{k+1} - {\bm{t}}^{k+1} + \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}}. \label{eq:proof_dual1
\end{align}
By \eqref{prof:bound1}, the RHS of this equation tends to $0$ as $k\to \infty$, whence we deduce that $\frac{\lVert \boldsymbol{\lambda}_{t}^{k+1} \rVert_{2}}{\sqrt{\beta_{t}^{k}}} \rightarrow 0 $. By monotonicity of the $(\beta^k_t)$ we then deduce that  $\frac{\lVert \boldsymbol{\lambda}_{t}^{k} \rVert_{2}}{\sqrt{\beta_{t}^{k}}} \rightarrow 0 $. Hence:
\begin{align*}
& \lVert \mathbf{D}_{h}{\bm{u}}^{k+1} - {\bm{t}}^{k+1} \rVert_{2} \leq  \dfrac{\lVert \boldsymbol{\lambda}_{t}^{k+1}\rVert_{2} + \lVert \boldsymbol{\lambda}_{t}^{k}\rVert_{2}}  {\sqrt{\beta_{t}^{k}}}, 
\end{align*}
where both quantities on the RHS tend to $0$ as $k\to\infty$. Therefore, ${\bm{t}}^{k}$ converges to ${\bm{t}}^{*}$ as well as $\mathbf{D}_{h}{\bm{u}}^{k}$  and ${\bm{s}}^{k}$ converges to ${\bm{s}}^{*}$ as well as $\mathbf{D}_{v}{\bm{u}}^{k}$, by uniqueness of the limit.

\paragraph*{Step 5}
We can now prove the convergence of ${\bm{u}}^{k}$ . For simplicity, let us define the quantities $\mathbf{A}:=\mathbf{S}\mathbf{H}$ and  $\mathbf{M}_{k}:=\dfrac{1}{\beta_t^{k}} \mathbf{A}^{T}\mathbf{A} + \mathbf{D}_{h}^{T}\mathbf{D}_{h} + \dfrac{\beta_{s}^{k}}{\beta_{t}^{k}}\mathbf{D}_{v}^{T}\mathbf{D}_{v}$, for every $k$. \textcolor{red}{Since \ref{cond:2_ADMManiso} holds}, we observe that the matrix $\mathbf{M}_k$ is invertible and symmetric for all $k$. The optimality condition of problem \eqref{eq:ADMM_aniso_3} reads:
\begin{equation}
\mathbf{M}_{k}{\bm{u}}^{k}= \mathbf{D}_{h}^{T} ({\bm{t}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{t}^{k}}{\beta^{k}}) + \dfrac{\beta_{s}^{k}}{\beta_{t}^{k}}\mathbf{D}_{v}^{T}({\bm{s}}^{k+1} - \dfrac{\boldsymbol{\lambda}_{s}^{k}}{\beta^{k}}) + \dfrac{1}{\beta_{t}^{k}}\mathbf{A}^{T}{\bm{g}}. \nonumber
\end{equation}

Since ${\bm{t}}^{k+1} \to {\bm{t}}^{*}$, ${\bm{s}}^{k+1} \to {\bm{s}}^{*}$,  $\frac{\boldsymbol{\lambda}_{t}^{k}}{\beta_{t}^{k}} \to 0$, $\frac{\boldsymbol{\lambda}_{s}^{k}}{\beta_{s}^{k}} \to 0$ by the previous steps, $\mathbf{D}_{h}^{T}$ and $\mathbf{D}_{v}^{T}$ are continuous operators, $\frac{\beta_{s}^{k}}{\beta_{t}^{k}} \to c$ by hypothesis, and  $\frac{1}{\beta_{t}^{k}}\mathbf{A}^{T}{\bm{g}} \to 0$ by the assumption $(\beta_t^k)\nearrow +\infty$, then the RHS converges pointwise to ${\bm{z}}^{*}=\mathbf{D}_{h}^{T}t^{*} + c \mathbf{D}_{v}^{T}s^{*} $. Therefore, $\mathbf{M}_{k}u^{k} \to {\bm{z}}^{*}$. Additionally, the sequence $\mathbf{M}_{k}^{-1}$ converges pointwise to a continuous operator $\mathbf{M}^{*}$. We thus have that $\mathbf{M}_{k}^{-1}\mathbf{M}_{k}{\bm{u}}^{k} \to \mathbf{M}^{*}{\bm{z}}^{*}$, or equivalently, ${\bm{u}}^{k} \to \mathbf{M}^{*}{\bm{z}}^{*}:={\bm{u}}^{*}$. We now want to show that ${\bm{t}}^*=\mathbf{D}_h {\bm{u}}^*$ and, similarly, that ${\bm{s}}^*=\mathbf{D}_v {\bm{u}}^*$. We show the details only for the former case. By the triangle inequality we get:
\begin{align}
\lVert {\bm{t}}^* -  \mathbf{D}_h {\bm{u}}^*\rVert_2  \leq & \lVert {\bm{t}}^* -  \mathbf{D}_h {\bm{u}}^k\rVert_2 + \lVert \mathbf{D}_h {\bm{u}}^k - \mathbf{D}_h {\bm{u}}^* \rVert_2 \nonumber \\
& \lVert {\bm{t}}^* -  \mathbf{D}_h {\bm{u}}^k\rVert_2 + \lVert \mathbf{D}_h \rVert \lVert {\bm{u}}^k - {\bm{u}}^* \rVert_2, \nonumber
\end{align}
where the both terms tend to $0$ by what shown in Step 3 and above. This concludes the proof.




\subsection{Proof of Theorem \eqref{th:conv_ADMMiso}}

The proof of Theorem \eqref{th:conv_ADMMiso} follows the same steps as the previous one. The only main difference in it is the definiton of $\mathbf{M}_{k}$, which reads in this case:
\begin{align}
    \mathbf{M}_{k}:= \dfrac{1}{\beta^{k}} \mathbf{A}^{T}\mathbf{A} + \mathbf{D}^{T}\mathbf{D} =
    \dfrac{1}{\beta^{k}} \mathbf{A}^{T}\mathbf{A}+
    \mathbf{D}_{h}^{T}\mathbf{D}_{h} + \mathbf{D}_{v}^{T}\mathbf{D}_{v}. \nonumber
\end{align}
By proceeding similarly as above the conclusion holds.






