# Image reconstructions using sparse dictionary representations and implicit,
non-negative mappings

Elizabeth Newman  Department of Mathematics, Emory University, Atlanta, GA,
USA $\big{\\{}$elizabeth.newman, jack.michael.solomon,
<EMAIL_ADDRESS>Jack Michael Solomon  Department of
Mathematics, Emory University, Atlanta, GA, USA $\big{\\{}$elizabeth.newman,
jack.michael.solomon<EMAIL_ADDRESS>Matthias Chung
Department of Mathematics, Emory University, Atlanta, GA, USA
$\big{\\{}$elizabeth.newman, jack.michael.solomon,
<EMAIL_ADDRESS>
###### Abstract

Many imaging science tasks can be modeled as a discrete linear inverse
problem. Solving linear inverse problems is often challenging, with ill-
conditioned operators and potentially non-unique solutions. Embedding prior
knowledge, such as smoothness, into the solution can overcome these
challenges. In this work, we encode prior knowledge using a non-negative patch
dictionary, which effectively learns a basis from a training set of natural
images. In this dictionary basis, we desire solutions that are non-negative
and sparse (i.e., contain many zero entries). With these constraints, standard
methods for solving discrete linear inverse problems are not directly
applicable. One such approach is the modified residual norm steepest descent
(MRNSD), which produces non-negative solutions but does not induce sparsity.
In this paper, we provide two methods based on MRNSD that promote sparsity. In
our first method, we add an $\ell_{1}$-regularization term with a new, optimal
step size. In our second method, we propose a new non-negative, sparsity-
promoting mapping of the solution. We compare the performance of our proposed
methods on a number of numerical experiments, including deblurring, image
completion, computer tomography, and superresolution. Our results show that
these methods effectively solve discrete linear inverse problems with non-
negativity and sparsity constraints.

##  1 Introduction

_Inverse problems_ are an important class of mathematical problems with the
aim to re- or uncover unobserved information from noisy data paired with an
underlying mathematical or statistical model. Inverse problems are ubiquitous
in modern science, from data analysis and machine learning to many imaging
applications, such as image deblurring, computer tomography, and super-
resolution [17].

In this work, we consider the common imaging problem that is given as a
discrete linear inverse problem

$\mathbf{b}=\mathbf{C}\mathbf{y}_{\rm true}+\text{\boldmath$\mathbf{\eta}$},$
(1)

where $\mathbf{y}_{\rm true}\in\mathcal{Y}\subset\mathbb{R}^{n}_{+}$ is the
desired solution (a vectorized representation of the image is some image
space), $\mathbf{C}\in\mathbb{R}^{m\times n}$ is a given forward imaging
process and $\mathbf{b}\in\mathbb{R}^{m}$ contains observed measurements. We
assume these observations $\mathbf{b}$ are corrupted by additive Gaussian
noise $\text{\boldmath$\mathbf{\eta}$}\sim\mathcal{N}({\bf
0},\sigma^{2}\mathbf{I}_{m})$. In most cases, the underlying operator
$\mathbf{C}$ is ill-conditioned (i.e., the solution is sensitive to small
perturbations in $\mathbf{b}$) or underdetermined ($m<n$) such that the
solution is not unique. In most of these common settings, prior knowledge is
required for selecting a meaningful solution. One approach to embedding prior
knowledge into the solution is to cast (1) as a _variational inverse problem_
of the form [10, 16]

$\min_{\mathbf{y}\in\mathcal{Y}}\quad\tfrac{1}{2}\left\|\mathbf{C}\mathbf{y}-\mathbf{b}\right\|_{2}^{2}+\mu\mathcal{R}(\mathbf{y}).$
(2)

The first term measures the data fidelity between the model prediction,
$\mathbf{C}\mathbf{y}$, and observations, $\mathbf{b}$, with discrepancy
measured in $\ell^{2}$-norm. For $\mu=0$, this corresponds to a maximum
likelihood estimator with Gaussian noise distribution. The second term,
$\mathcal{R}(\,\cdot\,)$, denotes an appropriate regularizer that promotes
desirably properties in the solution, $\mathbf{y}$, (e.g., small norm,
smoothness). The regularization parameter $\mu>0$ balances the weight of the
data fidelity and regularization terms, and is often chosen heuristically or
based on noise estimates [15].

Regularization for imaging problems is motivated by the fact that the _natural
image space_ $\mathcal{Y}\subset\mathbb{R}^{n}_{+}$ is minuscule compared to
the set of all images of dimension $\mathbb{R}^{n}_{+}$. For instance, drawing
a random sample from $\mathbb{R}^{n}_{+}$ will most likely look like noise
[4]. In other words, natural images exhibit a high degree of redundancy and
can be efficiently captured within a low-dimensional manifold. Hence, the
features and the structure of natural images must be included in an inverse
problem process.One common approach is to utilize an edge-preserving _total
variation_ (TV) regularization
$\mathcal{R}(\mathbf{y})=\left\|\mathbf{L}\mathbf{y}\right\|_{p}$. Here,
$\mathbf{L}\mathbf{y}$ is a finite difference approximation of the directional
pixel gradients, and $\left\|\mathbf{y}\right\|_{p}$ denotes the
$\ell^{p}$-norm [30]. This type of regularization has shown tremendous success
in various applications. However, disadvantages of this approach include the
potential over-smoothing of the reconstructed image (apart from the edges) and
potential inconsistent intensity shifts in $\mathbf{y}$, [8, 38, 13].

An entirely different approach to regularize is to represent the solution to
(2) in terms of a pre-trained image data “basis” or _dictionary_ [36].
Dictionaries of natural images or image elements provide an informative prior
and incorporate expected image features (e.g., edges, smoothness) into
$\mathbf{y}$ by directly. To be precise, we represent $\mathbf{y}$ as a linear
combination of vectors (atoms) $\mathbf{g}_{k}\in\mathbb{R}^{n}_{+}$ for
$k=1,\ldots,K$, i.e., $\mathbf{y}=\mathbf{G}\mathbf{x}$, where
$\mathbf{G}=[\mathbf{g}_{1},\ldots,\mathbf{g}_{K}]$ is referred to as the
_dictionary_ and $\mathbf{x}\in\mathbb{R}^{K}_{+}$ are appropriate non-
negative selection coefficients. If $K>n$ or even $K\gg n$, we say the
dictionary is _overcomplete_. As a result, the coefficient vector is not
unique and larger than the original image. In this setting, we desire
coefficients that exhibit sparsity, i.e., many elements in $\mathbf{x}$ are
being zero, to ensure we do not increase the storage of our solution.

To select such a suitable and sparse coefficient vector $\mathbf{x}$ (and
therefore a suitable $\mathbf{y}$), we ideally minimize

$\min_{\mathbf{x}\geq{\bf
0}}\quad\left\|\mathbf{x}\right\|_{0}\quad\mbox{subject
to}\quad\tfrac{1}{2}\left\|\mathbf{C}\mathbf{G}\mathbf{x}-\mathbf{b}\right\|_{2}^{2}\leq\delta^{2},$
(3)

where $\left\|\mathbf{x}\right\|_{0}$ is not a proper norm and is defined as
the cardinality of non-zero elements in $\mathbf{x}$, while $\delta$
represents the given noise level. Solving (3) is NP-hard and efficient
approximation approaches need to be utilized, [13, 40]. One approach in
tackling (3) is to incrementally generate non-zero elements in $\mathbf{x}$ by
iteratively selecting and updating the element in $\mathbf{x}$ minimizing the
remaining residual until a prescribed accuracy threshold is reached. Methods
following this approach include _matching pursuit_ and _orthogonal matching
pursuit_ , [25, 41]. A main disadvantage of matching pursuit approaches is
that we cannot ensure optimally. Further, each iteration requires the
computation of inner products with the remaining dictionary adding to the
computational complexity of the algorithm.

A common alternative approach to provide a non-negative and sparse
$\mathbf{x}$ is based on convex relaxation. We approximate (3) by

$\min_{\mathbf{x}\geq{\bf
0}}\quad\tfrac{1}{2}\left\|\mathbf{A}\mathbf{x}-\mathbf{b}\right\|_{2}^{2}+\lambda\left\|\mathbf{x}\right\|_{1},$
(4)

where $\mathbf{A}=\mathbf{C}\mathbf{G}$ and $\lambda>0$ is an appropriate
regularization parameter [6]. When $\mathbf{x}$ is unrestricted, the
optimization problem (4) is the so-called least absolute shrinkage and
selection operator (LASSO) problem. Methods such as _basis pursuit_ ,
_Alternating Direction Method of Multipliers_ (ADMM), and _Variable Projected
Augmented Lagrangian_ (VPAL) have been shown to efficiently solve an
unconstrained version of (4); see [9, 3, 11]. Additionally, Krylov subspace
methods are a major tool for large-scale least squares problems [34]. However,
none of these approaches are directly applicable in this setting because they
do not inherently enforce non-negativity or promote sparse solutions.
Consequently, alternative approaches have been developed to address these
constraints. Here, we follow another approach, we propose to use a method
designed for non-negative least squares and extend these approaches to also
provide sparse solutions.

Our work is structured as follows. In Section 2 we provide details on the
modified residual norm steepest descent methods which have been developed to
solve least squares problems with non-negativity constraints and we further
provide a background on patch dictionary learning. In Section 3 we discuss our
proposed methods to solve non-negative least squares problems with sparsity
priors. We demonstrate the effectiveness of our approaches with various
numerical experiments in Section 4 including deblurring (Section 4.2), image
completion (Section 4.3), and computer tomography (Section 4.4). We close up
our investigations with some conclusions and further outlook in Section 5.

##  2 Background

Because images are inherently non-negative, we are interested in the
constrained optimization problem

$\min_{\mathbf{x}\geq{\bf 0}}\quad
f(\mathbf{x})\equiv\tfrac{1}{2}\left\|\mathbf{A}\mathbf{x}-\mathbf{b}\right\|_{2}^{2}.$
(5)

While standard constrained optimization techniques may be employed for (5),
they often suffer from poor computational efficiency [29]. To address this
issue, specialized methods like the expectation-maximization (EM) algorithm or
modified residual norm steepest descent (_MRNSD_) have been proposed. These
methods have shown great success by directly handling the linear least squares
fitting term with non-negativity constraints [21, 26].

###  2.1 Modified Residual Norm Steepest Descent (MRNSD)

The main idea behind MRNSD is to enforce non-negativity by introducing an
implicit, element-wise exponential variable transformation and performing a
gradient descent update directly on $\mathbf{x}$ with appropriate step size
selection. Let us introduce the element-wise variable transformation
$x_{i}=\textnormal{e}^{z_{i}}$, for $i=1,\ldots,K$, short we write
$\mathbf{x}=\textnormal{e}^{\mathbf{z}}$. This mapping ensure
$\mathbf{x}\geq{\bf 0}$ and the optimization problem (5) reduces to

$\min_{\mathbf{z}}\quad\hat{f}(\mathbf{z})\equiv\tfrac{1}{2}\left\|\mathbf{A}\textnormal{e}^{\mathbf{z}}-\mathbf{b}\right\|_{2}^{2}.$
(6)

Note the gradient of $\hat{f}$, the objective function over
$\mathbf{z}$-space, is given as

$\nabla_{\mathbf{z}}\hat{f}(\mathbf{z})={\rm
diag\\!}\left(\textnormal{e}^{\mathbf{z}}\right)\mathbf{A}^{\top}(\mathbf{A}\textnormal{e}^{\mathbf{z}}-\mathbf{b})={\rm
diag\\!}\left(\mathbf{x}\right)\nabla_{\mathbf{x}}f(\mathbf{x}).$ (7)

MRNSD performs a steepest descent direction update with respect $\mathbf{x}$
using the gradient of $\widehat{f}$, i.e.,

$\displaystyle\mathbf{s}=-\nabla_{\mathbf{z}}\hat{f}(\mathbf{z})=-{\rm
diag\\!}\left(\mathbf{x}\right)\mathbf{g},$ (8)

where
$\mathbf{g}=\nabla_{\mathbf{x}}f(\mathbf{x})=\mathbf{A}^{\top}(\mathbf{A}\mathbf{x}-\mathbf{b})$
is the gradient of the $f$, the objective function in $\mathbf{x}$-space. We
note that $\mathbf{s}$ is a descent direction in $\mathbf{x}$-space, i.e.,
$\mathbf{s}^{\top}\mathbf{g}<0$, because of the assumed positivity of
$\mathbf{x}$. This ensures ${\rm diag\\!}\left(\mathbf{x}\right)$ is a
positive definite matrix, which in turn ensures we move in a descent
direction. Further note that an optimal step size can be computed, i.e.,

$\displaystyle\alpha=\frac{-\mathbf{s}^{\top}\mathbf{g}}{\mathbf{s}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{s}}.$
(9)

However, since MRNSD updates $\mathbf{x}$ directly (and not the surrogate
variable $\mathbf{z}$), there may exist a selected step length $\alpha>0$ for
which the update $\mathbf{x}+\alpha\mathbf{s}$ can have negative elements.
Hence, thresholding must be performed to maintain non-negativity in
$\mathbf{x}$ which is provided by

$\displaystyle\alpha=\min\left\\{-\mathbf{s}^{\top}\mathbf{g}/(\mathbf{s}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{s}),\min_{s_{i}<0}\\{-x_{i}/s_{i})\\}_{i=1}^{K}\right\\}.$
(10)

The MRNSD method is summarized in Algorithm 1.

Algorithm 1 Modified Residual Norm Steepest Descent (MRNSD) [26, 21]

1:Solve $\min_{\mathbf{x}\geq
0}\frac{1}{2}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}^{2}$

2:while not converged do

3: $\mathbf{g}=\mathbf{A}^{\top}(\mathbf{A}\mathbf{x}-\mathbf{b})$

4: $\mathbf{s}=-{\rm diag\\!}\left(\mathbf{x}\right)\mathbf{g}$
$\triangleright$ modified steepest descent direction

5:
$\alpha=\min\left\\{-\mathbf{s}^{\top}\mathbf{g}/(\mathbf{s}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{s}),\min_{s_{i}<0}\\{-x_{i}/s_{i})\\}_{i=1}^{K}\right\\}$
$\triangleright$ optimal step size with bounded line search

6: $\mathbf{x}\leftarrow\mathbf{x}+\alpha\mathbf{s}$

7:end while

Note that MRNSD only performs a surrogate variable transformation implicitly
but never constructs nor updates $\mathbf{z}$. While MRNSD has been shown to
effectively address non-negativity constraints, the utilized exponential
mapping prevents values of $\mathbf{x}$ from becoming zero and hence does not
naturally promote sparsity in the solution. Achieving sparsity requires ad hoc
threshholding in postprocessing, and motivates our exploration of efficient
methods that enforce both non-negativity and sparsity during optimization.

###  2.2 Dictionary Learning

Before we turn our attention to sparsity-promoting alternatives to MRNSD, we
next provide some details on patch dictionary learning and the resulting
sparse representations [36, 37]. Suppose we have a set of $t$ training
$M\times N$ images $\\{\mathbf{Y}_{1},\cdots,\mathbf{Y}_{t}\\}$ with ${\rm
vec}(\mathbf{Y}_{j})\in\mathcal{Y}\subset\mathbb{R}_{+}^{MN}$, where
$\operatorname{vec}(\cdot)$ vectorizes the image column-wise. We collect
(overlapping) patches of size $p\times q$ with $p<M$ and $q<N$, “patchify” the
images in our training set, and vectorize those patches and store them as
columns of a matrix $\mathbf{Y}^{\rm patch}\in\mathbb{R}_{+}^{pq\times tr}$
where $r$ is the number of patches per image; see Figure 1 for an
illustration.

$\mathbf{Y}$
---
$M\times N$
$\texttt{patchify}(\mathbf{Y})$
---
$pq\times r$
Figure 1: Illustration of image patchification. For simplicity, we illustrate
non-overlapping patches that exactly partition the image with $r=(M/p)(N/q)$
total patches. In Matlab, this operation is given by
im2col(Y,[p,q],’distinct’) and is reversible.

A goal in dictionary learning is to construct an appropriate non-negative
patch dictionary $\mathbf{D}\in\mathbb{R}_{+}^{pq\times s}$ with a
corresponding set of non-negative, column-wise sparse coefficients
$\mathbf{X}^{\rm patch}\in\mathbb{R}_{+}^{s\times r}$ such that
$\mathbf{Y}^{\rm patch}\approx\mathbf{D}\mathbf{X}^{\rm patch}$. The desired
properties of the decomposition ensure that each patch of our training set
will be constructed using only a few dictionary atoms, producing a “parts-of-
the-whole” type of basis [23].

Following the presentation in [36], we form this non-negative matrix
factorization by solving the optimization problem

$\displaystyle\min_{\mathbf{D}\in\mathcal{D},\mathbf{X}^{\rm
patch}\in\mathbb{R}_{+}^{s\times r}}\tfrac{1}{2}\left\|\mathbf{Y}^{\rm
patch}-\mathbf{D}\mathbf{X}^{\rm patch}\right\|_{{\rm
F}}^{2}+\beta\sum_{j=1}^{r}\|\mathbf{x}_{j}^{\rm patch}\|_{1},$ (11)

where $\mathbf{x}_{j}^{\rm patch}\in\mathbb{R}^{s}$ is the $j$-th column of
$\mathbf{X}^{\rm patch}$, $\beta>0$ is the regularization parameter for the
sparsity-promoting $\ell_{1}$-regularization, and $\mathcal{D}$ restricts the
values of the dictionary to lie between $0$ and $1$. In our experiments, we
use box constraints

$\displaystyle\mathcal{D}:=\left\\{\mathbf{D}\in\mathbb{R}^{pq\times s}\mid
0\leq d_{ij}\leq 1\text{ for }i=1,\dots,pq\text{ and }j=1,\dots,s\right\\}.$
(12)

The dictionary learning problem (11) is non-convex with respect to
$(\mathbf{D},\mathbf{X}^{\rm patch})$ jointly, but is convex with respect to
each variable separately. This motivates the use of an alternating
optimization strategy. To this end, we consider the alternating direction
method of multipliers (ADMM) approach proposed by [36], which reformulates
(11) as

$\displaystyle\min_{\mathbf{D},\mathbf{X}^{\rm patch},\mathbf{U},\mathbf{V}}$
$\displaystyle\quad\tfrac{1}{2}\left\|\mathbf{Y}^{\rm
patch}-\mathbf{U}\mathbf{V}\right\|_{{\rm
F}}^{2}+\beta\sum_{i=1}^{t}\left\|\mathbf{x}_{j}^{\rm
patch}\right\|_{1}+\chi_{\mathcal{D}}(\mathbf{D})+\chi_{\mathbb{R}_{+}^{s\times
r}}(\mathbf{X}^{\rm patch}),$ (13a) subject to
$\displaystyle\quad\mathbf{D}=\mathbf{U}\quad\text{ and }\quad\mathbf{X}^{\rm
patch}=\mathbf{V},$ (13b)

where $\chi_{\mathcal{Z}}$ is the characteristic function over the set
$\mathcal{Z}$ such that $\chi_{\mathcal{Z}}(z)=0$ if $z\in\mathcal{Z}$ and
$\chi_{\mathcal{Z}}(z)=+\infty$ if $z\not\in\mathcal{Z}$. Introducing the
auxiliary $\mathbf{U}$ and $\mathbf{V}$ ensures that the function is separable
over the optimization variables, a necessary condition for ADMM. We solve (13)
using the augmented Lagrangian approach [3]; for further details, see [36, 37,
28]. Examples of generated dictionary patches are shown in Figure 2.

The patch dictionary presentation is related to the original global dictionary
formulation in Section 1. Let $\mathbf{X}_{i}\in\mathbb{R}_{+}^{pq\times r}$
be the coefficients of the patches drawn from training image $\mathbf{Y}_{i}$
and assume the patches are non-overlapping. Then,

$\displaystyle\operatorname{vec}(\mathbf{Y}_{i})=\mathbf{P}\operatorname{vec}(\mathbf{D}\mathbf{X}_{i})=\underbrace{\mathbf{P}(\mathbf{I}\otimes\mathbf{D})}_{\mathbf{G}}\operatorname{vec}(\mathbf{X}_{i}),$
(14)

where $\otimes$ denotes the Kronecker product, $\mathbf{I}$ is the $r\times r$
identity matrix, and $\mathbf{P}$ is a permutation matrix mapping patches to
linear indices [32]. We require the permutation matrix because
$\operatorname{vec}(\mathbf{X}_{i})$ stacks vectorized patches column-wise,
but $\operatorname{vec}(\mathbf{Y}_{i})$ concatenates the columns of the image
vertically, which order the pixels differently.

(a)

(b)

Figure 2: Sample of $10$ patches from dictionaries generated from two
different sources. For training, we collect overlapping patches of size
$16\times 16$. The patches share similar binary-looking patterns because of
the box constraints (12).

##  3 Non-negativity and Sparsity-Promoting Mappings

We consider two alternative approaches based on MRNSD to promote sparse
solutions. In the first approach, we introduce an $\ell_{1}$-regularization
term following the presentation [28]. In the second approach, we modify the
implicit non-negativity mapping to allow for entries to be set equal to zero.
In both cases, choose an optimal step size based on the new iterations.

#### Sparsity-Promoting MRNSD

We consider a sparsity-promoting MRNSD that minimizes the non-negativity-
constrained lasso problem as introduced in (4). The $\ell_{1}$-norm promotes
sparsity when the solution switches from positive to negative entries.
However, because unregularized MRNSD imposes the implicit mapping that
enforces strict positivity of the solution, the $\ell_{1}$-regularization will
not promote sparsity. In response, we slightly modify the implicit mapping as
$\mathbf{x}=\textnormal{e}^{\mathbf{z}}-\text{\boldmath$\mathbf{\varepsilon}$}$
where $\text{\boldmath$\mathbf{\varepsilon}$}>{\bf 0}$ to allow for small
negative values. Following proximal operator theory [31], the resulting
sparsity-promoting MRNSD iteration is

$\displaystyle\mathbf{x}\leftarrow
S_{\alpha\lambda}(\mathbf{x}+\alpha\mathbf{s}),$ (15)

where the soft-thresholding operator $S_{\alpha\lambda}$ is applied entry-wise
and given by the piecewise function111In Matlab, we write softThreshold(y,c) =
sign(y) .* max(abs(y) - c, 0).

$\displaystyle(S_{\gamma}(\mathbf{x}))_{i}=\begin{cases}x_{i}+\gamma,&\mbox{for
}x_{i}<-\gamma,\\\ 0,&\mbox{for }|x_{i}|\leq\gamma,\\\ x_{i}-\gamma,&\mbox{for
}x_{i}>\gamma,\end{cases}\qquad\mbox{for }i=1,\ldots,K.$ (16)

In this work, we choose the optimal step size for the soft-threshold-updated
solution by solving

$\displaystyle\alpha^{*}\in\min_{\alpha\in[0,u]}\tfrac{1}{2}\|\mathbf{A}S_{\alpha\lambda}(\mathbf{x}+\alpha\mathbf{s})-\mathbf{b}\|_{2}^{2},\qquad\text{where}\quad
u=\min_{s_{i}<-\lambda}\\{-x_{i}/(s_{i}+\lambda)\\}_{i=1}^{K}$ (17)

at each iteration. Solving the scalar optimization problem (17) can be
performed efficiently utilizing nonlinear solvers for single variables such as
golden section search with parabolic interpolation [1]. We summarize sparsity-
promoting MRNSD (spMRNSD) in Algorithm 2. We note that in practice we can pre-
apply the operator $\mathbf{A}$ before solving (17), significantly reducing
the computational cost.

Note that the upper bound $u$ in (17) is slightly different from the boundary
for the original MRNSD step size (Algorithm 1, 5). Specifically, when we step
in $\mathbf{x}$-space in MRNSD, we adjust the step size $\alpha$ to ensure
that the step $x_{i}+\alpha s_{i}$ does not violate the non-negativity
constraint. Similarly, for spMRNSD, we must ensure that the smallest value in
$S_{\alpha\lambda}$ remains non-negative; that is, whenever $x_{i}+\alpha
s_{i}<-\alpha\lambda$ or rewritten $x_{i}+\alpha(s_{i}+\lambda)<0$ is
satisfied. Because $x_{i}$, $\lambda$ and $\alpha$ are assumed to be non-
negative, this condition arises when $s_{i}$ is negative and sufficiently
smaller than $-\lambda$. Thus, $\alpha$ cannot be larger than
$-x_{i}/(s_{i}+\lambda)$, which completes the upper bound of (17).

The original sparsity-promoting MRNSD from [28] used the optimal MRNSD step
size in Algorithm 1, 5. Instead of the optimal step size for the soft-
thresholding operator. We note that the MRNSD step size is acceptable for
spMRNSD in the sense that it will not introduce negative values in
$\mathbf{x}$, but may produce inferior updates. However, an optimal step size
based on the soft-thresholding operator could potentially be larger because it
allows for negative search directions up to the level of $-\lambda$.
Additionally, a larger step size will make the middle region of (16) wider,
potentially yielding sparser solutions. This being said, the MRNSD step size
may be larger or yield sparser solutions for certain problems, but these cases
are difficult to predict and further analysis will be considered for future
directions.

Algorithm 2 Sparsity-Promoting MRNSD (spMRNSD) [28]

1:Solve $\min_{\mathbf{x}\geq{\bf 0}}\
\frac{1}{2}\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}^{2}+\lambda\|\mathbf{x}\|_{1}$

2:while not converged do

3: $\mathbf{g}=\mathbf{A}^{\top}(\mathbf{A}\mathbf{x}-\mathbf{b})$

4: $\mathbf{s}=-{\rm diag\\!}\left(\mathbf{x}\right)\mathbf{g}$

5: compute $\alpha^{*}$ solving (17)

6: $\mathbf{x}\leftarrow
S_{\alpha^{*}\lambda}(\mathbf{x}+\alpha^{*}\mathbf{s})$

7:end while

#### Sparsity-Promoting Non-negative Mapping

As a competing approach to spMRNSD, we focus on adapting the non-negativity
MRNSD mapping. The core of MRNSD is to represent the solution via the
_implicit_ non-negative mapping $\mathbf{x}=\textnormal{e}^{\mathbf{z}}$.
Here, we consider a broader class of point-wise, non-negativity mappings
$w:\mathbb{R}^{K}\to\mathbb{R}_{+}^{K}$ and solve

$\displaystyle\min_{\mathbf{z}}\quad\tfrac{1}{2}\left\|\mathbf{A}w(\mathbf{z})-\mathbf{b}\right\|_{2}^{2}.$
(18)

There are two key differences between our approach and traditional MRNSD. The
first is that MRNSD performs all of the updates in $\mathbf{x}$-space,
potentially yielding small step sizes near the boundary of the non-negative
domain. In comparison, our approach will take steps in $\mathbf{z}$-space,
which is unconstrained and can allow for larger steps. The second difference
is the broader choice of mapping. MRNSD restricts itself to the exponential
mapping. However, any other non-negativity mapping may be viable as well,
e.g., $\mathbf{x}=w(\mathbf{z})=\mathbf{z}^{2}$, or
$\mathbf{x}=w(\mathbf{z})=\left|\mathbf{z}\right|$. Here, we propose the
following,

$\displaystyle
w_{a,c}(z)=\begin{cases}\textnormal{e}^{a(z-c)}-a(z-c)-1,&\mbox{for }z>c,\\\
0,&\mbox{for }z\leq c,\end{cases}$ (19)

where $a>0$ controls the steepness of the mapping and $c$ controls where the
mapping switches from exponential to zero. The mapping is continuously
differentiable and allows for the solution $\mathbf{x}=w(\mathbf{z})$ to have
entries equal to zero. An illustrative toy example of how different mappings
perform is provided in Figure 3.

$\widehat{\mathbf{x}}$$\mathbf{x}^{*}$$\mathbf{x}$$w_{1,0}(\mathbf{z})$$|\mathbf{z}|$$\mathbf{z}^{2}$$e^{\mathbf{z}}$
Figure 3: Consider the toy problem $\min_{\mathbf{x}\geq{\bf 0}}\
\left\|\mathbf{A}\mathbf{x}-\mathbf{b}\right\|_{2}^{2}$, with
$\mathbf{A}=\begin{bmatrix}20&5\\\ 5&20\end{bmatrix}$ and
$\mathbf{b}=\begin{bmatrix}2\\\ 23\end{bmatrix}$. The unconstrained minimizer
is at $\mathbf{x}^{*}=1/10\begin{bmatrix}-2&12\end{bmatrix}^{\top}$ (purple),
while the minimizer of the constrained optimization problem is obtained at
$\widehat{\mathbf{x}}=\begin{bmatrix}0&1.1059\end{bmatrix}^{\top}$ (dark red).
We illustrate the different steepest descent paths in the $\mathbf{x}$-space.
Let the current iteration start from
$\mathbf{x}=\begin{bmatrix}0.15&0.13\end{bmatrix}^{\top}$. We depict various
choices for MRNSD (line search in $\mathbf{x}$, dashed lines) and spNNGD (line
search in $\mathbf{z}$, solid lines), with four different functions for $w$,
i.e., $w(\mathbf{z})=\mathbf{z}^{2}$ (blue),
$w(\mathbf{z})=\textnormal{e}^{\mathbf{z}}$ (red),
$w(\mathbf{z})=\left|\mathbf{z}\right|$ (yellow),
$w(\mathbf{z})=w_{1,0}(\mathbf{z})$ from (19) (green). Note that with optimal
step size selection, spNNGD with the $w_{1,0}(\mathbf{z})$ would find the
minimizer in just one iteration.

This new method provides sparsity while ensuring non-negativity in
$\mathbf{x}$, hence we refer to this method as Gradient Descent with a
Sparsity-Promoting Non-negative Mapping (spNNGD). We summarize the proposed
approach in Algorithm 3.

Algorithm 3 Gradient Descent with Sparsity-Promoting Non-negative Mapping
(spNNGD)

1:Solve
$\min_{\mathbf{z}}\frac{1}{2}\|\mathbf{A}w(\mathbf{z})-\mathbf{b}\|_{2}^{2}$

2:while not converged do

3: $\mathbf{g}=\mathbf{A}^{\top}(\mathbf{A}w(\mathbf{z})-\mathbf{b})$

4: $\mathbf{s}=-{\rm diag\\!}\left(w^{\prime}(\mathbf{z})\right)\mathbf{g}$

5: choose an optimal step size  $\alpha^{*}\in\min_{\alpha>0}\
\tfrac{1}{2}\|\mathbf{A}w(\mathbf{z}+\alpha\mathbf{s})-\mathbf{b}\|_{2}^{2}\hskip
199.0pt$

6: $\mathbf{z}\leftarrow\mathbf{z}+\alpha^{*}\mathbf{s}$

7:end while

This approach shares some similarities with sparsity-promoting MRNSD. First,
the sparsity pattern will not change during the iterations; that is, once a
zero appears in the search direction, we do not change the value of
$\mathbf{z}$, so the zero remains throughout the iterations. This is because
$w^{\prime}(\mathbf{z})$ will have the same pattern of zeros as
$w(\mathbf{z})$. Second, the implicit vertical shift defined by
$\mathbf{\varepsilon}$ plays a similar role to $a$ and $c$ in that it enables
zero values to exist in the mapping.

##  4 Numerical Experiments

###  4.1 Experiment Setup

In the following, we discuss four experiments comparing the effectiveness of
our two proposed methods. Specifically, we investigate an image deblurring, an
image completion, a medical tomography, and a super-resolution problem. For
reproducible science, we made the experiments and code available at
https://github.com/elizabethnewman/DL4IP.

#### Metrics and Presented Results

We consider three metrics that measure (relative) accuracy and sparsity,
respectively, given by

rel. residual
$\displaystyle=\frac{\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_{2}}{\|\mathbf{b}\|_{2}},$
rel. error $\displaystyle=\frac{\|\mathbf{y}_{\rm
true}-\mathbf{G}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}},$ and rel.
sparsity
$\displaystyle=\frac{\texttt{nnz}(\mathbf{x})}{\texttt{nnz}(\mathbf{y}_{\rm
true})},$ (20)

where $\texttt{nnz}(\cdot)$ computes the number of nonzero entries for the
given array. For all metrics, lower values indicate better performance. While
the relative residual and the relative error are common measures of quality,
the relative sparsity metric, storing the cost of the coefficients
$\mathbf{x}$ relative to storing the original image, $\mathbf{y}_{\rm true}$,
is a critical measure in dictionary learning. The image’s storage cost is
reduced when the relative sparsity is below $1$. For every experiment and each
algorithm, we present the reconstruction results using the above metrics and
the convergence results, including the relative residual, a proxy for the
relative sparsity, and the optimal step sizes for each iteration.

#### Operators.

Each experiment in this section relies on a specific linear operator
$\mathbf{C}$. The shapes and sparsity patterns vary significantly for each
application. For completeness, we illustrate the sparsity patterns of the
operators in LABEL:fig:operators_C. In every problem, we incorporate a
generalized Tikhonov regularization on the approximation
$\mu\mathcal{R}(\mathbf{y})=\mu\left\|\mathbf{L}\mathbf{y}\right\|_{2}^{2}$.
We utilize different regularization operators $\mathbf{L}$ depending on the
problem. All of the operators we consider promote smoother approximations, and
hence are variations of two-dimensional finite difference operators. By a
“finite difference” operator, we mean that each row of $\mathbf{L}$ contains a
single $1$ and a single $-1$ to compute the difference between pixels. The
three finite differencing operators we consider are the following: the
difference between (1) all adjacent pixels, (2) adjacent pixels at the
boundary of patches, and (3) adjacent patches. The sparsity patterns of the
three operators are depicted in LABEL:fig:operators_L.

$\texttt{size}(\mathbf{C})=$$\frac{\texttt{nnz}(\mathbf{C})}{\texttt{numel}(\mathbf{C})}\approx$$65,\\!536\times
65,\\!536$$3.8\times 10^{-4}$
$\frac{\texttt{nnz}(\mathbf{C})}{\texttt{numel}(\mathbf{C})}\approx$Deblurring$6,\\!554\times
16,\\!384$$6.1\times 10^{-5}$
$\frac{\texttt{nnz}(\mathbf{C})}{\texttt{numel}(\mathbf{C})}\approx$Dead
Pixel$36,\\!200\times 65,\\!536$$3.5\times 10^{-3}$
$\frac{\texttt{nnz}(\mathbf{C})}{\texttt{numel}(\mathbf{C})}\approx$Tomography$40,\\!960\times
262,\\!144$$2.5\times 10^{-4}$
$\frac{\texttt{nnz}(\mathbf{C})}{\texttt{numel}(\mathbf{C})}\approx$Superresolution
(a)

Pixel Smoother$[(N-1)M+M(N-1)]\times MN$Patch Border
Smoother$[(\frac{M}{p}(\frac{N}{q}-1)+\frac{N}{q}(\frac{M}{p}-1)]\times
MN$Patch
Smoother$[pq((\frac{M}{p}(\frac{N}{q}-1)+\frac{N}{q}(\frac{M}{p}-1))]\times
MN$ (b)

Figure 4: Sparsity patterns of operators. The operators are illustrated with
the correct proportions based on their row and column dimensions. The width of
the operators is correctly scaled (e.g., the superresolution operator has
twice as many columns as the deblurring operator).

###  4.2 Deblurring

Capturing digital images can often introduce noise, blur, and other
distortions that hide the true representation. Image deblurring has been a
long-studied technique to remove unwanted artifacts in and enhance digital
images [19]. Here, we consider applying the Gaussian blurring operator to an
$N\times N$ image. To be precise, we construct the blurring operator
$\mathbf{C}\in\mathbb{R}^{N^{2}\times N^{2}}$ as a banded, Toeplitz matrix
with exponentially decaying bands222In Matlab, we write  C1 =
toeplitz(exp(-(0:b-1) / (2 * sigma * sigma)), zeros(1, n - b)); C =
kron(C1,C1);  to construct this blurring operator.. We use a bandwidth of
$b=3$ and an exponential decay rate of $\sigma=4$. Given a vectorized true
image $\mathbf{y}_{\rm true}\in\mathbb{R}_{+}^{N^{2}}$, we construct a
blurred, noisy image $\mathbf{b}=\mathbf{C}\mathbf{y}_{\rm
true}+\text{\boldmath$\mathbf{\eta}$}$ where each entry of the noise vector
$\mathbf{\eta}$ is drawn from a Gaussian distribution $\mathcal{N}({\bf
0},\beta\|\mathbf{C}\mathbf{y}_{\rm true}\|_{2}^{2}\mathbf{I}_{N^{2}})$ with
noise level $\beta=10^{-4}$. We include a pixel smoothing Tikhonov regularizer
term with regularization parameter $\mu=10^{-4}$. We train for $100$
iterations for both spNNGD and spMRNSD and report the results in Figure 5 and
the convergence behavior in Figure 6.

Original, $\mathbf{y}_{\rm true}$
---
$256\times 256$
Blurred, Noisy, $\mathbf{b}$
---
$256\times 256$
spMRNSD, $\mathbf{G}\mathbf{x}$
---
$\lambda=10^{-8}$
spNNGD, $\mathbf{G}\mathbf{x}$
---
$a=0.1$, $c=-0.75$
$0.011$$0$$0.007$$0$$\frac{\|\mathbf{y}_{\rm
true}-{\mathbf{G}\mathbf{x}}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$
$\bf 1.15\times
10^{-1}$$\frac{\texttt{nnz(coefficients)}}{\texttt{nnz(image)}}\approx$$\bf
0.15$$\frac{\|\mathbf{y}_{\rm
true}-{\mathbf{G}\mathbf{x}}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$
$1.57\times
10^{-1}$$\frac{\texttt{nnz(coefficients)}}{\texttt{nnz(image)}}\approx$$0.23$
absolute difference
---
$|\mathbf{y}_{\rm true}-\mathbf{G}\mathbf{x}|$
relative error, $\frac{\|\mathbf{y}_{\rm
true}-{\mathbf{G}\mathbf{x}}\|_{2}}{\|\mathbf{y}_{\rm
true}\|_{2}}\approx$relative sparsity,
$\frac{\texttt{nnz}(\mathbf{x})}{\texttt{nnz}(\mathbf{y}_{\rm
true})}\approx$matricized coefficients, $\mathbf{x}\equiv$ Figure 5:
Reconstructions for the image deblurring task. The top row (left-to-right)
contains the original image, the blurred, noisy image, the reconstruction
using spMRNSD, and the reconstruction using spNNGD. The middle row displays
the absolute differences between the original and the reconstructions with the
relative error and sparsity metrics below. Here, spMRNSD produces the best
approximation in terms of both metrics of interest. The bottom row depicts the
sparsity patterns of matricized solutions $\mathbf{x}$ of size $400\times 256$
where each column corresponds to the coefficients for one vectorized
dictionary patch.

In Figure 5, we observe that both spNNGD and spMRNSD produce qualitatively
good approximations to the original image for similar sparsity levels. Both
sets of coefficients require less than $25\%$ of the storage cost of the
original image. Both approximations produce some blocky artifacts due to the
dictionary patch representation; this is most noticeable in the spNNGD
reconstruction of small vegetables, such as the garlic in the bottom right. In
the bottom row, we observe similar sparsity patterns in the matricized
coefficients for spMRNSD and spNNGD. In particular, the two have some rows
with many nonzeros and some sparse rows, and the rows with high density are
similar for both solutions. It can be observed, that the solution obtained
from spMRNSD is slightly sparser than that obtained from spNNGD.

(a)

(b)

(c)

Figure 6: Convergence behavior for the deblurring problem. The relative
residual norm (left) converges at similar rates for both algorithms. The proxy
for relative sparsity (middle) converges faster for spNNGD. On average, the
step sizes (right) for spNNGD are larger and more consistent.

The convergence behavior in LABEL:fig:deblur_relerr shows that spNNGD and
spMRNSD converge to similar levels of accuracy at similar rates. In
comparison, the sparsity convergence plot LABEL:fig:deblur_sparsity, we see
that spNNGD converges quickly to a sparse solution and then stagnates, whereas
spMRNSD converges to a sparse solution more gradually. This behavior indicates
that spNNGD fixes a sparsity pattern early and learns the best coefficients
for that pattern while spMRNSD adjusts the sparsity pattern during iterations.
We note that both methods are constrained by the current sparsity pattern
because the search direction $\mathbf{s}={\rm
diag\\!}\left(w^{\prime}(\mathbf{z})\right)\mathbf{g}$ has the same zero
entries as the approximate solution $\mathbf{x}=w(\mathbf{z})$. The step size
plot in LABEL:fig:deblur_alpha supports our original hypothesis that we take
larger step sizes in $\mathbf{z}$-space (spNNGD) than in $\mathbf{x}$-space
(spMRNSD).

Relative ErrorRelative Sparsity$1$$0$ (a)

(b)

Figure 7: Parameter sensitivity exploration for the deblurring problem.

Examining the relative error plot, it becomes apparent that the relative error
of spNNGD is more susceptible to changes in the steepness parameter $a$ rather
than the shift parameter $c$. This observation is substantiated by the
relatively constant rows in the relative error plot. We also observe a
positive correlation between $a$ and $c$, implying that increasing both
parameters simultaneously results in small changes to the relative error.
However, while the relative error exhibits a degree of sensitivity to the
choice of parameters, the sparsity plot appears to be less sensitive to its
parameter selection.

For spMRNSD, we observe that the relative error remains mainly constant across
a wide range of regularization parameter $\lambda$ values, while the sparsity
undergoes drastic changes. This behavior is consistent with expectations, that
larger values for $\lambda$ will lead to poorer approximations and sparser
solutions.

There are several major takeaways from our exploration of the image deblurring
problem. First, the spMRNSD produced better results in terms of approximation
quality and sparsity, but spNNGD is competitive. In terms of convergence,
spNNGD obtains a sparser solution in fewer iterations than spMRNSD with larger
step sizes, but as a result, prescribes a sparsity pattern of the solution
early. In comparison, spMRNSD gradually increases sparsity which can bring
benefits in the long run. In terms of sensitivity to parameters, the spNNGD
parameters can drastically change the relative error, but the sparsity stays
more consistent. Conversely, the regularization parameter for spMRNSD can
drastically change the sparsity but does not impact the relative error as
dramatically.

###  4.3 Image Completion

Matrix completion arises in numerous fields [7, 24], such as collaborative
filtering to predict users preferences with incomplete information [5], pixels
[2], recovery of incomplete signals [12], and digital image inpainting to
restore damaged or missing information. There are similarly many different
approaches to completing matrices, including forming a low-rank approximation
and solving an inverse problem. In this section, we explore the inverse
problem approach for digital inpainting applications.

We construct the missing pixel operator $\mathbf{C}\in\mathbb{R}^{m\times
N^{2}}$ with $m\ll N^{2}$ that samples uniformly random pixels from an
$N\times N$ image. In practice, the rows of $\mathbf{C}$ are a subset of rows
from the $N^{2}\times N^{2}$ identity matrix. As a result,
$\mathbf{b}=\mathbf{C}\mathbf{y}$ contains the same as the original image with
a percentage of pixels randomly set to zero. In our experiment, we use
Johannes Vermeer’s _Girl with the Pearl Earring_
painting,333https://en.wikipedia.org/wiki/Girl_with_a_Pearl_Earring resized to
$N\times N=128\times 128$, and removed $60\%$ of the pixels in our
experiments. We report the representation results in Figure 8 and the
convergence behavior in Figure 9.

Original, $\mathbf{y}_{\rm true}$
---
$128\times 128$
Missing Pixels, $\mathbf{b}$
---
$60\%$ removed
spMRNSD, $\mathbf{G}\mathbf{x}$
---
$\lambda=10^{-3}$
spNNGD, $\mathbf{G}\mathbf{x}$
---
$a=1$, $c=-0.5$
absolute difference
---
$|\mathbf{y}_{\rm true}-\mathbf{G}\mathbf{x}|$
$1$$0.047$$0.5$$0$$\frac{\|\mathbf{y}_{\rm
true}-\mathbf{G}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$ $\bf
1.65\times
10^{-1}$$\frac{\texttt{nnz($\mathbf{x}$)}}{\texttt{nnz($\mathbf{y}_{\rm
true}$)}}\approx$ $\bf 0.09$$\frac{\|\mathbf{y}_{\rm
true}-\mathbf{G}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$
$1.83\times
10^{-1}$$\frac{\texttt{nnz($\mathbf{x}$)}}{\texttt{nnz($\mathbf{y}_{\rm
true}$)}}\approx$ $0.28$relative error, $\frac{\|\mathbf{y}_{\rm
true}-\mathbf{G}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$
relative sparsity,
$\frac{\texttt{nnz($\mathbf{x}$)}}{\texttt{nnz($\mathbf{y}_{\rm
true}$)}}\approx$matricized coefficients, $\mathbf{x}\equiv$ Figure 8:
Reconstructions for the image completion task. The top row (left-to-right)
contains the original image, the blurred, noisy image, the reconstruction
using spMRNSD, and the reconstruction using spNNGD. The middle row displays
the absolute differences between the original and the reconstructions with the
relative error and sparsity metrics below. Here, spMRNSD produces the best
approximation in terms of both metrics of interest. The bottom row depicts the
sparsity patterns of matricized solutions $\mathbf{x}$ of size $400\times 64$
where each column corresponds to the coefficients for one vectorized
dictionary patch.

In Figure 8, we observe that both spNNGD and spMRNSD produce qualitatively
good approximations of the original image. The approximation obtained from
spMRNSD is a sharper approximation to the original, particularly for the
facial features in the image. Additionally, the spMRNSD-generated coefficients
are notably sparser than those generated by spNNGD. The sparsity patterns of
the un-vectorized coefficients reflect these metrics.

(a)

(b)

(c)

Figure 9: Convergence behavior for the image completion problem. The relative
residual norm (left) converges at similar rates for both algorithms. The proxy
for relative sparsity (middle) converges faster for spNNGD. On average, the
step sizes (right) for spMRNSD are larger and the step sizes for spNNGD are
more consistent.

The convergence behavior in LABEL:fig:indicator_relerr and
LABEL:fig:indicator_sparsity matches that of Section 4.2. In particular, it
shows that both algorithms converge to similar levels of accuracy at similar
rates and spNNGD quickly reaches a sparse solution while spMRNSD converges
more gradually. Notably and surprisingly, the step size plot in
LABEL:fig:indicator_alpha reveals that larger steps are taken in
$\mathbf{x}$-space (spMRNSD) than in $\mathbf{z}$-space (spNNGD). The
difference in behavior in terms of both relative sparsity and step size can be
attributed to the image completion operator and data normalization.
Specifically, $\mathbf{C}$ has rows from an identity matrix, the dictionary
$\mathbf{G}$ contains entries between $0$ and $1$, and the right-hand side is
normalized so its maximum value is $1$. Thus, the approximation
$\mathbf{C}\mathbf{G}w_{a,c}(\mathbf{z})$ is more directly sensitive to values
in the coefficients. For spNNGD, a large update in $\mathbf{z}$ could
introduce unwanted zeros or could yield large values for
$w_{a,c}(\mathbf{z})$. For this reason, we take smaller step sizes. In
comparison, for spMRNSD, the soft thresholding operator reduces large values
and zeroes out small values, enabling slightly larger step sizes.

To summarize, the main takeaways from the image completion experiment are that
spNNGD and spMRNSD produce approximations with similar relative errors, and
spMRNSD produces a sparser solution. As before, spNNGD converges to a sparser
solution more quickly, but the gradual progress of spMRNSD produces better
metrics overall. The operator plays a role in determining the effectiveness of
the chosen method and may inform the preferred approach.

###  4.4 Tomography

Tomography is an imaging technique that reconstructs the characteristics of an
object’s interior by taking cross-sectional measurements from the exterior.
Tomography is widely used in a variety of applications ranging from but not
limited to geophysics, atmospheric science, physics, materials science, and
medical imaging. In the medical field, computer tomography (CT) is a common
non-invasive technique used for diagnostic purposes and surgical planning. CT
employs penetrating X-rays to pass through an object, with the amount of
energy absorbed dependent on the object’s properties along the X-ray path.
Detectors on the opposite side measure the intensity of the X-rays emitted by
the source. By rotating both the X-ray source and receiver around the object,
measurements are collected in a sinogram. CT is a classic example of an
inverse problem, where the goal is to deduce the object’s energy absorbency
from the sinogram. For illustration, we consider the 2D computer tomography
[27, 20]. The ill-posedness of this problem requires regularization and a
common regularization functional is total variation.

x-ray sourcedetector Figure 10: Illustration of parallel beam tomography,
[33].

In this experiment, we leverage the setup from IRTools [14] and AIRTools II
[18] to construct a parallel beam tomography problem; see Figure 10 for
illustration. We reconstruct the Shepp-Logan phantom of size $N\times
N=256\times 256$. The model matrix $\mathbf{C}$ is constructed using
$\lfloor\sqrt{2}N\rfloor=362$ rays and $100$ rotation angles equispaced
between $0$ and $179$ degrees. For the generalized Tikhonov regularization, we
use a patch smoothing operator with $\mu=1$. We train for $100$ iterations,
use the flower dictionary for representation, and report our results in Figure
11 and Figure 12.

Original, $\mathbf{y}_{\rm true}$
---
$256\times 256$
Sinogram, $\mathbf{b}$
---
$n_{\rm rays}\times n_{\rm proj}=362\times 100$
spMRNSD, $\mathbf{G}\mathbf{x}$
---
$\lambda=10^{0}$
spNNGD, $\mathbf{G}\mathbf{x}$
---
$a=0.5$, $c=-1$
absolute difference
---
$|\mathbf{y}_{\rm true}-\mathbf{G}\mathbf{x}|$
$1$$0$$1.08$$0$$\frac{\|\mathbf{y}_{\rm
true}-\mathbf{A}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$
$4.99\times
10^{-1}$$\frac{\texttt{nnz(coefficients)}}{\texttt{nnz(image)}}\approx$
$0.45$$\frac{\|\mathbf{y}_{\rm
true}-\mathbf{G}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$ $\bf
4.84\times
10^{-1}$$\frac{\texttt{nnz(coefficients)}}{\texttt{nnz(image)}}\approx$ $\bf
0.17$relative error, $\frac{\|\mathbf{y}_{\rm
true}-\mathbf{G}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$
relative sparsity,
$\frac{\texttt{nnz($\mathbf{x}$)}}{\texttt{nnz($\mathbf{y}_{\rm
true}$)}}\approx$matricized coefficients, $\mathbf{x}\equiv$ Figure 11:
Reconstructions for the computing tomography task. The top row (left-to-right)
contains the original image, the blurred, noisy image, the reconstruction
using spMRNSD, and the reconstruction using spNNGD. The middle row displays
the absolute differences between the original and the reconstructions with the
relative error and sparsity metrics below. Here, spMRNSD produces the best
approximation in terms of both metrics of interest. The bottom row depicts the
sparsity patterns of matricized solutions $\mathbf{x}$ of size $400\times 256$
where each column corresponds to the coefficients for one vectorized
dictionary patch.

In Figure 12, we observe spNNGD produces the best approximation results
combined with the sparsest solution. The approximation is particularly sharp
in the center of the Shepp-Logan phantom, while the reconstruction obtained
from spMRNSD retains many artifacts from the dictionary patches. The sparsity
patterns of the solutions are slightly different for the two algorithms as
well. The coefficients from spMRNSD are dense in the center while the spNNGD
contain more nonzeros at the edges. We note that for both algorithms, the
largest errors occur on the edges of the Shepp-Logan phantom. This is perhaps
unsurprising because the dictionary elements from the flower dictionary do not
contain similar curved features. A different dictionary prior could yield
improved performance on the edge features, but learning the dictionaries is
outside the scope of this work.

We also mention that the structure of the data motivations for using a patch
smoother rather than a pixel smoother or patch border smoother. We train a
dictionary on natural images, but the Shepp-Logan phantom is constructed
synthetically. The Shepp-Logan patches are closer to piecewise constant than
natural image patches. Using a patch smoother, we encourage patches of the
images to be more similar, which leads to larger regions of almost constant
pixel values.

(a)

(b)

(c)

Figure 12: Convergence behavior for the computed tomography problem. The
relative residual norm (left) converges faster for spNNGD and exhibits more
erratic behavior for spMNRSD. The proxy for relative sparsity (middle)
converges faster for spNNGD. On average, the step sizes (right) for spNNGD are
larger and more consistent.

As in the previous experiments, spNNGD converges more quickly to a smaller
relative error and a sparser solution and takes larger step sizes on average.
However, we observe starkly different convergence behavior in
LABEL:fig:tomography_relerr than in the previous experiments. The convergence
of spMRNSD is erratic after approximately $30$ iterations. This behavior is
due to the soft-thresholding operator and the mismatch of a search direction
in $\mathbf{z}$-space and an update in $\mathbf{x}$-space. These two factors
make it possible for the algorithm to step in a non-descent direction
$\mathbf{x}$-space. We see this occur when the optimal $\alpha$ is large and
as the solution becomes sparser. In contrast, spNNGD, which computes search
directions and steps in $\mathbf{z}$-space, consistently decays.

###  4.5 Superresolution

Superresolution is a technique to reconstruct a high-resolution image by a
sequence of low-resolution images to reveal image details that are
unobservable in lower resolution. This task arises in various applications,
from enlarging and restoring blurry or damaged photos to improving medical
imaging and enhancing surveillance footage [42]. For simplicity, let us assume
we desire to reconstruct a high-resolution image $\mathbf{x}_{\rm
true}\in\mathbb{R}^{N^{2}}$ from a sequence of low-resolution images
$\mathbf{b}_{1},\ldots,\mathbf{b}_{M}\in\mathbb{R}^{\nu^{2}}$ with $\nu^{2}\ll
N^{2}$. We assume that the (ridged) image transformation, such as shifts,
rotations, and scaling of the high-resolution image to each low-resolution
image is known and denoted by
$\mathbf{S}_{1},\ldots,\mathbf{S}_{M}\in\mathbb{R}^{N^{2}\times N^{2}}$. With
a downsampling/restriction operator $\mathbf{R}\in\mathbb{R}^{\nu^{2}\times
N^{2}}$ we state corresponding superresolution problem as

$\min_{\mathbf{x}\geq 0}\quad\left\|\begin{bmatrix}\mathbf{R}\mathbf{S}_{1}\\\
\vdots\\\ \mathbf{R}\mathbf{S}_{M}\\\ \end{bmatrix}\mathbf{y}_{\rm
true}-\begin{bmatrix}\mathbf{b}_{1}\\\ \vdots\\\ \mathbf{b}_{M}\\\
\end{bmatrix}\right\|_{2}^{2}+\mu\left\|\mathbf{L}\mathbf{y}_{\rm
true}\right\|_{2}^{2}$ (21)

where we selected a $\ell_{2}$ regularization term, [35].

In our experiments, we consider an resized $N\times N=512\times 512$ builtin
Matlab image of the moon.444img = imread(’moon.tif’); We form $\theta=10$ low-
resolution images, downsampled by a factor of eight. Thus, the low resolution
images $\mathbf{B}_{i}$ are of size $64\times 64$ and the restriction operator
is a matrix $\mathbf{R}_{+}\in\mathbb{R}^{(N/8)^{2}\times N^{2}}$. Poetically,
we reconstruct the moon using the earth dictionary (Figure 2). We use $100$
iterations and a patch neighbor smoothing regularizer with $\mu=10^{-2}$. We
report the results of the superresolution experiment in Figure 13 and Figure
15.

Original, $\mathbf{y}_{\rm true}$
---
$512\times 512$
Low Resolution, $\mathbf{b}_{i}$
---
$10$, $64\times 64$ images
spMRNSD, $\mathbf{G}\mathbf{x}$
---
$\lambda=10^{-4}$
spNNGD, $\mathbf{G}\mathbf{x}$
---
$a=0.2975$, $c=-0.05$
absolute difference
---
$|\mathbf{y}_{\rm true}-\mathbf{G}\mathbf{x}|$
$1$$0$$1$$0$$\frac{\|\mathbf{x}-\hat{\mathbf{x}}\|_{2}}{\|\mathbf{x}\|_{2}}\approx$
$\bf 1.10\times
10^{-1}$$\frac{\texttt{nnz(coefficients)}}{\texttt{nnz(image)}}\approx$ $\bf
0.19$$\frac{\|\mathbf{x}-\hat{\mathbf{x}}\|_{2}}{\|\mathbf{x}\|_{2}}\approx$
$1.35\times
10^{-1}$$\frac{\texttt{nnz(coefficients)}}{\texttt{nnz(image)}}\approx$
$0.98$relative error, $\frac{\|\mathbf{y}_{\rm
true}-\mathbf{G}\mathbf{x}\|_{2}}{\|\mathbf{y}_{\rm true}\|_{2}}\approx$
relative sparsity,
$\frac{\texttt{nnz($\mathbf{x}$)}}{\texttt{nnz($\mathbf{y}_{\rm
true}$)}}\approx$matricized coefficients, $\mathbf{x}\equiv$ Figure 13:
Reconstructions for the superresolution task. The top row (left-to-right)
contains the original image, the blurred, noisy image, the reconstruction
using spMRNSD, and the reconstruction using spNNGD. The middle row displays
the absolute differences between the original and the reconstructions with the
relative error and sparsity metrics below. Here, spMRNSD produces the best
approximation in terms of both metrics of interest. The bottom row depicts the
sparsity patterns of matricized solutions $\mathbf{x}$ of size $400\times
1024$ where each column corresponds to the coefficients for one vectorized
dictionary patch. Note that the lower right patch for both approximations has
significantly more error than the other patches due to the boundary conditions
of the operator.

In Figure 13, we observe that both algorithms produce similar quality
approximations, and both have the same artifact in the bottom right corner. We
conjecture that the cause of this artifact is twofold. First, the location of
this error may be due to the operator $\mathbf{C}$ which rotates the images
counterclockwise about the upper left corner. The second is that the initial
guess plays a crucial role in the approximation, and we see this artifact
appear in early iterations (see Figure 14). In this work, we opted for
sacrificing some approximation quality to promote more sparsity, and hence
focus on the results with the lower right corner error.

spNNGDspMRNSD$1$$0$iteration Figure 14: Superresolution approximations for
first five iterations. Approximations from both algorithms contain the lower
right corner due to the boundary conditions of the operator. We show the
approximations on the same color scale as the original image, which can be
misleading in the initial iteration. The first approximation contains values
larger than $1$, hence all of the pixels appear to be yellow.

We note in Figure 13 spNNGD does not produce any sparsity. We conjecture that
sparsity is harder to obtain in superresolution partially because the right-
hand sides already have many nonzeros (e.g.,
$\frac{\texttt{nnz}(\mathbf{b}_{1})}{\texttt{numel}(\mathbf{b}_{1})}\approx
56\%$ and overall
$\frac{\texttt{nnz}(\mathbf{b})}{\texttt{numel}(\mathbf{b})}\approx 75\%$) and
partially because we use small patches to reconstruct a higher resolution
image. Per the latter point, as shown in [28], the ratio between patch size
and image size can significantly affect the sparsity pattern. We have seen in
previous experiments that spNNGD settles on a sparsity pattern in early
iterations. Fixating on one sparsity pattern while maintaining a good
approximation leads to a non-sparse solution. In contrast, spMRNSD gradually
promotes sparsity, and performs better along this metric.

(a)

(b)

(c)

Figure 15: Convergence behavior for the superresolution problem. The relative
residual norm (left) converges at similar rates for both algorithms. The proxy
for relative sparsity (middle) reaches a lower level to spMRNSD. On average,
the step sizes (right) for spNNGD are larger and more consistent.

The first iteration of Figure 14 exhibits similar behavior to the previous
experiments with the two algorithms converging at similar rates in terms of
relative error and spNNGD promoting sparsity faster. However, we quickly see
spNNGD stagnate with minuscule step sizes, eventually halting progress at
about $40$ iterations. As expected, the early step sizes in
LABEL:fig:superresolution_alpha are larger for spNNGD. Interestingly, both
step sizes start to decay at similar rates for iterations $20$ through $30$,
which is the iteration at which spMNRSD stops becoming sparser in
LABEL:fig:superresolution_sparsity.

##  5 Conclusions

We have demonstrated success in solving inverse problems using patch
dictionary representations to obtain non-negative, sparse solutions in image
deblurring (Section 4.2), image completion (Section 4.3), computed tomography
(Section 4.4), and superresolution (Section 4.5). We proposed and compared two
different modifications to the modified residual norm steepest descent (MRNSD)
algorithm to produce solutions with the desired structure. Motivated by [28],
we first introduced a sparsity-promoting MRNSD (spMRNSD) by adding $\ell_{1}$
regularization and iterating with a soft-thresholding operator. Our key
contribution was using the correct optimal step size based on the soft-
thresholding operator update. We compared this soft-thresholding approach to a
new gradient descent approach with a sparsity-promoting, non-negative mapping
in spNNGD. We observed that typically spNNGD converges faster to a sparse
solution than spMRNSD and takes larger step sizes. However, the sparsity
pattern can be fixed after a few iterations, which can limit improvements
along that metric. In comparison, spMRNSD gradually increases sparsity via the
soft-thresholding operator. The performance of each algorithm can vary
significantly based on the task, depending on the given operator, the
algorithmic hyperparameters, and the dictionary used for the representation.
In general, it seems the gradual approach of spMRNSD offers more promise than
the sparsity-promoting, non-negative mapping of spNNGD, but the advantages of
the latter could lead to interesting future work.

The complex relationship between all of the parameters led to our comparison
of the two algorithms and motivates several exciting extensions. For spMRNSD,
we can we can explore other implicit mappings. Currently, we use the
exponential mapping to enable simple steps in the solution space. We could
perform similar implicit mappings and solution space steps with other
mappings. For example, with $\mathbf{x}=\mathbf{z}^{2}$, our search in
solution space (without forming $\mathbf{z}$ explicitly) would use search
directions $\mathbf{s}={\rm
diag\\!}\left(\tfrac{1}{2}\sqrt{\mathbf{x}}\right)\mathbf{g}$. For spNNGD, we
can consider adaptive approaches to modify the hyperparameters as we iterate.
For example, if we observe stagnation of the sparsity pattern, we may
gradually shift the mapping to promote more sparsity. Additionally, we could
consider accelerating the selection of the step size via linearization rather
than solving a one-dimensional optimization problem explicitly. Combining the
two algorithms, we could consider initializing with a few iterations of spNNGD
to reach a sparser solution, and then refine with spMRNSD.

Acknowledgement: This work was partially supported by the National Science
Foundation (NSF) under grants [DMS-2309751] and [DMS-2152661] for Newman and
Chung, respectively. Any opinions, findings, conclusions, or recommendations
expressed in this material are those of the authors and do not necessarily
reflect the views of the National Science Foundation.

## References

  * [1] Uri M Ascher and Chen Greif  “A first course on numerical methods”  SIAM, 2011
  * [2] Marcelo Bertalmio, Guillermo Sapiro, Vincent Caselles and Coloma Ballester  “Image Inpainting”  In _Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques_ , SIGGRAPH ’00  USA: ACM Press/Addison-Wesley Publishing Co., 2000, pp. 417–424  DOI: 10.1145/344779.344972
  * [3] Stephen Boyd et al.  “Distributed optimization and statistical learning via the alternating direction method of multipliers”  In _Foundations and Trends® in Machine learning_ 3.1 Now Publishers, Inc., 2011, pp. 1–122
  * [4] Steven L Brunton and J Nathan Kutz  “Data-driven science and engineering: Machine learning, dynamical systems, and control”  Cambridge University Press, 2019
  * [5] Emmanuel Candes and Yaniv Plan  “Matrix Completion With Noise”  In _Proceedings of the IEEE_ 98, 2010, pp. 925–936  DOI: 10.1109/JPROC.2009.2035722
  * [6] Emmanuel J Candes, Justin K Romberg and Terence Tao  “Stable signal recovery from incomplete and inaccurate measurements”  In _Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences_ 59.8 Wiley Online Library, 2006, pp. 1207–1223
  * [7] Emmanuel J. Candès and Benjamin Recht  “Exact Matrix Completion via Convex Optimization”  In _Foundations of Computational Mathematics_ 9.6, 2009, pp. 717–772  DOI: 10.1007/s10208-009-9045-5
  * [8] Tony F Chan and Jianhong Shen  “Image processing and analysis: variational, PDE, wavelet, and stochastic methods”  SIAM, 2005
  * [9] Scott Shaobing Chen, David L Donoho and Michael A Saunders  “Atomic decomposition by basis pursuit”  In _SIAM review_ 43.1 SIAM, 2001, pp. 129–159
  * [10] Julianne Chung and Silvia Gazzola  “Computational methods for large-scale inverse problems: a survey on hybrid projection methods”  In _Siam Review_ Society for IndustrialApplied Mathematics Publications, 2023
  * [11] Matthias Chung and Rosemary Renaut  “The Variable Projected Augmented Lagrangian Method”, 2022  arXiv:2207.08216 [math.NA]
  * [12] Rong Du, Cailian Chen, Bo Yang and Xinping Guan  “VANET based traffic estimation: A matrix completion approach”  In _2013 IEEE Global Communications Conference (GLOBECOM)_ , 2013, pp. 30–35  DOI: 10.1109/GLOCOM.2013.6831043
  * [13] Michael Elad  “Sparse and redundant representations: from theory to applications in signal and image processing”  Springer, 2010
  * [14] Silvia Gazzola, Per Christian Hansen and James G. Nagy  “IR Tools: a MATLAB package of iterative regularization methods and large-scale test problems”  In _Numerical Algorithms_ 81.3, 2019, pp. 773–811  DOI: 10.1007/s11075-018-0570-7
  * [15] Gene H. Golub, Michael Heath and Grace Wahba  “Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter”  In _Technometrics_ 21.2 Taylor & Francis, 1979, pp. 215–223  DOI: 10.1080/00401706.1979.10489751
  * [16] Per Christian Hansen  “Rank-Deficient and Discrete Ill-Posed Problems”  Society for IndustrialApplied Mathematics, 1998  DOI: 10.1137/1.9780898719697
  * [17] Per Christian Hansen  “Discrete inverse problems: insight and algorithms”  SIAM, 2010
  * [18] Per Christian Hansen and Jakob Sauer Jørgensen  “AIR Tools II: algebraic iterative reconstruction methods, improved implementation”  In _Numerical Algorithms_ 79.1, 2018, pp. 107–137  DOI: 10.1007/s11075-017-0430-x
  * [19] Per Christian Hansen, James G. Nagy and Dianne P. O’Leary  “Deblurring Images”  Society for IndustrialApplied Mathematics, 2006  DOI: 10.1137/1.9780898718874
  * [20] Avinash C Kak, Malcolm Slaney and Ge Wang  “Principles of Computerized Tomographic Imaging”  Wiley Online Library, 2002
  * [21] Linda Kaufman  “Maximum likelihood, least squares, and penalized least square for PET”  In _Medical Imaging, IEEE Transactions on_ 12, 1993, pp. 200–214  DOI: 10.1109/42.232249
  * [22] Dara Kerr  “Stunning high-resolution photo shows Earth’s many hues”, 2012  URL: https://www.cnet.com/science/stunning-high-resolution-photo-shows-earths-many-hues/
  * [23] Daniel D. Lee and H. Seung  “Learning the parts of objects by non-negative matrix factorization”  In _Nature_ 401.6755, 1999, pp. 788–791  DOI: 10.1038/44565
  * [24] Xiao Peng Li, Lei Huang, Hing Cheung So and Bo Zhao  “A Survey on Matrix Completion: Perspective of Signal Processing”, 2019  arXiv:1901.10885 [eess.SP]
  * [25] Stéphane G Mallat and Zhifeng Zhang  “Matching pursuits with time-frequency dictionaries”  In _IEEE Transactions on signal processing_ 41.12 IEEE, 1993, pp. 3397–3415
  * [26] James G. Nagy and Zdenek Strakos  “Enforcing nonnegativity in image reconstruction algorithms”  In _Mathematical Modeling, Estimation, and Imaging_ 4121 SPIE, 2000, pp. 182–190  International Society for OpticsPhotonics  DOI: 10.1117/12.402439
  * [27] Frank Natterer  “The Mathematics of Computerized Tomography”  SIAM, 2001
  * [28] Elizabeth Newman and Misha E. Kilmer  “Nonnegative Tensor Patch Dictionary Approaches for Image Compression and Deblurring Applications”  In _SIAM Journal on Imaging Sciences_ 13.3, 2020, pp. 1084–1112  DOI: 10.1137/19M1297026
  * [29] Jorge Nocedal and Stephen J Wright  “Numerical optimization”  Springer, 1999
  * [30] Stanley Osher, Andrés Solé and Luminita Vese  “Image Decomposition and Restoration Using Total Variation Minimization and the H1”  In _Multiscale Modeling & Simulation_ 1.3, 2003, pp. 349–370  DOI: 10.1137/S1540345902416247
  * [31] Neal Parikh and Stephen Boyd  “Proximal Algorithms”  In _Found. Trends Optim._ 1.3 Hanover, MA, USA: Now Publishers Inc., 2014, pp. 127–239  DOI: 10.1561/2400000003
  * [32] K.. Petersen and M.. Pedersen  “The Matrix Cookbook” Version 20121115  Technical University of Denmark, 2012  URL: http://www2.compute.dtu.dk/pubdb/pubs/3274-full.html
  * [33] Lars Ruthotto, Julianne Chung and Matthias Chung  “Optimal experimental design for inverse problems with state constraints”  In _SIAM Journal on Scientific Computing_ 40.4 SIAM, 2018, pp. B1080–B1100
  * [34] Yousef Saad  “Iterative methods for sparse linear systems”  SIAM, 2003
  * [35] J Tanner Slagel et al.  “Sampled Tikhonov regularization for large linear inverse problems”  In _Inverse Problems_ 35.11 IOP Publishing, 2019, pp. 114008
  * [36] Sara Soltani, Martin S. Andersen and Per Christian Hansen  “Tomographic Image Reconstruction using Dictionary Priors”  In _CoRR_ abs/1503.01993, 2015  arXiv: http://arxiv.org/abs/1503.01993
  * [37] Sara Soltani, Misha E. Kilmer and Per Christian Hansen  “A tensor-based dictionary learning approach to tomographic image reconstruction”  In _BIT Numerical Mathematics_ 56.4, 2016, pp. 1425–1454  DOI: 10.1007/s10543-016-0607-z
  * [38] David Strong and Tony Chan  “Edge-preserving and scale-dependent properties of total variation regularization”  In _Inverse problems_ 19.6 IOP Publishing, 2003, pp. S165
  * [39] The TensorFlow Team  “Flowers”, 2019  URL: http://download.tensorflow.org/example_images/flower_photos.tgz
  * [40] Ivana Tošić and Pascal Frossard  “Dictionary learning”  In _IEEE Signal Processing Magazine_ 28.2 IEEE, 2011, pp. 27–38
  * [41] Joel A Tropp  “Greed is good: Algorithmic results for sparse approximation”  In _IEEE Transactions on Information theory_ 50.10 IEEE, 2004, pp. 2231–2242
  * [42] Jianchao Yang, John Wright, Thomas S Huang and Yi Ma  “Image super-resolution via sparse representation”  In _IEEE transactions on image processing_ 19.11 IEEE, 2010, pp. 2861–2873