# Learning Multiscale Convolutional Dictionaries for Image Reconstruction

Tianlin Liu1, Anadi Chaman2, David Belius1 and Ivan Dokmanić12 This work was
supported by the European Research Council Starting Grant 852821—SWING. The
code to reproduce our experiments is available at
https://github.com/liutianlin0121/MUSC 1Department of Mathematics and Computer
Science, University of Basel.  2Department of Electrical and Computer
Engineering, University of Illinois at Urbana-Champaign.

###### Abstract

Convolutional neural networks (CNNs) have been tremendously successful in
solving imaging inverse problems. To understand their success, an effective
strategy is to construct simpler and mathematically more tractable
convolutional sparse coding (CSC) models that share essential ingredients with
CNNs. Existing CSC methods, however, underperform leading CNNs in challenging
inverse problems. We hypothesize that the performance gap may be attributed in
part to how they process images at different spatial scales: While many CNNs
use multiscale feature representations, existing CSC models mostly rely on
single-scale dictionaries. To close the performance gap, we thus propose a
multiscale convolutional dictionary structure. The proposed dictionary
structure is derived from the U-Net, arguably the most versatile and widely
used CNN for image-to-image learning problems. We show that incorporating the
proposed multiscale dictionary in an otherwise standard CSC framework yields
performance competitive with state-of-the-art CNNs across a range of
challenging inverse problems including CT and MRI reconstruction. Our work
thus demonstrates the effectiveness and scalability of the multiscale CSC
approach in solving challenging inverse problems.

##  I Introduction

Convolutional neural networks (CNNs) obtain state-of-the-art performance in
many image processing tasks. To understand their success, an active line of
recent research reduces CNNs into conceptually simpler and mathematically
better-understood building blocks. Examples of these simplified convolutional
models include _convolutional kernels_ [1, 2, 3], _convolutional scattering
transforms_ [4, 5, 6, 7], and _convolutional sparse coding_ [8, 9, 10]. In
addition to being mathematically tractable, these models have achieved
remarkable empirical success, sometimes matching state-of-the-art CNNs.

This work studies convolutional representations arising from the convolutional
sparse coding (CSC) paradigm, which provides a natural connection between
sparse representation models and CNNs. Indeed, many CNN instances can be
interpreted as optimizing a CSC objective through cascaded layers [8].
Moreover, CSC models compete favorably with state-of-the-art CNNs in several
image processing tasks including denoising, single image super-resolution, and
inpainting [11, 12, 10, 13, 14, 15, 16, 17, 18].

While these emerging results are promising, the successful applications of CSC
in imaging inverse problems are still confined to problems with relatively
simple forward operators, including Gaussian noise addition, blurring, and
uniformly random pixel removal. Common to these forward operators is their
_spatial locality_ – they introduce artifacts that are spatially correlated
only, if at all, within small pixel neighbourhoods. By contrast, a broad range
of imaging inverse problems involve forward models that mix distant parts of
the image and are _highly spatially heterogeneous_ ; examples include the
Radon transform for computed tomography, which computes line integrals along
radiating paths, and the Fourier transform for magnetic resonance imaging,
which computes inner products with globally-supported sinusoids. Working with
these forward models presents different challenges since they introduce
structured noise, such as streak artifacts, with long-range spatial
correlations. We thus ask a natural question: Can CSC models also yield strong
performance on such inverse problems with non-local operators?

To deal with spatially heterogeneous imagery data, one natural strategy is to
employ _multiscale_ dictionaries. Indeed, seminal works have shown that
multiscale dictionaries, either analytical or learned, are advantageous in
representing and processing images [19, 20, 21, 22, 23, 24]. Separating scales
is useful because it gives efficient descriptions of structural correlations
at different distances. Yet, these existing CSC models [10, 25, 26, 27, 28]
mostly employ single-scale dictionaries, whose dictionary atoms all have the
same size. While there exist proposals for multiscale CSC architectures, they
are tailored for specific tasks [29, 30]. In addition, CSC models do away with
flexible skip connections between non-consecutive layers, which are
nonetheless essential for many successful CNNs such as the U-Net and its
variants [31, 32, 33] to fuse features across scales. This challenge of
harnessing multiscale features in the CSC paradigm motivates our work.

To address the challenge, we introduce a multiscale convolutional dictionary
inspired by the highly successful U-Net [31]. We then apply the multiscale
convolutional dictionary to challenging, large-scale inverse problems in
imaging. The main contribution of this paper is twofold:

  * •

We propose a new convolutional dictionary, whose representation incorporates
atoms of different spatial scales. The proposed multiscale dictionary augments
standard, single-scale convolutional dictionaries to exploit the spatially-
heterogeneous properties of images.

  * •

We study the effectiveness of the multiscale convolutional dictionary through
experiments on large-scale datasets. We find that the performance of the
multiscale CSC approach is competitive with leading CNNs on datasets including
two major CT and MRI benchmarks. We additionally show that our model matches
(and slightly improves) the state-of-the-art performance on the deraining task
achieved by a deep neural network [34]. Notably, the single-scale CSC model
performs significantly worse on this task [27].

Overall, our work makes a step forward in closing the performance gap between
end-to-end CNNs and sparsity-driven dictionary models. At a meta level, it
(re)validates the fundamental role of sparsity in representations of images
and imaging operators [35, 36, 20].

The rest of this article is organized as follows. In Section II, we first
briefly review the sparse representation model and its relationship to CNNs.
Section III explains how we incorporate multiscale atoms in a dictionary
model; we also explain how to learn the multiscale dictionary from data under
the task-driven dictionary learning framework. Section IV reports experimental
results on tasks including CT reconstruction and MRI reconstruction.

##  II Background and Related Work

In this section, we briefly review the related work; a summary of notation is
given in Table I.

###  II-A Sparse representation models

Sparse representation has been extensively studied and widely used in imaging
inverse problems [37, 38, 39]. It is motivated by the idea that many signals,
images being a prime example, can be approximated by a linear combination of a
few elements from a suitable overcomplete basis. The sparse representation
framework posits that we can decompose a signal of interest111For simplicity
we write all signals as 1d vectors, but the formulation is valid in any
dimension. $\boldsymbol{z}\in\mathbb{R}^{d}$ as
$\boldsymbol{z}=\boldsymbol{D}\boldsymbol{\alpha}$, where
$\boldsymbol{D}\in\mathbb{R}^{d\times N}$ is an overcomplete dictionary of $N$
atoms ($N>d$) and $\boldsymbol{\alpha}\in\mathbb{R}^{N}$ is a sparse vector
with few non-zero entries. Learning a sparse representation model thus
comprises two sub-problems: (i) given a dictionary $\boldsymbol{D}$, encode
the signal $\boldsymbol{z}$ into a sparse vector $\boldsymbol{\alpha}$ (the
_sparse coding_ problem), and (ii) given a set of signals, learn an
appropriate dictionary $\boldsymbol{D}$ that sparsifies them (the _dictionary
learning_ problem). We briefly review these two problems and show how they are
related to neural network models such as CNNs.

###  II-B The sparse coding problem

The sparse coding problem is often formulated as basis pursuit denoising [40]
or Lasso regression [41]. Most relevant to our work is its formulation with
non-negative constraints on the sparse code $\boldsymbol{\alpha}$:

$\operatorname*{minimize}_{\boldsymbol{\alpha}\geq
0}~{}~{}\frac{1}{2}\left\lVert\boldsymbol{z}-\boldsymbol{D}\boldsymbol{\alpha}\right\rVert_{2}^{2}+\lambda\left\lVert\boldsymbol{\alpha}\right\rVert_{1}.$
(1)

Here, the first term
$\frac{1}{2}\left\lVert\boldsymbol{z}-\boldsymbol{D}\boldsymbol{\alpha}\right\rVert_{2}^{2}$
ensures that the code $\boldsymbol{\alpha}$ yields a faithful representation
of $\boldsymbol{z}$, the second term
$\lambda\left\lVert\boldsymbol{\alpha}\right\rVert_{1}$ controls the sparsity
of the code, and the two terms are balanced by a parameter $\lambda>0$. An
effective solver for the minimization problem (1) is the iterative shrinkage-
thresholding algorithm (ISTA) [42], which executes the following iteration

$\begin{split}\boldsymbol{\alpha}^{[k+1]}&\coloneqq\mathcal{S}(\boldsymbol{\alpha}^{[k]},\boldsymbol{z};\boldsymbol{D},\boldsymbol{\lambda})\\\
&\coloneqq\sigma(\boldsymbol{\alpha}^{[k]}+\eta\boldsymbol{D}^{\top}(\boldsymbol{z}-\boldsymbol{D}\boldsymbol{\alpha}^{[k]})-\eta\boldsymbol{\lambda}),\end{split}$
(2)

where the superscript $[k]$ denotes the iteration number, $\eta$ is a step-
size parameter, $\boldsymbol{\lambda}$ is a vector whose entries are all
$\lambda$, and $\sigma(x)\coloneqq\max(x,0)$ is a component-wise rectifier
function. For simplicity, we use
$\mathcal{S}(\boldsymbol{\alpha},\boldsymbol{z};\boldsymbol{D},\boldsymbol{\lambda})$
to denote one execution of ISTA with measurement $\boldsymbol{z}$, sparse code
$\boldsymbol{\alpha}$, dictionary $\boldsymbol{D}$, and threshold
$\boldsymbol{\lambda}$. The ISTA algorithm is a composition of such
executions; we write $\mathrm{ISTA}_{K}$ for the $K$-fold composition of
$\mathcal{S}$ with itself:

$\begin{split}&\mathrm{ISTA}_{K}(\boldsymbol{z};\boldsymbol{D},\boldsymbol{\lambda})\\\
&\coloneqq\underbrace{\Big{(}\mathcal{S}(\cdot,\boldsymbol{z};\boldsymbol{D},\boldsymbol{\lambda})\circ\cdots\circ\mathcal{S}(\cdot,\boldsymbol{z};\boldsymbol{D},\boldsymbol{\lambda})\Big{)}}_{K\text{
times}}(\boldsymbol{\alpha}^{[0]}),\end{split}$ (3)

where $\boldsymbol{\alpha}^{[0]}$ is the initial sparse code; throughout this
work, this initial code $\boldsymbol{\alpha}^{[0]}$ is assumed to contain zero
in all entries. We emphasize that ISTA is a nonlinear transform of its input
$\boldsymbol{z}$.

Table I: The notations used in this paper. $d$ | the dimension of an image
---|---
$N$ | the dimension of a sparse code
$M$ | the number of training samples
$K$ | the number of sparse coding steps
$\sigma(\cdot)$ | the ReLU non-linearlity
$\boldsymbol{z}$ | a noisy image to be processed
$\boldsymbol{D}$ | an overcomplete dictionary
$\boldsymbol{x}^{\dagger}$ | a ground-truth image
$\widehat{\boldsymbol{x}}$ | a predicted image
$\boldsymbol{\alpha}$ | a sparse code
$\boldsymbol{\lambda}$ | the thresholding parameters of ISTA

###  II-C The task-driven dictionary learning problem

We now briefly recall the task-driven dictionary learning framework [43].
Consider a supervised learning setting, in which we aim to identify a
parametric function that associates each input $\boldsymbol{z}$ (e.g., a
corrupted image) with its target $\boldsymbol{x}^{\dagger}$ (e.g., a clean
image) for all
$(\boldsymbol{z},\boldsymbol{x}^{\dagger})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$
drawn from some joint distribution. In the task-driven framework, we proceed
by first representing the signal $\boldsymbol{z}$ by a sparse code
$\boldsymbol{\alpha}_{\boldsymbol{z}}$ with respect to a dictionary
$\boldsymbol{D}$. One way to achieve this is to let

$\boldsymbol{\alpha}_{\boldsymbol{z}}\coloneqq\operatorname*{arg\,min}_{\boldsymbol{\alpha}\geq
0}\frac{1}{2}\left\lVert\boldsymbol{z}-\boldsymbol{D}\boldsymbol{\alpha}\right\rVert_{2}^{2}+\lambda\left\lVert\boldsymbol{\alpha}\right\rVert_{1},$
(4)

which can be approximated by $K$ iterations of ISTA as in (3). Next, we
approximate the desired target $\boldsymbol{x}^{\dagger}$ using the sparse
code $\boldsymbol{\alpha}_{\boldsymbol{z}}$ through a regression model
$f(\cdot,\boldsymbol{w})$ with learnable parameter $\boldsymbol{w}$. For
instance, $f(\cdot,\boldsymbol{w})$ could be a linear regression model with
weights and biases $\boldsymbol{w}$. The model output
$f(\boldsymbol{\alpha}_{\boldsymbol{z}},\boldsymbol{w})$ thus depends on the
regression model parameters $\boldsymbol{w}$ as well as the sparse code
$\boldsymbol{\alpha}_{\boldsymbol{z}}$, which in turn depends on the
dictionary $\boldsymbol{D}$ through the ISTA iterations. In this way, the
regression parameters $\boldsymbol{w}$ and dictionary $\boldsymbol{D}$ can be
_jointly_ optimized, for instance, with respect to the quadratic loss
objective evaluated on a dataset of $M$ input-target pairs
$\\{(\boldsymbol{z}_{i},\boldsymbol{x}^{\dagger}_{i})\\}_{i=1}^{M}$:

$\operatorname*{minimize}_{\boldsymbol{w},\boldsymbol{D}}~{}~{}\frac{1}{2M}\sum_{i=1}^{M}\|f(\boldsymbol{\alpha}_{\boldsymbol{z}_{i}},\boldsymbol{w})-\boldsymbol{x}^{\dagger}_{i}\|_{2}^{2}.$
(5)

Importantly, the task-driven objective in (5) implies that the dictionary
$\boldsymbol{D}$ is optimized to solve the supervised learning task and not
just to sparsely represent data.

Figure 1: Schematic illustration of the U-Net (left panel) and the dictionary
model considered in this work (right panel). (a): The U-Net processes input
images using convolution, scale separations, and skip connections in
conjunction with ReLU non-linearities and batch-normalization (BN) modules
indicated by colored arrows. The multi-channel feature maps produced by these
operations are illustrated as boxes with the channel numbers indicated on the
top of these boxes. Dark grey boxes indicate the feature maps produced by the
encoding branch of the U-Net, which are sent to the decoding branch either
through channel-wise concatenation (“skip connection”) or through the
bottleneck layer. (b): The dictionary considered in this work is a
simplification of the decoder branch of the U-Net: We retain convolution and
multiscale representation from the decoder branch of the U-Net but remove all
non-linearities, batch-normalization, and additive biases; additionally, we
remove a convolution at each spatial resolution level and halve the number of
convolutional channels for all convolutions. Grey boxes indicate the
multiscale sparse code
$\boldsymbol{\alpha}=(\boldsymbol{\alpha}_{0},\ldots,\boldsymbol{\alpha}_{4})$
that the dictionary takes as input. Dashed boxes indicate the position that
each $\boldsymbol{\alpha}_{i}$ feed into the dictionary. (c): The proposed as
a computational graph that uses multiscale dictionaries
$\boldsymbol{D}_{\textrm{enc}}$, $\widetilde{\boldsymbol{D}}_{\textrm{enc}}$,
and $\boldsymbol{D}_{\textrm{dec}}$; although each dictionary is linear, the
computational graph is nonlinear due to the thresholding operator.

###  II-D Convolutional sparse coding

Our work is inspired by the convolutional sparse coding (CSC) model [44, 45,
46, 47, 8], which bridges deep CNNs and sparse representation models.
Concretely, Papyan et al. [8] noticed that if the dictionary $\boldsymbol{D}$
has a convolutional structure and if the sparse code $\boldsymbol{\alpha}$ is
assumed to be non-negative, a single iteration of ISTA with
$\boldsymbol{\alpha}^{[0]}$ initialized as a zero vector and step-size
$\eta=1$ is equivalent to the forward pass of a single-layer convolutional
network

$\boldsymbol{\alpha}=\sigma(\boldsymbol{D}^{\top}\boldsymbol{z}+\boldsymbol{b}),$
(6)

where $\boldsymbol{b}$ is a vector whose components are $-\lambda$ (cf.
Equation (2)). This single-layer formulation can be extended to characterize a
deep CNN of multiple layers. Specifically, the forward-pass of a deep CNN of
$L$-layers can be interpreted to approximate the sparse codes
$\boldsymbol{\alpha}_{1},\cdots,\boldsymbol{\alpha}_{L}$ sequentially with
respect to different dictionaries
$\boldsymbol{D}_{1},\cdots,\boldsymbol{D}_{L}$; the back-propagation pass is
interpreted as an update to these dictionaries
$\\{\boldsymbol{D}_{i}\\}_{i=1}^{L}$ in a task-driven way.

###  II-E CNNs for solving inverse problems

Deep CNNs achieve state-of-the-art performance in many image processing tasks
[48, 49, 50, 51]. In particular, the U-Net [31] and its variants [32, 33, 52]
are among the most extensively used CNN architectures in solving image-to-
image learning tasks. U-Nets represent images via multiscale features computed
from measurements using an encoding (or downsampling) branch and a synthesized
into an estimated image in a decoding (or upsampling) branch (Figure 1a). In
the downsampling branch, the spatial resolutions of feature maps are reduced
while the number of feature maps is increased; in the decoding branch, these
features are recombined with previous high-resolution features via channel
concatenation (“skip connections”) and convolution. Heuristically, low-
resolution feature maps of a U-Net capture large-scale image properties,
whereas the high-resolution feature maps capture more fine-grained image
properties [52]. In a related line of work, Ye et al. [53, 54, 55] proposed to
use the framelets formalism [56] to study aspects of U-Net-like encoder-
decoder CNNs. A key observation they make is that a U-Net model is closely
related to convolutional framelets whose frame basis selection depends non-
linearly on input data.

##  III CSC with multiscale dictionaries

The structure of a convolutional dictionary is crucial to a CSC model since
the dictionary atoms characterize the signals that can be represented
sparsely. In the existing formulation of CSC, atoms of a convolutional
dictionary have a single scale, in the sense that they all share the same
spatial shape. However, many image classes and imaging artifacts exhibit
structured correlations over multiple scales. To exploit these correlations in
imaging inverse problems, we construct multiscale convolutional dictionaries.

Our construction is based on the U-Net reviewed in Section II. Indeed, the
tremendous success of U-Nets has in part been attributed to their ability to
represent images at multiple scales [33, 55], which is achieved by using up-
and downsampling operations together with skip connections as in Figure 1a.
Another property of the U-Net is its shared parameters across scales: Low-
resolution features (the grey boxes at the bottom of Figure 1a) and high-
resolution features (the grey boxes at the top of Figure 1a) undergo an
overlapping synthesizing path parameterized by shared weights. This weight-
sharing strategy has not been employed by existing proposals for multiscale
CSC dictionaries [29, 30]. In what follows, we describe the construction
process of a linear dictionary inspired by and closely following the standard
U-Net.

###  III-A Encoder–decoder dictionaries

We denote the encoding branch of the U-Net by
$f_{\mathrm{enc}}(\cdot,\boldsymbol{\theta}_{\mathrm{enc}}):\mathbb{R}^{d}\to\mathbb{R}^{N}$
with parameters $\boldsymbol{\theta}_{\mathrm{enc}}$; the encoding branch maps
the input $\boldsymbol{z}\in\mathbb{R}^{d}$ to convolutional feature maps
$\boldsymbol{\alpha}_{\boldsymbol{z}}=f_{\mathrm{enc}}(\boldsymbol{z},\boldsymbol{\theta}_{\mathrm{enc}})\in\mathbb{R}^{N}$,
illustrated as the dark grey boxes in Figure 1a. Note that, for a U-Net, the
intermediate feature map dimension $N$ (number of scalar coefficients in
$\boldsymbol{\alpha}$) is typically much greater than the image dimension $d$.
These feature maps are then fed into the decoding branch of the U-Net either
through skip connections or through the bottleneck layer. To describe this
process, we write the decoding branch of the U-Net as a function
$f_{\mathrm{dec}}(\cdot,\boldsymbol{\theta}_{\mathrm{dec}}):\mathbb{R}^{N}\to\mathbb{R}^{d}$
with parameters $\boldsymbol{\theta}_{\mathrm{dec}}$. That is, the function
$f_{\mathrm{dec}}(\cdot,\boldsymbol{\theta}_{\mathrm{dec}})$ takes the
convolutional feature maps produced by the encoding branch and transforms them
to produce the model output. We can thus write the output produced by a U-Net
as

$\widehat{\boldsymbol{z}}\coloneqq
f_{\mathrm{dec}}(\boldsymbol{\alpha}_{\boldsymbol{z}},\boldsymbol{\theta}_{\mathrm{dec}})=f_{\mathrm{dec}}\Big{(}f_{\mathrm{enc}}(\boldsymbol{z},\boldsymbol{\theta}_{\mathrm{enc}}),\boldsymbol{\theta}_{\mathrm{dec}}\Big{)}.$

We now focus on the image synthesis process of the U-Net, described by the
decoding function
$f_{\mathrm{dec}}(\cdot,\boldsymbol{\theta}_{\mathrm{dec}})$. This function
synthesizes convolutional feature maps at different spatial scales through
skip connections and upsampling. As such, the decoding branch of the U-Net
approximates an image $\boldsymbol{x}^{\dagger}\in\mathbb{R}^{d}$ using
multiscale feature maps
$\boldsymbol{\alpha}_{\boldsymbol{z}}\in\mathbb{R}^{N}$ of a much higher
dimension, so that $\boldsymbol{x}^{\dagger}\approx
f_{\mathrm{dec}}(\boldsymbol{\alpha}_{\boldsymbol{z}},\boldsymbol{\theta}_{\mathrm{dec}})$.
Conceptually, this representation is similar to the sparse and overcomplete
representation in a dictionary, except that the U-Net decoder is non-linear.

To construct a multiscale dictionary, we thus consider a stripped-down version
of the image synthesis process of U-Net by removing all non-linearities, batch
normalization, and additive biases from the function
$f_{\mathrm{dec}}(\cdot,\boldsymbol{\theta}_{\mathrm{dec}})$, as shown in
Figure 1b; to further simplify the architecture, at each spatial scale, we
additionally remove a convolution and halve the number of convolutional
channels for all convolutions. The resulting function is then simply a
_linear_ transformation

$\boldsymbol{\alpha}\coloneqq(\boldsymbol{\alpha}_{0},\ldots,\boldsymbol{\alpha}_{4})\mapsto\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha},$
(7)

where $\boldsymbol{\alpha}_{0},\ldots,\boldsymbol{\alpha}_{4}$ are sparse code
having different resolutions (visualized as the grey boxes in Figure 1b). This
dictionary shares the essential ingredients of convolution, multiscale
representation, and skip connections with the U-Net decoding branch and
therefore we refer to it as the decoder dictionary. A precise definition of
the decoder dictionary $\boldsymbol{D}_{\textrm{dec}}$ through convolution and
upsampling is provided in Appendix A.

###  III-B The dictionary-based sparsity prior

With a given decoder dictionary $\boldsymbol{D}_{\textrm{dec}}$ to describe
the image synthesis process, we next consider how to infer an associated
sparse code $\boldsymbol{\alpha}$, so that
$\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}$ is a good approximation of
the image we wish to model. In a supervised learning setting where the input
image $\boldsymbol{z}$ is given, it is natural to interpret
$\boldsymbol{\alpha}$ as an encoded representation of $\boldsymbol{z}$. Since
the encoding must produce a coefficient vector whose structure is compatible
with $\boldsymbol{\alpha}$, we endow an encoder dictionary
$\boldsymbol{D}_{\textrm{enc}}\in\mathbb{R}^{d\times N}$ with the same
structure of $\boldsymbol{D}_{\textrm{dec}}$ albeit with a different set of
atoms. This setup is analogous to U-Net’s encoding and decoding branches: the
encoder and decoder dictionaries $\boldsymbol{D}_{\textrm{enc}}$ and
$\boldsymbol{D}_{\textrm{dec}}$ are employed to process input signals and
produce output signals, respectively. The sparse code
$\boldsymbol{\alpha}_{\boldsymbol{z}}$ induced by an input $\boldsymbol{z}$
and the encoder dictionary $\boldsymbol{D}_{\textrm{enc}}$ then facilitate the
subsequent task for approximating the ground-truth image $\boldsymbol{x}$:

$\boldsymbol{z}\xrightarrow{\hskip 2.84526pt\textrm{ Sparse coding with
}\boldsymbol{D}_{\textrm{enc}}\hskip
2.84526pt}\boldsymbol{\alpha}_{\boldsymbol{z}}\xrightarrow{\hskip
2.84526pt\textrm{ Synthesis with }\boldsymbol{D}_{\textrm{dec}}\hskip
2.84526pt}\widehat{\boldsymbol{x}}\coloneqq\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}_{\boldsymbol{z}}.$
(8)

In what follows, we derive a supervised learning method that turns each
$\boldsymbol{z}$ into a prediction $\widehat{\boldsymbol{x}}$ using encoder
and decoder dictionaries.

###  III-C The task-driven dictionary learning objective

Under the task-driven framework introduced in Section II, we formulate a
supervised learning problem via sparse coding and dictionary learning. We
consider the following minimization problem over a dataset of $M$ input-target
pairs $\\{(\boldsymbol{z}_{i},\boldsymbol{x}^{\dagger}_{i})\\}_{i=1}^{M}$:

$\displaystyle\operatorname*{minimize}_{\\{\boldsymbol{D}_{\textrm{enc}},\boldsymbol{D}_{\textrm{dec}}\\},\boldsymbol{\lambda}>0}$
$\displaystyle\frac{1}{2M}\sum_{i=1}^{M}\|\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}_{\boldsymbol{z}_{i}}-\boldsymbol{x}^{\dagger}_{i}\|_{2}^{2}$
(9) where
$\displaystyle\boldsymbol{\alpha}_{\boldsymbol{z}_{i}}\coloneqq\mathrm{ISTA}_{K}(\boldsymbol{z}_{i};\boldsymbol{D}_{\textrm{enc}},\boldsymbol{\lambda}).$

The objective in (9) penalizes the discrepancy between the ground-truth signal
$\boldsymbol{x}^{\dagger}$ and the model prediction
$\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}_{\boldsymbol{z}}$, where the
latter is a signal synthesized from a sparse code
$\boldsymbol{\alpha}_{\boldsymbol{z}}$ via the decoder dictionary
$\boldsymbol{D}_{\textrm{dec}}$; the code
$\boldsymbol{\alpha}_{\boldsymbol{z}}$ is a sparse representation of the input
image $\boldsymbol{z}$ with respect to the encoder dictionary
$\boldsymbol{D}_{\textrm{enc}}$ by unrolling a fixed number $K$ of ISTA
iterations. The sparsity-controlling parameter $\boldsymbol{\lambda}$ is
multi-dimensional, weighting codes component-wise. The intuition behind this
choice is that the different convolutional features, especially those at
different scales, should be thresholded differently. The sparse code
$\boldsymbol{\alpha}$, illustrated as the grey boxes in Figure 1b, is a
collection of multi-dimensional tensors, each corresponds a spatial scale.

The task driven objective (9) defines a computational graph that transforms an
input image $\boldsymbol{z}$ into a prediction
$\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}_{\boldsymbol{z}}$. We term
this computational graph MUSC, since it involves multiscale U-Net-like sparse
coding. We note the MUSC is an instance of _optimization-driven networks_ [26]
derived by unrolling an optimization algorithm. It incorporates two modules
with meaningful objectives, one implementing sparse coding and the other
dictionary-based synthesis. This composition is arguably conceptually more
interpretable than end-to-end layerwise composition of deep networks.

While a traditional compressed sensing approach uses _a single dictionary_ for
reconstruction, our approach uses _two dictionaries_
$\boldsymbol{D}_{\textrm{enc}}$ and $\boldsymbol{D}_{\textrm{dec}}$ in the
task-driven learning objective (9). This discrepancy is due to different
assumptions in measurement-to-image reconstruction (the compressed sensing
approach) and image-to-image reconstruction (our approach). Consider an
inverse problem with a forward operator $\boldsymbol{A}$, a unknown ground-
truth signal $\boldsymbol{x}^{\dagger}$, and measurements
$\boldsymbol{y}\coloneqq\boldsymbol{A}\boldsymbol{x}^{\dagger}$; in CT
reconstruction, $\boldsymbol{A}$ is the Radon transform and $\boldsymbol{y}$
is the measured sinogram. The compressed sensing approach estimates
$\boldsymbol{x}^{\dagger}$ as $\boldsymbol{D}\boldsymbol{\alpha}^{\ast}$ for
some dictionary $\boldsymbol{D}$, where

$\boldsymbol{\alpha}^{\ast}=\operatorname*{arg\,min}_{\boldsymbol{\alpha}}\left\lVert\boldsymbol{A}\boldsymbol{D}\boldsymbol{\alpha}-\boldsymbol{y}\right\rVert_{2}+\lambda\left\lVert\boldsymbol{\alpha}\right\rVert_{1}$
(10)

is the inferred sparse code based on the dictionary $\boldsymbol{D}$. Note
that (10) and the synthesis $\boldsymbol{D}\boldsymbol{\alpha}^{*}$ require
only a single dictionary $\boldsymbol{D}$. However, this approach assumes that
we know the measurements $\boldsymbol{y}$ and the forward operator
$\boldsymbol{A}$.

If we were to apply a single dictionary
$\boldsymbol{D}\coloneqq\boldsymbol{D}_{\textrm{enc}}=\boldsymbol{D}_{\textrm{dec}}$
in our image-to-image learning approach in (9), we would find a sparse code
$\boldsymbol{\alpha}$ such that
$\boldsymbol{D}\boldsymbol{\alpha}\approx{\boldsymbol{x}^{\dagger}}$ and
$\boldsymbol{D}\boldsymbol{\alpha}\approx{\boldsymbol{A}^{+}\boldsymbol{A}\boldsymbol{x}^{\dagger}}$.
This is difficult when $\boldsymbol{A}^{+}\boldsymbol{A}$ significantly
differs from the identity operator as in the case of highly ill-posed
problems. On the other hand, using two dictionaries
$\boldsymbol{D}_{\textrm{enc}}$ and $\boldsymbol{D}_{\textrm{dec}}$ in (9)
requires finding a sparse code $\boldsymbol{\alpha}$ such that
$\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}\approx{\boldsymbol{x}^{\dagger}}$
and
$\boldsymbol{D}_{\textrm{enc}}\boldsymbol{\alpha}\approx{\boldsymbol{A}^{+}\boldsymbol{A}\boldsymbol{x}^{\dagger}}$,
a formulation that is more flexible when $\boldsymbol{A}^{+}\boldsymbol{A}$
substantially differs from the identity. Experiments in Section IV-E confirm
that allowing $\boldsymbol{D}_{\textrm{enc}}\neq\boldsymbol{D}_{\textrm{dec}}$
yields better performance. We note that our approach is morally related to
setting $\boldsymbol{D}_{\textrm{enc}}=\boldsymbol{A}\boldsymbol{D}$ in (10),
but since we do not know $\boldsymbol{A}$ we have to learn
$\boldsymbol{D}_{\textrm{enc}}$ from samples together with
$\boldsymbol{D}_{\textrm{dec}}$. Such a learned encoder dictionary captures
information about $\boldsymbol{A}$, entangled with information about the data
distribution.

###  III-D Relaxation on dictionaries

We now describe computational techniques that stabilize the gradient-descent-
based dictionary learning of MUSC. Following earlier work [57, 10, 26, 25, 6],
we untie the encoder dictionary from its adjoint during dictionary update.
That is, we replace the execution in (2) by

$\widetilde{\mathcal{S}}(\boldsymbol{\alpha},\boldsymbol{z};\boldsymbol{D}_{\textrm{enc}},\widetilde{\boldsymbol{D}}_{\textrm{enc}},\boldsymbol{\lambda})\coloneqq\sigma(\boldsymbol{\alpha}+\eta\widetilde{\boldsymbol{D}}_{\textrm{enc}}^{\top}(\boldsymbol{z}-\boldsymbol{D}_{\textrm{enc}}\boldsymbol{\alpha})-\eta\boldsymbol{\lambda}),$
(11)

where the dictionary $\widetilde{\boldsymbol{D}}_{\textrm{enc}}$ is
initialized to be identical to $\boldsymbol{D}_{\textrm{enc}}$ but is allowed
to evolve independently during training. Even though the theoretical effects
of this relaxation remain unclear, the dictionary
$\widetilde{\boldsymbol{D}}_{\textrm{enc}}$ can be interpreted as a learned
preconditioner that accelerates training [26, 25]; see also the investigation
in [58, 59, 6]. The _learned ISTA_ (LISTA) algorithm [57] corresponding to
(11) is written as

$\begin{split}&\mathrm{LISTA}_{K}(\boldsymbol{z};\boldsymbol{D}_{\textrm{enc}},\widetilde{\boldsymbol{D}}_{\textrm{enc}},\boldsymbol{\lambda})\\\
&\small\coloneqq\underbrace{\Big{(}\widetilde{\mathcal{S}}(\cdot,\boldsymbol{z};\boldsymbol{D}_{\textrm{enc}},\widetilde{\boldsymbol{D}}_{\textrm{enc}},\boldsymbol{\lambda}_{K})\circ\cdots\circ\widetilde{\mathcal{S}}(\cdot,\boldsymbol{z};\boldsymbol{D}_{\textrm{enc}},\widetilde{\boldsymbol{D}}_{\textrm{enc}},\boldsymbol{\lambda}_{1})\Big{)}}_{K\text{
times}}(\boldsymbol{\alpha}^{[0]}),\end{split}$ (12)

where $\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{K}$ are the soft-
thresholding parameters for each ISTA execution. Note that, in (12), the soft-
thresholding parameters $\\{\boldsymbol{\lambda}_{i}\\}_{i=1}^{K}$ depend on
the execution step. As shown in [6], incorporating step-dependent soft-
thresholding parameters can be beneficial. While [6] uses a homotopy
continuation strategy to adjust these parameters we treat them as learnable
parameters for simplicity. Taking these considerations into account, we define
a new regression loss:

$\displaystyle\mathcal{L}(\boldsymbol{D}_{\textrm{enc}},\boldsymbol{D}_{\textrm{dec}},\boldsymbol{\lambda})\coloneqq\frac{1}{2M}\sum_{i=1}^{M}\|\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}_{\boldsymbol{z}_{i}}-\boldsymbol{x}^{\dagger}_{i}\|_{2}^{2},$
(13)
$\displaystyle\text{where}\quad\boldsymbol{\alpha}_{\boldsymbol{z}_{i}}=\mathrm{LISTA}_{K}(\boldsymbol{z}_{i};\boldsymbol{D}_{\textrm{enc}},\widetilde{\boldsymbol{D}}_{\textrm{enc}},\boldsymbol{\lambda}).$

Unless mentioned otherwise, we use the loss (13) to train MUSC throughout our
paper. In Section IV-E, we compare the performance of trained model using (13)
and (9).

###  III-E Training the MUSC

Training the MUSC entails the following three steps:

  1. 1.

Dictionary initialization: We randomly initialize the dictionary
$\boldsymbol{D}_{\textrm{enc}}$ and initialize
$\boldsymbol{D}_{\textrm{dec}}$, and
$\widetilde{\boldsymbol{D}}_{\textrm{enc}}$ as identical copies of
$\boldsymbol{D}_{\textrm{enc}}$.

  2. 2.

Model forward pass: For each input image $\boldsymbol{z}_{i}$, we evaluate the
model prediction
$\boldsymbol{D}_{\textrm{dec}}\boldsymbol{\alpha}_{\boldsymbol{z}_{i}}$ as in
Equation (13). For ISTA executions, we initialize all sparse code
$\boldsymbol{\alpha}_{\boldsymbol{z}}^{[0]}$ as a collection of all-zero
tensors; the ISTA step-size parameter $\eta$ is initialized as the inverse of
the dominant eigenvalue of the matrix
$\boldsymbol{D}_{\textrm{enc}}^{\top}\boldsymbol{D}_{\textrm{enc}}$, which can
be approximated using by power iteration (Appendix C).

  3. 3.

Task-driven dictionary learning: For a mini-batch of input-target pairs, solve
the optimization problem in (13) with gradient descent.

Figure 2: Reconstructions of a test sample from the Rain12 dataset. | Rain12 | Rain100L
---|---|---
| PSNR | SSIM | PSNR | SSIM
LP [60] | 32.02 | 0.91 | 29.11 | 0.88
DDN [61] | 31.78 | 0.90 | 32.16 | 0.94
DSC [62] | 30.02 | 0.87 | 24.16 | 0.87
JORDER [63] | 33.92 | 0.95 | 36.61 | 0.97
PReNet [34] |  grey 36.69 | 0.96 |  grey 37.10 | 0.98
LGM (single-scale CSC; [27]) | 35.46 | 0.95 | 34.07 | 0.96
MUSC (multiscale CSC; ours) |  36.77 |  0.96 |  37.25 |  0.98
Table II: Performance on the deraining test set. Boldface indicates the best
performance; second-best results are highlighted in grey. All results are
collected from [34] and [27] except MUSC.

##  IV Experiments

We report the performance of MUSC on deraining, CT reconstruction, and MRI
reconstruction tasks. The motivations for choosing these tasks are as follows.
First, we note that single-scale CSC models have recently been applied to the
deraining task, achieving performance slightly worse than state-of-the-art
deep networks [27]; we thus aim to test the capability of our multiscale
approach on the same task. We additionally choose CT and MRI reconstruction
tasks as there exist challenging, large-scale, and up-to-date benchmark
datasets for these tasks. Two such datasets that we use are the LoDoPaB-CT
[64] and the fastMRI [65]. An additional strength of these two datasets is
that the model evaluation process is carefully controlled: The evaluation on
the challenge fold (for LoDoPaB-CT) or the test fold (for fastMRI) is
restricted through an online submission portal with the ground truth hidden
from the public. As a result, overfitting to these evaluation folds is
difficult and quantitative comparisons are transparent.

Throughout this section, we use the MUSC architecture whose encoder and
decoder dictionaries are displayed in Figure 1b and mathematically defined in
Appendix A. Hyper-parameter choices for the experiments are provided in
Appendix D. For each task, we use well-known CNN models as baselines. We note
that, for the CT and MRI reconstruction tasks, there are two major approaches
to employ CNNs. In the first, _model-based_ approach, one applies neural
networks on raw measurement data (sinogram data in CT and k-space data in MRI)
by embedding a task-dependent forward operator (the Radon transform for CT and
the Fourier transform for MRI) into multiple layers or iterations of the
network. Learning methods of this approach can be highly performant at the
cost of being computationally expensive, especially during training, since one
needs to apply the forward operator (and the adjoint of its derivative)
repeatedly [49]. In the second, _model-free_ approach, the (pseudoinverse of
the) forward operator is used at most _once_ during data preprocessing and is
never used during subsequent supervised training. These preprocessed images
contain imaging artifacts. During supervised learning, one applies a CNN
directly on these preprocessed images. The proposed MUSC is in this sense a
model-free approach and we compare it to model-free baselines. We note that in
this case one does not need to know the forward operator at all. The leading
model-free baseline CNN methods in this approach are typically U-Net variants
tuned to the task at hand. For a more thorough comparison, we also implemented
the original U-Net architecture proposed in [31] (schematically illustrated in
Figure 1a) in these tasks as additional baselines.

While model-free approaches perform somewhat worse than model-based ones, our
purpose here is to show that a general-purpose multiscale convolutional model
can perform as well as convolutional neural networks ceteris paribus, rather
than to propose state-of-the-art reconstruction algorithms for specific
problems. This general-purpose approach further allows us to tackle structured
denoising problems such as deraining where the forward operator is simply the
identity.

Figure 3: Reconstructions of a test sample from the LoDoPaB-CT dataset. | Number of parameters | PSNR | PSNR-FR | SSIM | SSIM-FR
---|---|---|---|---|---
FBP | - | 30.19 $\pm$ 2.55 | 34.46 $\pm$ 2.18 | 0.727 $\pm$ 0.127 | 0.836 $\pm$ 0.085
TV | - | 33.36 $\pm$ 2.74 | 37.63 $\pm$ 2.70 | 0.830 $\pm$ 0.121 | 0.903 $\pm$ 0.082
CINN | 6.43M | 35.54 $\pm$ 3.51 | 39.81 $\pm$ 3.48 | 0.854 $\pm$ 0.122 | 0.919 $\pm$ 0.081
U-Net++ | 9.17M | 35.37 $\pm$ 3.36 | 39.64 $\pm$ 3.40 | 0.861 $\pm$ 0.119 | 0.923 $\pm$ 0.080
MS-D-CNN | 181.31K | 35.85 $\pm$ 3.60 | 40.12 $\pm$ 3.56 | 0.858 $\pm$ 0.122 | 0.921 $\pm$ 0.082
U-Net | 31.04M | 35.87 $\pm$ 3.59 | 40.14 $\pm$ 3.57 | 0.859 $\pm$ 0.121 | 0.922 $\pm$ 0.081
LoDoPaB U-Net | 613.32K |  grey 36.00 $\pm$ 3.63 |  grey 40.28 $\pm$ 3.59 |  grey 0.862 $\pm$ 0.119 |  grey 0.923 $\pm$ 0.079
MUSC (ours) | 13.87M |  36.08 $\pm$ 3.68 |  40.35 $\pm$ 3.64 |  0.863 $\pm$ 0.119 |  0.924 $\pm$ $0.080$
Table III: Performance on the LoDoPaB-CT challenge set. All values are taken
from the official challenge leaderboard.

###  IV-A Deraining

Image deraining aims to remove rain streaks from an image. Formally, a rainy
image $\boldsymbol{z}$ is expressed as
$\boldsymbol{z}=\boldsymbol{x}^{\dagger}+\boldsymbol{s}$, where
$\boldsymbol{x}^{\dagger}$ is a clean image and $\boldsymbol{s}$ is the rain
streaks component. The goal is to reconstruct the clean image
$\boldsymbol{x}^{\dagger}$ based on the rainy image $\boldsymbol{z}$.
Recently, single-scale CSC models have been applied to the draining task [27].
Despite theoretical progress, these single-scale CSC models still fall short
competing with leading deep learning models, as remarked in [27]. In this
section, we demonstrate that our multiscale CSC model closes this performance
gap.

Throughout this subsection, we follow the experiment setup of [27]. We use 200
clean and rainy image pairs as the training dataset. A rainy image is created
by adding synthesized rain streaks to its clean version. We use two test sets,
Rain12 [60] and Rain100L [63], to benchmark our results. Similar to [27], we
train our model to restore rain streaks based on rainy images; a derained
image is then the difference between a rainy image and the restored rain
streaks. To be consistent with [27, 63, 34], the evaluation result is
calculated after transforming the image into the luma component in the YCbCr
domain using the software provided by [34]. Additional details of the
experiment are provided in Appendix D.

We report the reconstruction performance in Table II and visualize the
reconstruction results in Figure 2. Our multiscale convolutional dictionary
approach matches or outperforms baseline methods. Notably, it improves upon
the LGM method (the single-scale CSC approach of [27]) by a non-trivial
margin.

Figure 4: Reconstructions of a test sample from the fastMRI single-coil knee dataset. | Number of parameters | NMSE | PSNR | SSIM
---|---|---|---|---
|  | PD | PDFS | PD | PDFS | PD | PDFS
TV | - | 0.0287 | 0.0900 | 31.4 | 27.7 | 0.645 | 0.494
U-Net | 31.04M | 0.0161 | 0.0538 | 33.8 | 29.9 | 0.809 |  grey 0.631
fastMRI U-Net-256 | 214.16M | 0.0154 | 0.0525 | 34.0 | 30.0 | 0.815 | 0.636
MUSC (ours) | 13.87M |  grey 0.0156 |  grey 0.0526 | 34.0 |  30.0 |  grey 0.811 |  grey 0.631
Table IV:  Performance on fastMRI single-coil knee validation data. Results are collected from [65] except U-Net and MUSC. The fastMRI U-Net-32 model refers to a U-Net variant defined in [65] whose output after the first convolution has 32 channels. Other models are defined similarly. PDFS and PD correspond to two MRI acquisition protocols with fat suppression (PDFS) and without fat suppression (PD) [65]. | Number of parameters | NMSE | PSNR | SSIM
---|---|---|---|---
|  | PD | PDFS | PD | PDFS | PD | PDFS
Zero-filled IFFT | - | 0.0266 | 0.0597 | 32.0 | 29.2 | 0.7765 | 0.6045
TV | - | 0.0221 | 0.0716 | 33.1 | 28.5 | 0.7036 | 0.5096
U-Net | 31.04M | 0.0130 | 0.0430 | 35.3 | 30.6 |  $0$.8513 | 0.6626
fastMRI U-Net-256 | 214.16M | 0.0115 | 0.0414 | 35.9 | 30.8 | 0.8618 | 0.6680
MUSC (ours) | 13.87M |  grey 0.0126 |  grey 0.0421 |  grey 35.5 |  grey 30.7 |  grey 0.8537 |  grey 0.6633
Table V: Performance on single-coil knee test data. Results are collected from
the fastMRI public leaderboard.

###  IV-B CT reconstruction

Computed tomography (CT) aims to recover images from their sparse-view
sinograms. We use the LoDoPaB-CT dataset [64] to benchmark our results. This
dataset contains more than 40000 pairs of human chest CT images and their
simulated low photon count measurements. The ground truth images of this
dataset are human chest CT scans corresponding to the LIDC/IDRI dataset [66],
cropped to 362 × 362 pixels. The low-dose projections are simulated using the
default setting of [64].

To train our MUSC, we use the default dataset split as recommended in [64]:
The dataset is divided into 35820 training samples, 3522 validation samples,
3553 test samples, and 3678 challenge samples. Here, the ground-truth samples
from the challenge dataset are hidden from the public; the evaluation on this
fold is performed through the online submission system of the LoDoPaB-CT
challenge222https://lodopab.grand-challenge.org/challenge/.

We compare the reconstruction performance of MUSCs with five modern CNN
baselines, namely CINN [67], U-Net++ [68], MS-D-CNN [69], U-Net [31], and
LoDoPaB U-Net [64]; the LoDoPaB U-Net refers to a modified U-Net architecture
tailored to the LoDoPaB-CT task. Figure 3 shows the reconstruction results of
a test sample. In Table III, we quantitatively compare MUSC with two classic
methods (FBP and TV) together with five CNN baseline methods mentioned above.
As shown in Table III, MUSC outperforms all baselines. The metrics PSNR and
PSNR-FR are taken from [49]: For a ground-truth signal
$\boldsymbol{x}^{\dagger}$ and its approximation $\widehat{\boldsymbol{x}}$,
we define

$\operatorname{PSNR}\left(\widehat{\boldsymbol{x}},\boldsymbol{x}^{\dagger}\right)\coloneqq
10\log_{10}\left(\frac{\max(\boldsymbol{x}^{\dagger})-\min(\boldsymbol{x}^{\dagger})}{\operatorname{MSE}\left(\widehat{\boldsymbol{x}},\boldsymbol{x}^{\dagger}\right)}\right),$
$\operatorname{PSNR-
FR}\left(\widehat{\boldsymbol{x}},\boldsymbol{x}^{\dagger}\right)\coloneqq
10\log_{10}\left(\frac{1}{\operatorname{MSE}\left(\widehat{\boldsymbol{x}},\boldsymbol{x}^{\dagger}\right)}\right).$

Figure 5: Ablation study on how different choices of model components affect
the overall performance. Left panel shows the PSNR (evaluated on validation
samples) of six trained models as the training progresses; right panel shows
the configuration of each trained model, where Case 2 corresponds to the usual
setting used in other sections of this paper. For training, we used first 10%
of training samples of the LoDoPab-CT dataset; the validation samples are 50
samples in the validation fold of the dataset.

Figure 6: Atoms in a randomly initialized (panel a) and learned decoder
dictionary based on the derain dataset (panel b), LoDoPab-CT (panel c) and
fastMRI (panel d) dataset. For all panels, each row corresponds to a support
size (denoted by $s$) of atoms. Top rows are atoms that have a small support
size; bottom rows are atoms that have a large support size. In each row, atoms
are displayed in a sorted order of decreasing $\ell_{2}$ norm; for the
visualization purpose, they are normalized into the range $[-1,1]$.

###  IV-C MRI reconstruction

We further considered the task of accelerated magnetic resonance imaging (MRI)
reconstruction using the fastMRI dataset [65] procured by Facebook and NYU.
Specifically, we used the single-coil knee dataset with a 4-fold acceleration
factor. This dataset contains 973 volumes or 34742 slices in the training set,
199 volumes or 7135 slices in the validation set, and 108 volumes or 3903
slices in the test set. The ground-truth images in the test set are not
provided to the public and the evaluation must be made through the fastMRI
online submission system333https://fastmri.org/.

Following the training protocol of [65], we first transformed the undersampled
k-space measurements into the image space using zero-filled Inverse Fast
Fourier Transform (IFFT); we use the transformed images as input to MUSC and
other CNN baselines. Consistent with previous work [65], we found that U-Net
variants deliver exceptional performance on validation samples (Table IV).
Remarkably, MUSC performs on-par with U-Net variants, yielding visually
indistinguishable results (Figure 4). We next evaluate the U-Net and the MUSC
on test samples through the fastMRI submission system. On the test data, the
proposed MUSC produces results comparable to the best-performing U-Net result
(fastMRI U-Net-256) provided by the fastMRI challenge organizer while having
an order of magnitude fewer parameters (Table V).

Figure 7: Sparsity of dictionary-induced convolutional features maps. Each bar
corresponds to the sparsity level of a feature map tensor from the “deepest”
activations corresponding to large-support atoms (“Middle”) to the
“shallowest’ activations corresponding to small-support atoms (“Up-4”).

###  IV-D Single-image super-resolution

We have additionally tested the MUSC on a standard super-resolution task,
whose results are deferred to Appendix G. The goal of this task is to recover
a high-resolution image from its degraded, low-resolution version. Unlike
tasks such as CT and MRI reconstruction, in which the image degradation
processes introduce long-range spatially correlated noise like streak
artifacts, the blurring process in the super-resolution task is spatially
local.

In this case, we do not observe a performance gain of using a multiscale model
– either U-Net or MUSC – over state-of-the-art single-scale CSC models.
Interestingly, MUSC outperforms the U-Net, but is up to 0.5 dB worse than
single-scale CSC.

In subsection IV-F, we study this phenomenon by analyzing the sparse code
yielded by MUSC. In the super-resolution tasks, the nonzeros in sparse codes
are confined to high-resolution channels, or, equivalently, small filter
supports which only leverage local information. This is well aligned with the
intuition that the blurring forward operator mixes information only locally.
It suggests that the right strategy is to use a large number of small-support
filters just like CSC does, instead of “wasting” trainable parameters on
unused large scales. We similarly find that a single-scale CSC model works
better than MUSC on a denoising (Gaussian noise removal) task. Together, these
findings suggest that multiscale features are no panacea for imaging inverse
problems; the configuration of scales needs to resonate with the task-
dependent forward operator that we aim to invert.

###  IV-E Ablations on the choices of model components

In Figure 5, we show ablation experiments that demonstrate how different
choices of model components affect the overall performance. There, Case 2 is
the off-the-shelf setup we have used in all other sections of this paper; this
option has the fastest learning speed and highest end-point accuracy.
Consistent with findings in [26, 25, 6], we find it advantageous to use untied
adjoints as described in (11): Untied dictionaries (Cases 1, 2, and 4 in
Figure 5) in general perform much better than tied dictionaries with
$\boldsymbol{D}_{\textrm{enc}}=\boldsymbol{D}_{\textrm{enc}}=\boldsymbol{D}_{\textrm{dec}}$
(Cases 5 and Case 6). What is more, we find that learnable threshold
$\boldsymbol{\lambda}$ gives better results than fixed threshold. The non-
negative constraint of sparse code $\boldsymbol{\alpha}\geq 0$ does not
greatly influence the end-point performance of models, although with the
constraint the model learns slightly faster (Case 2) than without (Case 1).

###  IV-F Probing multiscale dictionary-based representations

Thus far, we have shown that our proposed multiscale CSC approach, dubbed
MUSC, performs comparably to state-of-the-art CNNs in a range of imaging
inverse problems. This is noteworthy, as the strong performance is achieved
simply by employing a multiscale dictionary – as opposed to a single-scale one
– in an otherwise standard CSC paradigm. The strong performance suggests the
usefulness of the multiscale representation. We now analyze our learned
dictionaries and their induced sparse representations.

#### Visualizing dictionary atoms

We visualize dictionary atoms of the MUSC. To extract a dictionary atom from a
dictionary $\boldsymbol{D}$, we first prepare an _indicator code_
$\boldsymbol{\delta}$, which is a collection of multichannel tensors that
takes a value $1$ at a certain entry and $0$ elsewhere; a dictionary atom
corresponds to that entry is computed as $\boldsymbol{D}\boldsymbol{\delta}$.
Note that, different positions of the nonzero entry may give rise to atoms of
different support sizes. This can be seen in Figure 1b: The indicator code is
illustrated as the grey boxes; depending on the grey box the nonzero entry
resides in, the sparse code activates different receptive fields under
composite convolutions and transposed convolutions. If the nonzero entry
resides in the top-most box, then the support of the atom is 3 as it undergoes
only a single 3 × 3 convolution; if the nonzero entry is in one of the lower
boxes, the support of the atom is larger as the code undergoes multiple
convolutions and one or more transposed convolutions.

In Figure 6, we show samples of multiscale atoms in
$\boldsymbol{D}_{\textrm{dec}}$ of varying sizes – we crop these atoms to only
show their nonzero support regions. As can be seen in Figure 6b-d, the learned
dictionaries contain Gabor-like or curvelet-like atoms with different spatial
widths, resolutions, and orientations. Thus the learned dictionaries indeed
exploit multiscale features. For comparison, we also show a randomly
initialized dictionary (Figure 6a). Unlike a learned dictionary, a random
dictionary does not exhibit structures in atoms. We also visualize atoms of
encoder dictionaries $\boldsymbol{D}_{\textrm{enc}}$ and
$\widetilde{\boldsymbol{D}}_{\textrm{enc}}$ in Appendix E. Using a similar
technique, we also probe the multiscale representations learned by U-Nets in
Appendix F.

#### Sparsity levels of representations

We anticipate that the trained dictionaries induce different sparsity levels
at different resolution levels in a task-dependent manner: More non-zeros
associated with large-support atoms are useful when imaging artifacts have
long-range correlations (e.g., streak artifacts in CT) than when the artifacts
are localized (e.g., deraining or super-resolution).

Figure 7 shows the sparsity levels across tasks, both before and after
dictionary learning. We observe that, prior to any learning, the sparsity
levels induced by randomly initialized dictionaries (grey bars) are
approximately uniform across scales. After learning, the sparsity levels of
feature maps differentiate in a task-dependent way (orange bars in all
panels). This task-dependent differentiation suggests the usefulness of
multiscale representations – the learned sparsity levels are neither collapsed
to a single scale nor remain uniform across spatial scales; instead, they are
weighted and combined across scales in a problem-dependent way. A curious
effect of multiscale learning arises in super-resolution (panel d): the
activations are nonzero _only in high-resolution features_ (“Up-2”, “Up-3”,
and “Up-4”), corroborating the intuition that low-resolution features are not
important for this task. Additionally, comparing the “Middle” bars across
panels, we see that CT and MRI reconstruction tasks indeed use more nonzero
coefficients on large-support atoms than tasks such as deraining and super-
resolution.

##  V Discussion

The CSC paradigm provides a natural connection between sparse modeling and
CNNs. Despite being mathematical principled, existing CSC models still fall
short competing with CNNs in terms of empirical performance on challenging
inverse problems. In this work, we report one simple and effective way to
close the performance gap between CSC and CNN models: incorporating a
multiscale structure in the CSC dictionaries. Crucial to our approach is the
structure of our constructed multiscale dictionary: It takes inspiration from
and closely follows the highly successful U-Net model. We show that the
constructed multiscale dictionary performs on par with leading CNNs in major
imaging inverse problems. These results suggest a strong link between
dictionary learning and CNNs– in both cases, multiscale structures are
essential ingredients.

Beyond empirical performance, we believe that the interpretability of the
proposed MUSC is showing the way towards an interpretable deep learning model.
An interpretable model consists of components whose objectives and
functionality have nominal values. The MUSC fulfills this desideratum by
incorporating two modules with well-understood objectives, one implementing
sparse coding and the other dictionary-based synthesis.

Overall, our work demonstrates the effectiveness and scalability of CSC models
on imaging inverse problems. While deep neural networks are profoundly
influencing image reconstruction, our work shows promise in a different
direction: the principles of sparsity and multiscale representation developed
decades ago are still useful in designing performant, parameter-efficient
(compared to mainstream CNNs), and interpretable architectures that push the
current limits of machine learning for imaging inverse problems.

## Acknowledgments

We thank Johannes Leuschner for helpful discussions regarding the LoDoPaB-CT
dataset. Numerical experiments were partly performed at the sciCORE
(http://scicore.unibas.ch/) scientific computing center at the University of
Basel. We thank all anonymous reviewers whose comments helped improve this
manuscript.

## References

  * [1] J. Mairal, P. Koniusz, Z. Harchaoui, and C. Schmid, “Convolutional kernel networks,” in _Proceedings of Neural Information Processing Systems (NeurIPS)_ , 2020.
  * [2] J. Mairal, “End-to-end kernel learning with supervised convolutional kernel networks,” in _Proceedings of Neural Information Processing Systems (NeurIPS)_ , 2016.
  * [3] A. Bietti and J. Mairal, “Group invariance, stability to deformations, and complexity of deep convolutional representations,” _Journal of Machine Learning Research_ , vol. 20, no. 25, pp. 1–49, 2019.
  * [4] S. Mallat, “Group invariant scattering,” _Communications on Pure and Applied Mathematics_ , vol. 65, no. 10, pp. 1331–1398, 2012.
  * [5] J. Bruna and S. Mallat, “Invariant scattering convolution networks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 35, no. 8, pp. 1872–1886, 2013.
  * [6] J. Zarka, L. Thiry, T. Angles, and S. Mallat, “Deep network classification by scattering and homotopy dictionary learning,” in _International Conference on Learning Representations (ICLR)_ , 2020.
  * [7] J. Zarka, F. Guth, and S. Mallat, “Separation and concentration in deep networks,” in _International Conference on Learning Representations (ICLR)_ , 2021.
  * [8] V. Papyan, Y. Romano, and M. Elad, “Convolutional neural networks analyzed via convolutional sparse coding,” _Journal of Machine Learning Research_ , vol. 18, no. 83, pp. 1–52, 2017.
  * [9] V. Papyan, Y. Romano, J. Sulam, and M. Elad, “Theoretical foundations of deep learning via sparse representations: A multilayer sparse model and its connection to convolutional neural networks,” _IEEE Signal Processing Magazine_ , vol. 35, no. 4, pp. 72–89, 2018.
  * [10] D. Simon and M. Elad, “Rethinking the CSC model for natural images,” in _Proceedings of Neural Information Processing Systems (NeurIPS)_ , vol. 32, 2019.
  * [11] K. Degraux, U. S. Kamilov, P. T. Boufounos, and D. Liu, “Online convolutional dictionary learning for multimodal imaging,” in _2017 IEEE International Conference on Image Processing (ICIP)_ , 2017.
  * [12] I. Y. Chun and J. A. Fessler, “Convolutional dictionary learning: Acceleration and convergence,” _IEEE Transactions on Image Processing_ , vol. 27, no. 4, pp. 1697–1712, 2018.
  * [13] E. Zisselman, J. Sulam, and M. Elad, “A local block coordinate descent algorithm for the CSC model,” in _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019, pp. 8200–8209.
  * [14] I. Rey-Otero, J. Sulam, and M. Elad, “Variations on the convolutional sparse coding model,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 519–528, 2020.
  * [15] N. Chodosh and S. Lucey, “When to use convolutional neural networks for inverse problems,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020, pp. 8223–8232.
  * [16] B. Tolooshams, A. H. Song, S. Temereanca, and D. E. Ba, “Convolutional dictionary learning based auto-encoders for natural exponential-family distributions,” in _Proceedings of the 37th International Conference on Machine Learning (ICML)_ , 2020.
  * [17] P. V. Arun, K. M. Buddhiraju, A. Porwal, and J. Chanussot, “CNN-based super-resolution of hyperspectral images,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 58, no. 9, pp. 6106–6121, 2020.
  * [18] J.-J. Huang and P. L. Dragotti, “Learning deep analysis dictionaries for image super-resolution,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 6633–6648, 2020.
  * [19] S. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 11, no. 7, pp. 674–693, 1989.
  * [20] E. J. Candès and D. L. Donoho, “New tight frames of curvelets and optimal representations of objects with piecewise ${C}^{2}$ singularities,” _Communications on Pure and Applied Mathematics_ , vol. 57, no. 2, pp. 219–266, 2004.
  * [21] J.-L. Starck, E. Candes, and D. Donoho, “The curvelet transform for image denoising,” _IEEE Transactions on Image Processing_ , vol. 11, no. 6, pp. 670–684, 2002.
  * [22] J. Mairal, G. Sapiro, and M. Elad, “Learning multiscale sparse representations for image and video restoration,” _Multiscale Modeling & Simulation_, vol. 7, no. 1, pp. 214–241, 2008.
  * [23] B. Ophir, M. Lustig, and M. Elad, “Multi-scale dictionary learning using wavelets,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 5, no. 5, pp. 1014–1024, 2011.
  * [24] M. Maggioni, S. Minsker, and N. Strawn, “Multiscale dictionary learning: Non-asymptotic bounds and robustness,” _Journal of Machine Learning Research_ , vol. 17, no. 2, pp. 1–51, 2016.
  * [25] B. Lecouat, J. Ponce, and J. Mairal, “Fully trainable and interpretable non-local sparse models for image restoration,” in _European Conference on Computer Vision (ECCV) 2020_ , 2020.
  * [26] ——, “A flexible framework for designing trainable priors with adaptive smoothing and game encoding,” in _Proceedings of Neural Information Processing Systems (NeurIPS)_ , 2020.
  * [27] R. Khatib, D. Simon, and M. Elad, “Learned greedy method (LGM): A novel neural architecture for sparse coding and beyond,” _Journal of Visual Communication and Image Representation_ , vol. 77, p. 103095, 2021.
  * [28] M. Scetbon, M. Elad, and P. Milanfar, “Deep K-SVD denoising,” _IEEE Transactions on Image Processing_ , vol. 30, pp. 5944–5955, 2021.
  * [29] M. Li, Q. Xie, Q. Zhao, W. Wei, S. Gu, J. Tao, and D. Meng, “Video rain streak removal by multiscale convolutional sparse coding,” in _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 6644–6653.
  * [30] H. Chang, J. Han, C. Zhong, A. M. Snijders, and J.-H. Mao, “Unsupervised transfer learning via multi-scale convolutional sparse coding for biomedical applications,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 40, no. 5, pp. 1182–1194, 2018.
  * [31] 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_ , N. Navab, J. Hornegger, W. M. Wells, and A. F. Frangi, Eds. Springer International Publishing, 2015, pp. 234–241.
  * [32] T. Falk, D. Mai, R. Bensch, Ö. Çiçek, A. Abdulkadir, Y. Marrakchi, A. Böhm, J. Deubner, Z. Jäckel, K. Seiwald, A. Dovzhenko, O. Tietz, C. Dal Bosco, S. Walsh, D. Saltukoglu, T. L. Tay, M. Prinz, K. Palme, M. Simons, I. Diester, T. Brox, and O. Ronneberger, “U-net: deep learning for cell counting, detection, and morphometry,” _Nature Methods_ , vol. 16, no. 1, pp. 67–70, Jan 2019.
  * [33] Z. Zhou, M. M. R. Siddiquee, N. Tajbakhsh, and J. Liang, “Unet++: Redesigning skip connections to exploit multiscale features in image segmentation,” _IEEE Transactions on Medical Imaging_ , vol. 39, no. 6, pp. 1856–1867, June 2020.
  * [34] D. Ren, W. Zuo, Q. Hu, P. Zhu, and D. Meng, “Progressive image deraining networks: A better and simpler baseline,” in _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019, pp. 3932–3941.
  * [35] G. Beylkin, R. Coifman, and V. Rokhlin, “Fast wavelet transforms and numerical algorithms i,” _Communications on Pure and Applied Mathematics_ , vol. 44, pp. 141–183, 1991.
  * [36] A. M. Bruckstein, D. L. Donoho, and M. Elad, “From sparse solutions of systems of equations to sparse modeling of signals and images,” _SIAM review_ , vol. 51, no. 1, pp. 34–81, 2009.
  * [37] M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing MRI,” _IEEE signal processing magazine_ , vol. 25, no. 2, pp. 72–82, 2008\.
  * [38] J. Charlety, S. Voronin, G. Nolet, I. Loris, F. J. Simons, K. Sigloch, and I. C. Daubechies, “Global seismic tomography with sparsity constraints: Comparison with smoothing and damping regularization,” _Journal of Geophysical Research: Solid Earth_ , vol. 118, no. 9, pp. 4887–4899, 2013.
  * [39] H. Liebgott, R. Prost, and D. Friboulet, “Pre-beamformed RF signal reconstruction in medical ultrasound using compressive sensing,” _Ultrasonics_ , vol. 53, no. 2, pp. 525 – 533, 2013.
  * [40] C. Chen and D. Donoho, “Basis pursuit,” in _Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers_ , vol. 1, 1994, pp. 41–44 vol.1.
  * [41] R. Tibshirani, “Regression shrinkage and selection via the Lasso,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [42] I. Daubechies, M. Defrise, and C. D. Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” _Communications on Pure and Applied Mathematics_ , vol. 57, no. 11, pp. 1413–1457, 2004.
  * [43] J. Mairal, F. Bach, and J. Ponce, “Task-driven dictionary learning,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 34, no. 4, p. 791–804, April 2012.
  * [44] H. Bristow, A. Eriksson, and S. Lucey, “Fast convolutional sparse coding,” in _2013 IEEE Conference on Computer Vision and Pattern Recognition_ , 2013, pp. 391–398.
  * [45] M. Šorel and F. Šroubek, “Fast convolutional sparse coding using matrix inversion lemma,” _Digital Signal Processing_ , vol. 55, pp. 44–51, 2016.
  * [46] J. Liu, C. Garcia-Cardona, B. Wohlberg, and W. Yin, “Online convolutional dictionary learning,” in _2017 IEEE International Conference on Image Processing (ICIP)_ , 2017.
  * [47] H. Sreter and R. Giryes, “Learned convolutional sparse coding,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 2191–2195.
  * [48] D. O. Baguer, J. Leuschner, and M. Schmidt, “Computed tomography reconstruction using deep image prior and learned reconstruction methods,” _Inverse Problems_ , vol. 36, no. 9, p. 094004, Sep 2020.
  * [49] J. Leuschner, M. Schmidt, P. S. Ganguly, V. Andriiashen, S. B. Coban, A. Denker, D. Bauer, A. Hadjifaradji, K. J. Batenburg, P. Maass, and M. van Eijnatten, “Quantitative comparison of deep learning-based image reconstruction methods for low-dose and sparse-angle CT applications,” _Journal of Imaging_ , vol. 7, no. 3, 2021.
  * [50] T. Huang, W. Dong, X. Yuan, J. Wu, and G. Shi, “Deep gaussian scale mixture prior for spectral compressive imaging,” in _2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2021, pp. 16 211–16 220.
  * [51] W. Dong, P. Wang, W. Yin, G. Shi, F. Wu, and X. Lu, “Denoising prior driven deep neural network for image restoration,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 41, no. 10, pp. 2305–2318, 2019\.
  * [52] C. Etmann, R. Ke, and C.-B. Schönlieb, “iunets: Learnable invertible up- and downsampling for large-scale inverse problems,” in _2020 IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP)_ , 2020, pp. 1–6.
  * [53] J. C. Ye, Y. Han, and E. Cha, “Deep convolutional framelets: A general deep learning framework for inverse problems,” _SIAM Journal on Imaging Sciences_ , vol. 11, no. 2, pp. 991–1048, 2018.
  * [54] Y. Han and J. C. Ye, “Framing U-Net via deep convolutional framelets: Application to sparse-view CT,” _IEEE Transactions on Medical Imaging_ , vol. 37, no. 6, pp. 1418–1429, 2018.
  * [55] J. C. Ye and W. K. Sung, “Understanding geometry of encoder-decoder CNNs,” in _Proceedings of the 36th International Conference on Machine Learning, ICML_ , C. Chaudhuri and S. Ruslan, Eds., 2019.
  * [56] R. Yin, T. Gao, Y. M. Lu, and I. Daubechies, “A tale of two bases: Local-nonlocal regularization on image patches with convolution framelets,” _SIAM Journal on Imaging Sciences_ , vol. 10, no. 2, pp. 711–750, 2017.
  * [57] K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” in _Proceedings of the 27th International Conference on Machine Learning (ICML)_ , J. Fürnkranz and T. Joachims, Eds., 2010, pp. 399–406.
  * [58] T. Moreau and J. Bruna, “Understanding trainable sparse coding with matrix factorization,” in _International Conference on Learning Representations (ICLR)_ , 2017.
  * [59] J. Liu, X. Chen, Z. Wang, and W. Yin, “ALISTA: Analytic weights are as good as learned weights in LISTA,” in _International Conference on Learning Representations (ICLR)_ , 2019.
  * [60] Y. Li, R. T. Tan, X. Guo, J. Lu, and M. S. Brown, “Rain streak removal using layer priors,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 2736–2744.
  * [61] X. Fu, J. Huang, D. Zeng, Y. Huang, X. Ding, and J. Paisley, “Removing rain from single images via a deep detail network,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 1715–1723.
  * [62] Y. Luo, Y. Xu, and H. Ji, “Removing rain from a single image via discriminative sparse coding,” in _2015 IEEE International Conference on Computer Vision (ICCV)_ , 2015, pp. 3397–3405.
  * [63] W. Yang, R. T. Tan, J. Feng, Z. Guo, S. Yan, and J. Liu, “Joint rain detection and removal from a single image with contextualized deep networks,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 42, no. 6, pp. 1377–1393, 2020.
  * [64] J. Leuschner, M. Schmidt, D. O. Baguer, and P. Maass, “LoDoPaB-CT, a benchmark dataset for low-dose computed tomography reconstruction,” _Scientific Data_ , vol. 8, no. 1, p. 109, Apr 2021.
  * [65] J. Zbontar, F. Knoll, A. Sriram, M. J. Muckley, M. Bruno, A. Defazio, M. Parente, K. J. Geras, J. Katsnelson, H. Chandarana, Z. Zhang, M. Drozdzal, A. R., M. Rabbat, P. Vincent, J. Pinkerton, D. Wang, N. Yakubova, E. Owens, C. L. Zitnick, M. P. Recht, D. K. Sodickson, and Y. W. Lui, “fastMRI: An open dataset and benchmarks for accelerated MRI,” _ArXiv e-prints_ , 2018\.
  * [66] S. G. Armato, G. McLennan, L. Bidaut, M. F. McNitt-Gray, C. R. Meyer, A. P. Reeves, B. Zhao, D. R. Aberle, C. I. Henschke, E. A. Hoffman, E. A. Kazerooni, H. MacMahon, E. J. Van Beeke, D. Yankelevitz, A. M. Biancardi, P. H. Bland, M. S. Brown, R. M. Engelmann, G. E. Laderach, D. Max, R. C. Pais, D. P. Qing, R. Y. Roberts, A. R. Smith, A. Starkey, P. Batrah, P. Caligiuri, A. Farooqi, G. W. Gladish, C. M. Jude, R. F. Munden, I. Petkovska, L. E. Quint, L. H. Schwartz, B. Sundaram, L. E. Dodd, C. Fenimore, D. Gur, N. Petrick, J. Freymann, J. Kirby, B. Hughes, A. V. Casteele, S. Gupte, M. Sallamm, M. D. Heath, M. H. Kuhn, E. Dharaiya, R. Burns, D. S. Fryd, M. Salganicoff, V. Anand, U. Shreter, S. Vastagh, and B. Y. Croft, “Med PhysThe Lung Image Database Consortium (LIDC) and Image Database Resource Initiative (IDRI): a completed reference database of lung nodules on CT scans,” _Med Phys_ , vol. 38, no. 2, pp. 915–931, Feb 2011.
  * [67] A. Denker, M. Schmidt, J. Leuschner, P. Maass, and J. Behrmann, “Conditional normalizing flows for low-dose computed tomography image reconstruction,” _ArXiv e-prints_ , 2020.
  * [68] Z. Zhou, M. M. Rahman Siddiquee, N. Tajbakhsh, and J. Liang, “UNet++: A nested U-Net architecture for medical image segmentation,” in _Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support_ , D. Stoyanov, Z. Taylor, G. Carneiro, T. Syeda-Mahmood, A. Martel, L. Maier-Hein, J. M. R. Tavares, A. Bradley, J. P. Papa, V. Belagiannis, J. C. Nascimento, Z. Lu, S. Conjeti, M. Moradi, H. Greenspan, and A. Madabhushi, Eds. Cham: Springer International Publishing, 2018, pp. 3–11.
  * [69] D. M. Pelt, K. J. Batenburg, and J. A. Sethian, “Improving tomographic reconstruction from limited data using mixed-scale dense convolutional neural networks,” _Journal of Imaging_ , vol. 4, no. 11, 2018.
  * [70] T. Salimans and D. P. Kingma, “Weight normalization: A simple reparameterization to accelerate training of deep neural networks,” in _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , ser. NIPS’16. Red Hook, NY, USA: Curran Associates Inc., 2016, p. 901–909.
  * [71] S. Gu, W. Zuo, Q. Xie, D. Meng, X. Feng, and L. Zhang, “Convolutional sparse coding for image super-resolution,” in _2015 IEEE International Conference on Computer Vision (ICCV)_ , 2015, pp. 1823–1831.
  * [72] J. He, L. Yu, Z. Liu, and W. Yang, “Image super-resolution by learning weighted convolutional sparse coding,” _Signal, Image and Video Processing_ , vol. 15, no. 5, pp. 967–975, Jul 2021.
  * [73] S. Schulter, C. Leistner, and H. Bischof, “Fast and accurate image upscaling with super-resolution forests,” in _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015, pp. 3791–3799.
  * [74] J. Kim, J. K. Lee, and K. M. Lee, “Accurate image super-resolution using very deep convolutional networks,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 1646–1654.
  * [75] W. Lai, J. Huang, N. Ahuja, and M. Yang, “Deep laplacian pyramid networks for fast and accurate super-resolution,” in _IEEE Conference on Computer Vision and Pattern Recognition_ , 2017.
  * [76] E. Agustsson and R. Timofte, “Ntire 2017 challenge on single image super-resolution: Dataset and study,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) Workshops_ , July 2017.
  * [77] R. Zeyde, M. Elad, and M. Protter, “On single image scale-up using sparse-representations,” in _Proceedings of the 7th International Conference on Curves and Surfaces_. Berlin, Heidelberg: Springer-Verlag, 2010, p. 711–730.
  * [78] M. Bevilacqua, A. Roumy, C. Guillemot, and M. Morel, “Low-complexity single-image super-resolution based on nonnegative neighbor embedding,” in _Proceedings of the British Machine Vision Conference_. BMVA Press, 2012, pp. 135.1–135.10.
  * [79] D. Martin, C. Fowlkes, D. Tal, and J. Malik, “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,” in _Proceedings of the 8th International Conference on Computer Vision_ , vol. 2, July 2001, pp. 416–423.
  * [80] J.-B. Huang, A. Singh, and N. Ahuja, “Single image super-resolution from transformed self-exemplars,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2015, pp. 5197–5206.
  * [81] J. Yang, J. Wright, T. S. Huang, and Y. Ma, “Image super-resolution via sparse representation,” _IEEE Transactions on Image Processing_ , vol. 19, no. 11, pp. 2861–2873, 2010.
  * [82] R. Timofte, V. D. Smet, and L. V. Gool, “A+: Adjusted anchored neighborhood regression for fast super-resolution,” in _Computer Vision – ACCV 2014_ , D. Cremers, I. Reid, H. Saito, and M. Yang, Eds. Springer International Publishing, 2015, pp. 111–126.

##  Appendix A The definition of the MUSC dictionary

In the main text, we illustrate the architecture of a U-Net (Figure 1a) and
the corresponding MUSC decoder dictionary (Figure 1b). Loosely speaking, the
decoder dictionary is the decoding branch of a standard U-Net with all ReLU
activations, batch normalization, and some convolution operations removed.
Here we provide a formal definition of the decoder dictionary. While in the
main text we assume for simplicty that signals are 1D vectors, in this section
we represent RGB images as multichannel tensors. To that end, we first
consider single-input single-output (SISO) operations whose input and output
are 2D signals having a single channel. We follow by considering multiple-
input multiple-output (MIMO) operations whose input and output are 3D tensors
having multiple channels.

#### SISO convolution and transposed convolution

Let $\boldsymbol{\xi}\in\mathbb{R}^{H\times W}$ be a 2D signal with height $H$
and width $W$. We regard $\boldsymbol{\xi}$ as a function defined on the
discrete domain $\\{1,\ldots,H\\}\times\\{1,\ldots,W\\}$. With standard zero
padding, we extend the domain of $\boldsymbol{\xi}$ to
$\mathbb{Z}\times\mathbb{Z}$. The notation $\boldsymbol{\xi}[i,j]$ represents
the value of the function $\boldsymbol{\xi}$ at the coordinate $(i,j)$.

#### SISO convolution

Given a 2D signal $\boldsymbol{\xi}\in\mathbb{R}^{H\times W}$ and parameters
(filter weights, filter impulse response)
$\boldsymbol{w}\in\mathbb{R}^{3\times 3}$, the convolution of
$\boldsymbol{\xi}$ and $\boldsymbol{w}$ is defined as

$(\boldsymbol{w}\ast\boldsymbol{\xi})[i,j]=\sum_{p=-1}^{1}\sum_{q=-1}^{1}\boldsymbol{\xi}[i+p,j+q]\cdot\boldsymbol{w}[p,q].$
(14)

#### SISO transposed convolution

A transposed convolution in the U-Net consists of a bed-of-nails upsampling by
a factor of two followed by a 2-by-2 convolution with an interpolating filter.
The bed-of-nails upsampling interleaves zeros between samples, which can be
written as $\boldsymbol{\xi}\otimes\begin{bmatrix}0&0\\\ 0&1\end{bmatrix}$
with $\otimes$ denoting the Kronecker product. The 2-by-2 convolution between
the filter $\boldsymbol{v}\in\mathbb{R}^{2\times 2}$ and signal
$\boldsymbol{\xi}\in\mathbb{R}^{H\times W}$ is written as

$(\boldsymbol{v}\circledast\boldsymbol{\xi})[i,j]=\sum_{p=0}^{1}\sum_{q=0}^{1}\boldsymbol{\xi}[i+p,j+q]\cdot\boldsymbol{v}[p,q].$
(15)

Putting these together, we have the following definition for a transposed
convolution.

######  Definition 1.

A transposed convolution $\mathcal{U}_{\boldsymbol{v}}$ parameterized by the
convolutional filter $\boldsymbol{v}\in\mathbb{R}^{2\times 2}$ is a map
defined by

$\displaystyle\mathcal{U}_{\boldsymbol{v}}\colon\mathbb{R}^{n\times n}$
$\displaystyle\to\mathbb{R}^{2n\times 2n}$ $\displaystyle\boldsymbol{\xi}$
$\displaystyle\mapsto\boldsymbol{v}\circledast(\boldsymbol{\xi}\otimes\begin{bmatrix}0&0\\\
0&1\end{bmatrix}).$

The following technical result simplifies the upsampling calculation.

######  Proposition 1.

Let $\boldsymbol{v}\in\mathbb{R}^{2\times 2}$ be a matrix and let
$\bar{\boldsymbol{v}}$ be its vertically and horizontally flipped version,
i.e., $\bar{\boldsymbol{v}}[i,j]=\boldsymbol{v}[3-i,3-j]$ for
$(i,j)\in\\{1,2\\}\times\\{1,2\\}$. Let
$\boldsymbol{\xi}\in\mathbb{R}^{n\times n}$ be a matrix with dimension $n\geq
2$. Then the following equation holds:

$\mathcal{U}_{\boldsymbol{v}}(\boldsymbol{\xi})=\boldsymbol{\xi}\otimes\bar{\boldsymbol{v}}.$
(16)

###### Proof.

Observe that the result of LHS and RHS are both $2n$-by-$2n$ matrices, each
composed of $n\times n$ number of 2-by-2 submatrices. It is enough to show
that each $2$-by-$2$ sub-matrix of the LHS equals that of the RHS if they have
the same position indices.

We first expand the RHS. By definition of Kronecker product, we have

$\boldsymbol{\xi}\otimes\bar{\boldsymbol{v}}={\begin{bmatrix}\boldsymbol{\xi}[1,1]\bar{\boldsymbol{v}}&\cdots&\boldsymbol{\xi}[1,n]\bar{\boldsymbol{v}}\\\
\vdots&\ddots&\vdots\\\
\boldsymbol{\xi}[n,1]\bar{\boldsymbol{v}}&\cdots&\boldsymbol{\xi}[n,n]\bar{\boldsymbol{v}}\end{bmatrix}}\in\mathbb{R}^{2n\times
2n},$

where the $(i,j)$-th 2-by-2 sub-matrix of
$\boldsymbol{\xi}\otimes\bar{\boldsymbol{v}}$ has the form

$\boldsymbol{\xi}[i,j]\bar{\boldsymbol{v}}=\begin{bmatrix}\boldsymbol{\xi}[i,j]\cdot{\boldsymbol{v}}[2,2]&\boldsymbol{\xi}[i,j]\cdot{\boldsymbol{v}}[2,1]\\\
\boldsymbol{\xi}[i,j]\cdot{\boldsymbol{v}}[1,2]&\boldsymbol{\xi}[i,j]\cdot{\boldsymbol{v}}[1,1]\end{bmatrix},$
(17)

for each $(i,j)\in\\{1,\cdots,n\\}\times\\{1,\cdots,n\\}$.

For the LHS, we have
$\mathcal{U}_{\boldsymbol{v}}(\boldsymbol{\xi})=\boldsymbol{v}\circledast(\boldsymbol{\xi}\otimes\begin{bmatrix}0&0\\\
0&1\end{bmatrix})$ by definition. Writing out the Kronecker product, we see
that the $(i,j)$-th 2-by-2 submatrix of the LHS has the form

$\boldsymbol{v}\circledast\begin{bmatrix}0&0&0\\\ 0&\boldsymbol{\xi}[i,j]&0\\\
0&0&0\end{bmatrix}=\begin{bmatrix}\boldsymbol{\xi}[i,j]\cdot\boldsymbol{v}[2,2]&\boldsymbol{\xi}[i,j]\cdot\boldsymbol{v}[2,1]\\\
\boldsymbol{\xi}[i,j]\cdot\boldsymbol{v}[1,2]&\boldsymbol{\xi}[i,j]\cdot\boldsymbol{v}[1,1]\end{bmatrix},$
(18)

for each $(i,j)\in\\{1,\cdots,n\\}\times\\{1,\cdots,n\\}$. Comparing Equation
(17) and Equation (18), we see the each 2-by-2 sub-matrix of the LHS and RHS
of Equation (16) are identical as claimed. ∎

#### MIMO convolution and transposed convolution

To extend our formulation to the MIMO case, we denote by
$\boldsymbol{x}=(\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{C})\in\mathbb{R}^{C\times
H\times W}$ a multi-channel signal.

#### MIMO convolution

We let
$\boldsymbol{W}=(\boldsymbol{w}_{11},\ldots,\boldsymbol{w}_{1C},\ldots,\boldsymbol{w}_{M1},\ldots,\boldsymbol{w}_{MC})\in\mathbb{R}^{M\times
C\times 3\times 3}$ be a multi-channel convolutional kernel of 3-by-3 filters,
where $C$ denotes the number of input channels and $M$ denotes the number of
output channels. The MIMO convolution is defined as

$\text{Conv}(\boldsymbol{W},\boldsymbol{x})=\sum_{c=1}^{C}\left(\boldsymbol{w}_{1c}\ast\boldsymbol{\xi}_{c},\ldots,\boldsymbol{w}_{Mc}\ast\boldsymbol{\xi}_{c}\right)\in\mathbb{R}^{M\times
H\times W}.$

#### MIMO transposed convolution

Let
$\boldsymbol{x}=(\boldsymbol{\xi}_{1},\ldots,\boldsymbol{\xi}_{C})\in\mathbb{R}^{C\times
H\times W}$ be a multi-channel signal. We let
$\boldsymbol{V}=(\boldsymbol{v}_{11},\ldots,\boldsymbol{v}_{1C},\ldots,\boldsymbol{v}_{M1},\ldots,\boldsymbol{v}_{MC})\in\mathbb{R}^{M\times
C\times 2\times 2}$ be a multi-channel convolutional kernel of 2-by-2 filters.
The MIMO transposed convolution is defined as

$\text{TConv}(\boldsymbol{V},\boldsymbol{x})=\sum_{c=1}^{C}\left(\boldsymbol{\xi}_{c}\otimes\bar{\boldsymbol{v}}_{1c},\ldots,\boldsymbol{\xi}_{c}\otimes\bar{\boldsymbol{v}}_{Mc}\right)\in\mathbb{R}^{M\times
2H\times 2W}.$

#### Up-block

An up-block receives multichannel input from (i) low-resolution features in
the shape of $(2C,H,W)$ and (ii) skip-connection features in the shape of
$(C,2H,2W)$. With these input features, an up-block applies a transposed
convolution to low-resolution features that halves the number of channels,
concatenates the resulting features with skip-connection features, and then
submits the concatenated results for another convolution. This can be written
as

$\displaystyle\text{Up}_{\boldsymbol{W},\boldsymbol{V}}:\mathbb{R}^{C\times
2H\times 2W}\times\mathbb{R}^{2C\times H\times W}$
$\displaystyle\to\mathbb{R}^{C\times 2H\times 2W}$
$\displaystyle(\boldsymbol{\beta},\boldsymbol{x})$
$\displaystyle\mapsto\text{Conv}([\boldsymbol{\beta};\text{TConv}(\boldsymbol{x},\boldsymbol{V})],\boldsymbol{W}),$

where $\boldsymbol{W}\in\mathbb{R}^{2C\times C\times 3\times 3}$ and
$\boldsymbol{V}\in\mathbb{R}^{2C\times C\times 2\times 2}$ are parameters and
$[\cdot;\cdot]$ denotes channel-wise concatenation.

#### The multiscale dictionary

We now define the multiscale dictionary

$\displaystyle\boldsymbol{D}:\mathbb{R}^{C\times H\times
W}\times\cdots\times\mathbb{R}^{C/2^{4}\times 2^{4}\cdot H\times 2^{4}\cdot
W}$ $\displaystyle\to\mathbb{R}^{C^{\textrm{out}}\times 2^{4}\cdot H\times
2^{4}\cdot W}$
$\displaystyle(\boldsymbol{\alpha}_{0},\ldots,\boldsymbol{\alpha}_{4})$
$\displaystyle\mapsto\boldsymbol{D}(\boldsymbol{\alpha}_{0},\ldots,\boldsymbol{\alpha}_{4})$
(19)

that is used as encoders and decoders in our work. Here,
$C^{\textrm{out}}=1\textrm{~{}or~{}}3$ for grayscale and RGB images. We
construct the dictionary $\boldsymbol{D}$ as a linear transform
$\boldsymbol{R}$ that is independent of $C^{\textrm{out}}$ followed by a one-
by-one convolution that produce $C^{\textrm{out}}$ channels:

$\boldsymbol{D}=\texttt{Conv1x1}\circ\boldsymbol{R}.$ (20)

The linear map $\boldsymbol{R}$ is defined by cascading $4$ Up-blocks with
parameters $\boldsymbol{W}_{i}$ and $\boldsymbol{V}_{i},i=1,\ldots,4$:

$\displaystyle\boldsymbol{R}:\mathbb{R}^{C\times H\times
W}\times\cdots\times\mathbb{R}^{C/2^{4}\times 2^{4}\cdot H\times 2^{4}\cdot
W}$ $\displaystyle\to\mathbb{R}^{C/2^{4}\times 2^{4}\cdot H\times 2^{4}\cdot
W}$ $\displaystyle(\boldsymbol{\alpha}_{0},\cdots,\boldsymbol{\alpha}_{4})$
$\displaystyle\mapsto\boldsymbol{\xi}_{4},$ (21)

where

$\boldsymbol{\xi}_{i}=\text{Up}_{\boldsymbol{W}_{i},\boldsymbol{V}_{i}}(\boldsymbol{\alpha}_{i},\boldsymbol{\xi}_{i-1})$

for $i=1,\ldots 4$ with $\boldsymbol{\xi}_{0}=\boldsymbol{\alpha}_{0}$. That
is, the function $\boldsymbol{R}$ transforms multiscale sparse code
$(\boldsymbol{\alpha}_{0},\ldots,\boldsymbol{\alpha}_{4})$ to a tensor
$\boldsymbol{\xi}_{4}$ of shape $\mathbb{R}^{C/2^{4}\times 2^{4}\cdot H\times
2^{4}\cdot W}$. The $1\times 1$ convolution operator $\boldsymbol{C}$ next
synthesizes these features into a tensor with channel number 1 or 3:

$\texttt{Conv1x1}:\mathbb{R}^{C/2^{4}\times 2^{4}\cdot H\times 2^{4}\cdot
W}\to\mathbb{R}^{C^{\textrm{out}}\times 2^{4}\cdot H\times 2^{4}\cdot W}.$
(22)

##  Appendix B Can a learned sparse coder yield dense outputs?

In the main text, we introduced the $K$-fold ISTA algoritm for sparse coding.
The sparse coding result reads
$\mathrm{ISTA}_{K}(\boldsymbol{z};\boldsymbol{D}_{\textrm{enc}},\boldsymbol{\lambda})$,
where $\boldsymbol{z}$ is a given input, $\boldsymbol{D}_{\textrm{enc}}$ is an
encoding dictionary, and $\boldsymbol{\lambda}$ is a threshold tensor. Since
we also learn a $\boldsymbol{\lambda}$ from data during training, one concern
is whether a learned $\boldsymbol{\lambda}$ be zero. Indeed, if this were to
occur, the sparse coding results would be dense and it would defeat the
purpose of sparse coding.

In this section, we argue that $\boldsymbol{\lambda}=0$ is (i) not optimal for
many inverse problems and (ii) unlikely to be learned from data. To see the
non-optimality of $\boldsymbol{\lambda}=0$, we observe that with the iteration
number $K$ is large enough, with $\boldsymbol{\lambda}=0$, and with zero
initialization of sparse code, ISTA yields the smallest $\ell^{2}$-norm
solution which can be written via the Moore–Penrose pseudoinverse,

$\lim_{K\to\infty}\mathrm{ISTA}_{K}(\boldsymbol{z};\boldsymbol{D}_{\textrm{enc}},\boldsymbol{0})=\boldsymbol{D}_{\textrm{enc}}^{+}\boldsymbol{z}.$

Our task-driven dictionary learning problem thus has the form

$\operatorname*{minimize}_{\boldsymbol{D}_{\textrm{enc}},\boldsymbol{D}_{\textrm{dec}}}\quad\frac{1}{2M}\sum_{i=1}^{M}\|\boldsymbol{D}_{\textrm{dec}}\boldsymbol{D}_{\textrm{enc}}^{+}\boldsymbol{\alpha}_{\boldsymbol{z}_{i}}-\boldsymbol{x}^{\dagger}_{i}\|_{2}^{2}.$

This implies that $\lambda=0$ results in a linear reconstruction method which
comes with all the known drawbacks of linear methods. In particular, it cannot
do better than the linear minimum mean square error (LMMSE) estimator (a
generalized Wiener filter). Since regularization in ill-posed inverse problems
entails the use of data models, and most useful data models are nonlinear
(e.g., natural and medical images are known to be sparse or compressible in
wavelet frames, but they do not belong to any linear subspace), these problems
demand $\boldsymbol{\lambda}>0$.

One may still wonder whether our learning procedure will overfit finite
datasets with $\boldsymbol{\lambda}=0$. To see that this will not happen, note
that $\boldsymbol{D}_{\textrm{dec}}$ and $\boldsymbol{D}_{\textrm{enc}}$ are
constrained to have a specific structure: They are multiscale variants of
block-Toeplitz matrices. Additionally, filters at high resolutions are
convolutions of filters at lower resolutions which induces a rather
complicated algebraic structure. As a result, the set of valid dictionaries in
our model has a much smaller dimension than the set of all possible
dictionaries of correct size, and solving (9) cannot be reduced to finding two
generic overcomplete dictionaries that “overfit” to training data to achieve
zero loss. In fact, forcing $\boldsymbol{\lambda}=0$ typically incurs a large
loss in (9), so $\boldsymbol{\lambda}=0$ is unlikely to be learned from data.
This is the magic of multiscale convolutional sparsity.

##  Appendix C Power iteration

We describe how to approximate the dominant eigenvalue of the matrix
$\boldsymbol{D}_{\textrm{enc}}^{\top}\boldsymbol{D}_{\textrm{enc}}$ using by
power iteration. We achieve this by first estimating the eigenvector
associated with the dominant eigenvalue using the power iteration method, by
recursively calculating

$\boldsymbol{b}_{k+1}=\frac{\boldsymbol{D}_{\textrm{enc}}^{\top}\boldsymbol{D}_{\textrm{enc}}\boldsymbol{b}_{k}}{\left\|\boldsymbol{D}_{\textrm{enc}}^{\top}\boldsymbol{D}_{\textrm{enc}}\boldsymbol{b}_{k}\right\|_{2}}$
(23)

up to some step $K$ with $\boldsymbol{b}_{0}$ being an all-zero vector. The
estimated dominant eigenvalue can then be derived from

$\lambda_{K}=\frac{\boldsymbol{b}_{K}^{\top}\boldsymbol{D}_{\textrm{enc}}^{\top}\boldsymbol{D}_{\textrm{enc}}\boldsymbol{b}_{K}}{\boldsymbol{b}_{K}^{\top}\boldsymbol{b}_{K}}.$
(24)

##  Appendix D Details of the experimental setup

| Conv $3\times 3$ | Trans-Conv $2\times 2$ | Trans-Conv $1\times 1$
---|---|---|---
| in-channels | out-channels | stride | in-channels | out-channels | stride | in-channels | out-channels | stride
Scale 1 | 512 | 512 | 1 | 512 | 256 | 2 | – | – | –
Scale 2 | 512 | 256 | 1 | 256 | 128 | 2 | – | – | –
Scale 3 | 256 | 128 | 1 | 128 | 64 | 2 | – | – | –
Scale 4 | 128 | 64 | 1 | 64 | 32 | 2 | – | – | –
Scale 5 | 64 | 32 | 1 | – | – | – | 32 | 1 or 3 | 1
Table VI:  Parameters used at each scale of the encoder and decoder
dictionaries. Scale 1 corresponds to the low-resolution scale (the bottom-most
gray box in Figure 1b) and Scale 5 correspond to the high-resolution scale
(the top-most gray box in Figure 1b).

In Table VI, we summarize parameters used in each scale of our dictionaries.
In Table VII and Table VIII, we summarize the hyperparameters used for
training MUSC and the U-Net baseline.

Task |  Epoch number |  Batch size |  Learning rate |  Num ISTA steps |  Lasso parameter
---|---|---|---|---|---
Derain | 20000 | 16 | 5e-4 | 5 | 0.001
LoDoPaB-CT | 70 | 2 | 2e-4 | 5 | 0.001
fastMRI | 70 | 2 | 5e-5 | 5 | 0.001
Super-resolution | 1500 | 2 | 5e-5 | 5 | 0.001
Table VII: Hyperparameters used for the MUSC.

For Derain and super-resolution task, we randomly crop images into $128\times
128$ patches during gradient-descent training.

Task |  Epoch ber |  Batch size |  Learning rate
---|---|---|---
LoDoPaB-CT | 70 | 32 | 0.001
fastMRI | 70 | 8 | 0.001
Table VIII: Hyperparameters used for the U-Net.

Note that, for CT and MRI tasks, U-Nets and MUSCs are trained on full images;
in deraining and super-resolution tasks, U-Nets and MUSCs are trained on
cropped images, where sizes of the cropped images are in Table VIII. Since in
the deraining and LoDoPaB-CT tasks the range of the target images is non-
negative, we clip the negative values of the synthesized image. During model
training, we use weight normalization [70], a reparametrization trick that
decouples the magnitude of a convolutional filter from its direction. To
enforce the positivity of the ISTA parameter $\boldsymbol{\lambda},$ we
reparametrize
$\boldsymbol{\lambda}=\textrm{ReLU}(\tilde{\boldsymbol{\lambda}})+$1e-5 and
perform gradient-based learning on $\tilde{\boldsymbol{\lambda}}$ instead.

##  Appendix E Visualizing the atoms in encoder dictionaries

Consistent to how we visualize decoder dictionary atoms in Section IV-F, we
visualize atoms from/ encoder dictionaries $\boldsymbol{D}_{\textrm{enc}}$ and
$\widetilde{\boldsymbol{D}}_{\textrm{enc}}$ in Figure 8 and Figure 9.

Figure 8: Atoms in a randomly initialized (panel a) and learned
$\boldsymbol{D}_{\textrm{enc}}$ based on the derain dataset (panel b),
LoDoPab-CT (panel c) and fastMRI (panel d) dataset. The visualization setup is
identical to Figure 6 in the main text.

Figure 9: Atoms in a randomly initialized (panel a) and learned
$\widetilde{\boldsymbol{D}}_{\textrm{enc}}$ based on the derain dataset (panel
b), LoDoPab-CT (panel c) and fastMRI (panel d) dataset. The visualization
setup is identical to Figure 6 in the main text.

##  Appendix F Visualizing the representation of U-Nets

Similar to how we visualize dictionary atoms in Section IV-F, we visualize
prototypical images that U-Nets synthesize through its decoder branch
$f^{\text{dec}}(\cdot,\boldsymbol{\gamma})$. Concretely, we first prepare a
set of indicator codes corresponding to different spatial resolutions as
described in Section IV-F. We then feed each indicator code
$\boldsymbol{\delta}$ into the decoder branch of a U-Net to yield
$f^{\text{dec}}(\boldsymbol{\delta},\boldsymbol{\gamma})$. Due to additive
biases and batchnorm modules of the U-Net, the synthesized output
$f^{\text{dec}}(\boldsymbol{\delta},\boldsymbol{\gamma})$ has the same support
as the full image. To focus on the region influenced by the indicator code, we
thus display the support of
$f^{\text{dec}}(\boldsymbol{\delta},\boldsymbol{\gamma})-f^{\text{dec}}(\mathbf{0},\boldsymbol{\gamma})$,
where $\mathbf{0}$ is an all-zero tensor; the purpose of this subtract is to
offset those image values solely influenced by batchnorm and additive biases
but not by the indicator code. These synthesized results are visualized in
Figure 10.

As it can be seen, compared to the randomly initialized U-Net (Figure 10a),
the representations of learned U-Nets (Figure 10b and c) are organized in a
more structured way at each scale. Compared to the (linear) representations
learned by the MUSC, the (nonlinear) U-Net atoms much less resemble the
classical oriented multiresolution systems such as curvelets.

Figure 10: Indicator-code-induced images of U-Nets with randomly initialized
(left panel) and learned parameters based on the LoDoPab-CT (middle panel) and
fastMRI (right panel). For all panels, each row corresponds to a support size
(denoted by $s$) of the receptive field of a convolution layer. For the
visualization purpose, these images are normalized into the range $[-1,1]$.

##  Appendix G Super-resolution

CSC models achieve competitive performance in image super-resolution [71, 72].
We train an out-of-the-box MUSC for this task to study the sparsity patterns
of its learned representation in Section IV-F. We follow the protocol in
earlier work [73, 74, 75] to train MUSC on the DIV2K dataset [76]. Low-
resolution images are prepared by downscaling high-resolution images by a
factor of four. We used bicubic interpolated images up-scaled from low-
resolution images as model inputs and high-resolution images as model targets.
The trained models were evaluated on standard datasets including Set-14 [77],
Set-5 [78], B-100 [79], and Urban-100 [80]. Table IX shows the performance of
the trained models. We observe that single-scale CSC has an edge of all other
models, indicating the limited usefulness of large filters and the benefits of
trading them off for a large number of small-support convolutional channels,
as also discussed in the main text.

| Set-5 | Set-14 | B-100 | Urban-100
---|---|---|---|---
ScSR [81] | 29.07 | 26.40 | 26.61 | 24.02
A+ [82] | 30.17 | 26.94 | 26.81 | 24.29
CSC [71] (results taken from [72]) | 31.82 | 28.29 | 27.44 | 25.59
U-Net | 30.52 | 27.43 | 27.04 | 24.68
MUSC (ours) | 31.30 | 28.01 | 27.20 | 25.08
Table IX: Performance on the single-image super-resolution task with a scaling
factor of four.