remarkRemark hypothesisHypothesis claimClaim exampleExample
$\texttt{graphLa+}\Psi$ Bianchi et al.

# A data-dependent regularization method based on the graph Laplacian

Davide Bianchi School of Mathematics (Zhuhai), Sun Yat-sen University, Zhuhai,
519082, China<EMAIL_ADDRESS>Davide Evangelista Department of
Computer Science and Engineering, University of Bologna, Bologna, 40126, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>Stefano Aleotti
Department of Science and High Technology, University of Insubria, Como,
22100, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>Marco
Donatelli33footnotemark: 3 Elena Loli Piccolomini22footnotemark: 2 Wenbin Li
School of Science, Harbin Institute of Technology, Shenzhen, Shenzhen, 518055,
China
(). Davide Bianchi is supported by NSFC (Grant No. 12250410253). Marco
Donatelli is partially supported by MIUR - PRIN 2022 N.2022ANC8HL and GNCS-
INDAM. Wenbin Li is supported by Natural Science Foundation of Shenzhen (Grant
No. JCYJ20190806144005645) and NSFC (Grant No. 41804096<EMAIL_ADDRESS>
###### Abstract

We investigate a variational method for ill-posed problems, named
$\texttt{graphLa+}\Psi$, which embeds a graph Laplacian operator in the
regularization term. The novelty of this method lies in constructing the graph
Laplacian based on a preliminary approximation of the solution, which is
obtained using any existing reconstruction method $\Psi$ from the literature.
As a result, the regularization term is both dependent on and adaptive to the
observed data and noise. We demonstrate that $\texttt{graphLa+}\Psi$ is a
regularization method and rigorously establish both its convergence and
stability properties.

We present selected numerical experiments in 2D computerized tomography,
wherein we integrate the $\texttt{graphLa+}\Psi$ method with various
reconstruction techniques $\Psi$, including Filter Back Projection
(graphLa+FBP), standard Tikhonov (graphLa+Tik), Total Variation (graphLa+TV),
and a trained deep neural network (graphLa+Net). The $\texttt{graphLa+}\Psi$
approach significantly enhances the quality of the approximated solutions for
each method $\Psi$. Notably, graphLa+Net is outperforming, offering a robust
and stable application of deep neural networks in solving inverse problems.

###### keywords:

ill-posed problems; non-local operators; graph Laplacian; deep learning; deep
neural networks; medical imaging

47A52; 05C90; 68T07; 92C55

##  1 Introduction

We consider the model equation

(ME) $K{\boldsymbol{x}}={\boldsymbol{y}}^{\delta},$

where $K\colon X\simeq\mathbb{R}^{n}\to Y\simeq\mathbb{R}^{m}$ represents a
discretized version of a linear operator that is inherently ill-posed. Given a
fixed ${\boldsymbol{x}}_{\textnormal{gt}}\in X$ and ${\boldsymbol{y}}\coloneqq
K{\boldsymbol{x}}_{\textnormal{gt}}$, we want to recover a good approximation
of the ground-truth ${\boldsymbol{x}}_{\textnormal{gt}}$ from a noisy
observation ${\boldsymbol{y}}^{\delta}$ of ${\boldsymbol{y}}$, that is

(1)
${\boldsymbol{y}}^{\delta}\coloneqq{\boldsymbol{y}}+\boldsymbol{\eta},\qquad\|{\boldsymbol{y}}^{\delta}-{\boldsymbol{y}}\|\leq\delta,$

where $\boldsymbol{\eta}$ is a random perturbation, typically unavoidable in
practical scenarios, and $\delta>0$ is the noise intensity.

The ill-posedness of $K$ and the presence of noise force the introduction of a
regularization strategy to solve our model equation. A standard variational
method for Eq. ME reads

(2) ${\boldsymbol{x}}^{\delta}_{\alpha}\in\underset{{\boldsymbol{x}}\in
X}{\operatorname{argmin}}\left\\{\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}+\alpha\|L{\boldsymbol{x}}\|_{1}\right\\},$

where $L$ is a linear mapping, characterized by the property that
$\ker(L)\cap\ker(K)=\\{\mathbf{0}\\}$, see [21, Chapter 8] and [43]. In this
context, the term $\frac{1}{2}\|K\mathbf{x}-\mathbf{y}^{\delta}\|_{2}^{2}$
quantifies the fidelity of the reconstruction, the $\|L\mathbf{x}\|_{1}$
component serves as regularization term, and $\alpha>0$ balances the trade-off
between data fidelity and the regularization effect. Formulation (2) is a
special case of generalized Tikhonov.

Typical choices for $L$ include linear differential operators, especially for
denoising applications involving signals that are nearly piecewise-constant,
such as those encountered in imaging [29].

In the majority of the cases, $L$ is a discretization of Euclidean first or
second order linear differential operators, but the last decade has been seen
an increasing interest for models and techniques from graph theory. In [25,
26], the authors introduced for the first time graph-based non-local operators
for image segmentation and denoising problems; in [10, 14, 7, 39, 51] this
approach is further investigated in the context of image deblurring and
computerized tomography problems, and in [45, 4] it is applied to other
different image processing inverse problems.

In essence, embedding a graph-based operator $L$ in the regularization term
acts as a guiding mechanism for the overall regularization process. This
operator helps in identifying the “correct” neighborhood to concentrate the
reconstruction efforts on, by capturing specific features that can be inferred
from the observed signal ${\boldsymbol{y}}^{\delta}$.

Initially, a graph structure $G$ is constructed from a discretized signal to
incorporate features such as interfaces and discontinuities. This discrete
space, which heavily depends on the signal itself, can provide more insights
into the neighborhood where ${\boldsymbol{x}}_{\textnormal{gt}}$ resides than
a flat manifold like the Euclidean space.

Then, by choosing a suitable graph operator $L$, the optimization process in
Equation (2) is oriented towards the (supposed) neighborhood of
${\boldsymbol{x}}_{\textnormal{gt}}$.

The key point is to construct the graph from a signal that closely
approximates the primary features of ${\boldsymbol{x}}_{\textnormal{gt}}$. In
[39], it was observed that generating the graph $G$ directly from the observed
and noisy data ${\boldsymbol{y}}^{\delta}$ results in poor outcomes for
imaging tasks such as deblurring or tomographic reconstruction. This is
because ${\boldsymbol{y}}^{\delta}$ exists in a different domain compared to
${\boldsymbol{x}}_{\textnormal{gt}}$. To address this, the authors suggested
an initial preprocessing step, transforming ${\boldsymbol{y}}^{\delta}$ to
$\Psi({\boldsymbol{y}}^{\delta})$, and subsequently constructing a graph from
$\Psi({\boldsymbol{y}}^{\delta})$. This preprocessing step involves a
reconstruction map $\Psi\colon Y\to X$, from the space of observations $Y$ to
the domain $X$ where ${\boldsymbol{x}}_{\textnormal{gt}}$ lives. This can be
achieved, for example, by employing a standard Tikhonov filter [21] or the
Filter Back Projection (FBP) method [34], depending on the inverse problem to
handle. In the same spirit and independently, iterative schemes incorporating
updates to graph weights were introduced in [45, 51, 4].

Taking inspiration from [39], we introduce the novel method
$\texttt{graphLa+}\Psi$ . The initial preprocessing step involves selecting a
family of reconstructor maps

$\Psi_{\Theta}\colon Y\to X,$

where $\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$ is a family of
parameters that can depend on $\delta$ and ${\boldsymbol{y}}^{\delta}$. It is
important to note that the pair $(\Psi_{\Theta},\Theta)$ is very general and
may not be a convergent regularizing method.

Then, a graph $G$ is built from
$\Psi^{\delta}_{\Theta}\coloneqq\Psi_{\Theta}({\boldsymbol{y}}^{\delta})\in
X$, as well as its associated graph Laplacian operator
$\Delta_{\Psi^{\delta}_{\Theta}}$ (see Section 2.1 for a proper definition).
Finally, we find a minimizer of (2) with $L=\Delta_{\Psi^{\delta}_{\Theta}}$,
that is,

(3)
${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}\in\underset{{\boldsymbol{x}}\in
X}{\operatorname{argmin}}\left\\{\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}\|_{1}\right\\}.$

Let us remark that the regularization term in (3) intrinsically relies on both
the noise intensity $\delta$ and the observed data
${\boldsymbol{y}}^{\delta}$. This dependency introduces a layer of complexity
into the analysis, making it highly non-standard.

$X$$Y$${\boldsymbol{x}}_{\textnormal{gt}}$${\boldsymbol{x}}_{0}$${\boldsymbol{x}}_{\textnormal{sol}}$$\Psi^{\delta}_{\Theta}$${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}$$y$$y^{\delta}$$K$$\eta$$\Psi$graphLa+$\Psi$$\delta\to
0$$\delta\to 0$ Figure 1: A schematic representation of the
$\texttt{graphLa+}\Psi$ method. The reconstructors $\Psi_{\Theta}$ do not
necessarily need to be a regularization method, and this is represented by the
piecewise linear path of $\Psi^{\delta}_{\Theta}$ as $\delta$ goes to $0$.
However, when combined with the graph Laplacian in the Tikhonov method (2), it
generates a convergent and stable regularization operator, that is,
$\texttt{graphLa+}\Psi$, which is represented by the smooth red path. See
Sections 2 and 3 for more details on the notation.

In this study, we demonstrate that under certain, albeit very weak,
hypotheses, $\texttt{graphLa+}\Psi$ is a convergent and stable regularization
method in the limit as $\delta\to 0$. At the best of our knowledge, this is
the first time that such a theoretical analysis has been carried out.

Properties of the graph Laplacian are the key ingredients. Informally
speaking, it helps to “chain” the reconstructors $\Psi_{\Theta}$ in the
original and well-established regularization method Eq. 2. See Figure 1 for a
visual representation of the method.

For these reasons, and due to the minimal assumptions about $\Psi_{\Theta}$,
it becomes feasible to select reconstructors that may not be convergent
regularizing methods, or those whose regularization properties lack rigorous
proof yet show empirical effectiveness in certain applications. Owing to the
influence of the graph Laplacian, the overall $\texttt{graphLa+}\Psi$ method
maintains regularization and stability regardless.

In this context, a very interesting choice for $\Psi_{\Theta}$ is represented
by the class of Deep Neural Networks (DNNs). They are highly nonlinear
operators usually characterized by a paramount number of trainable parameters.
Generally speaking, the set of parameters is trained by feeding the DNNs with
a large number of data and minimizing a loss function. Thanks to the
exponential increase of dedicated computational power, DNNs achieved state-of-
the art performances in various applications. Particularly in recent years,
they have been extensively studied in the field of inverse problems, see [5]
and references therein.

However, applying DNNs to ill-posed inverse problems presents notable
drawbacks, most importantly instability and their “black-box” nature. Firstly,
DNNs are often sensitive to data perturbations and have a tendency to produce
hallucinations, i.e. false yet realistic-looking artifacts, see [3, 20].
Secondly, the complexity of their internal mechanisms, which involve millions
of parameters and nonlinear mappings, makes them challenging to understand or
explain [49].

These shortcomings contribute to a general mistrust in DNNs, particularly in
real-world situations where accuracy, stability, and reliability are
paramount, such as in medical applications. To address these concerns, we
propose the integration of the $\texttt{graphLa+}\Psi$ method with a DNN,
named graphLa+Net. This combination leverages the regularization properties of
$\texttt{graphLa+}\Psi$ , resulting in a stable, convergent, and
mathematically robust regularization method that incorporates a DNN,
effectively countering the aforementioned flaws. An initial demonstration of
the potential of graphLa+Net is presented in [11].

Our strategy with graphLa+Net follows other studies seeking to embed DNNs
within a rigorous mathematical framework, see [38, 48, 44, 12]. This approach
is crucial for the reliable application of DNNs to ill-posed inverse problems.

While the methodology and theory we develop herein apply broadly and are not
limited solely to ill-posed inverse problems in imaging, to keep the paper
self-contained we will focus on 2D computerized tomography (CT) applications.

The manuscript is organized as follows. Section 2 introduces the notation and
provides preliminaries on graph and regularization theories relevant to our
model setting. This section particularly focuses on a general approach to
generating a graph from an image. Section 3 presents the
$\texttt{graphLa+}\Psi$ method, covering its theoretical basis. In Section 4,
we specialize in the integration of the $\texttt{graphLa+}\Psi$ method with a
DNN, namely graphLa+Net. Sections 6 and 5 showcase numerical experiments
within the CT framework, employing both synthetic and real data. These
experiments demonstrate that the $\texttt{graphLa+}\Psi$ approach notably
improves the quality of approximated solutions for each method $\Psi$, with
graphLa+Net exhibiting exceptional performance, ensuring a robust and stable
implementation of DNNs for solving inverse problems. In particular, thanks to
the regularization property granted by the $\texttt{graphLa+}\Psi$ approach,
we were able to avoid any estimate of the noise distribution during the DNN
training, focusing only on the best accuracy performance of the DNN on an
unperturbed dataset.

The manuscript concludes in Section 7, where we summarize our findings and
suggest future research directions. To maintain the readability and flow of
the manuscript, the technical proofs involved in Section 3 are moved in
Appendix A.

##  2 The model setting

Hereafter, we indicate in bold any finite dimensional (column) vector, i.e.
${\boldsymbol{x}}\in X\simeq\mathbb{R}^{n}$, and for the $p$-component of
${\boldsymbol{x}}$ we write ${\boldsymbol{x}}(p)$, for $p=1,\ldots,n$. We use
the standard notation $\|\cdot\|_{r}$ for an $\ell^{r}$-norm,
$r\in[1,\infty]$.

As a baseline assumption we say that the unperturbed observation
${\boldsymbol{y}}\in Y\simeq\mathbb{R}^{m}$, which is given in Eq. 1 for
$\delta=0$, is the realization of the action of $K$ on an element
${\boldsymbol{x}}_{\textnormal{gt}}\in X$. That is, There exists
${\boldsymbol{x}}_{\textnormal{gt}}$ such that
$K{\boldsymbol{x}}_{\textnormal{gt}}={\boldsymbol{y}}$.

On the one hand, the goal is to obtain a reliable approximation of the ground
truth signal ${\boldsymbol{x}}_{\textnormal{gt}}$ when we have access only to
an observation ${\boldsymbol{y}}^{\delta}$ corrupted by unknown noise. On the
other hand, the unperturbed system Eq. ME for $\delta=0$ can be
underdetermined, the observed data ${\boldsymbol{y}}^{\delta}$ can not belong
to the range of $K$, or the (pseudo) inverse of $K$ can be very sensitive even
to small levels of $\delta>0$. For those reasons, we say that the inverse
problem of finding a solution for Eq. ME is ill-posed, and need to be
regularized. For an overview of ill-posed inverse problems and regularization
techniques, we refer the reader to [21, 28, 47].

Our analysis will heavily depend on graph structures, and so in Section 2.1 we
provide the groundwork of Graph Theory, and in Section 2.2 we explore the
relationship between images and graphs. This relationship will be exploited in
the numerical experiments of Section 6.

###  2.1 Graph Theory

Fix a finite set $P$. We indicate by $C(P)$ the set of real-valued functions
on $P$, that is,

$C(P)\coloneqq\\{{\boldsymbol{x}}\mid{\boldsymbol{x}}:P\to\mathbb{R}\\}.$

Assuming an implicit ordering of the elements in $P$, it is evident that if
the cardinality of the set $P$ is $n$ then $C(P)\simeq\mathbb{R}^{n}\simeq X$.
Hence, we use the same notation for finite-dimensional vectors, even when
referring to elements ${\boldsymbol{x}}\in C(P)$. From Section 3 we will often
make the identifications $C(P)=\mathbb{R}^{n}=X$.

For a detailed introduction to the graph setting as presented here, see [35].

######  Definition 2.1.

A _graph_ over a finite set $P$ is a pair $G=(w,\mu)$ given by:

  * •

A nonnegative _edge-weight_ function $w\colon P\times P\to[0,\infty)$ that
satisfies

    1. i)

Symmetry: $w(p,q)=w(q,p)$ for every $p,q\in P$;

    2. ii)

No loops: $w(p,p)=0$ for every $p\in P$.

  * •

A positive _node measure_ $\mu\colon P\to(0,\infty)$.

Two nodes are connected if $w(p,q)>0$, and in that case we write $p\sim q$. A
finite _walk_ is a finite sequence of nodes $\\{p_{i}\\}_{i=0}^{k}$ such that
$w(p_{i},p_{i+1})>0$ for $i=0,\ldots,k-1$. A subset $Q\subseteq P$ is
_connected_ if for every pair of nodes $p,q\in Q$ there is a finite walk such
that $p_{0}=p$, $p_{k}=q$, and each $p_{i}$ belongs to $Q$. A connected subset
$Q\subseteq P$ is a _connected component_ of $P$ if it is maximal with respect
to the ordering of inclusion.

######  Definition 2.2 (graph Laplacian).

The _graph Laplacian_ $\Delta\colon C(P)\to C(P)$ associated to the graph
$G=(w,\mu)$ is defined by the action

$\displaystyle\Delta{\boldsymbol{x}}(p)$
$\displaystyle\coloneqq\frac{1}{\mu(p)}\sum_{q\in
P}w(p,q)\left({\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)\right).$

Notice that $\Delta$ does not depend on the ordering (or labeling) of the
elements in $P$.

####  2.1.1 Construction of a graph from an element of $C(P)$

Given a finite set $P=\\{p\mid p=1,\ldots,n\\}$ , fix a distance
$\operatorname{dist}(\cdot,\cdot)$ on the set $P$ and a nonnegative function
$h_{\textnormal{d}}\colon\mathbb{R}\to[0,+\infty)$ such that
$h_{\textnormal{d}}(0)=0$. Then it follows easily that

$w_{\textnormal{d}}(p,q)\coloneqq
h_{\textnormal{d}}(\operatorname{dist}(p,q))$

is an edge-weight function on $P$, since $w_{\textnormal{d}}$ satisfies Items
i and ii from Definition 2.1. It is based on the geometric properties of $P$
induced by the distance $\operatorname{dist}(\cdot,\cdot)$. The magnitude of
the connections between two nodes $p$ and $q$ is then regulated by
$h_{\textnormal{d}}$.

Fix now an element ${\boldsymbol{x}}\in C(P)$, another nonnegative function
$h_{\textnormal{i}}\colon\mathbb{R}\to[0,+\infty)$ and finally define

(4)
$w_{{\boldsymbol{x}}}(p,q)\coloneqq\underbrace{w_{\textnormal{d}}(p,q)}_{\text{geometry}}\cdot\underbrace{h_{\textnormal{i}}(|{\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)|)}_{{\boldsymbol{x}}\text{
intensity}}.$

The edge-weight function Eq. 4 between two nodes depends on both the
“physical” distance, thanks to $w_{\textnormal{d}}$, and the “variation of
intensity” of ${\boldsymbol{x}}$, thanks to
$h_{\textnormal{i}}(|{\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)|)$. This dual
dependence has a twofold effect: it can separate nodes that reside in
different and unrelated regions of the space, and it weights the magnitude of
the connections depending on the difference of intensities of
${\boldsymbol{x}}$. For an example of application, see Section 2.2.

Now choose any strictly positive function, which may depend on
${\boldsymbol{x}}$,

(5) $\mu_{{\boldsymbol{x}}}(p)>0\quad\forall\;p\in P.$

Therefore, following Definition 2.1, for any fixed ${\boldsymbol{x}}\in C(P)$,
the pair Eq. 4-Eq. 5 is a graph $G$, on the set $P$, induced by
${\boldsymbol{x}}$. We will write then $\Delta_{{\boldsymbol{x}}}$ to indicate
the associated graph Laplacian, as per Definition 2.2.

###  2.2 Images and graphs

We briefly describe here how images and graphs can be related, following the
strategy and notation described in the previous Section 2.1.1, and provide a
practical example.

First, any image is made by the union of several pixels $p\in P$ disposed on a
grid of dimension $d_{1}\times d_{2}$, where $d_{1},d_{2}\in\mathbb{N}$. It is
then natural to identify them as ordered pairs $p=(i,j)\in\mathbb{Z}^{2}$ with
$i=1,\ldots,d_{1}$ and $j=1,\ldots,d_{2}$.

As a consequence, a reasonable choice to connect two pixels by their
“physical” distance is to set

$\operatorname{dist}(p,q)\coloneqq\|p-q\|_{\infty}\quad\mbox{and}\quad
h_{\textnormal{d}}(t)\coloneqq\mathds{1}_{(0,R]}(t),$

where $\mathds{1}_{(0,R]}$ is the indicator function of the set $(0,R]$ and
$R$ is a parameter of control which tells the maximum distance allowed for two
pixels to be neighbors. If $0<\|p-q\|_{\infty}\leq R$, then $p$ and $q$ are
connected with an edge of magnitude $1$, viz. $w_{\textnormal{d}}(p,q)=1$.

Figure 2: Simple outline of how to build a graph from an image
${\boldsymbol{x}}$. To be read from left to right. Left: a $7\times 7$ pixels
image made by orange-like and blue-like square pixels. The color intensity of
each pixel is given by the pixel-wise evaluation of a function
${\boldsymbol{x}}$. Center-left: each pixel corresponds one node, represented
by a black circle. Since the pixels are disposed on a grid, each node can be
associated to an ordered pair in $\mathbb{Z}^{2}$. Center-right: the geometric
edge-weight function $w_{\textnormal{d}}$ in Eq. 4 is given by
$\mathds{1}_{(0,1]}(\|p-q\|_{\infty})$, that is, two nodes $p,q$ are connected
if and only if $\|p-q\|_{\infty}=1$, and in that case the magnitude of the
connection is one. Right: the magnitude of an edge between two nodes is then
weighted by
$h_{\textnormal{i}}(\|{\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)\|)\in(0,1]$,
where $h_{\textnormal{i}}$ is the Gaussian function, see Eq. 6. The role of
$h_{\textnormal{i}}$ is to measure the difference of intensity between two
adjacent pixels, and it is close to zero when two pixels have very different
color intensities. This is represented by the different thickness of the edges
connecting two adjacent pixels, where a thick edge means very similar color
intensity and a thin edge means a very different color intensity.

Second, a gray-scale image is given by the light intensities of its pixels.
That is, a gray-scale image can be represented by a function
${\boldsymbol{x}}\in C(P)$ such that ${\boldsymbol{x}}(p)\in[0,1]$, where $0$
means black and $1$ means white. A common choice to weight the connection of
two different pixels by their light intensities is to use the Gaussian
function, that is

$h_{\textnormal{i}}(t)\coloneqq\rm{e}^{-\frac{t^{2}}{\sigma^{2}}},\qquad\sigma>0.$

The reason lies on the relationship between the heat kernel and the
discretization of the Laplacian on a manifold. See [13] for discussions about
the Laplacian on discretized manifolds, and [17] for another application of
the Gaussian function to define the weights of the graph Laplacian.

We arrive then at the definition of the edge-weight function in (4), applied
to our case, that is

(6)
$w_{{\boldsymbol{x}}}(p,q)=\mathds{1}_{(0,R]}(\|p-q\|_{\infty})\rm{e}^{-\frac{|{\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)|^{2}}{\sigma^{2}}}.$

For colored images, it is a little bit different. For example, in the RGB
representation, the color of a pixel is given by the combination of the light
intensities of three channels R(ed), G(reen) and B(lue). Therefore, in
principle a colored image should be regarded as a function
${\boldsymbol{x}}:P\to\mathbb{R}^{3}$, where
${\boldsymbol{x}}(p)=({\boldsymbol{x}}_{R}(p),{\boldsymbol{x}}_{G}(p),{\boldsymbol{x}}_{B}(p))$
is a vector-valued function whose elements represent the light intensities for
each channel. However, it is common to assume that the ill-posed operator $K$
acts independently on each channel, and therefore the regularization is made
on each channel separately. If for some reasons this assumption can not be
made, then we can simply modify the definition of (4) in the following way:

$w_{{\boldsymbol{x}}}(p,q)\coloneqq w_{\textnormal{d}}(p,q)\cdot
h_{\textnormal{i}}(\|{\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)\|),$

where $\|\cdot\|$ is any appropriate norm in $\mathbb{R}^{3}$. See Figure 2
for a simple example of building a graph from an image.

##  3 The $\texttt{graphLa+}\Psi$ regularization method

In this section we introduce the $\texttt{graphLa+}\Psi$ method and we show
that it is a regularization method.

Let us begin by considering a family of operators
$\\{\Psi_{\Theta}:Y\rightarrow X\\}$ which we generally refer to as
_reconstructors_. $\Psi_{\Theta}$ s not necessarily linear, and
$\Theta\in\mathbb{R}^{k}$ denotes its parameters. This family of
reconstructors is at the core of our $\texttt{graphLa+}\Psi$ method. It has
the very important role to give a first approximation of
${\boldsymbol{x}}_{\textnormal{gt}}$, since upon this approximation the
reconstruction step of our method is built .

We need to enforce some (weak) regularity on $\Psi_{\Theta}$. Assume that
there exists an element ${\boldsymbol{x}}_{0}\in X$, a nonnegative function
$g=g(\delta,{\boldsymbol{x}}_{0})$, and a parameter choice rule
$\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$ such that

$\|\Psi_{\Theta(\delta,{\boldsymbol{y}}^{\delta})}({\boldsymbol{y}}^{\delta})-{\boldsymbol{x}}_{0}\|_{2}\leq
g(\delta,{\boldsymbol{x}}_{0})\to 0\qquad\mbox{as }\delta\to 0.$

The above Section 3 will guarantee that the $\texttt{graphLa+}\Psi$ method,
which we will introduce properly in Eq. 7, is a convergent regularization
method. Let now ${\boldsymbol{y}}^{\delta}$ and
$\Theta\coloneqq\Theta(\delta,{\boldsymbol{y}}^{\delta})$ be fixed, and
$\\{\delta_{k}\\}$ and $\\{{\boldsymbol{y}}^{\delta_{k}}\\}$ be sequences such
that $\delta_{k}\to\delta$ and
${\boldsymbol{y}}^{\delta_{k}}\rightarrow{\boldsymbol{y}}^{\delta}$ for
$k\rightarrow\infty$. Write
$\Theta_{k}\coloneqq\Theta(\delta_{k},{\boldsymbol{y}}^{\delta_{k}})$. The
next hypothesis will guarantee stability for Eq. 7.

$\Theta_{k}\to\Theta\quad\mbox{and}\quad\|\Psi_{\Theta_{k}}({\boldsymbol{y}}^{\delta_{k}})-\Psi_{\Theta}({\boldsymbol{y}}^{\delta})\|_{2}\leq
f(k)\to 0\;\mbox{for }k\to\infty.$

In the next two examples, we make clear that the above assumptions are pretty
weak and they can be satisfied by several large classes of reconstructors. In
particular, the pair $(\Psi_{\Theta},\Theta)$ does not need to be a convergent
regularization method.

######  Example 3.1.

A simple example of a family of reconstructors that satisfies Section 3 is the
case when we identify them with a single, (locally) Lipschitz continuous
operator. That is, fix $\Theta\equiv\hat{\Theta}$ for every $\delta$ and
${\boldsymbol{y}}^{\delta}$, and choose an operator $\Psi_{\hat{\Theta}}$ such
that
$\|\Psi_{\hat{\Theta}}({\boldsymbol{y}}_{1})-\Psi_{\hat{\Theta}}({\boldsymbol{y}}_{2})\|_{2}\leq
c\|{\boldsymbol{y}}_{1}-{\boldsymbol{y}}_{2}\|_{2}$. Thanks to Eq. 1 and the
Lipschitz condition, Section 3 is then trivially verified with
${\boldsymbol{x}}_{0}\coloneqq\Psi_{\hat{\Theta}}({\boldsymbol{y}})$ and
$g(\delta,{\boldsymbol{x}}_{0})=c\delta$. In the same way, Section 3 is
verified by
$f(k)=c\|{\boldsymbol{y}}^{\delta_{k}}-{\boldsymbol{y}}^{\delta}\|$.

The situation of the preceding example will occur later, when
$\Psi_{\hat{\Theta}}$ is implemented as a trained DNN, as discussed in
Sections 4 and 5.3. Let us emphasize that DNNs are typically non-regularizing
methods.

######  Example 3.2.

A less trivial family of reconstructors that satisfy Section 3 is the one
given by any typical regularization operators, as per [21, Definition 3.1],
for example. Let $\Theta\in(0,\infty)$ and $\Psi_{\Theta}$ be continuous, and
$K^{\dagger}$ be the usual Moore-Penrose pseudo-inverse of $K$. Fix
${\boldsymbol{x}}_{0}\coloneqq K^{\dagger}{\boldsymbol{y}}$ and observe that
$K{\boldsymbol{x}}_{0}={\boldsymbol{y}}$, since
${\boldsymbol{y}}\in\operatorname{range}(K)$ by Section 2. By definition of
regularization operator, there exists a parameter choice rule
$\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$ such that

$g(\delta,{\boldsymbol{x}}_{0})\coloneqq\sup\\{\|\Psi_{\Theta(\delta,{\boldsymbol{y}}^{\delta})}({\boldsymbol{y}}^{\delta})-{\boldsymbol{x}}_{0}\|_{2}\mid{\boldsymbol{y}}^{\delta}\in
Y,\;\|{\boldsymbol{y}}^{\delta}-K{\boldsymbol{x}}_{0}\|_{2}\leq\delta\\}\to
0\quad\mbox{as }\delta\to 0,$

and Section 3 is then verified. About Section 3, this is a bit more involved
and depends on the regularization method itself and the parameter choice rule.
We will show a specific example in Appendix A.

The class of convergent regularization methods which fits Example 3.2 is vast,
and not necessarily restricted to the Hilbert setting. We mention Tikhonov-
Phillips methods (both ordinary and iterative variants) [31, 36, 8, 9],
popular $\ell^{p}$-$\ell^{q}$ methods with general regularization term [37,
19], and framelets [15, 32, 16].

As detailed in Section 2.1.1, we can define a graph Laplacian induced by

$\Psi_{\Theta}^{\delta}\coloneqq\Psi_{\Theta}({\boldsymbol{y}}^{\delta})\in
X\simeq C(P),$

and we will denoted it by $\Delta_{\Psi^{\delta}_{\Theta}}$. Consider then the
following Tikhonov-type variational method for Eq. ME,

(7)
${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}\in\underset{{\boldsymbol{x}}\in
X}{\operatorname{argmin}}\left\\{\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}\|_{1}\right\\}.$

We call Eq. 7 the $\texttt{graphLa+}\Psi$ method.

Let us focus our attention on the regularization term
$\mathcal{R}({\boldsymbol{x}},{\boldsymbol{y}}^{\delta})\coloneqq\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}\|_{1}$.
As a first remark, notice that, unlike typical regularization methods,
$\mathcal{R}$ depends not only on ${\boldsymbol{x}}$ but also on the data
${\boldsymbol{y}}^{\delta}$. This complicates any attempt at studying the
convergence of Eq. 7 for $\delta\to 0$.

As a second remark, indicating with $w_{\Psi_{\Theta}^{\delta}}$ the edge-
weight function Eq. 4 with ${\boldsymbol{x}}$ replaced by
$\Psi_{\Theta}^{\delta}$, we have that minimizing $\mathcal{R}$ means to force
${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}$ to be constant on
the regions where $w_{\Psi_{\Theta}^{\delta}}(p,q)$ is “large”. This is a
direct consequence of Definition 2.2,

$|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}(p)|=\frac{1}{\mu_{\Psi_{\Theta}^{\delta}}(p)}\left|\sum_{q\in
P}w_{\Psi_{\Theta}^{\delta}}(p,q)({\boldsymbol{x}}(p)-{\boldsymbol{x}}(q))\right|.$

For example, setting $w_{\Psi_{\Theta}^{\delta}}$ as in Eq. 6, then
$w_{\Psi_{\Theta}^{\delta}}(p,q)$ attains its maximum value of $1$ if and only
if $p\sim q$ and $\Psi_{\Theta}^{\delta}(p)=\Psi_{\Theta}^{\delta}(q)$.
Therefore, $\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}\|_{1}\approx 0$
if ${\boldsymbol{x}}$ is constant on the regions where
$\Psi_{\Theta}^{\delta}$ is constant.

Put simply, we can say that $\mathcal{R}$ helps to distinguish the regions of
uniformity from the regions of interfaces. In intuition, once
$\Psi_{\Theta}^{\delta}$ reconstructs the region of interfaces in
${\boldsymbol{x}}_{\textnormal{gt}}$, the final solution
${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}$ will be a good
approximation of the ground truth ${\boldsymbol{x}}_{\textnormal{gt}}$, even
though $\Psi_{\Theta}^{\delta}$ deviates from
${\boldsymbol{x}}_{\textnormal{gt}}$ in other aspects.

As a last comment, the $\ell^{1}$-norm in $\mathcal{R}$ is introduced to
enforce sparsity on the approximated solution
${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}$. This is mainly in
view of the imaging applications of Section 6. All the theory developed here
works with the $\ell^{1}$-norm replaced by any $\ell^{r}$-norm for $r>1$.

###  3.1 Convergence and stability results

In this subsection we study the regularizing properties of
$\texttt{graphLa+}\Psi$ . Since the proofs involved are quite technical, to
help the readability of this subsection we move them to Appendix A.

First, we need a definition of solution for Eq. ME. Under Section 3, let
${\boldsymbol{x}}_{0}$ and $\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$
be such that

(8) ${\boldsymbol{x}}_{0}\coloneqq\lim_{\delta\to
0}\Psi_{\Theta(\delta,{\boldsymbol{y}}^{\delta})}({\boldsymbol{y}}^{\delta}).$

######  Definition 3.3.

We call ${\boldsymbol{x}}_{\textnormal{sol}}$ a _graph-minimizing solution_
with respect to ${\boldsymbol{x}}_{0}$, defined in Eq. 8, if
$K{\boldsymbol{x}}_{\textnormal{sol}}={\boldsymbol{y}}$ and

(9)
$\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{\textnormal{sol}}\|_{1}=\min\\{\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}\|_{1}\mid{\boldsymbol{x}}\in
X,\;K{\boldsymbol{x}}={\boldsymbol{y}}\\}.$

Let us indicate with $w_{\Psi_{\Theta}^{\delta}}$ and $w_{0}$ the edge-weight
functions in Eq. 4 induced by $\Psi_{\Theta}^{\delta}$ and
${\boldsymbol{x}}_{0}$, respectively. To make our analysis work, we need three
last hypotheses which are related to properties of
$w_{\Psi_{\Theta}^{\delta}}$ and $\Delta_{\Psi^{\delta}_{\Theta}}$. As we will
discuss here, those hypotheses are easy to be satisfied.

For every $p,q\in P$, $w_{\Psi_{\Theta}^{\delta}}>0$ if and only if
$w_{0}(p,q)>0$. The following lemma is an immediate consequence of the above
hypothesis.

######  Lemma 3.4.

Under Section 3.1, there is an invariant subspace $V\subseteq C(P)$ such that
$\ker(\Delta_{\Psi^{\delta}_{\Theta}})=\ker(\Delta_{{\boldsymbol{x}}_{0}})=V$
for every $\Psi_{\Theta}^{\delta}$.

######  Proof 3.5.

See Appendix A.

The invariant subspace $V$ replaces $\ker(L)$ in the typical null-space
condition $\ker(K)\cap\ker(L)=\\{\mathbf{0}\\}$, invoked for functionals of
the form Eq. 2.

$\ker(K)\cap V=\\{\boldsymbol{0}\\}$.

The above hypotheses are not difficult to check in practice. Indeed, the null
space of the graph Laplacian is given by the constant functions on the
connected components of $P$. So, for example, from Eq. 4, if $P$ is connected
with respect to $w_{\textnormal{d}}$ and $h_{\textnormal{i}}$ is strictly
positive, then $P$ is connected with respect to the edge-weight function
$w_{{\boldsymbol{x}}}$ for any ${\boldsymbol{x}}\in C(P)$. As a consequence,
$\ker(\Delta_{\Psi^{\delta}_{\Theta}})=\\{t\boldsymbol{1}\mid
t\in\mathbb{R}\\}$ for any $\Psi_{\Theta}^{\delta}\in C(P)$, where we indicate
with $\boldsymbol{1}$ the vector whose elements are all ones.

The last assumption we need is on the function $h_{\textnormal{i}}$ in Eq. 4
and the node measure $\mu_{\boldsymbol{x}}$ in Eq. 5.

The function $h_{\textnormal{i}}$ and the map
${\boldsymbol{x}}\mapsto\mu_{{\boldsymbol{x}}}(p)$, for every fixed $p\in P$,
are locally Lipschitz.

In Section 6, we will provide specific choices of Eqs. 4 and 5 for our
numerical experiments, where we will show that Sections 3.1, 3.1, and 3.1 are
satisfied.

Throughout the remainder of this section, we will assume the validity of all
the hypotheses introduced thus far: Sections 2, 3.1, 3.1, 3, and 3.1.

The first three results are about existence and uniqueness of the graph-
minimizing solution and the well-posedness of Eq. 7.

######  Proposition 3.6.

There exists a graph-minimizing solution
${\boldsymbol{x}}_{\textnormal{sol}}$.

######  Proof 3.7.

See Appendix A.

######  Proposition 3.8.

For every fixed $\delta,\alpha>0$ and ${\boldsymbol{y}}^{\delta}\in Y$, there
exists a solution ${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}$
for the variational problem Eq. 7.

######  Proof 3.9.

See Appendix A.

######  Corollary 3.10.

If $K$ is injective, then ${\boldsymbol{x}}_{\textnormal{sol}}$ and
${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}$ are unique.

######  Proof 3.11.

See Appendix A.

######  Remark 3.12.

Without the injectivity property, uniqueness can fail. The main culprit is the
$\ell^{1}$-norm in the regularization term. However, it is possible to achieve
uniqueness in a less stringent manner, namely, for every
${\boldsymbol{y}}^{\delta}$ outside a set of negligible measure. The approach
should be in line with [18, 2], but adapted to this specific context. That
being said, relaxing the assumptions to regain uniqueness of the solutions
falls beyond the scope of the current work.

The next theorems provide convergence and stability results.

######  Theorem 3.13.

Assume that $\alpha\colon(0,+\infty)\to(0,+\infty)$ satisfies

(10a) $\lim_{\delta\to 0}\alpha(\delta)=0,$ (10b) $\lim_{\delta\to
0}\frac{\delta^{2}}{\alpha(\delta)}=0.$

Fix a sequence $\\{\delta_{k}\\}$ such that

$\lim_{k\to\infty}\delta_{k}=0,\qquad\|{\boldsymbol{y}}^{\delta_{k}}-{\boldsymbol{y}}\|_{2}\leq\delta_{k},$

and set $\alpha_{k}\coloneqq\alpha(\delta_{k})$. Let
$\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$ be the parameter choice rule
as in Section 3, and set
$\Theta_{k}\coloneqq\Theta(\delta_{k},{\boldsymbol{y}}^{\delta_{k}})$. Then
every sequence $\\{{\boldsymbol{x}}_{k}\\}$ of elements that minimize the
functional Eq. 7, with $\delta_{k}$ and $\Theta_{k}$, has a convergent
subsequence. The limit ${\boldsymbol{x}}_{\textnormal{sol}}$ of the convergent
subsequence $\\{{\boldsymbol{x}}_{k^{\prime}}\\}$is a graph-minimizing
solution with respect to ${\boldsymbol{x}}_{0}$, and

$\lim_{k^{\prime}}\|\Delta_{\Psi^{\delta_{k^{\prime}}}_{\Theta_{k^{\prime}}}}{\boldsymbol{x}}_{k^{\prime}}\|_{1}=\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{\textnormal{sol}}\|_{1}.$

If ${\boldsymbol{x}}_{\textnormal{sol}}$ is unique, then
${\boldsymbol{x}}_{k}\to{\boldsymbol{x}}_{\textnormal{sol}}$.

######  Proof 3.14.

See Appendix A.

######  Theorem 3.15.

Let now ${\boldsymbol{y}}^{\delta}$ be fixed and $\\{\delta_{k}\\}$ and
$\\{{\boldsymbol{y}}^{\delta_{k}}\\}$ be sequences such that
$\delta_{k}\to\delta$ and
${\boldsymbol{y}}^{\delta_{k}}\to{\boldsymbol{y}}^{\delta}$ for $k\to\infty$.
Then every sequence $\\{{\boldsymbol{x}}_{k}\\}$ with

${\boldsymbol{x}}_{k}\in\underset{{\boldsymbol{x}}\in
X}{\operatorname{argmin}}\left\\{\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta_{k}}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{x}}\|_{1}\right\\},$

has a converging subsequence $\\{{\boldsymbol{x}}_{k^{{}^{\prime}}}\\}$ such
that

$\lim_{k^{{}^{\prime}}}{\boldsymbol{x}}_{k^{{}^{\prime}}}\in\underset{{\boldsymbol{x}}\in
X}{\operatorname{argmin}}\left\\{\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}\|_{1}\right\\}.$

######  Proof 3.16.

See Appendix A.

##  4 graphLa+Net

Among the possible choices for reconstructors $\Psi_{\Theta}$, we propose to
use a Deep Neural Network (DNN). In this case, the set of parameters $\Theta$
contains matrices and vectors, which are the building blocks of DNNs. For a
mathematical introduction to neural networks we refer to [30]. Informally
speaking, a DNN is a long chain of compositions of affine operators and
nonlinear activation functions. The set of parameters $\Theta$ is then trained
by minimizing a loss function over a large number of data, see Section 5.3.
When considering $\Psi_{\Theta}$ as a DNN, we specify it by calling the method
graphLa+Net. A first overview of graphLa+Net was proposed in [11].

Note that, in principle, it is possible to make the network parameters
$\Theta$ dependent on the noise level $\delta$ by training it multiple times
for different values of $\delta$. However, this is rarely done in practice due
to the significant amount of time and energy consumption it would require.
This limitation has always been a crucial challenge in employing DNNs for
regularizing ill-posed problems, as it necessitates the estimate of an optimal
noise level $\delta$ which is suitable for diverse applications. Additionally,
opting for $\delta=0$ is generally not a good choice due to the typical high
sensitivity of DNNs to the noise.

Nonetheless, the regularizing property of $\texttt{graphLa+}\Psi$ effectively
addresses this issue. In the subsequent discussion, we will consider DNNs as
reconstructors with a fixed $\hat{\Theta}$, where $\hat{\Theta}$ has been
trained over a noiseless dataset. As will be shown in Sections 6.1 and 6.2,
the resulting graphLa+Net is not only regularizing and stable but also
significantly superior in performance, despite the inherent instability of the
original DNN.

###  4.1 The architecture

In our DNN model, we utilize a modified version of the U-net, as detailed in
[46]. The structure of this modified U-net is illustrated in Figure 3. U-net
is a widely recognized multi-scale Convolutional Neural Network architecture,
known for its effectiveness in processing images with global artifacts. This
fully convolutional network features a symmetrical encoder-decoder structure,
employing strided convolutions to expand its receptive field. The encoder
layers’ strides create distinct levels of resolution within the network. Each
level comprises a set number of blocks, where a block consists of a
convolutional layer with a fixed number of channels, followed by batch
normalization and a ReLU activation function. The number of convolutional
channels is doubled at each successive level, starting from a baseline number
in the first layer. Specifically, our network is designed with four levels and
a baseline of 64 convolutional channels.

As mentioned, the decoder mirrors the encoder but uses upsampling
convolutional layers in place of strided convolutions. Furthermore, to
preserve high-frequency details, skip connections link the final layer of each
encoder level to the corresponding first layer of the decoder.

The network we consider in this paper is called _Residual U-net (ResU-net)_
and it was proposed in [23]. This network is a variation of the original U-net
architecture, modified through two key changes. First, we have reconfigured
the skip connections to work as additions rather than concatenations, a
strategy aimed at reducing the total number of parameters. Second, we
introduce a residual connection that links the input and output layers
directly, which implies that the network learns the residual mapping between
the input and the expected output. The importance of the residual connection
has been observed in [27], where the authors proved that the residual manifold
containing the artifacts is easier to learn than the true image manifold.

Figure 3: A diagram of the ResU-net architecture used in the experiments.

##  5 Experimental setup

In this section we describe the setting of the $\texttt{graphLa+}\Psi$ method
in our numerical experiment.

###  5.1 Sparse view Computed Tomography

X-rays Computed Tomography (CT) is a widespread imaging system particularly
useful to detect injuries, tumors, internal bleeding, bone fractures, and
various medical conditions. It is essentially constituted by a source rotating
around an object and emitting X-rays beams from a fixed number of angles along
its arc trajectory. Passing through the interior of the object, a quantity of
radiation proportional to the density of the tissues gets absorbed and the
resulting rays are measured by a detector. The collection of all the
measurements at different projection angles is called the sinogram.

In this case, the matrix $K\in\mathbb{R}^{m\times n}$ represents the discrete
Radon operator, ${\boldsymbol{x}}\in\mathbb{R}^{n}$ is the 2D object
discretized in $n$ pixels (voxels) and
${\boldsymbol{y}}^{\delta}\in\mathbb{R}^{m}$ is the sinogram. Here, the
dimension $m$ depends on the number of projection angles $n_{\alpha}$. In
medical applications, sparse view CT denotes the case where only a small
number of angles $n_{\alpha}$ is considered, in order to reduce the radiation
dose absorbed by the patient. In this particular setting, the matrix $K$ is
under-determined (i.e. $m<n$). To faithfully simulate the real medical
scenario in our experiments we model $K$ as the fan beam Radon transform,
modeling the spread of X-rays with a fan trajectory.

###  5.2 The datasets and the test problem

We evaluate the $\texttt{graphLa+}\Psi$ algorithms using two distinct image
datasets. The first is the COULE dataset, which comprises synthetic images
with a resolution of 256x256 pixels. These images feature ellipses and lines
of varying gray intensities against a dark background. This dataset is
publicly available on Kaggle [22]. The second dataset is a subsampled version
of the AAPM Low Dose CT Grand Challenge dataset, provided by the Mayo Clinic
[40]. It contains real chest CT image acquisitions, each also at a resolution
of 256x256 pixels.

To simulate the sinogram ${\boldsymbol{y}}^{\delta}$ we consider $n_{\alpha}$
different angles evenly distributed within the closed interval $[0,179]$,
where $n_{\alpha}=120$ in the experiments with COULE dataset and
$n_{\alpha}=180$ in the experiments with Mayo dataset. The sinograms are
generated using the IRtools toolbox [24].

(a) True image ${\boldsymbol{x}}_{\textnormal{gt}}$ (b) Sinogram

Figure 4: (a): Example of a ${\boldsymbol{x}}_{\textnormal{gt}}$ image from
COULE dataset. (b): The resulting sinogram

(a) True image ${\boldsymbol{x}}_{\textnormal{gt}}$ (b) Sinogram

Figure 5: (a): Example of a ${\boldsymbol{x}}_{\textnormal{gt}}$ image from
Mayo dataset. (b): The resulting sinogram

In Figures 4 and 5, we present an example of the true image and its resulting
sinogram of size $n_{d}\times n_{\alpha}$, where
$n_{d}=\lfloor\sqrt{2}n\rfloor$ is the number of pixels of the detector. To
simulate real-world conditions, we add white Gaussian noise $\boldsymbol{\xi}$
to the sinogram at an intensity level of $\delta$, indicating that the norm of
the noise is $\delta$ times the norm of the sinogram. In particular, we
compute ${\boldsymbol{y}}^{\delta}$ as:

${\boldsymbol{y}}^{\delta}={\boldsymbol{y}}+\delta\|{\boldsymbol{y}}\|\frac{\boldsymbol{\xi}}{\|\boldsymbol{\xi}\|}.$

###  5.3 Neural Network training

We randomly select 400 pairs of images from COULE and 3,305 pairs of images
from Mayo as training sets, all of the form
$({\boldsymbol{x}}_{\textnormal{gt}},{\boldsymbol{y}}^{\delta})$. Following
the discussion in Section 4, we train the DNN over the training sets in a
supervised manner without extra noise, specifically by setting $\delta=0$. The
process involves finding $\hat{\Theta}$ that minimizes the Mean Squared Error
between the predicted reconstruction $\Psi_{\Theta}({\boldsymbol{y}})$ and the
ground-truth solution ${\boldsymbol{x}}_{\textnormal{gt}}$. Once the training
process is completed, we set $\Theta\equiv\hat{\Theta}$.

Since the ResU-net architecture is fully convolutional, the input required by
the network is an image. Hence, the input sinogram ${\boldsymbol{y}}$ has to
be pre-processed through a fast algorithm mapping the sinogram to a coarse
reconstructed image, such as FBP [41, 33] or a few iterations of a
regularizing algorithm [42, 23]. For those experiments, we choose the FBP.

The networks has been trained on an NVIDIA RTX A4000 GPU card with 16Gb of
VRAM, for a total of 50 epochs and a batch size of 10. We use Adam optimizer
with a learning rate of 0.001, $\beta_{1}=0.9$, and $\beta_{2}=0.999$, in all
the experiments.

##  6 Numerical experiments

In the numerical example, to compute the graph Laplacian we consider the edge-
weight function $\omega_{\boldsymbol{x}}$ defined as in Eq. 6, and the node
function $\mu_{\boldsymbol{x}}$ in the following form

(11) $\mu_{\boldsymbol{x}}(p)\coloneqq\sqrt{\sum_{p,q\in
P}w_{\boldsymbol{x}}^{2}(p,q)}.$

It is straightforward to check that
$h_{\textnormal{i}}=\textrm{e}^{-\frac{t^{2}}{\sigma^{2}}}$ used in Eq. 6 is
Lipschitz, and so ${\boldsymbol{x}}\mapsto\mu_{\boldsymbol{x}}(p)$ is locally
Liptschitz for every $p\in P$. Therefore, the choices Eqs. 6 and 11 satisfy
Section 3.1.

To find approximate solution of our $\ell^{2}-\ell^{1}$ model Eq. 7, we use a
Majorization–Minimization strategy combined with a Generalized Krilov Subspace
approach; we refer readers to [37] for the details of the algorithm.

This section is divided into two parts, with each part concentrating on tests
for a specific dataset. In both scenarios, we rigorously examined the
performance of the $\texttt{graphLa+}\Psi$ method to provide a comprehensive
analysis of the method’s robustness and adaptability.

More in details, we consider a wide range of reconstructors $\Psi$, including
Filter Back Projection (graphLa+FBP), standard Tikhonov (graphLa+Tik), Total
Variation (graphLa+TV), and the trained DNN (graphLa+Net) described in the
previous Section 4. In all cases, the ground truth images
${\boldsymbol{x}}_{\textnormal{gt}}$ are drawn outside of the training sets
defined for the DNN in Section 5.3.

The regularization parameters that appear in the Tikhonov and TV methods are
calculated using Generalized Cross Validation and the discrepancy principle,
respectively.

For comparison, we will include the reconstruction achieved by our method
using the ground truth image ${\boldsymbol{x}}_{\textnormal{gt}}$ as a first
approximation, labeled as
$\texttt{graphLa+}{\boldsymbol{x}}_{\textnormal{gt}}$. This serves as an upper
bound reference for the effectiveness of all the $\texttt{graphLa+}\Psi$
methods.

The quantitative results of our experiments will be measured by the
reconstruction relative error (RRE) and peak signal-to-noise ratio (PSNR),
where

$\operatorname{RRE}({\boldsymbol{x}})\coloneqq\frac{\|{\boldsymbol{x}}_{\textnormal{gt}}-{\boldsymbol{x}}\|^{2}}{\|{\boldsymbol{x}}_{\textnormal{gt}}\|^{2}},\qquad\operatorname{PSNR}({\boldsymbol{x}})\coloneqq
20\log_{10}\left(\frac{255}{\|{\boldsymbol{x}}_{\textnormal{gt}}-{\boldsymbol{x}}\|}\right),$

and by the structural similarity index (SSIM) [50]. Finally, all the numerical
tests are replicable and the codes can be downloaded from [1].

###  6.1 Example 1: COULE

In this first example, we tested our proposal on an image of the COULE test
set corrupted with some white Gaussian noise with level intensity of $2\%$.
Regarding the parameter selection for the edge-weight function in Eq. 6, we
chose $R=5$ and $\sigma=10^{-3}$.

Initial reconstructors $\Psi$ | RRE | SSIM | PSNR
---|---|---|---
FBP | 0.0761 | 0.2219 | 22.3671
Tik | 0.0478 | 0.4607 | 24.4129
TV | 0.0328 | 0.8275 | 29.6926
Net | 0.0164 | 0.9461 | 35.6828
$\texttt{graphLa+}\Psi$ |
graphLa+FBP | 0.0279 | 0.7594 | 31.0908
graphLa+Tik | 0.0272 | 0.9380 | 31.3000
graphLa+TV | 0.0162 | 0.9822 | 35.8047
graphLa+Net | 0.0098 | 0.9844 | 40.1793
$\texttt{graphLa+}{\boldsymbol{x}}_{\textnormal{gt}}$ | 0.0043 | 0.9866 | 47.2474
Table 1: Quality of initial and final reconstruction for different $\Psi$ for
the COULE dataset.

The quality of the reconstructions achieved with different operators $\Psi$ is
presented in Table 1. Specifically, the upper part of the table displays the
values of the initial reconstructions while the lower part displays the values
obtained by $\texttt{graphLa+}\Psi$ combined with the corresponding initial
reconstructor $\Psi$.

Notably, the $\texttt{graphLa+}\Psi$ method results in a substantially greater
improvement across all metrics for all initial reconstructors $\Psi$, with the
highest performance attained by graphLa+Net.

As further confirmation, Figure 6 displays the reconstructions obtained for
different $\Psi$. The level of details and sharpness in the graphLa+Net image
is incomparable with all the other methods.

(a) FBP (b) Tik (c) TV (d) Net

(e) graphLa+FBP (f) graphLa+Tik (g) graphLa+TV (h) graphLa+Net

Figure 6: Initial and final reconstructions using $\texttt{graphLa+}\Psi$
method for different $\Psi$.

As a final investigation into the capabilities of our proposal, we tested the
stability of the method against varying noise intensity. In Figure 7, we
present the PSNR and SSIM values for various levels of noise. Similar to the
previous analysis, the purple line represents the reconstruction obtained by
utilizing the true image ${\boldsymbol{x}}_{\textnormal{gt}}$ to compute the
graph Laplacian. Notably, even though our neural network was trained with a
$0\%$ noise level, the graphLa+Net method consistently outperforms all other
cases across all noise levels. Additionally, Figure 7 shows that the
integration of the graph Laplacian with the DNN serves as an effective
regularization method, in contrast to the standalone application of the DNN.

(a) PSNR

(b) SSIM

Figure 7: PSNR and SSIM for different level of noises and different
reconstructors $\Psi$ for the COULE dataset.

###  6.2 Example 2: Mayo

In this second example, we tested our proposal on an image of the Mayo test
set corrupted with some white Gaussian noise. Regarding the parameter
selection for the edge-weight function in Eq. 6, we chose $R=5$ and
$\sigma=2\times 10^{-4}$, except for the graphLa+Net method for which we used
$R=3$ and $\sigma=10^{-3}$.

Given that the Mayo dataset reflects a real-world scenario, the ground truth
images ${\boldsymbol{x}}_{\textnormal{gt}}$, used for generating sinograms and
for comparison, are not the actual true images. Instead, they are
reconstructions themselves, inherently containing some level of noise.
Consequently, comparing metrics for a fixed level of additional noise, as done
in Table 1, is a bit less informative. Instead, it is still interesting to
evaluate the $\texttt{graphLa+}\Psi$ method across different reconstructors
$\Psi$ and various levels of noise intensity using both PSNR and SSIM metrics,
see Figure 8. As previously noted in the COULE example, the graphLa+Net method
consistently outperforms all other cases, even if the neural network was
trained with a $0\%$ noise level.

(a) PSNR

(b) SSIM

Figure 8: PSNR and SSIM for different level of noises and different
reconstructors $\Psi$ for the Mayo dataset

In Figure 9, we present a visual inspection of some of the reconstructions for
$\delta=1\%$. Notably, the graphLa+Net image exhibits the sharpest quality
compared to all other cases, and being very close to the upper limit given by
$\texttt{graphLa+}{\boldsymbol{x}}_{\textnormal{gt}}$.

(a) True image ${\boldsymbol{x}}_{\textnormal{gt}}$ (b) Tik (c) TV (d) Net

(e) $\texttt{graphLa+}{\boldsymbol{x}}_{\textnormal{gt}}$ (f) graphLa+Tik (g)
graphLa+TV (h) graphLa+Net

Figure 9: Initial and final reconstructions using $\texttt{graphLa+}\Psi$
method for different $\Psi$.

As additional confirmation, in Figure 10 we zoom on the central part of the
considered image. In this way, we can clearly note that the graphLa+Net
approach achieves also an extraordinary quality of detail in the
reconstruction.

(a) True image ${\boldsymbol{x}}_{\textnormal{gt}}$ (b) graphLa+Tik (c)
graphLa+TV (d) graphLa+Net

Figure 10: Zoom in of the central lower part for different methods.

##  7 Conclusions

In this work we introduced and analyzed a novel regularization method which
uses a graph Laplacian operator built upon a first approximation of the
solution through a reconstructor $\Psi$. We demonstrated that under certain,
albeit very weak, hypotheses on the recontructor $\Psi$,
$\texttt{graphLa+}\Psi$ is a convergent and stable regularization method. In
all the proposed numerical examples, $\texttt{graphLa+}\Psi$ greatly enhanced
the quality of the reconstructions for every initial reconstructor $\Psi$.

Moreover, leveraging the regularization property of $\texttt{graphLa+}\Psi$,
we suggested to use a DNN as first reconstructor $\Psi$. This new hybrid
method, called graphLa+Net, merges the regularity of a standard variational
method and the accuracy of a DNN, giving as a result a stable regularization
method with very high accuracy.

In forthcoming works we will focus more on sharpening the choices for the
edge-weight function $w_{\boldsymbol{x}}$ such as building an automatic rule
for estimating its parameters.

##  Appendix A Theoretical analysis

In this appendix we collect the proofs of the results presented in Section
3.1. We will assume that all the hypotheses introduced so far are satisfied,
i.e. Sections 3.1, 2, 3.1, 3.1, and 3.

We begin with the invariance property of the null space of the graph
Laplacians built upon $\Psi^{\delta}_{\Theta}$.

Lemma 3.4.

######  Proof A.1.

The null space of a generic graph Laplacian, as per Definition 2.2, is given
by the subspace of functions which are constant on the connected components of
the node set $P$. By Section 3.1, it is easy to check that a sequence
$\\{p_{i}\\}_{i=0}^{k}$ is a walk with respect to $w_{\Psi^{\delta}_{\Theta}}$
if and only if it is a walk with respect to $w_{0}$. Therefore, all the
connected components of $P$, identified by $w_{\Psi^{\delta}_{\Theta}}$ and
$w_{0}$, are invariant. This concludes the proof.

###  A.1 Existence of solutions and well-posedness of $\texttt{graphLa+}\Psi$

Define

(12)
$\Gamma({\boldsymbol{x}})\coloneqq\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}\|_{1}.$

We recall that a (nonnegative) functional $\Gamma\colon X\to[0,\infty)$ is
_coercive_ if $\Gamma({\boldsymbol{x}})\to\infty$ for
$\|{\boldsymbol{x}}\|\to\infty$, where $\|\cdot\|$ can be any norm on
$X\simeq\mathbb{R}^{n}$.

Thanks to Section 3.1, the functional Eq. 12 is coercive, for every fixed
$\alpha>0$ and $\delta\geq 0$. Although this should certainly be well-known,
for the convenience of the reader we include a short proof in the next lemma.

######  Lemma A.2.

$\Gamma$ is coercive for every fixed $\alpha>0$ and $\delta\geq 0$.

######  Proof A.3.

Let $V$ be the invariant null space of $\Delta_{\Psi^{\delta}_{\Theta}}$ from
Lemma 3.4, and let us indicate with $\pi$ and $\pi_{\perp}$ the projection
into $V$ and $V^{\perp}$, respectively. In general, it holds that

(13) $\qquad\inf_{\begin{subarray}{c}{\boldsymbol{u}}\in V^{\perp}\\\
{\boldsymbol{u}}\neq\boldsymbol{0}\end{subarray}}\frac{\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{u}}\|_{1}}{\|{\boldsymbol{u}}\|_{1}}\geq\gamma_{1}>0.$

Since $\ker(K)\cap V=\\{\mathbf{0}\\}$ by Section 3.1, then it holds that too

(14) $\qquad\inf_{\begin{subarray}{c}{\boldsymbol{v}}\in V\\\
{\boldsymbol{v}}\neq\boldsymbol{0}\end{subarray}}\frac{\|K{\boldsymbol{v}}\|_{2}}{\|{\boldsymbol{v}}\|_{2}}\geq\gamma_{2}>0.$

Fix a sequence $\\{{\boldsymbol{x}}_{j}\\}$ such that
$\|{\boldsymbol{x}}_{j}\|_{2}\to\infty$. We want to show that
$\Gamma({\boldsymbol{x}}_{j})\to\infty$. Clearly, for $j\to\infty$

$\|{\boldsymbol{x}}_{j}\|_{2}\to\infty\quad\mbox{if and only
if}\quad\|\pi{\boldsymbol{x}}_{j}\|_{2}^{2}+\|\pi_{\perp}{\boldsymbol{x}}_{j}\|_{1}\to\infty.$

There are two cases:

  1. (i)

$\lim_{j}\|\pi_{\perp}{\boldsymbol{x}}_{j}\|_{1}=\infty$;

  2. (ii)

$\liminf_{j}\|\pi_{\perp}{\boldsymbol{x}}_{j}\|_{1}\leq c<\infty$ and
$\lim_{j}\|\pi{\boldsymbol{x}}_{j}\|_{2}^{2}=\infty$.

If we are in i, then by Eq. 13

$\displaystyle\alpha\gamma_{1}\|\pi_{\perp}{\boldsymbol{x}}_{j}\|_{1}$
$\displaystyle\leq\alpha\|\Delta_{\Psi^{\delta}_{\Theta}}\pi_{\perp}{\boldsymbol{x}}_{j}\|_{1}$
(15)
$\displaystyle\leq\frac{1}{2}\|K{\boldsymbol{x}}_{j}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}_{j}\|_{1}=\Gamma({\boldsymbol{x}}_{j}).$

If we are in ii, then by Eq. 14

$\displaystyle\frac{\gamma_{2}}{4}\|\pi{\boldsymbol{x}}_{j}\|_{2}^{2}\leq\frac{1}{4}\|K\pi{\boldsymbol{x}}_{j}\|_{2}^{2}$
$\displaystyle\leq\frac{1}{2}\|K{\boldsymbol{x}}_{j}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}_{j}\|_{1}+\frac{1}{2}\|K\pi_{\perp}{\boldsymbol{x}}_{j}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}$
(16)
$\displaystyle=\Gamma({\boldsymbol{x}}_{j})+\frac{1}{2}\|K\pi_{\perp}{\boldsymbol{x}}_{j}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}.$

Notice that from ii it follows that
$\liminf_{j}\|K\pi_{\perp}{\boldsymbol{x}}_{j}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}$
is bounded.

Passing to the $\liminf$ in both Eqs. 15 and 16 we conclude.

######  Remark A.4.

The same result as in Lemma A.2 applies when replacing
$\Psi^{\delta}_{\Theta}$ with ${\boldsymbol{x}}_{0}$ in $\Gamma$, and the
proof remains unchanged.

Proposition 3.6.

######  Proof A.5.

The existence relies on standard topological arguments, but we provide a
comprehensive exposition for the convenience of the reader.

Let

$c\coloneqq\inf\\{\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}\|_{1}\mid{\boldsymbol{x}}\in
X,\;K{\boldsymbol{x}}={\boldsymbol{y}}\\},$

which is well-defined thanks to Section 2. Therefore there exists a sequence
$\\{{\boldsymbol{x}}_{j}\\}$ such that
$K{\boldsymbol{x}}_{j}={\boldsymbol{y}}$ for every $j$ and
$\lim_{j}\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{j}\|_{1}=c$. In
particular, there exists $c_{1}>0$ such that

(17) $\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{j}\|_{1}\leq
c_{1}\quad\mbox{for every }j.$

There are two possible cases:

  1. (i)

$\|{\boldsymbol{x}}_{j}\|_{2}\leq c_{2}$ for some $c_{2}>0$, for every $j$;

  2. (ii)

There exists a subsequence $\\{{\boldsymbol{x}}_{j^{\prime}}\\}$, such that
$\lim_{j^{\prime}}\|{\boldsymbol{x}}_{j^{\prime}}\|_{2}=\infty$.

If we are in case i, then by compactness and continuity arguments we can
conclude that there exists ${\boldsymbol{x}}_{\textnormal{sol}}$ such that

$\lim_{j^{\prime}}{\boldsymbol{x}}_{j^{\prime}}={\boldsymbol{x}}_{\textnormal{sol}},\quad\mbox{and}\quad\begin{cases}K{\boldsymbol{x}}_{\textnormal{sol}}={\boldsymbol{y}},\\\
\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{\textnormal{sol}}\|_{1}=c,\end{cases}$

that is, ${\boldsymbol{x}}_{\textnormal{sol}}$ is a graph-minimizing solution
with respect to ${\boldsymbol{x}}_{0}$.

Suppose now we are in case ii. By Lemma A.2 and Remark A.4,
$\Gamma({\boldsymbol{x}}_{j^{\prime}})\to\infty$ for any fixed
$\alpha,\delta$. Since
$\|K{\boldsymbol{x}}_{j^{\prime}}-{\boldsymbol{y}}^{\delta}\|^{2}_{2}=\|{\boldsymbol{y}}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}\leq\delta^{2}$
for every $j^{\prime}$, it follows necessarily that
$\lim_{j^{\prime}}\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{j^{\prime}}\|_{1}=\infty$.
This leads to an absurdity in light of (17).

Proposition 3.8.

######  Proof A.6.

From Lemma A.2, the nonnegative functional $\Gamma$ is coercive on a finite
dimensional vector space. By standard theory, there exists a minimizer, see
for example [6, Proposition 11.15].

Corollary 3.10.

######  Proof A.7.

The uniqueness of ${\boldsymbol{x}}_{\textnormal{sol}}$ is straightforward. On
the other hand, the uniqueness of
${\boldsymbol{x}}^{\delta}_{\Psi^{\delta}_{\Theta},\alpha}$ is derived from
the fact that if $K$ is injective then the functional
${\boldsymbol{x}}\mapsto\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta}\|_{2}^{2}$
is strongly convex. This property leads to the strong (and therefore strict)
convexity of $\Gamma$. According to [6, Corollary 11.9], the desired result
follows.

###  A.2 Convergence and stability analysis

In this subsection we prove Theorems 3.13 and 3.15. The main difficulty is
given by the regularization term
$\mathcal{R}({\boldsymbol{x}},{\boldsymbol{y}}^{\delta})$, that depends on the
observed data ${\boldsymbol{y}}^{\delta}$. Therefore, all standard techniques
can not be applied straightforwardly.

A crucial role will be played by the following two lemmas, which guarantee
uniform convergence of $\Delta_{\Psi^{\delta}_{\Theta}}$, and a special
uniform coercivity property for $\Gamma$ in Eq. 12.

######  Lemma A.8.

Let $\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$, ${\boldsymbol{x}}_{0}$
and $g$ be defined as in Section 3. For every ${\boldsymbol{x}}\in X$ it holds
that

$\|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}-\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}\|_{1}\leq
c\|{\boldsymbol{x}}\|_{\infty}g(\delta,{\boldsymbol{x}}_{0})\to 0\quad\mbox{as
}\delta\to 0,$

where $c$ is a positive constant independent of ${\boldsymbol{x}}$. In
particular, the convergence is uniform on every compact subset.

######  Proof A.9.

We use the notations introduced in Sections 2.1.1 and 3.1. In particular,
$w_{\Psi_{\Theta}^{\delta}}$ and $w_{0}$ are the edge-weight functions in Eq.
4 induced by
$\Psi^{\delta}_{\Theta}\coloneqq\Psi_{\Theta}({\boldsymbol{y}}^{\delta})$ and
${\boldsymbol{x}}_{0}$, respectively; $\mu_{\Psi_{\Theta}^{\delta}}$ and
$\mu_{0}$ are the node measures in Eq. 5 induced by $\Psi^{\delta}_{\Theta}$
and ${\boldsymbol{x}}_{0}$, respectively. Let us observe that, with this
notation, Section 3 now reads

$\|\Psi^{\delta}_{\Theta}-{\boldsymbol{x}}_{0}\|_{2}\leq
g(\delta,{\boldsymbol{x}}_{0})\quad\mbox{as }\delta\to 0.$

Define

$t_{\delta,p,q}\coloneqq|\Psi^{\delta}_{\Theta}(p)-\Psi^{\delta}_{\Theta}(q)|,\quad
t_{0,p,q}\coloneqq|{\boldsymbol{x}}_{0}(p)-{\boldsymbol{x}}_{0}(q)|,\quad
d\coloneqq\max_{p,q}\\{w_{\textnormal{d}}(p,q)\\}.$

Then,

$\displaystyle|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}(p)-\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}(p)|$
$\displaystyle\leq\sum_{q\in
P}\left|\frac{w_{\Psi_{\Theta}^{\delta}}(p,q)}{\mu_{\Psi_{\Theta}^{\delta}}(p)}-\frac{w_{0}(p,q)}{\mu_{0}(p)}\right||{\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)|$
$\displaystyle=\sum_{q\in
P}w_{\textnormal{d}}(p,q)\left|\frac{h_{\textnormal{i}}(t_{\delta,p,q})}{\mu_{\Psi_{\Theta}^{\delta}}(p)}-\frac{h_{\textnormal{i}}(t_{0,p,q})}{\mu_{0}(p)}\right||{\boldsymbol{x}}(p)-{\boldsymbol{x}}(q)|$
(18)
$\displaystyle\leq\frac{2d\|{\boldsymbol{x}}\|_{\infty}}{\mu_{\Psi_{\Theta}^{\delta}}(p)\mu_{0}(p)}\sum_{q\in
P}\left|\mu_{0}(p)h_{\textnormal{i}}(t_{\delta,p,q})-\mu_{\Psi_{\Theta}^{\delta}}(p)h_{\textnormal{i}}(t_{0,p,q})\right|.$

We need to bound the right-hand side of the last inequality. Because of the
locally Lipschitz continuity of the maps $h_{\textnormal{i}}$ and
$\Psi_{\Theta}^{\delta}\mapsto\mu_{\Psi_{\Theta}^{\delta}}(p)$, granted by
Section 3.1, then for $\delta$ small enough there exists a positive constant
$\bar{c}$, independent of $p$, such that

$\frac{2d}{\mu_{\Psi_{\Theta}^{\delta}}(p)\mu_{0}(p)}\leq\bar{c},\quad|\mu_{\Psi_{\Theta}^{\delta}}(p)-\mu_{0}(p)|\leq\bar{c}\|\Psi_{\Theta}^{\delta}-{\boldsymbol{x}}_{0}\|_{2},\quad\mbox{and}\quad|h_{\textnormal{i}}(t_{\delta,p,q})-h_{\textnormal{i}}(t_{0,p,q})|\leq\bar{c}|t_{\delta,p,q}-t_{0,p,q}|.$

Therefore, by adding and subtracting the auxiliary term
$\mu_{0}(p)h_{\textnormal{i}}(t_{0,p,q})$, it holds that

(19)

Let us bound now $|t_{\delta,p,q}-t_{0,p,q}|$:

$\displaystyle|t_{\delta,p,q}-t_{0,p,q}|$
$\displaystyle=\left||\Psi^{\delta}_{\Theta}(p)-\Psi^{\delta}_{\Theta}(q)|-|{\boldsymbol{x}}_{0}(p)-{\boldsymbol{x}}_{0}(q)|\right|$
$\displaystyle\leq\left|\left(\Psi^{\delta}_{\Theta}(p)-{\boldsymbol{x}}_{0}(p)\right)+\left({\boldsymbol{x}}_{0}(q)-\Psi^{\delta}_{\Theta}(q)\right)\right|$
$\displaystyle\leq\left|\Psi^{\delta}_{\Theta}(p)-{\boldsymbol{x}}_{0}(p)\right|+\left|\Psi^{\delta}_{\Theta}(q)-{\boldsymbol{x}}_{0}(q)\right|$
(20) $\displaystyle\leq
2\sqrt{n}\|\Psi^{\delta}_{\Theta}-{\boldsymbol{x}}_{0}\|_{2}\leq
2\sqrt{n}g(\delta,{\boldsymbol{x}}_{0}).$

Combining Eqs. 18, 19, and 20, we get that

$|\Delta_{\Psi^{\delta}_{\Theta}}{\boldsymbol{x}}(p)-\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}(p)|\leq
n\bar{c}^{2}\left(2\sqrt{n}\max_{p}\\{\mu_{0}(p)\\}+\max_{p,q}\\{h_{\textnormal{i}}(t_{0,p,q})\\}\right)\|{\boldsymbol{x}}\|_{\infty}g(\delta,{\boldsymbol{x}}_{0}),$

and then the thesis easily follows.

Let us point out that the constants that appears in the proof depends on the
dimension $n$. Appropriate choices of $w_{\textnormal{d}}$,
$h_{\textnormal{i}}$ and $\mu_{\boldsymbol{x}}$, with respect to $n$, could
make $c$ independent of $n$ itself. However, there is still the issue of the
convergence of the (sequence of) graph Laplacians $\Delta_{\boldsymbol{x}}$ as
$n\to\infty$, and this is beyond the scope of this work.

######  Lemma A.10.

Let $\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$ be the parameter choice
rule as in Section 3 and let $\delta_{k}\to 0$. Write
$\Theta_{k}\coloneqq\Theta(\delta_{k},{\boldsymbol{y}}^{\delta_{k}})$ and fix
$\alpha>0$. For any sequence $\\{{\boldsymbol{x}}_{k}\\}$ such that
$\limsup_{k}\|{\boldsymbol{x}}_{k}\|_{2}=\infty$, then

$\limsup_{k}\frac{1}{2}\|K{\boldsymbol{x}}_{k}-{\boldsymbol{y}}\|_{2}^{2}+\alpha\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{x}}_{k}\|_{1}=\infty.$

######  Proof A.11.

Let $V$ be the invariant null space from Lemma 3.4. Define

$\gamma_{k}\coloneqq\inf_{\begin{subarray}{c}{\boldsymbol{u}}\in V^{\perp}\\\
\|{\boldsymbol{u}}\|_{1}=1\end{subarray}}\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{u}}\|_{1}\quad\mbox{and}\quad\gamma_{0}\coloneqq\inf_{\begin{subarray}{c}{\boldsymbol{u}}\in
V^{\perp}\\\
\|{\boldsymbol{u}}\|_{1}=1\end{subarray}}\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{u}}\|_{1}>0.$

By Lemma A.8 and uniform convergence on compact sets, it holds that

$\forall\,{\boldsymbol{u}}\in V^{\perp}\mbox{ s.t.
}\|{\boldsymbol{u}}\|_{1}=1,\quad\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{u}}\|_{1}-\frac{\gamma_{0}}{2}\leq\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{u}}\|_{1}\qquad\mbox{for
}k\geq N=N(\gamma_{0}).$

Therefore,

$\inf_{\begin{subarray}{c}{\boldsymbol{u}}\in V^{\perp}\\\
\|{\boldsymbol{u}}\|_{1}=1\end{subarray}}\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{u}}\|_{1}\geq\frac{\gamma_{0}}{2}\qquad\forall
k\geq N(\gamma_{0}).$

In particular, there exists $\gamma_{1}>0$ such that

$\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{x}}\|_{1}\geq\gamma_{1}\|\pi_{\perp}{\boldsymbol{x}}\|_{1}\qquad\forall\,{\boldsymbol{x}},\;\forall\,k.$

The rest of the proof follows like in Lemma A.2.

The next theorem presents a convergence result for the $\texttt{graphLa+}\Psi$
method Eq. 7. The overall proof follows a fairly standard approach, as for
example in [47, Theorem 3.26]. However, since the regularizing term in Eq. 7
depends on the data ${\boldsymbol{y}}^{\delta}$ as well, it involves a few
nontrivial technical aspects.

We will use a slight modification of the notation introduced in Eq. 12.
Specifically,

$\Gamma_{k}({\boldsymbol{x}})\coloneqq\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta_{k}}\|_{2}^{2}+\alpha_{k}\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{x}}\|_{1}.$

Theorem 3.13.

######  Proof A.12.

The sequence $\\{{\boldsymbol{x}}_{k}\\}$ is well-posed thanks to Proposition
3.8. Fix a graph-minimizing solution ${\boldsymbol{x}}_{\textnormal{sol}}$ as
in Definition 3.3, which exists because of Proposition 3.6. Then, by
definition, it holds that

(21) $\displaystyle\Gamma_{k}({\boldsymbol{x}}_{k})$
$\displaystyle\leq\Gamma_{k}({\boldsymbol{x}}_{\textnormal{sol}})\leq\frac{\delta_{k}^{2}}{2}+\alpha_{k}\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{\textnormal{sol}}\|_{1}+\alpha_{k}c\|{\boldsymbol{x}}_{\textnormal{sol}}\|_{\infty}g(\delta_{k},{\boldsymbol{x}}_{0})\to
0$

as $\delta\to 0$, where the last inequality comes from Lemma A.8, and the
convergence to zero is granted by Eq. 10a. Therefore,
$\|K{\boldsymbol{x}}_{k}-{\boldsymbol{y}}^{\delta_{k}}\|_{2}^{2}\to 0$, and in
particular

(22)
$\|K{\boldsymbol{x}}_{k}-{\boldsymbol{y}}\|_{2}\leq\|K{\boldsymbol{x}}_{k}-{\boldsymbol{y}}^{\delta_{k}}\|_{2}+\delta_{k}\to
0\quad\mbox{as }\delta_{k}\to 0.$

Since

$\alpha_{k}\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{x}}_{k}\|_{1}\leq\Gamma_{k}({\boldsymbol{x}}_{k}),$

then by Eqs. 10b and 21, we get that

(23)
$\limsup_{k}\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{x}}_{k}\|_{1}\leq\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{\textnormal{sol}}\|_{1}.$

Let $\alpha^{+}\coloneqq\max\\{\alpha_{k}\mid k\in\mathbb{N}\\}$. Then,
combining Eqs. 22 and 23,

$\frac{1}{2}\|K{\boldsymbol{x}}_{k}-{\boldsymbol{y}}\|_{2}^{2}+\alpha^{+}\|\Delta_{\Psi^{\delta_{k}}_{\Theta_{k}}}{\boldsymbol{x}}_{k}\|_{1}\leq
c<\infty,$

and from Lemma A.10 we deduce that $\\{{\boldsymbol{x}}_{k}\\}$ is bounded.
Therefore, there exists a convergent subsequence
$\\{{\boldsymbol{x}}_{k^{\prime}}\\}$ which converges to a point
${\boldsymbol{x}}^{*}$. From Eq. 22, it holds that
$K{\boldsymbol{x}}^{*}={\boldsymbol{y}}$. Moreover, by the boundedness of
$\\{{\boldsymbol{x}}_{k^{\prime}}\\}$ and the uniform convergence granted by
Lemma A.8, we infer that

$\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}^{*}\|_{1}=\lim_{k^{\prime}}\|\Delta_{\Psi^{\delta_{k^{\prime}}}_{\Theta_{k^{\prime}}}}{\boldsymbol{x}}_{k^{\prime}}\|_{1}.$

Applying Eq. 23, we finally get

$\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}^{*}\|_{1}=\lim_{k^{\prime}}\|\Delta_{\Psi^{\delta_{k^{\prime}}}_{\Theta_{k^{\prime}}}}{\boldsymbol{x}}_{k^{\prime}}\|_{1}\leq\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}_{\textnormal{sol}}\|_{1}\leq\|\Delta_{{\boldsymbol{x}}_{0}}{\boldsymbol{x}}^{*}\|_{1}.$

That is, ${\boldsymbol{x}}^{*}$ is a graph-minimizing solution. If the graph-
minimization solution is unique, then we have just proven that every
subsequence of $\\{{\boldsymbol{x}}_{k}\\}$ has a subsequence converging to
${\boldsymbol{x}}_{\textnormal{sol}}$, and therefore
${\boldsymbol{x}}_{k}\to{\boldsymbol{x}}_{\textnormal{sol}}$ by a standard
topological argument.

We proceed now to prove the stability result in Theorem 3.15. However, as
first step we demonstrate that the standard Tikhonov reconstruction method
coupled with the discrepancy principle satisfies Section 3.

######  Example A.13.

Consider a standard Tikhonov reconstruction method, that is

$\Psi_{\Theta}({\boldsymbol{y}})\coloneqq\left(K^{T}K+\Theta
I\right)^{-1}K^{T}{\boldsymbol{y}},$

where $\Theta\in(0,\infty)$. Let
$\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})$ defined by the discrepancy
principle as per [21, Equation (4.57)]. Then it is easy to check that
$\Theta_{k}\to\Theta=\Theta(\delta,{\boldsymbol{y}}^{\delta})\ \mathrm{for}\
k\rightarrow\infty$. Setting

$\mathcal{T}^{k}=(K^{T}K+\Theta_{k}I)^{-1}\quad\mbox{and}\quad\mathcal{T}=(K^{T}K+\Theta
I)^{-1},$

it holds that

$\displaystyle\|\Psi_{\Theta_{k}}({\boldsymbol{y}}^{\delta_{k}})-\Psi_{\Theta}({\boldsymbol{y}}^{\delta})\|_{2}\leq\underbrace{\|(\mathcal{T}^{k}-\mathcal{T})K^{T}{\boldsymbol{y}}^{\delta_{k}}\|_{2}}_{\mathbf{I}}+\underbrace{\|\mathcal{T}K^{T}({\boldsymbol{y}}^{\delta}-{\boldsymbol{y}}^{\delta_{k}})\|_{2}}_{\mathbf{II}}.$

Now,

$\displaystyle\mathbf{I}\leq|\Theta-\Theta_{k}|\|\mathcal{T}^{k}\|\|\mathcal{T}K^{T}\|\|{\boldsymbol{y}}^{\delta_{k}}\|_{2}\leq\frac{|\Theta-\Theta_{k}|}{2\Theta_{k}\sqrt{\Theta}}\|{\boldsymbol{y}}^{\delta_{k}}\|_{2}=f_{1}(k)\to
0,$

and

$\displaystyle\mathbf{II}\leq\|(K^{T}K+\Theta
I)^{-1}K^{T}\|\|{\boldsymbol{y}}^{\delta}-{\boldsymbol{y}}^{\delta_{k}}\|_{2}\leq\frac{1}{2\sqrt{\Theta}}\|{\boldsymbol{y}}^{\delta}-{\boldsymbol{y}}^{\delta_{k}}\|=f_{2}(k)\to
0.$

Now we need a couple of preliminary lemmas.

######  Lemma A.14.

For all ${\boldsymbol{x}}\in X$ and
${\boldsymbol{y}}^{\delta_{k_{1}}},{\boldsymbol{y}}^{\delta_{k_{2}}}\in Y$, we
have

$\Gamma_{k_{1}}({\boldsymbol{x}})\leq
2\Gamma_{k_{2}}({\boldsymbol{x}})+\|{\boldsymbol{y}}^{\delta_{k_{1}}}-{\boldsymbol{y}}^{\delta_{k_{2}}}\|+\|(\Delta_{\Psi_{\Theta_{k_{1}}}^{\delta_{k_{1}}}}-\Delta_{\Psi_{\Theta_{k_{2}}}^{\delta_{k_{2}}}}){\boldsymbol{x}}\|_{1}$

######  Proof A.15.

By standard $p$-norm inequalities

$\displaystyle\Gamma_{k_{1}}({\boldsymbol{x}})$
$\displaystyle=\frac{1}{2}\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta_{k_{1}}}\|_{2}^{2}+\alpha\|\Delta_{\Psi_{\Theta_{k_{1}}}^{\delta_{k_{1}}}}{\boldsymbol{x}}\|_{1}$
$\displaystyle\leq\|K{\boldsymbol{x}}-{\boldsymbol{y}}^{\delta_{k_{2}}}\|_{2}^{2}+\|{\boldsymbol{y}}^{\delta_{k_{1}}}-{\boldsymbol{y}}^{\delta_{k_{2}}}\|_{2}^{2}+\alpha\|\Delta_{\Psi_{\Theta_{k_{1}}}^{\delta_{k_{1}}}}{\boldsymbol{x}}\|_{1}$
$\displaystyle\leq
2\Gamma_{k_{2}}({\boldsymbol{x}})+\|{\boldsymbol{y}}^{\delta_{k_{1}}}-{\boldsymbol{y}}^{\delta_{k_{2}}}\|+\|(\Delta_{\Psi_{\Theta_{k_{1}}}^{\delta_{k_{1}}}}-\Delta_{\Psi_{\Theta_{k_{2}}}^{\delta_{k_{2}}}}){\boldsymbol{x}}\|_{1}$

######  Lemma A.16.

Let ${\delta_{k}}$, $\Theta_{k}$, $\Theta$ and $f$ defined as in Section 3. It
holds that

$\|\Delta_{\Psi_{\Theta_{k}}^{\delta_{k}}}{\boldsymbol{x}}-\Delta_{\Psi_{\Theta}^{\delta}}{\boldsymbol{x}}\|_{1}\leq
c\|{\boldsymbol{x}}\|_{\infty}f(k)\to 0\quad\mbox{as }k\to\infty,$

where $c$ is a positive constant independent of ${\boldsymbol{x}}$.

######  Proof A.17.

The proof is similar to Lemma A.8 using hypothesis 3.

Theorem 3.15.

######  Proof A.18.

Because ${\boldsymbol{x}}_{k}$ is a minimizer of $\Gamma_{k}$, we have

(24)
$\Gamma_{k}({\boldsymbol{x}}_{k})\leq\Gamma_{k}({\boldsymbol{x}}),\quad\forall\,{\boldsymbol{x}}\in
X.$

Choose now a vector $\bar{{\boldsymbol{x}}}\in X$. By applying the previous
equation to ${\boldsymbol{x}}=\bar{{\boldsymbol{x}}}$ and using twice Lemma
A.14, it follows that

$\displaystyle\Gamma({\boldsymbol{x}}_{k})$ $\displaystyle\leq
2\Gamma_{k}({\boldsymbol{x}}_{k})+\|{\boldsymbol{y}}^{\delta}-{\boldsymbol{y}}^{\delta_{k}}\|+\|(\Delta_{\Psi_{\Theta}^{\delta}}-\Delta_{\Psi_{\Theta_{k}}^{\delta_{k}}}){\boldsymbol{x}}_{k}\|_{1}$
$\displaystyle\leq
2\Gamma_{k}(\bar{{\boldsymbol{x}}})+\|{\boldsymbol{y}}^{\delta}-{\boldsymbol{y}}^{\delta_{k}}\|+\|(\Delta_{\Psi_{\Theta}^{\delta}}-\Delta_{\Psi_{\Theta_{k}}^{\delta_{k}}}){\boldsymbol{x}}_{k}\|_{1}$
$\displaystyle\leq
4\Gamma(\bar{{\boldsymbol{x}}})+2\|{\boldsymbol{y}}^{\delta}-{\boldsymbol{y}}^{\delta_{k}}\|+2\|(\Delta_{\Psi_{\Theta}^{\delta}}-\Delta_{\Psi_{\Theta_{k}}^{\delta_{k}}}){\boldsymbol{x}}_{k}\|_{1}.$

From Lemma A.16, since ${\boldsymbol{y}}^{\delta_{k}}$ converges to
${\boldsymbol{y}}^{\delta}$ with respect to the norm topology, then there
exist $k_{0}\in\mathbb{N}$ such that

$M:=2\Gamma(\bar{{\boldsymbol{x}}})+1\geq\Gamma({\boldsymbol{x}}_{k}),\quad\forall\,k\geq
k_{0}.$

The thesis now follows easily by adapting the proof in [47, Theorem 3.4].

## References

  * [1] S. Aleotti, D. Bianchi, D. Evangelista, M. Donatelli, and E. Loli Piccolomini, Official GitHub repository for the $\textnormal{{graphla+}}\psi$ codes, 2023, https://github.com/devangelista2/GraphLaPlus.  Online; accessed 22-December-2023.
  * [2] A. Ali and R. J. Tibshirani, The Generalized Lasso Problem and Uniqueness, Electronic Journal of Statistics, 13 (2019), pp. 2307–2347.
  * [3] V. Antun, F. Renna, C. Poon, B. Adcock, and A. C. Hansen, On instabilities of deep learning in image reconstruction and the potential costs of AI, in Proceedings of the National Academy of Sciences, vol. 117, 2020, pp. 30088–30095.
  * [4] P. Arias, V. Caselles, and G. Sapiro, A variational framework for non-local image inpainting, in International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition, Springer, 2009, pp. 345–358.
  * [5] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb, Solving inverse problems using data-driven models, Acta Numerica, 28 (2019), pp. 1–174.
  * [6] H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, CMS books in mathematics, Springer Cham, 2 ed., 2017.
  * [7] D. Bianchi, A. Buccini, M. Donatelli, and E. Randazzo, Graph Laplacian for image deblurring, ETNA, 55 (2022), pp. 169–186.
  * [8] D. Bianchi, A. Buccini, M. Donatelli, and S. Serra-Capizzano, Iterated fractional Tikhonov regularization, Inverse Problems, 31 (2015), p. 055005.
  * [9] D. Bianchi and M. Donatelli, On generalized iterated Tikhonov regularization with operator-dependent seminorms, Electronic Transactions on Numerical Analysis, 47 (2017), pp. 73–99.
  * [10] D. Bianchi and M. Donatelli, Graph approximation and generalized Tikhonov regularization for signal deblurring, in 2021 21st International Conference on Computational Science and Its Applications (ICCSA), IEEE, 2021, pp. 93–100.
  * [11] D. Bianchi, M. Donatelli, D. Evangelista, W. Li, and E. L. Piccolomini, Graph Laplacian and Neural Networks for Inverse Problems in Imaging: GraphLaNet, in International Conference on Scale Space and Variational Methods in Computer Vision, 2023, pp. 175–186.
  * [12] D. Bianchi, G. Lai, and W. Li, Uniformly convex neural networks and non-stationary iterated network Tikhonov (iNETT) method, Inverse Problems, 39 (2023), p. 055002.
  * [13] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, Geometric Deep Learning: Going beyond Euclidean data, IEEE Signal Processing Magazine, 34 (2017), pp. 18–42, https://doi.org/10.1109/MSP.2017.2693418.
  * [14] A. Buccini and M. Donatelli, Graph Laplacian in $\ell^{2}-\ell^{q}$ regularization for image reconstruction, in 2021 21st International Conference on Computational Science and Its Applications (ICCSA), IEEE, 2021, pp. 29–38.
  * [15] J.-F. Cai, R. H. Chan, and Z. Shen, A framelet-based image inpainting algorithm, Applied and Computational Harmonic Analysis, 24 (2008), pp. 131–149.
  * [16] Y. Cai, M. Donatelli, D. Bianchi, and T.-Z. Huang, Regularization preconditioners for frame-based image deblurring with reduced boundary artifacts, SIAM Journal on Scientific Computing, 38 (2016), pp. B164–B189.
  * [17] L. Calatroni, Y. van Gennip, C.-B. Schönlieb, H. M. Rowland, and A. Flenner, Graph clustering, variational image segmentation methods and Hough transform scale detection for object measurement in images, Journal of Mathematical Imaging and Vision, 57 (2017), pp. 269–291.
  * [18] E. Candes and B. Recht, Simple bounds for recovering low-complexity models, Mathematical Programming, 141 (2013), pp. 577–589.
  * [19] J. Chung and S. Gazzola, Flexible krylov methods for $\ell_{p}$ regularization, SIAM Journal on Scientific Computing, 41 (2019), pp. S149–S171.
  * [20] M. J. Colbrook, V. Antun, and A. C. Hansen, The difficulty of computing stable and accurate neural networks: On the barriers of deep learning and Smale’s 18th problem, in Proceedings of the National Academy of Sciences, vol. 119, 2022, p. e2107151119.
  * [21] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of inverse problems, Springer Science & Business Media, 1996.
  * [22] D. Evangelista, E. Morotti, and E. Loli Piccolomini, COULE dataset.  https://www.kaggle.com/datasets/loiboresearchgroup/coule-dataset, 2023.  Accessed on 01/12/2023.
  * [23] D. Evangelista, E. Morotti, and E. L. Piccolomini, Rising: A new framework for model-based few-view ct image reconstruction with deep learning, Computerized Medical Imaging and Graphics, 103 (2023), p. 102156.
  * [24] S. Gazzola, P. C. Hansen, and J. G. Nagy, Ir tools: a matlab package of iterative regularization methods and large-scale test problems, Numerical Algorithms, 81 (2019), pp. 773–811.
  * [25] G. Gilboa and S. Osher, Nonlocal linear image regularization and supervised segmentation, Multiscale Modeling & Simulation, 6 (2007), pp. 595–630.
  * [26] G. Gilboa and S. Osher, Nonlocal operators with applications to image processing, Multiscale Modeling & Simulation, 7 (2009), pp. 1005–1028.
  * [27] Y. S. Han, J. Yoo, and J. C. Ye, Deep residual learning for compressed sensing CT reconstruction via persistent homology analysis, arXiv preprint arXiv:1611.06391, (2016).
  * [28] P. C. Hansen, Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion, SIAM, 1998.
  * [29] P. C. Hansen, J. G. Nagy, and D. P. O’leary, Deblurring images: matrices, spectra, and filtering, SIAM, 2006.
  * [30] C. F. Higham and D. J. Higham, Deep learning: An introduction for applied mathematicians, Siam review, 61 (2019), pp. 860–891.
  * [31] M. E. Hochstenbach and L. Reichel, Fractional Tikhonov regularization for linear discrete ill-posed problems, BIT Numerical Mathematics, 51 (2011), pp. 197–215.
  * [32] J. Huang, M. Donatelli, and R. H. Chan, Nonstationary iterated thresholding algorithms for image deblurring, Inverse Probl. Imaging, 7 (2013), pp. 717–736.
  * [33] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser, Deep convolutional neural network for inverse problems in imaging, IEEE transactions on image processing, 26 (2017), pp. 4509–4522.
  * [34] A. C. Kak and M. Slaney, Principles of computerized tomographic imaging, SIAM, 2001.
  * [35] M. Keller, D. Lenz, and R. K. Wojciechowski, Graphs and Discrete Dirichlet Spaces, Grundlehren der mathematischen Wissenschaften, Springer, Cham, 2021.
  * [36] E. Klann and R. Ramlau, Regularization by fractional filter methods and data smoothing, Inverse Problems, 24 (2008), p. 025018.
  * [37] A. Lanza, S. Morigi, L. Reichel, and F. Sgallari, A generalized krylov subspace method for $\ell_{p}-\ell_{q}$ minimization, SIAM Journal on Scientific Computing, 37 (2015), pp. S30–S50.
  * [38] H. Li, J. Schwab, S. Antholzer, and M. Haltmeier, NETT: Solving inverse problems with deep neural networks, Inverse Problems, 36 (2020), p. 065005.
  * [39] Y. Lou, X. Zhang, S. Osher, and A. Bertozzi, Image recovery via nonlocal operators, Journal of Scientific Computing, 42 (2010), pp. 185–197.
  * [40] T. R. Moen, B. Chen, D. R. Holmes III, X. Duan, Z. Yu, L. Yu, S. Leng, J. G. Fletcher, and C. H. McCollough, Low-dose CT image and projection dataset, Medical physics, 48 (2021), pp. 902–911.
  * [41] E. Morotti, D. Evangelista, and E. Loli Piccolomini, A green prospective for learned post-processing in sparse-view tomographic reconstruction, Journal of Imaging, 7 (2021), p. 139.
  * [42] E. Morotti, D. Evangelista, and E. L. Piccolomini, Increasing noise robustness of deep learning-based image processing with model-based approaches, Numerical Computations: Theory and Algorithms NUMTA 2023, (2023), p. 155.
  * [43] V. A. Morozov, Methods for solving incorrectly posed problems, Springer New York, NY, 1984.
  * [44] S. Mukherjee, S. Dittmer, Z. Shumaylov, S. Lunz, O. Öktem, and C.-B. Schönlieb, Learned Convex Regularizers for Inverse Problems, 2020, https://arxiv.org/abs/arXiv:2008.02839.
  * [45] G. Peyré, S. Bougleux, and L. Cohen, Non-local regularization of inverse problems, in Computer Vision – ECCV 2008, D. Forsyth, P. Torr, and A. Zisserman, eds., Springer Berlin Heidelberg, 2008, pp. 57–68.
  * [46] O. Ronneberger, P. Fischer, and T. Brox, U-net: Convolutional networks for biomedical image segmentation, in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, Springer, 2015, pp. 234–241.
  * [47] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen, Variational methods in imaging, Springer, 2009.
  * [48] J. Schwab, S. Antholzer, and M. Haltmeier, Deep null space learning for inverse problems: convergence analysis and rates, Inverse Problems, 35 (2019), p. 025008.
  * [49] E. Tjoa and C. Guan, A survey on explainable artificial intelligence (xai): Toward medical xai, IEEE transactions on neural networks and learning systems, 32 (2020), pp. 4793–4813.
  * [50] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE transactions on image processing, 13 (2004), pp. 600–612.
  * [51] X. Zhang, M. Burger, X. Bresson, and S. Osher, Bregmanized nonlocal regularization for deconvolution and sparse reconstruction, SIAM Journal on Imaging Sciences, 3 (2010), pp. 253–276.