A signal processing interpretation of noise-reduction convolutional
  neural networks
Encoding-decoding CNNs play a central role in data-driven noise reduction and
can be found within numerous deep-learning algorithms. However, the development
of these CNN architectures is often done in ad-hoc fashion and theoretical
underpinnings for important design choices is generally lacking. Up to this
moment there are different existing relevant works that strive to explain the
internal operation of these CNNs. Still, these ideas are either scattered
and/or may require significant expertise to be accessible for a bigger
audience. In order to open up this exciting field, this article builds
intuition on the theory of deep convolutional framelets and explains diverse ED
CNN architectures in a unified theoretical framework. By connecting basic
principles from signal processing to the field of deep learning, this
self-contained material offers significant guidance for designing robust and
efficient novel CNN architectures.
Luis A. Zavala-Mondragón, Peter H. N. de With, Fons van der Sommen
# A signal processing interpretation of noise-reduction convolutional neural networks

###### Abstract

Encoding-decoding CNNs play a central role in data-driven noise reduction and can be found within numerous deep-learning algorithms. However, the development of these CNN architectures is often done in _ad-hoc_ fashion and theoretical underpinnings for important design choices is generally lacking. Up to this moment there are different existing relevant works that strive to explain the internal operation of these CNNs. Still, these ideas are either scattered and/or may require significant expertise to be accessible for a bigger audience. In order to open up this exciting field, this article builds intuition on the theory of deep convolutional framelets and explains diverse ED CNN architectures in a unified theoretical framework. By connecting basic principles from signal processing to the field of deep learning, this self-contained material offers significant guidance for designing robust and efficient novel CNN architectures.

 Denoising, convolutional neural network, encoding-decoding. +
Footnote †: publication: 0000-0000/00000.0 © 2021 IEEE

## I Introduction

A well-known image-processing application is noise/artifact reduction of images, which is consists in estimating a noise/artifact-free signal out of a noisy observation. In order to achieve this, conventional signal processing algorithms often employ explicit assumptions on the signal and noise characteristics, which has resulted in well-known algorithms such as wavelet shrinkage [1], sparse dictionaries [2], total-variation minimization [3] and low-rank approximation [4]. With the advent of deep learning techniques, signal processing algorithms applied to image denoising, have been regularly outperformed and increasingly replaced by encoding-decoding convolutional neural networks (CNNs).

In this article, rather than conventional signal processing algorithms, we focus on the so-called encoding-decoding CNNs. These models contain an _encoder_ which maps the input to a multi-channel/redundant representations and a _decoder_, which maps the encoded signal back to the original domain. In both, the encoder and decoder, sparsifying non-linearities which suppress parts of the signal are applied. In contrast with conventional signal processing algorithms, encoding-decoding CNNs are often presented as a solution which does not make explicit assumptions on the signal and noise. For example, in supervised algorithms, an encoding-decoding CNN _learns_ the optimal parameters to filter the signal from a set of paired examples of noise/artifact-free images and images contaminated with noise/artifacts [5, 6, 7], which highly simplifies the solution of the noise-reduction problems, since this circumvents the use of explicit modeling of the signal and noise. Furthermore, the good performance and simple use of encoder-decoder CNNs/autoencoders have enabled additional data-driven noise-reduction algorithms, where CNNs are embedded as part of a larger system. Examples of such approaches are unsupervised noise-reduction [8], denoising based on generative adversarial networks [9]. Besides this, smoothness in signals can be obtained also by advanced regularization using CNNs, e.g. by exploiting data-driven model-based iterative reconstruction [10].

Despite of the impressive noise-reduction performance and flexibility of encoding-decoding convolutional neural networks, these models have also downsides that should be considered. First, the complexity and heuristic nature of such designs often offers restricted understanding of the internal operation of such architectures [11]. Second, the training and deployment of CNNs requires specialized hardware and the use of significant computational resources. Third and final, the restricted understanding of the signal modeling in encoding-decoding CNNs does not clearly reveal the limitations of such models and, consequently, it is not obvious how to overcome these problems.

In order to overcome the limitations of encoding-decoding CNNs, new research has tackled the lack of explainability of these models by acknowledging the similarity of the building blocks of encoding-decoding CNNs applied to image noise reduction and the elements of well-known signal processing algorithms, such as wavelet decomposition, low-rank approximation [12, 13, 14], variational methods [15], lower-dimensional manifolds [8] and convolutional sparse coding [16]. Furthermore, practical works based on shrinkage-based CNNs inspired in well-established wavelet shrinkage algorithms has further deepened the connections between signal processing and CNNs [17, 18]. This unified treatment of signal processing-inspired CNNs has resulted in more explainable [6, 8], better performing [6] and more memory-efficient designs [19].

This article has three main objectives. First, to summarize the diverse explanations of the components of encoding-decoding convolutional neural networks applied to image noise reduction based on the concept of deep convolutional framelets [12] and on elementary signal processing concepts.
Both aspects are considered with the aim of achieving an in-depth understanding of the internal operation of encoding-decoding CNNs and to show that the design choices have _implicit_ assumptions about the signal behavior inside the CNN. A second objective is to offer practitioners tools for optimizing their CNN designs with signal processing concepts. Third and final, the aim is to show practical use cases, where existing CNNs are analyzed in a unified framework, thereby enabling a better comparison of different designs by making their internal operation explicitly visible. Our analysis are based on existing works [12, 6, 20], who analyzed CNNs where the non-linearities are ignored. In this article, we overcome this limitation and present a complete analysis including the non-linear activations, which reveals important assumptions implicit in the analyzed models.

The structure of this article is as follows. Section II introduces the notation that is used in this text. Section III describes the signal model and the architecture of encoding-decoding networks. Afterwards, Section IV addresses fundamental aspects of signal processing, such as singular value decomposition, low-rank approximation, framelets, as well as the estimation of signals in the framelet domain. All the concepts of Sections III and IV converge in Section V, where the encoding-decoding CNNs are interpreted in terms of a data-driven low-rank approximation and of wavelet shrinkage. Afterwards, based on the learnings from Section V, Section VI shows the analysis of diverse architectures from a signal processing perspective and under a set of explicit assumptions. Afterwards, Section VII explores if some of the theoretical properties exposed here are related to trained models. Based on the diverse described models and theoretical operation of CNNs, Section VIII addresses a design criterion which can be used to design or choose new models and briefly describes the state-of-the art for noise reduction with CNNs. Finally, Section IX elaborates concluding remarks and discusses diverse elements that have not yet been (widely) explored by current CNN designs.

## II **Notation**

Convolutional neural networks (CNNs) are composed by basic elements, such as convolution, activation and down/up-sampling layers. In order to achieve better clarity in the explanations given in this paper, we define the mathematical notation to represent the basic operations of CNNs. Part of the definitions presented here are based on the work of Zavala _et al._[19].

In the following, a scalar is represented by a lower-case letter (e.g. \(a\)), while a vector is represented by an underlined lower-case letter (e.g. \(\underline{a}\)). Furthermore, a matrix, such as an image or convolution mask, is represented by a boldface lowercase letter (e.g. variables \(\mathbf{x}\) and \(\mathbf{y}\)). Finally, a tensor is defined by a boldface uppercase letter. For example, the two arbitrary tensors \(\mathbf{A}\) and \(\mathbf{Q}\) are defined by

\[\mathbf{A}=\begin{pmatrix}\mathbf{a}_{0}^{0}&\dots&\mathbf{a}_{N_{\text{C}}- 1}^{0}\\ \vdots&\ddots&\vdots\\ \mathbf{a}_{0}^{N_{\text{R}}-1}&\dots&\mathbf{a}_{N_{\text{C}}-1}^{N_{\text{R} }-1}\end{pmatrix},\mathbf{Q}=\begin{pmatrix}\mathbf{q}^{0}\\ \vdots\\ \mathbf{q}^{N_{\text{R}}-1}\end{pmatrix}. \tag{1}\]

Here, entries \(\mathbf{a}_{c}^{r}\) and \(\mathbf{q}^{r}\) represent two-dimensional arrays (matrices). Since the defined tensors are used in the context of CNNs, matrices \(\mathbf{a}_{c}^{r}\) and \(\mathbf{q}^{r}\) are learned filters, which have dimensions \((N_{\text{V}}\times N_{\text{H}})\), where \(N_{\text{V}}\) and \(N_{\text{H}}\) denote the filter dimensions in the vertical and horizontal directions, respectively. Finally, we define the total tensor dimension of \(\mathbf{A}\) and \(\mathbf{Q}\) by \((N_{\text{C}}\times N_{\text{R}}\times N_{\text{V}}\times N_{\text{H}})\) and \((N_{\text{R}}\times 1\times N_{\text{V}}\times N_{\text{H}})\), where \(N_{\text{R}}\) and \(N_{\text{C}}\) are the number of row and column entries, respectively. If the tensor \(\mathbf{A}\) contains the convolution weights in a CNN, the row-entry dimensions represent the input number of channels to a layer, while the number of column elements denotes the number of output channels.

Having defined the notation for the variables, we focus on a few relevant operators. First, the transpose of a tensor \((\cdot)^{\intercal}\), expressed by

\[\mathbf{Q}^{\intercal}=\begin{pmatrix}\mathbf{q}^{0}&\dots&\mathbf{q}^{N_{ \text{R}}-1}\end{pmatrix}. \tag{2}\]

Furthermore, the convolution of two tensors is written as \(\mathbf{AQ}\) and specified by

\[\mathbf{AQ}=\begin{pmatrix}\sum_{r=0}^{N_{\text{R}}-1}\mathbf{a}_{r}^{0}* \mathbf{q}^{r}\\ \vdots\\ \sum_{r=0}^{N_{\text{R}}-1}\mathbf{a}_{r}^{N_{\text{R}}-1}*\mathbf{q}^{r} \end{pmatrix}. \tag{3}\]

Here, the symbol \(*\) defines the convolution between two matrices (images).

In this paper, images which are 2D arrays (matrices), are often convolved with 4D tensors. When this operation is performed, images are considered to have dimensions \((1\times 1\times N_{\text{V}}\times N_{\text{H}})\). In addition, in this paper matrix \(\mathbf{I}\) is the identity signal for the convolution operator, which for a 2D image is the Kronecker delta/discrete impulse (an image with a single non-zero pixel with unity amplitude at the center of the image). Furthermore, we indicate that variables in the _decoding path_ of a CNN are distinguished with a tilde (e.g. \(\tilde{\mathbf{K}}\), \(\tilde{\underline{b}}\)).

Additional symbols that will be used throughout the article are the down- and up-sampling operations by a factor \(s\), which are denoted by \(f_{(s\downarrow)}(\cdot)\) for down-sampling and for up-sampling \(f_{(s\uparrow)}(\cdot)\). In this paper, both operations are defined in the same way as in multi-rate filter banks. For example, consider the signal

\[\underline{x}=\begin{pmatrix}1,&2,&3,&4,&5,&6,&7,&8,&9,&10\end{pmatrix}. \tag{4}\]

If we apply the down-sampling operator to \(\underline{x}\) by a factor 2, this results into

\[\underline{z}=f_{(2\downarrow)}(\underline{x})=\begin{pmatrix}1,&3,&5,&7,&9 \end{pmatrix}, \tag{5}\]

where \(\underline{z}\) is the down-sampled version of \(\underline{x}\). Conversely, the result of applying the up-sample operator \(f_{(2\uparrow)}(\cdot)\) gives as result

\[f_{(2\uparrow)}(\underline{z})=\begin{pmatrix}1,&0,&3,&0,&5,&0,&7,&0,&9,&0 \end{pmatrix}. \tag{6}\]

Additional operators used in the article are the _rectified linear unit_ (ReLU), the _shrinkage/thresholding_ and the _clipping_, which are represented by \((\cdot)_{+}\), \(\tau_{(\cdot)}(\cdot)\) and \(\mathcal{C}_{(\cdot)}(\cdot)\), respectively.

For better clarity, the most important symbols used in this article are summarized in Table I. In addition, the graphical
representations of some of the symbols that will be used to graphically describe CNNs are shown in Fig. 1.

## III **Encoding-decoding CNNs**

### **Signal model and noise reduction configurations**

In noise-reduction applications, the common additive signal model is defined by

\[\mathbf{y}=\mathbf{x}+\boldsymbol{\eta}, \tag{7}\]

where the observed signal \(\mathbf{y}\) is the result of contaminating a noiseless image \(\mathbf{x}\) with additive noise \(\boldsymbol{\eta}\). Assume that the noiseless signal \(\mathbf{x}\) is to be estimated from the noisy observation \(\mathbf{y}\). In deep learning applications, this is often achieved by models with the form

\[\hat{\mathbf{x}}=G(\mathbf{y}). \tag{8}\]

Here, \(G(\cdot)\) is a generic encoding-decoding CNN. We refer to this form of noise reduction as _non-residual_. Alternatively, it is possible to find \(\hat{\mathbf{x}}\) by training \(G(\cdot)\) to estimate the noise component \(\hat{\boldsymbol{\eta}}\), and subtract it from the noisy image \(\mathbf{y}\) to estimate the noiseless image \(\hat{\mathbf{x}}\), or equivalently

\[\hat{\mathbf{x}}=\mathbf{y}-G(\mathbf{y}). \tag{9}\]

This model is referred to as _residual_[5, 21, 7], because the output of the network is subtracted from its input. For reference, Fig. 2 portrays the difference of the placement of the encoding-decoding structure in residual and non-residual configurations.

### **Encoding-decoding CNNs**

Encoding-decoding (convolutional) neural networks are rooted in the techniques for data-dimensionality reduction and unsupervised feature extraction, where a given signal is mapped to an alternative space via a non-linear transformation. This space should have properties which are somehow attractive for the considered task. For example, for dimensionality reduction, the alternative space should be lower-dimensional than the original input. In this article, we are interested in models that are useful for noise-reduction applications. Specifically, this manuscript addresses models that are referred to as encoding-decoding _convolutional_ neural networks such as the model by Ranzato _et al._[22], in which the encoder uses convolution filters to produce multi-channel/redundant representations, in which sparsifying non-linearities are applied. The sparsified signal is later mapped back to the original representation. It should be noted that despite the fact that the origins of the encoding-decoding CNNs are linked to feature extraction, this type of architecture quickly showed to be useful for other applications such as noise reduction, which is the topic of this article. For the rest of this manuscript whenever we mention an encoding-decoding CNN, we are referring to a design which follows the same basic principles of Ranzato's design.

It can be observed that encoding-decoding CNNs are constituted of three main parts. (1) The **encoder**, which maps the incoming image to a representation with more image channels with a convolution layer. Every channel of the resulting redundant representation contains a fraction of the content of the original signal. It should be noted that the encoder often (but not necessarily) decreases the resolution of the higher dimensional representation, to enable multi-resolution processing and to decrease the memory requirements of the design. (2) The **decoder**, which maps the multi-channel representation back to the original space. (3) The **non-linearities**, which suppress specific parts of the signal. In summary, the most basic encoding-decoding step in a CNN \(G(\cdot)\) is expressed by

\[\boxed{G(\mathbf{y})=G_{\text{dec}}(G_{\text{enc}}(\mathbf{y}))}\enspace, \tag{10}\]

Fig. 1: Symbols used for the schematic representations of the CNNs addressed in this article.

Fig. 2: Residual and non-residual network configurations. Note that the main difference between both designs is the global skip connection occurring in the residual structure. Still, it can be observed that the network \(G(\cdot)\) may contain skip connections _internally_.
where \(G_{\text{enc}}(\cdot)\) is the encoder, which is generally defined by

\[\begin{split}\mathbf{C}_{0}=& E_{0}(\mathbf{y}),\\ \mathbf{C}_{1}=& E_{1}(\mathbf{C}_{0}),\\ \mathbf{C}_{2}=& E_{2}(\mathbf{C}_{1}),\\ \vdots&\\ \mathbf{C}_{N-1}=& E_{N-1}(\mathbf{C}_{N-2}),\\ G_{\text{enc}}(\mathbf{y})=&\mathbf{C}_{N-1}. \end{split} \tag{11}\]

Here, \(\mathbf{C}_{n}\) represents the code generated by the \(n\)-th encoding \(E_{n}(\cdot)\), which can be expressed by the equation

\[\boxed{\mathbf{C}_{n}=E_{n}(\mathbf{C}_{n-1})=f_{(s\downarrow)}\big{(}A_{( \underline{b}_{n-1})}(\mathbf{K}_{n-1}\mathbf{C}_{n-1})\big{)}}\enspace. \tag{12}\]

Here, the function \(A(\cdot)_{(\cdot)}\) is a generic activation used in the encoder and \(f_{(s\downarrow)}(\cdot)\) is a down-sampling function by factor \(s\). Complementary to the encoder, the decoder network maps the multi-channel sparse signal back to the original domain. Here, we define the decoder by

\[\begin{split}\tilde{\mathbf{C}}_{N-2}=& D_{N-1}( \mathbf{C}_{N-1}),\\ \vdots&\\ \tilde{\mathbf{C}}_{1}=& D_{2}(\tilde{\mathbf{C}}_{2 }),\\ \tilde{\mathbf{C}}_{0}=& D_{1}(\tilde{\mathbf{C}}_{1 }),\\ G(\mathbf{y})=& D_{0}(\tilde{\mathbf{C}}_{0}),\end{split} \tag{13}\]

where \(\tilde{\mathbf{C}}_{n}\) is the \(n-\)th decoded signal, which is produced by the \(n\)-th decoder layer, yielding the general expression:

\[\boxed{\tilde{\mathbf{C}}_{n-1}=D_{n}(\tilde{\mathbf{C}}_{n})=\tilde{A}_{( \underline{\underline{\mathbf{j}}})}\big{(}\tilde{\mathbf{K}}_{n}^{\intercal }f_{(s\uparrow)}(\tilde{\mathbf{C}}_{n})\big{)}}\enspace. \tag{14}\]

In the above, \(\tilde{A}(\cdot)_{(\cdot)}\) is the activation function used in the decoder and \(f_{(s\uparrow)}(\cdot)\) is an up-sampling function of factor \(s\).

An important remark is that the encoder-decoder CNN does not always contain down/up-sampling layers in which case, the decimation factor \(s\) is unity, which causes \(f_{(1\uparrow)}(\mathbf{x})=f_{(1\downarrow)}(\mathbf{x})=\mathbf{x}\) for any matrix \(\mathbf{x}\). Furthermore, it should be noted also that we assume that the number of channels of the code \(\mathbf{C}_{N}\) is always larger than the previous one \(\mathbf{C}_{N-1}\). Furthermore, it should be noted that a single encoder layer \(E_{n}(\cdot)\) and its corresponding decoder layer \(D_{n}(\cdot)\) can be considered a single-layer encoder-decoder network/pair.

For this article, the encoding convolution filter for a given layer \(\mathbf{K}\) has dimensions \((N_{\text{o}}\times N_{\text{i}}\times N_{\text{h}}\times N_{\text{v}})\), where \(N_{\text{i}}\) and \(N_{\text{o}}\) are the number of input and output channels for a convolution layer, respectively. Similarly, \(N_{\text{h}}\) and \(N_{\text{v}}\) are the number of elements in the horizontal and vertical directions, respectively. Note that the encoder increases the number of channels of the signal (e.g. \(N_{\text{o}}>N_{\text{i}}\)), akin to Ranzatto's design [22]. Furthermore, it is assumed that the _decoder_ is symmetric in the number of channels to the _encoder_, therefore, the dimensions of the decoding convolution kernel \(\mathbf{K}^{\intercal}\) are \((N_{\text{i}}\times N_{\text{o}}\times N_{\text{h}}\times N_{\text{v}})\). The motivation of this symmetry is to emphasize the similarity between the signal processing and the CNN elements.

## IV **Signal processing fundamentals**

As shown by Ye _et al._[12], within encoding-decoding CNNs, the signal is treated akin to well-known sparse representations, where the coefficients used for the transformation are directly learned from the training data. Prior to addressing this important concept in more detail, relevant supporting concepts such as sparsity, sparse transformations and non-linear signal estimation in the wavelet domain are explained.

### **Sparsity**

A sparse image is a signal where most coefficients are small and the relatively few large coefficients capture most of the information [23]. This characteristic allows to discard low-amplitude components with relatively small perceptual changes. Hereby, the use of sparse signals is attractive for applications such as image compression, denoising and suppression of artifacts.

Despite the convenient characteristics of sparse signals, natural images are often non-sparse. Still, there are numerous transformations, which allow to map the signal to a sparse domain and which are analogous to the internal operations of CNNs. For example, _singular value decomposition_ factorizes the image in terms of two sets of orthogonal bases of which few basis pairs contain most of the energy of the image. An alternative transformation is based on _framelets_, where an image is decomposed in a multi-channel representation, whereby each resulting channel contains a fragment of the Fourier spectrum. In the remainder of this section we will address all of these representations in more detail.

### **Sparse signal representations**

#### Iv-B1 **Singular value decomposition (SVD) and low-rank approximation**

Assume that an image (patch) is represented by a matrix \(\mathbf{y}\) with dimensions \((N_{\text{r}}\times N_{\text{c}})\), where \(N_{\text{r}}\) and \(N_{\text{c}}\) are the number of rows and columns, respectively. Then, the singular value decomposition factorizes \(\mathbf{y}\) as

\[\mathbf{y}=\sum_{n=0}^{N_{\text{N}}-1}(\underline{u}_{n}\underline{v}_{n}^{ \intercal})\cdot\underline{\sigma}[n], \tag{15}\]

in which \(N_{\text{SV}}\) is the number of singular values, \(n\) is a scalar index, while \(\underline{u}_{n}\) and \(\underline{v}_{n}\) are the \(n^{\text{th}}\) left- and right-singular vectors, respectively. Furthermore, vector \(\underline{\sigma}\) contains the singular values and each of its entries \(\underline{\sigma}[n]\) is the weight assigned to every basis pair \(\underline{u}_{n}\), \(\underline{v}_{n}\). This means that the product \((\underline{u}_{n}\underline{v}_{n}^{\intercal})\) contributes more to the image content for higher values of \(\underline{\sigma}[n]\). It is customary that the singular values are ranked in descending order and the amplitudes of the singular values \(\underline{\sigma}\) are sparse, therefore \(\underline{\sigma}[0]\gg\underline{\sigma}[N_{\text{SV}}-1]\). The reason for this sparsity is because image (patches) intrinsically have high correlation. For example, many images contain repetitive patterns (e.g. a wall with bricks, a fence, the tiles of a rooftrop or the stripes of a zebra) or uniform regions (for example, the sky, the skin of a person). This means an image patch may contain only a few linearly independent vectors that describe most of the image content. Consequently, a higher weight is assigned to such image bases.
Given that the amplitudes of the singular values of \(\mathbf{y}\) in SVD are sparse, it is possible approximate \(\hat{\mathbf{y}}\) with only a few bases \((\underline{u}_{n}\underline{v}_{n}^{\intercal})\). Note that this procedure reduces the rank of signal \(\mathbf{y}\) and hence it is known as _low-rank approximation_. This process is equivalent to

\[\hat{\mathbf{y}}=\sum_{n=0}^{N_{\mathrm{LR}}-1}(\underline{u}_{n}\underline{v }_{n}^{\intercal})\cdot\sigma[n], \tag{16}\]

where \(N_{\mathrm{SV}}>N_{\mathrm{LR}}\). Note that this effectively cancels the product \((\underline{u}_{n}\underline{v}_{n}^{\intercal})\) where the weight given by \(\sigma[n]\) is low. Alternatively, it is possible to assign a weight of zero to the product \((\underline{u}_{n}\underline{v}_{n}^{\intercal})\) for \(n\geq N_{\mathrm{LR}}\).

The low-rank representation of a matrix is desirable for diverse applications among which we can find image denoising. The motivation for using low-rank approximation for this application results from the fact that -as mentioned earlier-natural images are considered low-rank due to the strong spatial correlation between pixels, whereas noise is high-rank (it is spatially uncorrelated). In consequence, reducing the rank/number of singular values decreases the presence of noise, while still providing a good approximation of the noise-free signal, as exemplified in Fig. 3.

#### Iii-B2 **Framelets**

Just as SVD, framelets are also commonly used for image processing. In a nutshell, a framelet transform is a signal representation that factorizes/decomposes an arbitrary signal into multiple bands/channels. Each of these channels contain a segment of the energy of the original signal. In image and signal processing, the framelet bands are the result of convolving the analyzed signal with a group of discrete filters that have finite length/support. In this article, the most important characteristic that the filters of the framelet transform should comply with, is that the bands they generate capture _all_ the energy contained on the input to the decomposition. This is important to avoid the loss of information of the decomposed signal. In this text, we refer to framelets that comply with the previous characteristics as _tight_ framelets and the following paragraphs will describe this property in more detail.

In its decimated version, the framelet decomposition for _tight_ frames is represented by

\[\mathbf{Y}_{\text{frame}}=f_{(2\downarrow)}(\mathbf{F}\mathbf{y}), \tag{17}\]

in which \(\mathbf{Y}_{\text{frame}}\) is the decomposed signal and \(\mathbf{F}\) is the framelet basis (tensor). Note that the signal \(\mathbf{Y}_{\text{frame}}\) has more channels than \(\mathbf{y}\). Furthermore, the original signal \(\mathbf{y}\) is recovered from \(\mathbf{Y}_{\text{frame}}\) by

\[\mathbf{y}=\tilde{\mathbf{F}}^{\intercal}f_{(2\uparrow)}(\mathbf{Y}_{\text{ frame}})\cdot c. \tag{18}\]

Here, \(\tilde{\mathbf{F}}\) is the filter of the inverse framelet transform and \(c\) denotes an arbitrary constant. If \(c=1\) the framelet is _normalized_. Finally, note that the framelet transform can also be _undecimated_. This means that in undecimated representations, the down-sampling and up-sampling layers \(f_{(2\downarrow)}(\cdot)\) and \(f_{(2\uparrow)}(\cdot)\) are not used. An important property of the undecimated representation is that it is less prone to aliasing than its decimated counterpart, but more computationally expensive. Therefore, for efficiency reasons, the decimated framelet decomposition is often preferred over the undecimated representation. In summary, the decomposition and synthesis of the decimated framelet decomposition is represented by

\[\boxed{\mathbf{y}=\tilde{\mathbf{F}}^{\intercal}f_{(2\uparrow)}\big{(}f_{(2 \downarrow)}(\mathbf{F}\mathbf{y})\big{)}\cdot c}\enspace, \tag{19}\]

while for the undecimated framelet it holds that

\[\boxed{\mathbf{y}=\tilde{\mathbf{F}}^{\intercal}(\mathbf{F}\mathbf{y})\cdot c }\enspace. \tag{20}\]

A notable normalized framelet is the _discrete wavelet transform_ (DWT), where variables \(\mathbf{F}\) and \(\tilde{\mathbf{F}}\) are replaced by tensors \(\mathbf{W}=\big{(}\mathbf{w}_{\mathrm{LL}},\mathbf{w}_{\mathrm{LH}},\mathbf{ w}_{\mathrm{HH}},\mathbf{w}_{\mathrm{HH}}\big{)}\) and \(\tilde{\mathbf{W}}=\big{(}\tilde{\mathbf{w}}_{\mathrm{LL}},\tilde{\mathbf{w}}_ {\mathrm{LH}},\tilde{\mathbf{w}}_{\mathrm{HL}},\tilde{\mathbf{w}}_{\mathrm{HH}}\big{)}\), respectively. Here, \(\mathbf{w}_{\mathrm{LL}}\) is the filter for the low-frequency band, while \(\mathbf{w}_{\mathrm{LH}}\), \(\mathbf{w}_{\mathrm{HL}}\), \(\mathbf{w}_{\mathrm{HH}}\) are the filters used to extract the detail in the horizontal, vertical and diagonal directions, respectively. Finally, \(\tilde{\mathbf{w}}_{\mathrm{LH}}\,\tilde{\mathbf{w}}_{\mathrm{LH}},\tilde{ \mathbf{w}}_{\mathrm{HL}},\tilde{\mathbf{w}}_{\mathrm{HH}}\) are the filters of the inverse decimated DWT.

In order to understand the DWT more intuitively, Fig. 4 shows the decimated framelet decomposition using the filters of the discrete wavelet transform. Note that the convolution \(\mathbf{W}\mathbf{y}\) results in a four-channel signal, where each channel contains only a fraction of the spectrum of image \(\mathbf{y}\). This allows to down-sample each channel with minimal aliasing. Furthermore, to recover the original signal, each individual channel is up-sampled, thereby introducing aliasing, which is then removed by the filters of the inverse transform. Finally, all the channels are added and the original signal is recovered.

Analogous to the low-rank approximation, in framelets, the reduction of noise is achieved by setting noisy components to zero. These components are typically assumed to have low-amplitude compared with the amplitude of the sparse signal, as expressed by

\[\hat{\mathbf{y}}=\tilde{\mathbf{F}}^{\intercal}f_{(2\uparrow)}\big{(}\tau_{( \underline{t})}\big{(}f_{(2\downarrow)}(\mathbf{F}\mathbf{y})\big{)}\big{)} \cdot c, \tag{21}\]

where \(\tau_{\underline{t}}(\cdot)\) is a generic thresholding/shrinkage function, which sets each of the pixels in \(f_{(2\downarrow)}(\mathbf{F}\mathbf{y})\) to zero when values are lower than the threshold level \(\underline{t}\).

Fig. 3: SVD reconstruction of clean and corrupted images with a different number of singular values. Note that the reconstruction of the clean image with 8 or 32 singular values (\(N_{\mathrm{SV}}=8\) or \(N_{\mathrm{SV}}=32\), respectively) yields to reconstructions indistinguishable from the original image. This contrasts with their noisy counterparts, where \(N_{\mathrm{SV}}=8\) reconstructs a smoother image in which the noise is attenuated, while \(N_{\mathrm{SV}}=32\) reconstructs the noise texture perfectly.
### _Nonlinear signal estimation in the framelet domain_

As mentioned in Section IV-B2, framelets decompose a given image \(\mathbf{y}\) by convolving it with a tensor \(\mathbf{F}\). Note that many of the filters that compose \(\mathbf{F}\) have a high-pass nature. Images often contain approximately uniform regions in which the variation is low, therefore, convolving a signal \(\mathbf{y}\) with a high-pass filter \(\mathbf{f}_{\mathbf{h}}\) -where \(\mathbf{f}_{\mathbf{h}}\in\mathbf{F}\)- produces the sparse detail band \(\mathbf{d}=\mathbf{f}_{\mathbf{h}}*\mathbf{y}\) in which uniform regions have low amplitudes, while transitions i.e. edges contain most of the energy of the bands.

Assuming a model in which a single pixel \(d\in\mathbf{d}\) is observed, which is contaminated with additive noise \(\eta\). Then, the resulting observed pixel \(z\) is defined by

\[z=d+\eta. \tag{22}\]

In order to recover the noiseless pixel \(d\) from observation \(z\), it is possible to use the point-_maximum a posteriori_ (MAP) estimate [1, 24], defined by the maximization problem

\[\hat{d}=\operatorname*{arg\,max}_{d}\bigl{[}\ln\big{(}P(d|z)\big{)}\bigr{]}. \tag{23}\]

Here, the log-posterior \(\ln\big{(}P(d|z)\big{)}\) is defined by

\[\ln\big{(}P(d|z)\big{)}=\ln\big{(}P(z|d)\big{)}+\ln\big{(}P(d)\big{)}, \tag{24}\]

where the conditional probability density function (PDF) \(P(z|d)\) expresses the noise distribution, which is often assumed Gaussian and is defined by

\[P(z|d)\propto\exp\bigg{(}-\frac{(z-d)^{2}}{2\sigma_{\eta}^{2}}\bigg{)}. \tag{25}\]

Here, \(\sigma_{\eta}^{2}\) is the noise variance. Furthermore, as prior probability, it is assumed that the distribution of \(P(d)\) corresponds to a Laplacian distribution, which has been used in wavelet-based denoising [1]. Therefore, \(P(d)\) is mathematically described by

\[P(d)\propto\exp\bigg{(}-\frac{|d|}{\sigma_{d}}\bigg{)}, \tag{26}\]

where \(\sigma_{d}\) is the dispersion measure of the Laplace distribution. For reference, Fig. 5 portrays an example of a Gaussian and a Laplacian PDF. Note that the Laplacian distribution has a higher probability of zero elements to occur than the Gaussian distribution for the same standard deviation.

Finally, substituting Eq. (25) and Eq. (26) in Eq (24) results in

\[\ln\big{(}P(d|z)\big{)}\propto-\frac{(z-d)^{2}}{2\sigma_{\eta}^{2}}-\frac{|d |}{\sigma_{d}}. \tag{27}\]

Fig. 4: 2D spectrum analysis of the decimated discrete framelet decomposition and reconstruction. In the figure, function \(\mathcal{F}\{\cdot\}\) stands for the amplitude Fourier spectrum of the input argument. The yellow squares indicate a region in the low-frequency area of the Fourier spectrum, while the orange, purple and blue squares indicate the high-pass/detail bands. For these images, ideal orthogonal bases are assumed. Note that the forward transform is composed by two steps. First, the signal is convolved with the wavelet basis \((\mathbf{W}\mathbf{y})\). Afterwards, down-sampling is applied to the signal (\(f_{(2i)}(\mathbf{W}\mathbf{y})\)). During the inverse transformation, the signal is up-sampled by inserting zeros between each sample (\(f_{(21)}(f_{(2i)}(\mathbf{W}\mathbf{y}))\)), which causes spatial aliasing (dashed blocks). Finally, the spatial aliasing is removed by the inverse transform filter \(\tilde{\mathbf{W}}\) and all the channels are added \((\mathbf{W}^{\intercal}f_{(21)}(f_{(2i)}(\mathbf{W}\mathbf{y})))\).

Fig. 5: Probability density function for Gaussian (left) and Laplacian (right) distributions.
In the above, maximizing \(d\) in \(\ln\left(P(d|z)\right)\) with the first derivative criterion -in an (un) constrained way- leads to two common activations in noise-reduction CNNs: the ReLU and the soft-shrinkage function. Furthermore, the solution also can be used to derive the so-called _clipping_ function, which is useful in residual networks.

For reference and further understanding, Fig. (6) portrays the elements composing the noise model of Eq. (22), the signal transfer characteristics of the ReLU, soft-shrinkage and clipping functions, and the effect that these functions have on the signal of the observed noisy detail band \(z\).

#### Iii-B1 **Rectified linear unit (ReLU)**

If Eq. (27) is solved for \(d\) while constraining the estimator to be positive, the noiseless estimate \(\hat{d}\) becomes

\[\boxed{\hat{d}=(z-t)_{+}}\enspace, \tag{28}\]

which is also expressed by

\[(z-t)_{+}=\begin{cases}z-t,&\text{if }z\geq t,\\ \hskip 14.226378pt0,&\text{if }t>z.\end{cases} \tag{29}\]

Here, the threshold level is defined by

\[t=\sigma_{\eta}^{2}/\sigma_{d}. \tag{30}\]

Note that this estimator cancels the negative and low-amplitude elements of \(d\) lower than the magnitude of the threshold level \(t\). For example, if the signal content on the feature map is low, then \(\sigma_{d}\to 0\). In such case, \(t\rightarrow+\infty\) and consequently \(\hat{d}\to 0\). This means that the channel is suppressed. Alternatively, if the feature map has strong signal presence i.e. \(\sigma_{d}\rightarrow\infty\), consequently \(t\to 0\) and then \(\hat{d}\rightarrow(z)_{+}\).

A final remark is made on the modeling of functions of a CNN. It should be noted that the estimator of Eq. (28) is analogous to the activation function of a CNN, known as _rectified linear unit_ (ReLU). However, in a CNN the value of \(t\) would be the bias \(b\) learned from the training data.

#### Iii-B2 **Soft-shrinkage/thresholding**

If Eq. (27) is maximized in an unconstrained way, the estimate \(\hat{d}\) is

\[\boxed{\hat{d}=\tau_{(t)}^{\text{Soft}}(z)=(z-t)_{+}-(z-t)_{+}}\enspace. \tag{31}\]

Here, \(\tau_{(t)}^{\text{Soft}}(\cdot)\) denotes the soft-shrinkage/-thresholding function, which is often also written in the form

\[\tau_{(t)}^{\text{Soft}}(z)=\begin{cases}z+t,&\text{if }z\geq t,\\ \hskip 14.226378pt0,&\text{if }t>z\geq-t,\\ z-t,&\text{if }-t>z.\end{cases} \tag{32}\]

It can be observed that the soft-threshold enforces the low-amplitude components whose magnitude is lower than the magnitude threshold level \(t\) to zero. In this case, \(t\) is also defined by Eq. (30). It should be noted that the soft-shrinkage estimator can also be obtained from a variational perspective [25]. Finally, it can be observed that the soft-shrinkage is the superposition of two ReLU functions, which has been pointed out by Fan _et al._[18].

#### Iii-B3 **Soft clipping**

In Section IV-C1 and Section IV-C2, the estimate \(\hat{d}\) is obtained directly from the noisy observation \(z\). Alternatively, it is possible to estimate the noise \(\eta\) and subtract it from \(z\) akin to the residual CNNs represented by Eq. (9). This can be achieved by solving the model

\[\boxed{\hat{\eta}=z-\hat{d}=z-\tau_{(t)}^{\text{Soft}}(z)}\enspace, \tag{33}\]

which is equivalent to

\[\hat{\eta}=\mathcal{C}_{(t)}^{\text{Soft}}(z)=z-((z-t)_{+}-(-z-t)_{+}), \tag{34}\]

where \(\mathcal{C}_{(t)}^{\text{Soft}}(\cdot)\) is the _soft clipping_ function. Note that this function also can be expressed by

\[\mathcal{C}_{(t)}^{\text{Soft}}(z)=\begin{cases}\enspace t,&\text{if }z\geq t,\\ \enspace z,&\text{if }t\geq z>-t,\\ -t,&\text{if }-t\geq z.\end{cases} \tag{35}\]

#### Iii-B4 **Other thresholding layers**

One of the main drawbacks of the soft-threshold activation is that it is a biased estimator. This limitation has been addressed by the hard and semi-hard thresholds, which are (asymptotically) unbiased estimators for large input values. In this section, we focus solely on the semi-hard threshold and avoid the hard variant, because is discontinuous and, therefore, not suited for models that rely on gradient-based optimization, such as CNNs.

Among the semi-hard thresholds, two notable examples are the _garrote shrink_ and the shrinkage functions generated by derivatives of Gaussians (DoG) [26, 19]. The garrote shrink function \(\tau_{(\cdot)}^{\text{Gar}}(\cdot)\) is defined by

\[\tau_{(t)}^{\text{Gar}}(z)=\frac{(z^{2}-t^{2})_{+}}{z}. \tag{36}\]

Fig. 6: Signals involved in the additive noise model, input/output transfer characteristics of activation layers and estimates produced by the activation layers when applied to the noise-contaminated signal. The first row shows the signals involved in the additive noise model. The second row depicts the output amplitude of activation functions with respect to the input amplitude. Finally, the last row depicts the application of the activation functions to the noisy observation \(z\).
Furthermore, an example of a shrinkage function based on the derivative of Gaussians is given by

\[\tau_{(t)}^{\text{DoG}}(z)=z-\mathcal{C}_{(t)}^{\text{DoG}}(z), \tag{37}\]

where the semi-hard clipping function with the derivative of Gaussians \(\mathcal{C}_{(\cdot)}^{\text{DoG}}(\cdot)\) is given by

\[\mathcal{C}_{(t)}^{\text{DoG}}(z)=z\cdot\exp\bigg{(}-\frac{z^{p}}{t^{p}}\bigg{)}, \tag{38}\]

in which \(p\) is an even number.

The garrote and semi-hard DoG shrinkage function are shown in Fig. 7, as well as their clipping counterparts. Note that the shrinkage functions approximate unity for \(|z|\gg t\). Therefore, they are asymptotically unbiased for large signal values.

The final thresholding function addressed in this section is the linear expansion of thresholds proposed by Blu and Luisier [26]. This technique combines multiple thresholding functions to improve the performance. This approach is known as _linear expansion of thresholds_ (LET) and it is defined by

\[\tau_{(\mathbf{\hat{L}})}^{\text{LET}}(z)=\sum_{n=0}^{N_{T}-1}a_{n}\cdot\tau_ {(t_{n})}(z), \tag{39}\]

where \(a_{n}\) is the weighting factor assigned to each threshold, where all weighting factor should add up to unity.

V **Bridging the gap between signal processing and CNNs: Deep convolutional framelets and shrinkage-based CNNs**

This section addresses the theoretical operation of noise-reduction convolutional neural networks based on ReLUs and shrinkage/thresholding functions. The first part of this section describes the theory of deep convolutional framelets [12], which is the most extensive study on the operation of encoding-decoding ReLU-based CNNs up to this moment. Afterwards, the section concentrates on the operation of networks which use shrinkage functions instead of ReLUs [17, 18, 19], with the aim of mimicking well-established denoising algorithms [1]. Finally, the last part of this section addresses the connections between both methods and additional links between convolutional neural networks and signal processing.

### **Theory of deep convolutional framelets**

The theory of deep convolutional framelets [12] describes the operation of encoding-decoding ReLU-based CNNs. Its most relevant contributions are as follows. (1) To establish the equivalence of framelets and the convolutional layers of CNNs. (2) The theory of deep convolutional framelets provides the conditions to preserve the signal integrity within a ReLU CNN. (3) Explain how ReLU and convolution layers reduce noise within an encoding-decoding CNN.

The similarity between framelets and the encoding and decoding convolutional filters can be observed when comparing Eqs. (12), (14) with Eqs. (17), (18), where it becomes visible that the convolution structure of encoding-decoding CNNs is analogous to the forward and inverse framelet decomposition.

Regarding the signal reconstruction characteristics, the theory of deep convolutional framelets [12] states the following. First, in order to be able to recover an arbitrary signal \(\mathbf{y}\in\mathbb{R}^{N}\), the number of output channels of a convolution layer with ReLU activation should _at least_ duplicate the number of input channels. Second, the encoding convolution kernel \(\mathbf{K}\) should be composed of pairs of filters with opposite phase. These two requirements ensure that any negative and positive values propagate through the network. Under these conditions, the encoding and decoding convolution filters \(\mathbf{K}\) and \(\tilde{\mathbf{K}}\) should comply with

\[\boxed{\mathbf{y}=\tilde{\mathbf{K}}^{\intercal}(\mathbf{K}\mathbf{y})_{+} \cdot c}\enspace. \tag{40}\]

It can be noticed that Eq. (40) is an extension of Eq. (20), which describes the reconstruction characteristics of tight framelets. From this point, we refer to convolutional kernels compliant with Eq (40) as _phase-complementary tight framelets_. As a final remark, it should be noted that a common practice in CNN designs is also to use ReLU non-linearities in the decoder, in such case the phase-complementary tight-framelet condition can still be met as long as the pixels \(y\in\mathbf{y}\) comply with \(y\geq 0\), which is equivalent to

\[\boxed{\mathbf{y}=(\mathbf{y})_{+}=\left(\tilde{\mathbf{K}}^{\intercal}( \mathbf{K}\mathbf{y})_{+}\cdot c\right)_{+}}\enspace. \tag{41}\]

It can be observed that the relevance of the properties defined in Eqs. (40) and (41) is that they ensure that a CNN can propagate any arbitrary signal, which is important to avoid any distortions (such as image blur) in the processed images.

An additional element of the theory of deep convolutional framelets regarding the reconstruction of the signal, is to show that conventional pooling layers (e.g. average pooling) discard high-frequency information of the signal, which effectively blurs the processed signals. Furthermore, Ye _et al._[12] have demonstrated that this can be fixed by replacing the conventional up/down-sampling layers by reversible operations, such

Fig. 7: Transfer characteristics of the semi-hard thresholds based on the difference of Gaussians and of the garrote shrink, as well as their clipping counterparts. Note that in contrast with the soft-shrinkage and clipping functions shown in Fig. 6, in the semi-hard thresholds tend to unity for large values, while the semi-hard clipping functions tend to zero for large signal intensities.
as the discrete wavelet transform. To exemplify this property, we refer to Fig. 4. If only an average pooling layer followed by an up-sampling stage would be applied, the treatment of the signal would be equivalent to the low-frequency branch of the DWT. Consequently, only the low-frequency spectrum of the signal would be recovered and images processed with that structure would become blurred. In contrast, if the full forward and inverse wavelet transform of Fig. 4 is used for up- and down-sampling, it is possible to reconstruct any signal, irrespective of its frequency content.

The ultimate key contribution of the theory of deep convolutional framelets is the explanation of the operation of ReLU-based noise-reduction CNNs. For the non-residual configuration, ReLU CNNs perform the following operations. (1) The convolution filters decompose the incoming signal into a sparse multi-channel representation. (2) The feature maps which are uncorrelated to the signal, contain mainly noise. In this case, the bias and the ReLU activation cancel the noisy feature maps in a process analogous to the MAP estimate shown in Section IV-C1. (3) The decoder reconstructs the
filtered image. Note that this process is analogous to the low-rank decomposition described in Section IV-B1. In the case of _residual_ networks, the CNN learns to estimate the noise, which means that in that configuration the ReLU non-linearities suppress the channels with high activation.

A visual example of the low-rank approximation in ReLU CNNs is shown in Fig. 8, which illustrates the operation of an idealized single-layer encoding-decoding ReLU CNN operating both, in residual and non-residual way. It can be noted The ReLU activation suppresses specific channels in the sparse decomposition provided by the encoder, thereby preserving the low-rank structures in the non-residual network. Alternatively, in the residual example, the ReLUs eliminate the feature maps with high activation, which results in a noise estimate that is subtracted from the input to estimate the noiseless signal.

### _Shrinkage and clipping-based CNNs_

Just as ReLU networks, the encoder of shrinkage networks [17, 18, 19] separates the input signal in a multi-channel representation. As a second processing stage, the shrinkage networks estimate the noiseless encoded signal by cancelling the low-amplitude pixels in the feature maps in a process akin to the MAP estimate of Section IV-C2. As final step, the encoder reconstructs the estimated noiseless image. Note that the use of shrinkage functions reduces the number of channels required by ReLU counterparts to achieve perfect signal reconstruction, because the shrinkage activation preserves positive and negative values, while ReLU only preserves the positive part of the signal.

Fig. 8: Operation of a simplified denoising (non-) residual ReLU CNN according to the theory of deep convolutional framelets (TDCF). In the figure, the noisy observation \(\mathbf{y}\) is composed by two vertical bars plus uncorrelated Gaussian noise. Furthermore, for this example, the encoding and decoding convolution filters (\(\mathbf{K}\) and \(\bar{\mathbf{K}}\), respectively) are the Haar basis of the 2D discrete wavelet transform and its phase-inverted counterparts. Given the content of the image, the image in the decomposed domain \(\mathbf{K}\mathbf{y}\) produces only a weak activation for the vertical and diagonal filters (\(\mathbf{w_{LH}}\) and \(\mathbf{w_{HH}}\), respectively) and those feature maps contain mainly noise. In the case of the non-residual network, the ReLUs and biases suppress the channels with low activation (see column \((\mathbf{K}\mathbf{y}+\underline{b})_{+}\)), which is akin to the low-rank approximation. In contrast, in the residual example, the channels with image content are suppressed, while preserving the uncorrelated noise. Finally, the decoding section reconstructs the noise-free estimate \(\bar{\mathbf{x}}\) for the non-residual network or the noise estimate \(\hat{\boldsymbol{\eta}}\) for the residual example, where it is subtracted from \(\mathbf{y}\) to compute the noiseless estimate \(\hat{\mathbf{x}}\).
As shown in Section III-A, in residual learning, a given encoding-decoding network estimates the noise signal \(\mathbf{\eta}\), so that it can be subtracted from the noisy observation \(\mathbf{y}\) to generate the noiseless estimate \(\hat{\mathbf{x}}\). As shown in Section IV-C3, in the framelet domain this is achieved by preserving the low-amplitude values of the feature maps by clipping the signal. Therefore in residual networks, the shrinkage functions can be explicitly replaced by clipping activations.

Visual examples of the operation of a single-layer shrinkage and clipping networks are presented in Fig. 8, where it can be noted that the operation of shrinkage and clipping networks is analogous to their ReLU counterparts, with the main difference that shrinkage and clipping networks do not require phase-complements in the encoding and decoding layers as ReLU-based CNNs do.

### _Shrinkage and clipping in ReLU networks_

As addressed in Section IV-C, the soft-threshold function is the superposition of two ReLU activations. As a consequence, it is feasible that in ReLU CNNs shrinkage behavior could arise in addition to the low-rankness enforcement mentioned in Section V-A. It should be noted that this only can happen if the number of channels of the encoder and decoder complies with the redundancy constraints of theory of deep convolutional framelets and if the decoder is linear. To prove this, Eq. (31) is reparametrized as

\[\hat{\mathbf{d}}=\tilde{\mathbf{K}}^{\intercal}(\mathbf{K}\mathbf{y}+\underline {b})_{+}, \tag{50}\]

where convolution filters \(\mathbf{K}\) and \(\tilde{\mathbf{K}}^{\intercal}\) are defined by \(\mathbf{K}=\begin{pmatrix}\left(\mathbf{I}&-\mathbf{I}\right)\end{pmatrix}^{\intercal}\) and \(\tilde{\mathbf{K}}^{\intercal}=\begin{pmatrix}\left(\mathbf{I}&-\mathbf{I} \right)\end{pmatrix}\), respectively, and \(\underline{b}=\begin{pmatrix}-t&-t\end{pmatrix}^{\intercal}\) represents the threshold value.

In addition to the soft-shrinkage, note that the clipping function described by Eq. (34) also can be expressed by Eq. (50) if \(\mathbf{K}=\begin{pmatrix}\left(\mathbf{I}&-\mathbf{I}&\mathbf{I}&-\mathbf{I} \right)\end{pmatrix}^{\intercal}\), \(\tilde{\mathbf{K}}^{\intercal}=\begin{pmatrix}\left(\mathbf{I}&-\mathbf{I}&- \mathbf{I}&\mathbf{I}\right)\end{pmatrix}\) and \(\underline{b}=\begin{pmatrix}0&0&-t\end{pmatrix}^{\intercal}\). It can be noted that representing the clipping function in convolutional form requires _four_ times more channels than the original input signal.

It should be noted that the ability of ReLUs to approximate other signals has also been noticed Daubechies _et al._[31], who have proven that deep ReLU CNNs are universal function approximators. In addition, Ye and Sung [13] have demonstrated that the ReLU function is the main source of the high-approximation power of CNNs.

Fig. 9: Operation of denoising in shrinkage and clipping networks. In the non-residual configuration, the noisy signal \(\mathbf{y}\) is decomposed by a set of convolution filters, which for this example are the 2D Haar basis functions of the discrete wavelet transform (\(\mathbf{K}\mathbf{y}\)). As a second step, the semi-hard shrinkage produces a MAP estimate of the noiseless detail bands/feature maps (\(\tau_{(\underline{0})}^{\text{Doc}}(\mathbf{K}\mathbf{y})\)). As third and final step, the decoder maps the estimated noiseless encoded signal to the original image domain. In the residual network, the behavior is similar, but the activation layer is a clipping function which performs a MAP estimate of the noise in the feature maps, which is reconstructed by the decoder to generate the noise estimate \(\hat{\mathbf{\eta}}\). After reconstruction, the noise estimate is subtracted from the noisy observation \(\mathbf{y}\) to generate the noise-free estimate \(\hat{\mathbf{x}}\).
### _Additional links between encoding-decoding CNNs and existing signal processing techniques_

Up to this moment, it has been assumed that the operation of the encoding and decoding convolution filters is limited to map the input image to a multi-channel representation and to reconstruct it (i.e. \(\mathbf{K}\) and \(\bar{\mathbf{K}}^{\intercal}\) comply with \(\bar{\mathbf{K}}^{\intercal}(\mathbf{K})_{+}=\mathbf{I}\cdot c\)). Still, it is possible that -besides decomposition and synthesis tasks- the encoding-decoding structure also filters/colors the signal in a way that improves the image estimates. It should be noted that this implies that the perfect reconstruction encoding-decoding structure is no longer preserved. For example, considering the following linear encoding-decoding structure

\[\hat{\mathbf{x}}=\bar{\mathbf{K}}^{\intercal}(\mathbf{K}\mathbf{y}), \tag{51}\]

which can be reduced to

\[\hat{\mathbf{x}}=\mathbf{k}*\mathbf{y}. \tag{52}\]

Here, \(\mathbf{k}=\bar{\mathbf{K}}^{\intercal}\mathbf{K}\), is optimized to reduce the distance between \(\mathbf{y}\) and the ground-truth \(\mathbf{x}\). Consequently, the equivalent filter \(\mathbf{k}\) can be considered a _Wiener filter_. It should be noted that this text is not the first in addressing the potential Wiener-like behavior of a CNN. For example, Mohan _et al._[14] suggested that by eliminating the bias of the convolution layers, the CNN could behave more akin to the Wiener filter and to be able to generalize better to unseen noise levels. It should be noted that by doing that the CNN can also behave akin to the switching behavior described by the theory of deep convolutional framelets, which can be described by the
equation

\[(z)_{+}=\begin{cases}z,&\text{if }z\geq 0,\\ 0&\text{if }z<0,\end{cases}, \tag{53}\]

where \(z\) is a pixel which belongs to the signal \(\mathbf{z}=\mathbf{k}*\mathbf{x}\). It can be observed that in contrast with the low-rank behavior described in Section V-A, in this case the switching behavior is only dependent on the correlation between signal \(\mathbf{x}\) and the filter \(\mathbf{k}\). Consequently, if the value of \(z\) is positive, its value is preserved. On the contrary, if the correlation between \(\mathbf{x}\) and \(\mathbf{k}\) is negative, then the value of \(z\) is cancelled. Consequently, the noise reduction becomes independent/invariant of the noise level. It can be observed that this effect can be considered a non-linear extension of the so-called signal annihilation filters [32].

It should be noticed that besides the low-rank approximation interpretation of ReLU-based CNNs, additional links to other techniques can be derived. For example, the decomposition and synthesis provided by the encoding-decoding structure is also akin to the non-negative matrix factorization (NMF) [33], in which a signal is factorized as a weighted sum of positive bases. In this conception, the feature maps are the bases, which are constrained to be positive by the ReLU function. Furthermore, an additional interpretation of encoding-decoding CNNs can be obtained by analyzing them from a low-dimensional manifold representation perspective [8]. Here, the convolution layers of CNNs are interpreted as two operations. On one hand, they can provide a Hankel representation. On the other hand, they provide a bottleneck which reduces the dimensionality of the manifold of image patches. It should be noticed that the Hankel-like structure that is attributed to the convolution layers of CNNs, has also been noticed by the theory of the deep convolutional framelets [12]. Two final connection with signal-processing and CNNs is the variational formulation combined with kernel-based methods [15] and the convolutional sparse coding interpretation of CNNs by Papyan _et al._[16].

## VI **Mathematical analysis of relevant designs**

In order to demonstrate the application of the principles summarized in Sections IV and V, this section analyzes relevant designs of ReLU and shrinkage CNNs. The analyses focus on three main aspects, which are: (1) the overall descriptions of the network architecture, (2) the signal reconstruction characteristics provided by the convolutional layers of the encoder and decoder sub-networks, and (3) the number operations \(\mathcal{O}(\cdot)\) executed by the trainable parts of the network, since this will give insight on the computational requirements to execute each network and its overall complexity.

The signal reconstruction analysis provides a theoretical indication that a given CNN design can propagate _any_ arbitrary signal when considering the use of ideal filters (i.e. they provide perfect reconstruction and are maximally sparse). In other words, for a fixed network architecture, there exists a selection of parameters (weights and biases) that make the neural network equal to the identity function. This result is important, because a design that cannot propagate arbitrary signals under ideal conditions, will potentially distort the signals that propagate through it _by design_. Consequently, this cannot be fixed by training with large datasets and/or with the application of any special loss term. In order to understand better the signal reconstruction analysis, we provide a brief example, where it is a non-residual CNN \(G(\cdot)\), where we propagate a noiseless signal \(\mathbf{x}\) contaminated with noise \(\boldsymbol{\eta}\), so that

\[\mathbf{x}\approx G(\mathbf{x}+\boldsymbol{\eta}). \tag{54}\]

Here, an ideal CNN allows to propagate _any_\(\mathbf{x}\), while cancelling the noise component \(\boldsymbol{\eta}\), irrespective of the content of \(\mathbf{x}\). If we switch our focus to an ideal _residual CNN_\(R(\cdot)\), it is possible to observe that

\[\hat{\mathbf{x}}\approx R(\mathbf{y})=\mathbf{y}-G(\mathbf{y}). \tag{55}\]

Here, \(G(\cdot)\) is the encoding-decoding section of the residual network \(R(\cdot)\). Consequently, it is desirable that the network \(G(\cdot)\) is able to propagate the noise \(\boldsymbol{\eta}\), while suppressing the noiseless signal \(\mathbf{x}\), which is equivalent to

\[\boldsymbol{\eta}\approx G(\mathbf{x}+\boldsymbol{\eta}). \tag{56}\]

It should be noted that in both residual and non-residual cases, there are two behaviors. On one hand, there is a signal which the network decomposes and reconstructs (almost) perfectly, and on the other hand a signal is suppressed. The signal reconstruction analysis focuses on the signals that the network can propagate or reconstruct, rather than the signal cancellation behavior. In consequence, we focus on the linear part of \(G(\cdot)\) (i.e. its convolution structure), of which, according to

Fig. 10: Simplified structure of encoding-decoding ReLU CNNs. The displayed networks are the U-Net/filtered back-projection network the encoder-decoder residual CNN (RED) and finally, the learned wavelet-frame shrinkage network. Note that for all the designs, the encoding-decoding structures are indicated by dashed blocks. It should be heard in mind that the drawings are simplified, they do not contain normalization layers, are shallow, commonly appearing dual convolutions are drawn as one layer.
Section V-A, we assume that it handles the decomposition and reconstruction of the signal within the CNN. It should be noted that the idealized model assumed here, is only considered for analysis purposes, since practical implementations do no guarantee that this exact behavior is factually obtained. The reader is referred to Section V-D and the text boxes _Fitting low-rank approximation in ReLU CNNs_ and _Network depth_.

In order to test the perfect reconstruction in non-residual CNNs, we propose the following procedure. (1) We assume an _idealized_ model \(G(\cdot)\), where its convolution filters \(\mathbf{K}_{n}\) and \(\tilde{\mathbf{K}}_{n}\) comply with the phase-complementary tight framelet condition and where the biases and non-linearities suppress low-amplitude (and negative for ReLU activations) samples from the feature maps. (2) The biases/thresholds of ReLU/shrinkage CNNs are set to zero (or to infinity for clipping activations). It can be observed that this condition prevents the low-rank (or high-rank for residual models) approximation behavior of the idealized CNN. Under this circumstance, it should be possible to prove that the analyzed CNN can perfectly reconstruct _any_ signal. (3) The last step involves simplifying the mathematical description of the resulting model of the previous point. The mathematical simplification of the model should lead to the identity function if the model complies with the perfect reconstruction property.

To conclude the explanation on the perfect reconstruction analysis, we provide two relevant considerations. First, it can be claimed that residual networks, such as the model \(R(\mathbf{y})=\mathbf{y}-G(\mathbf{y})\) discussed in Eq. (55), is able to reconstruct any signal when \(G(\mathbf{y})=0\) for any \(\mathbf{y}=\mathbf{x}+\mathbf{\eta}\). Still, this does not convey information about the behavior of the encoding-decoding network \(G(\cdot)\), which should be able to perform perfect decomposition and reconstruction of the noise signal \(\mathbf{\eta}\), as discussed in Eq. (55). To avoid this trivial solution, instead of analyzing the network \(R(\cdot)\), the analysis described for non-residual models is applied to the encoding-decoding structure \(G(\cdot)\), which means that the residual connection is excluded from the analysis.

The second concluding remark is that in order to distinguish the equations of the perfect signal reconstruction analysis from other models, we specify the analyzed designs of the perfect reconstruction models, in which the low-rank approximation behavior is avoided by setting the bias values to zero, with a special operator \(\mathcal{P}\{\cdot\}\).

For the analyses regarding the total number of operations of the trainable parameters, it is assumed that the tensors \(\mathbf{K}_{0}\), \(\tilde{\mathbf{K}}_{0}^{\mathsf{T}}\), \((\tilde{\mathbf{K}}_{0}^{\mathsf{u}})^{\mathsf{T}}\), \((\tilde{\mathbf{K}}_{0}^{\mathsf{d}})^{\mathsf{T}}\), \(\mathbf{K}_{1}\) and \(\tilde{\mathbf{K}}_{1}^{\mathsf{T}}\) shown in Fig. 10 have dimensions \((C_{0}\times 1\times N_{\mathsf{f}}\times N_{\mathsf{f}})\), \((1\times C_{0}\times N_{\mathsf{f}}\times N_{\mathsf{f}})\), \((1\times C_{0}\times N_{\mathsf{f}}\times N_{\mathsf{f}})\), \((1\times C_{0}\times N_{\mathsf{f}}\times N_{\mathsf{f}})\), \((C_{1}\times C_{0}\times N_{\mathsf{f}}\times N_{\mathsf{f}})\), \((C_{0}\times C_{1}\times N_{\mathsf{f}}\times N_{\mathsf{f}})\), \(N_{\mathsf{f}}\)), respectively. Here, \(C_{0}\) and \(C_{1}\) represent the number of channels after the first and second convolution layer, all the convolution filters are assumed to be square with \(N_{\mathsf{f}}\times N_{\mathsf{f}}\) pixels. Furthermore, the input signal \(\mathbf{x}\) has dimensions \((1\times 1\times N_{\mathsf{f}}\times N_{\mathsf{c}})\), where \(N_{\mathsf{r}}\) and \(N_{\mathsf{c}}\) denote the numbers of rows and columns, respectively.

The analyses shown for the different networks in this article have the following limitations. (1) The analyzed networks have only enough decomposition levels and convolution layers to understand their basic operation. The motivation for this simplification is to keep the analyses short and clear. Moreover, the same principles can be extended to deeper networks. Since the same single-decomposition CNNs would be recursively embedded within the given architectures. (2) The normalization layers are not considered, because they are linear operators which provide mean shifts and amplitude scaling. Consequently, for analysis purposes it can be assumed that they are embedded in the convolution weights. (3) For every encoder convolution kernel it is assumed that there is at least one decoder filter. (4) No co-adaptation between the filters of the encoder and decoder layers is considered.

The remainder of this section shows analyses of a selection of a few representative designs. Specifically, the chosen designs are the _U-Net_[34] and its residual counterpart, the _filtered back-projection network_[21]1. Additional designs analyzed here are the _the residual encoder-decoder CNN_[5]2, as well as the _learned wavelet-shrinkage network_ (LWFSN) 3. For reference, all the designs are portrayed in Fig. 10.

Footnote 1: Matlab implementation by their authors available at [https://github.com/panakinfo/FBPConvNet](https://github.com/panakinfo/FBPConvNet)

Footnote 2: Pytorch implementation by their authors available at [https://github.com/Sinyu/REDC-CNN](https://github.com/Sinyu/REDC-CNN)

Footnote 3: Pytorch implementation available at [https://github.com/LuisAlbertZM/demo_LWFSN_TMI](https://github.com/LuisAlbertZM/demo_LWFSN_TMI) and interactive demo available at IEEE’s code ocean [https://codeocean.com/capsule/9027829/tree/v1](https://codeocean.com/capsule/9027829/tree/v1). The demo also includes as reference pytorch implementations of FBPComNet and the tight-frame U-Net.

### **U-Net/filtered back-projection network**

#### V-A1 **U-Net - overview of the design**

The first networks analyzed are the U-Net and filtered back-projection networks both of which share the encoding-decoding structure \(U(\cdot)\). However, they differ in the fact that the U-Net is non-residual, while the filtered back-projection network operates in residual configuration. Therefore, the estimate of the noiseless signal \(\hat{\mathbf{x}}\) from the noisy observation \(\mathbf{y}\) in the conventional U-Net is achieved by

\[\hat{\mathbf{x}}=U(\mathbf{y}), \tag{57}\]

whereas in the filtered back-projection network, \(U(\cdot)\) is used in residual configuration, which is equivalent to

\[\hat{\mathbf{x}}=\mathbf{y}-U(\mathbf{y}). \tag{58}\]

If we now switch focus to the encoding-decoding structure of the U-Net \(U(\mathbf{y})\), it can be shown that it is described by

\[U(\mathbf{y})=U^{\mathsf{u}}(\mathbf{y})+U^{\mathsf{d}}(\mathbf{y}), \tag{59}\]

where \(U^{\mathsf{u}}(\mathbf{y})\) corresponds to the _undeimated path_ and is defined by

\[U^{\mathsf{u}}(\mathbf{y})=\big{(}\tilde{\mathbf{K}}_{0}^{\mathsf{u}}\big{)}^{ \mathsf{T}}\big{(}\mathbf{K}_{0}\mathbf{y}+\underline{b}_{0}\big{)}_{+}, \tag{60}\]

while the decimated path is

\[U^{\mathsf{d}}(\mathbf{y})=\big{(}\tilde{\mathbf{K}}_{0}^{\mathsf{d}}\big{)}^{ \mathsf{T}}\tilde{\mathbf{W}}_{\mathsf{L}}^{\mathsf{T}}f_{(2\uparrow)}\Big{(} \big{(}\tilde{\mathbf{K}}_{1}^{\mathsf{T}}(\mathbf{K}_{1}\mathbf{Z}+ \underline{b}_{1})_{+}+\tilde{b}_{1}\big{)}_{+}\Big{)}. \tag{61}\]

Here, signal \(\mathbf{Z}\) is defined by

\[\mathbf{Z}=f_{(2\downarrow)}\big{(}\mathbf{W}_{\mathsf{L}}(\mathbf{K}_{0} \mathbf{y}+\underline{b}_{0})_{+}\big{)}. \tag{62}\]

Note the decimated path contains two nested encoding-decoding architectures, as observed by Jin _et al._[21], who
has acknowledged that the nested filtering structure is akin to the (learned) iterative soft thresholding algorithm [27, 28].

#### V-A2 **U-Net - signal reconstruction analysis**

To prove if the U-Net can perfectly reconstruct _any_ signal, we assume that the biases are equal to zero, on this condition the network \(\mathcal{P}\{U\}(\mathbf{y})\) is defined by

\[\mathcal{P}\{U\}(\mathbf{y})=\mathcal{P}\{U^{\mathrm{u}}\}(\mathbf{y})+ \mathcal{P}\{U^{\mathrm{d}}\}(\mathbf{y}), \tag{63}\]

where sub-network \(\mathcal{P}\{U^{\mathrm{u}}\}(\cdot)\) is defined by

\[\mathcal{P}\{U^{\mathrm{u}}\}(\mathbf{y})=\big{(}\tilde{\mathbf{K}}_{0}^{ \mathrm{u}}\big{)}^{\mathsf{T}}(\mathbf{K}_{0}\mathbf{y})_{+}. \tag{64}\]

Assuming that \(\big{(}\mathbf{K}_{0},\tilde{\mathbf{K}}_{0}^{\mathrm{u}}\big{)}\) is a complementary-phase tight-framelet pair then \(\mathcal{P}\{U^{\mathrm{u}}\}(\mathbf{y})\) is simplified to

\[\mathcal{P}\{U^{\mathrm{u}}\}(\mathbf{y})=\mathbf{y}\cdot c_{0}. \tag{65}\]

Furthermore, the low-frequency path is

\[U^{\mathrm{d}}(\mathbf{y})=\big{(}\tilde{\mathbf{K}}_{0}^{\mathrm{d}}\big{)}^ {\mathsf{T}}\tilde{\mathbf{W}}_{\mathrm{L}}^{\mathsf{T}}f_{(2\uparrow)}\Big{(} \big{(}\tilde{\mathbf{K}}_{1}^{\mathsf{T}}(\mathbf{K}_{1}\mathcal{P}\{ \mathbf{Z}\})_{+}\big{)}_{+}\Big{)}, \tag{66}\]

Where \(\mathcal{P}\{\mathbf{Z}\}\) is defined by

\[\mathcal{P}\{\mathbf{Z}\}=f_{(2\downarrow)}\big{(}\mathbf{W}_{\mathrm{L}}( \mathbf{K}_{0}\mathbf{y})_{+}\big{)}. \tag{67}\]

If \(\mathbf{K}_{1}\) is a phase-complementary tight frame, we know that \(\tilde{\mathbf{K}}_{1}^{\mathsf{T}}(\mathbf{K}_{1}\mathbf{Z})_{+}=\mathbf{Z} \cdot c_{1}\). Consequently, Eq. (66) becomes

\[\mathcal{P}\{U^{\mathrm{d}}\}(\mathbf{y})=\big{(}\tilde{\mathbf{K}}_{0}^{ \mathrm{d}}\big{)}^{\mathsf{T}}\tilde{\mathbf{W}}_{\mathrm{L}}^{\mathsf{T}}f_{ (2\uparrow)}\big{(}(\mathbf{W}_{\mathrm{L}}\mathbf{K}_{0}\mathbf{y})_{+}\big{)} \big{)}\cdot c_{1}. \tag{68}\]

Here, it can be noticed that if \(\mathbf{K}_{0}\) is a phase-complementary tight framelet, then \(\mathcal{P}\{U^{\mathrm{d}}\}(\mathbf{y})\) approximates a low-pass version of \(\mathbf{y}\), or equivalently

\[\mathcal{P}\{U^{\mathrm{d}}\}(\mathbf{y})\approx\boldsymbol{\mathcal{W}}_{ \mathrm{L}}\mathbf{y}\cdot c_{1}, \tag{69}\]

where \(\boldsymbol{\mathcal{W}}_{\mathrm{L}}\) is a low-pass filter. Finally, substituting Eq. (65) and Eq. (69) in Eq. (63) results in

\[\mathcal{P}\{U\}(\mathbf{y})\approx(\mathbf{I}\cdot c_{0}+\boldsymbol{ \mathcal{W}}_{\mathrm{L}}\cdot c_{1})\mathbf{y}. \tag{70}\]

This result proves that the design of the U-Net cannot evenly reconstruct all the frequency of \(\mathbf{y}\) unless \(c_{1}=0\), in which case, the whole low-frequency branch of the network is ignored. Note that this limitation is inherent to its design and cannot be circumvented by training with large datasets and/or with any loss function.

#### V-A3 **U-Net - number of operations**

It can be noted that encoding filter \(\mathbf{K}_{0}\), convolves \(\mathbf{x}\) at its original resolution and maps it to a tensor with \(C_{0}\) channels. Therefore, the number of operations \(\mathcal{O}(\cdot)\) for kernel \(\mathbf{K}_{0}\) is \(\mathcal{O}(\mathbf{K}_{0})=C_{0}\cdot N_{t}\cdot N_{\mathrm{c}}\cdot N_{f}^{2}\) [FLOPS] (floating-point operations). Conversely, due to the symmetry between encoder and decoder filters, \(\mathcal{O}(\tilde{\mathbf{K}}_{0}^{\mathrm{u}})=\mathcal{O}(\tilde{\mathbf{K }}_{0}^{\mathrm{d}})=\mathcal{O}(\mathbf{K}_{0})\). Furthermore, for this design, filter \(\mathbf{K}_{1}\) processes the signal encoded by \(\mathbf{K}_{0}\), which is down-sampled by a factor of one half and maps it from \(C_{0}\) to \(C_{1}\) channels, this results in the estimated operation cost \(\mathcal{O}(\mathbf{K}_{1})=\mathcal{O}(\tilde{\mathbf{K}}_{1})=C_{0}\cdot C _{1}\cdot N_{\mathrm{r}}\cdot N_{\mathrm{c}}\cdot N_{\mathrm{f}}^{2}\cdot(2)^ {-2}\) [FLOPS]. Finally, adding the contributions of filters \(\mathbf{K}_{0}\), \(\tilde{\mathbf{K}}_{0}^{\mathrm{u}}\), \(\tilde{\mathbf{K}}_{0}^{\mathrm{d}}\), \(\mathbf{K}_{1}\) and \(\tilde{\mathbf{K}}_{1}\) results in

\[\mathcal{O}(U)=(3+2^{-1}\cdot C_{1})\cdot C_{0}\cdot N_{\mathrm{r}}\cdot N_{ \mathrm{c}}\cdot N_{\mathrm{f}}^{2}\]

[FLOPS] (71)

#### V-A4 **U-Net - Concluding remarks**

The U-Net/FBPConvNet is a flexible multi-resolution architecture. Still, as it has been shown, the pooling structure of this CNN may be sub-optimal for noise reduction applications because this configuration does not allow to recover the frequency information of the signal evenly. This has been noted and fixed by Han and Ye [6], who introduced the so-called tight-frame U-Net. In which the down/up-sampling structure is replaced by the discrete wavelet transform and its inverse. This simple modification overcomes the limitations of the U-Net and improved its performance for artifact removal in compressed sensing imaging.

### **Residual encoder-decoder CNN**

#### V-B1 **Residual encoder-decoder CNN - overview of the design**

The residual encoder-decoder CNN shown in Fig. 10 consists of nested single-layer residual encoding-decoding networks. For example, in the network showcased in Fig. 10 it is visible that network \(R_{1}(\cdot)\) is nested into \(R_{0}(\cdot)\). Furthermore, for this case the image estimate is given by

\[\hat{\mathbf{x}}=(\mathbf{y}+R_{0}(\mathbf{y})+\tilde{\underline{b}}_{0})_{+}, \tag{72}\]

in which \(R_{0}(\cdot)\) is the outer residual network and \(\tilde{\underline{b}}_{0}\) is the bias for the output layer. Note that the ReLU placed at the output layer intrinsically assumes that the estimated signal \(\hat{\mathbf{x}}\) is positive.

From Eq. (72), the output of the sub-network \(R_{0}(\cdot)\) is defined by

\[\mathbf{Z}=R_{0}^{\text{dec}}(\hat{\mathbf{Q}}). \tag{73}\]

Here, the decoder \(R_{0}^{\text{dec}}(\cdot)\) is defined by

\[R_{0}^{\text{dec}}(\hat{\mathbf{Q}})=\tilde{\mathbf{K}}_{0}^{\mathsf{T}}\hat{ \mathbf{Q}}. \tag{74}\]

In the above, \(\hat{\mathbf{Q}}\) is the noiseless estimate of the intermediate signal \(\mathbf{Q}\) and it is defined by

\[\hat{\mathbf{Q}}=(\mathbf{Q}+R_{1}(\mathbf{Q})+\tilde{\underline{b}}_{1})_{+}, \tag{75}\]

where the network \(R_{1}(\cdot)\) is

\[R_{1}(\mathbf{Q})=\tilde{\mathbf{K}}_{1}^{\mathsf{T}}\big{(}\mathbf{K}_{1} \mathbf{Q}+\underline{b}_{1})_{+}. \tag{76}\]

Furthermore, \(\mathbf{Q}\) represents the signal encoded by \(R_{0}(\cdot)\), or equivalently

\[\mathbf{Q}=R_{0}^{\text{enc}}(\mathbf{y}), \tag{77}\]

where \(R_{0}^{\text{enc}}(\cdot)\) is defined by

\[R_{0}^{\text{enc}}(\mathbf{y})=\mathbf{K}_{0}\mathbf{y}. \tag{78}\]

#### V-B2 **Residual encoder-decoder CNN - signal reconstruction analysis**

As mentioned earlier, the residual encoder-decoder CNN is composed by nested residual blocks, which are independently analyzed to study the reconstruction characteristics of this network. First, block \(R_{1}(\cdot)\), is given by

\[\mathcal{P}\{R_{1}\}(\mathbf{Q})=\tilde{\mathbf{K}}_{1}^{\mathsf{T}}(\mathbf{K }_{1}\mathbf{Q})_{+}. \tag{79}\]

Under complementary-phase tight-frame assumptions for the pair \((\mathbf{K}_{1},\tilde{\mathbf{K}}_{1})\), Eq. (79) reduces to

\[\mathcal{P}\{R_{1}\}(\mathbf{Q})=\mathbf{Q}, \tag{80}\]
which shows that the encoder and decoder \(R_{1}(\cdot)\) can approximately reconstruct any signal. Now switching to \(R_{0}\) it can be observed that the linear part is

\[\mathcal{P}\{R_{0}\}(\mathbf{y})=\tilde{\mathbf{K}}_{0}^{\intercal}(\mathbf{K}_{ 0}\mathbf{y})_{+}. \tag{81}\]

Just as with \(R_{1}(\cdot)\), it is assumed that the convolution kernels are tight-framelets. Therefore Eq. (81) becomes

\[\mathcal{P}\{R_{0}\}(\mathbf{y})=\mathbf{y}. \tag{82}\]

Consequently, \(R_{0}(\cdot)\) and \(R_{1}(\cdot)\) can reconstruct any arbitrary signal under complementary-phase tight-frame assumptions.

#### V-B3 **Residual encoder-decoder CNN - number of operations**

In this case, all the convolution layers operate at the original resolution of image \(\mathbf{x}\). Therefore, the number of operations \(O(\cdot)\) for kernel \(\mathbf{K}_{0}\) and \(\tilde{\mathbf{K}}_{0}\) is \(O(\mathbf{K}_{0})=O(\tilde{\mathbf{K}}_{0})=C_{0}\cdot N_{\text{r}}\cdot N_{ \text{c}}\cdot N_{\text{f}}^{2}\) [FLOPs], while \(\mathbf{K}_{1}\) and \(\tilde{\mathbf{K}}_{1}\) are requiring \(O(\mathbf{K}_{1})=O(\tilde{\mathbf{K}}_{1})=C_{0}\cdot C_{1}\cdot N_{\text{r} }\cdot N_{\text{c}}\cdot N_{\text{f}}^{2}\) [FLOPs]. By adding the contributions of both encoding-decoding pairs, the total operations for the residual encoder-decoder become

\[\mathcal{O}(R)=2\cdot(1+C_{1})\cdot C_{0}\cdot N_{\text{r}}\cdot N_{\text{c}} \cdot N_{\text{f}}^{2}\ [\text{FLOPS}]. \tag{83}\]

#### V-B4 **Residual encoder-decoder CNN - concluding remarks**

The residual encoder-decoder network consists on a set of nested single-resolution residual encoding-decoding CNNs. The single-resolution design increases its computation cost with respect to multi-resolution designs such as the U-Net. In addition, it should be noted that the use of a ReLU as output layer of the encoder-decoder residual network forces the signal estimates to be positive, but this is not always convenient. For example in computed tomography imaging, it is common that images contain positive and negative values.

### **Learned wavelet-frame shrinkage network - description of the architecture**

The learned wavelet-frame shrinkage network is a multi-resolution architecture in which the discrete wavelet transform is used for down/up-sampling and also as part of the decomposition where shrinkage is applied. In this CNN, the noiseless estimates are produced by

\[\hat{\mathbf{x}}=L(\mathbf{y}), \tag{84}\]

where \(L(\cdot)\) represents the encoding-decoding structure of the learned wavelet-frame shrinkage network and the encoding-decoding network \(L(\cdot)\) is

\[L(\mathbf{y})=L_{\text{L}}(\mathbf{y})+L_{\text{H}}(\mathbf{y}). \tag{85}\]

Here, the high-frequency path is given by

\[L_{\text{H}}(\mathbf{y})=\tilde{\mathbf{K}}_{0}^{\intercal}\tilde{\mathbf{W}} _{\text{H}}^{\intercal}f_{(2\uparrow)}\big{(}\tau_{(\downarrow)}^{\text{LET} }\big{)}\big{(}f_{(2\downarrow)}\big{(}\mathbf{W}_{\text{H}}\mathbf{K}_{0} \mathbf{y}\big{)}\big{)}\big{)}. \tag{86}\]

Note that in this design the encoder leverages the filter \(\mathbf{W}_{\text{H}}\) for generating a sparse signal prior to the shrinkage stage, i.e. \(\tau_{(\downarrow)}^{\text{LET}}\big{(}f_{(2\downarrow)}\big{(}\mathbf{W}_{ \text{H}}\mathbf{K}_{0}\mathbf{y}\big{)}\big{)}\). Meanwhile, the low-frequency path \(L_{\text{L}}(\cdot)\) is

\[L_{\text{L}}(\mathbf{y})=\tilde{\mathbf{K}}_{0}^{\intercal}\tilde{\mathbf{W}} _{\text{L}}^{\intercal}f_{(2\uparrow)}\big{(}f_{(2\downarrow)}\big{(}\mathbf{ W}_{\text{L}}\mathbf{K}_{0}\mathbf{y}\big{)}\big{)}. \tag{87}\]

#### V-B2 **Learned wavelet-frame shrinkage network - signal reconstruction analysis**

When analyzing the signal propagation of the learned wavelet-frame shrinkage network, we set the threshold level \(\underline{t}_{0}=0\). This turns Eq. (85) turns into

\[\mathcal{P}\{L\}(\mathbf{y})=\mathcal{P}\{L_{\text{L}}\}(\mathbf{y})+\mathcal{ P}\{L_{\text{H}}\}(\mathbf{y}). \tag{88}\]

Here, \(\mathcal{P}\{L_{\text{H}}\}(\cdot)\) is defined by

\[\mathcal{P}\{L_{\text{H}}\}(\mathbf{y})=\tilde{\mathbf{K}}_{0}^{\intercal} \tilde{\mathbf{W}}_{\text{H}}^{\intercal}f_{(2\uparrow)}\big{(}f_{(2 \downarrow)}\big{(}\mathbf{W}_{\text{H}}\mathbf{K}_{0}\mathbf{y}\big{)}\big{)}, \tag{89}\]

while the low-frequency path \(\mathcal{P}\{L_{\text{L}}\}(\cdot)\) is mathematically described by

\[\mathcal{P}\{L_{\text{L}}\}(\mathbf{y})=\tilde{\mathbf{K}}_{0}^{\intercal} \tilde{\mathbf{W}}_{\text{L}}^{\intercal}f_{(2\uparrow)}\big{(}f_{(2 \downarrow)}\big{(}\mathbf{W}_{\text{L}}\mathbf{K}_{0}\mathbf{y}\big{)}\big{)}. \tag{90}\]

Substituting Eq. (89) and Eq. (90) in Eq. (88) results in the equation

\[\mathcal{P}\{L\}(\mathbf{y})=\tilde{\mathbf{K}}_{0}^{\intercal}\tilde{\mathbf{ W}}^{\intercal}f_{(2\uparrow)}\big{(}f_{(2\downarrow)}\big{(}\mathbf{W} \mathbf{K}_{0}\mathbf{y}\big{)}\big{)}. \tag{91}\]

For the discrete wavelet transform, it holds that \(\mathbf{Q}=\tilde{\mathbf{W}}^{\intercal}f_{(2\uparrow)}\big{(}f_{(2\downarrow )}\big{(}\mathbf{W}\mathbf{Q}\big{)}\big{)}\). Consequently, Eq. (91) is simplified to

\[\mathcal{P}\{L\}(\mathbf{y})=\tilde{\mathbf{K}}_{0}^{\intercal}\mathbf{K}_{0} \mathbf{y}. \tag{92}\]

Assuming that \(\mathbf{K}_{0}\) is a tight framelet i.e. \(\tilde{\mathbf{K}}_{0}^{\intercal}\mathbf{K}_{0}=\mathbf{I}\cdot c\), with \(c=1\), then

\[\mathcal{P}\{L\}(\mathbf{y})=\mathbf{y}. \tag{93}\]

This proves that the encoding-decoding section of the learned wavelet-frame shrinkage network allows for perfect signal reconstruction.

#### V-B3 **Learned wavelet-frame shrinkage network - number of operations**

The learned wavelet-frame shrinkage network contains a simpler convolution structure than the networks reviewed up to this moment. Therefore, for a single-level decomposition architecture, the total number of operations is

\[\mathcal{O}(L)=2\cdot C_{0}\cdot N_{\text{r}}\cdot N_{\text{c}}\cdot N_{\text{ f}}^{2}\ [\text{FLOPS}]. \tag{94}\]

#### V-B4 **Learned wavelet-frame shrinkage network - residual variant**

To illustrate the use of clipping activations in residual noise reduction, the residual version of the learned wavelet-frame shrinkage network is included. Note that there are two main differences with the conventional learned wavelet-frame shrinkage network. First, the shrinkage functions are replaced by clipping activations. Second, the low-frequency signal is suppressed. This is performed because the original design of the learned wavelet-frame shrinkage network does

Fig. 11: Residual version of the learned wavelet-frame shrinkage network. It can be noticed that the low-frequency branch of the network is nulled. In deeper networks it would further decomposed and the nulling would be activated at the deepest level (lowest resolution).
not have any non-linearities in that section. This is akin to the low-frequency nulling proposed by Kwon and Ye [35]. The modified learned wavelet-frame shrinkage network is shown in Fig. 11. It can be observed that by setting to zero the low-frequency branch of the design the model is inherently assuming that the noise is high-pass.

#### Vi-B5 (Residual) Learned wavelet-frame shrinkage network - concluding remarks

The (residual) learned wavelet-frame shrinkage network is a design which explicitly mimics wavelet-shrinkage algorithms. It can be observed that the (r)LWFSN inherently assume that noise is high-frequency and explicitly avoid non-linear processing in the low-frequency band. Follow-up experiments also included non-linearities in the low-frequency band of the learned wavelet-frame shrinkage network [36] and obtained similar results to the original design.

## VII **What happens in trained models?**

### **Properties of convolution kernels and low-rank approximation**

The assumption that the convolution filters of a CNN behave as (complementary-phase) tight framelets is useful for analyzing the theoretical ability of a CNN to propagate signals. Albeit, it is difficult to prove that trained models comply with this assumption, because there are diverse elements affecting the optimization of the model, e.g. the initialization of the network, the data presented to the model, the optimization algorithm as well as its parameters. In addition, in real CNNs, there may be co-adaptation between the diverse CNN layers, which may prevent that the individual filters of the CNN behave as tight framelets, since the decomposition and filtering performed by one layer is not independent from the rest [37].

To test the behavior if the filters of trained CNN can converge to complementary-phase tight framelets, at least, on a simplified environment, we propose to train a toy model, as displayed in Fig. 13. If the trained filters of an encoder-decoder pair of the toy model (\(\mathbf{K}_{l}\), \(\mathbf{\tilde{K}}_{l}\)), (where \(l\) denotes one of the decomposition levels) behave as a complementary-phase tight framelet, then the pair (\(\mathbf{K}_{l}\), \(\mathbf{\tilde{K}}_{l}\)) approximately complies with the condition presented in Eq. (40), which for identity input \(\mathbf{I}\) simplifies to

\[\mathbf{\tilde{K}}_{n}^{\intercal}(\mathbf{K}_{n})_{+}=\mathbf{I}\cdot c_{n}\,, \tag{95}\]

in which \(c_{n}\) is an arbitrary constant.

The toy model is trained on images that contain multiple randomly-generated overlapping triangles. All the images were scaled to the range [0,1]. For this experiment, the input to the images is the noise-contaminated image and the objective/desired output is the noiseless image. For training the CNNs, normally-distributed noise with standard deviation of 0.1 was added to the ground truth. For every epoch, a batch of 192 training images is generated. As validation and test images we use the "astronaut" and "cameraman" images included in the software package _Scipy_. The model is optimized with Adam for 25 epochs with a linearly decreasing learning rate. The initial learning rate for the optimizer is set to 10\({}^{-3}\) and the batch size is set to 1 sample. The convolution kernels were initialized with Xavier initialization using a uniform distribution (see Glorot and Bengio [38]). The code is available at IEEE's code ocean [https://codeocean.com/capsule/7845737/tree](https://codeocean.com/capsule/7845737/tree).

With the described settings, we have trained the toy model, and have tested if the phase-complementary tight-framelet property holds for the filters of the deepest level \(l\)=2. The results for the operation \(\mathbf{\tilde{K}}_{2}^{\intercal}(\mathbf{K}_{2})_{+}\) are displayed in Fig. 12 (left), which shows that when the weights of the encoder and decoder have different initial values, the kernel pair (\(\mathbf{K}_{2}\), \(\mathbf{\tilde{K}}_{2}\)) are not complementary-phase tight framelets. We have observed that the forward and inverse filters of wavelets/framelets are often the same or at least very similar. Based on this reasoning, we have initialized the toy model with the same initial values of the kernel pairs (\(\mathbf{K}_{n}\), \(\mathbf{\tilde{K}}_{n}\)). As shown by Fig. 12 (right), with the proposed initialization, the filters of the CNN converge to tensors with properties reminiscent of complementary-phase tight-framelets. This suggests that the initialization of the CNN has an important influence on the convergence of the model to a specific solution.

Fig. 14 displays a test image processed with two toy models, one trained trained with different and one trained with the same initial values for the encoding-decoding pairs. It can be observed that there are no significant differences between the images produced by both models. In the same figure (lower row), we have set the bias of both networks to zero. In this case, it is expected that the networks reconstruct the noisy input, as confirmed by the figure, where both CNNs partly reconstruct the original noisy signal. This result suggests that the ReLU plus bias pairs operate akin to the low-rank approx

Fig. 12: Phase-complementary tight-framelet test for the trained toy network, initialized with random weights. The left figure shows the product \(\mathbf{\tilde{K}}_{l}^{\intercal}(\mathbf{K}_{2})_{+}\), where the initialization of \(\mathbf{K}_{2}\) and \(\mathbf{\tilde{K}}_{2}\) is different. It can be seen that the pair \((\mathbf{K}_{2},\mathbf{\tilde{K}}_{2})\) does not comply with the complementary-phase framelel criterion of Eq. (95). This contrasts with the right result, which displays the result of the product \(\mathbf{\tilde{K}}_{l}^{\intercal}(\mathbf{K}_{2})_{+}\), for the same CNN, but where the initial values of \(\mathbf{\tilde{K}}_{2}\) and \(\mathbf{K}_{2}\) are identical. For this initialization, the filters approximate the complementary-phase tight-framelet criterion.

Fig. 13: Toy model used for experiment on the properties of the filters of a trained CNN. The dimensions for tensors \(\mathbf{K}_{0}\), \(\mathbf{K}_{1}\) and \(\mathbf{K}_{2}\) are \((6\times 1\times 3\times 3)\), \((12\times 6\times 3\times 3)\) and \((24\times 12\times 3\times 3)\). The network is symmetric and the filter dimensions for the decoder convolution kernels \(\mathbf{\tilde{K}}_{n}\) are the same as there corresponding encoding kernel \(\mathbf{K}_{n}\).
imation mechanism proposed the theory of deep convolutional framelets.

The following conclusions can be drawn from this experiment. First, the filters of the CNN may not necessarily converge to complementary-phase tight framelets. This is possibly caused by the initialization of the network and/or the interaction/co-adaptation between the multiple encoder/decoder layers. Second, we confirm that for our experimental setting, the low-rank approximation behavior in the CNN can be observed. For example, when setting the biases and thresholds to zero, part of the noise texture (high-rank information) is recovered. Third, it is possible that linear filtering happens in the network as well, which may explain why noise texture is not fully recovered when setting the biases to zero. Fourth and final, we have observed that the behavior of the trained models change drastically depending on factors such as the learning rate and the initialization values of the model. For this reason, we consider this experiment and its outcome more as a proof of concept, where further investigation is needed.

### **Generalization**

From the explanations in Section IV-C, it can be noted that the bias/threshold used in CNNs can modulate how much of the signal is suppressed by the nonlinearities. In addition, Section V-D established that there are additional mechanisms for noise reduction within the CNN, such as the Wiener-like behavior observed by Mohan _et al._[14]. This raises the question how robust conventional CNNs are to noise-level changes different from the level that the model has been trained with. To perform such experiment, we have trained two variants of the toy model. The first variant is inspired by the multi-scale sparse coding network by Mentl _et al._[17], where the biases of each of the nonlinearities (ReLU in this case) are multiplied by an estimate of the standard deviation of the noise. In the design of this example, the noise estimate \(\hat{\sigma}_{\eta}\), which in accordance to Chang _et al._[1] is defined by

\[\hat{\sigma}_{\eta}=1.4826\cdot\text{Median}(|\mathbf{f}_{\text{HH}}*\mathbf{ x}|). \tag{96}\]

Here, variable \(\mathbf{f}_{\text{HH}}\) is the diagonal convolution filter of the discrete wavelet transform with Haar basis. For comparison purposes, we will refer to this model as _adaptive_ toy model. The second variant of the toy model being tested, examines the case where the convolution layers of the model do not add bias to the signal. This model is based in the so-called bias-free CNNs proposed by Mohan _et al._, in which the bias of every convolution filter is set to zero during training. This setting has the purpose of achieving better generalization on the model, since it is claimed that this modification causes the model to behave independently of the noise level.

We have trained the described variants of the toy models with the same settings of the experiment in Section VII-A. The three models are evaluated on the test image with varying noise levels \(\sigma_{n}\in[0.100,0.150,0.175,0.200,0.225]\), the result for this evaluation is displayed in Fig. 15. These results confirm that the performance of the original toy model degrades for higher noise levels. In contrast, the adaptive and bias-free toy

Fig. 14: Processed “cameraman” image for (in)dependently sampled initialization for the encoding and decoding filters. The top-left picture represents the noise-contaminated input (\(\sigma_{\eta}=0.1\)) and the bottom-left, the noiseless reference. The middle-column images are the processed noisy image with the toy model trained with different initialization for its convolution filters, while the right-column images are processed with the model where the same initial values are used for the encoding and decoding filters. The top-middle and top-right images are nearly identical in terms of quality and SNR, so that initialization has no effect. The middle and bottom-right images are the same model presented that processed the middle and top-right figures, but where its bias is set to zero. As expected, the noise is partly reconstructed.
model perform better than the original toy model for most noise levels.

The results of this experiment confirm the diverse noise-reduction mechanisms within a CNN, as well as showing that the CNNs have certain modeling limitations. For example, noise invariance, which can be addressed by further incorporating prior knowledge to the model, such as the case of the adaptive model, or by forcing the model to have a more Wiener-like behavior such as the case of the bias-free model. In the case of the bias-free model, note that, theoretically, it should be possible to obtain exactly the same behavior with the original toy model if the biases of the model would have converged to zero. This reasoning suggests that the large amount of free parameters and non-linear behavior of the model can potentially prevent to find the optimal/robust solution, in which case the incorporation of prior knowledge can help to improve the model.

## VIII **Which network fits to my problem?**

### **Design elements**

When choosing or designing a CNN for a specific noise-reduction application, multiple choices and design elements should be considered. For example, the required performance, the memory required to train/deploy models, if certain signal preservation characteristics are required, the target execution time for the model, the characteristics of the images being processed, etc. Based on these requirements diverse design elements of CNNs can be more or less desirable, for example, the activation functions, the use of single/multi-resolution models, the need for skip connections, and so forth. This section briefly discusses such elements by focusing on the impact that such elements have in terms of performance

Fig. 15: Comparison of the baseline (original) toy model against its adaptive and bias-free variants. The models are evaluated in the cameraman picture with increasing noise levels. The top row displays the noisy input. The second top-row represents the images processed with the original toy model. Meanwhile, the third row are the results of the adaptive toy model. Finally, the bottom-row are the results corresponding to the bias-free model. It can be observed that the performance original toy model degrades as the noise level increases, while the performance adaptive and bias-free model degrade less with increased noise levels, resulting in pictures with lower noise levels.
and potential computational cost. A summary of the main conclusions of these elements is included in Table II.

#### V-A1 **Nonlinearity**

In literature, the most common activation function in CNNs is the ReLU. There are two main advantages of the ReLU with respect to other activations. First, ReLUs potentially enforce more sparsity in the feature maps than -for example- the soft shrinkage, because ReLUs cancel not only small values of the feature maps like the shrinkage functions do, but also all the negative values. The second advantage of the ReLU is its capacity to approximate other functions (see Section V-C). Note that the high-capacity of the ReLU to represent other functions [31, 13] (often referred to as _expressivity_) may also be one of the reasons why these models are prone to overfitting.

The better expressivity of the ReLU-CNNs may be the reason why -at the time of writing this manuscript- ReLU-based CNNs perform marginally better than the shrinkage-based models in terms of metrics such as signal-to-noise ratio or the structural similarity index metric [39, 40, 19]. Despite this small benefit, the visual characteristics of estimates produced by ReLU and shrinkage-based networks are very similar. Furthermore, the computational cost of ReLU-based designs is potentially higher those with shrinkage functions, because ReLUs require more feature maps to preserve the signal integrity. For example, the LWFSN shown in Section VI-C achieves a performance very close to the FBPConvNet and the tight-frame U-Net for noise reduction in computed tomography, but only with a small fraction of the total trainable parameters, which allows for a faster and less computation-expensive model [19].

As concluding remark, it can be noted that regardless of the expressivity of the ReLU activation, it is not entirely clear if this means that ReLU activations outperform other functions such as the soft threshold _in general_. Because we could not find articles, which are specifically focused on comparing the performance of ReLU/shrinkage-based models. In spite of this, there are some works that compare shrinkage-based CNNs with other (architecturally different) models based on ReLUs that indicate that the compared ReLU-based designs slightly outperform the shrinkage-based ones. For example, Herbetaeu and Kevran [40] proposed the so-called DCT2-Net, which is a shrinkage-based CNN, which despite of its good performance, it is still outperformed by the ReLU-based DnCNN [7] CNN. A similar behavior was observed by Zavala _et al._[19], where their shrinkage-based LWFSN could not outperform the ReLU-based FBPConvNet [21] and the tight-frame U-Net [6]. Another similar case is the deep K-SVD network [39], which achieves performance close (but slightly less good) than the ReLU-based DnCNN. Among the few examples were we found that a ReLU CNN performed better than shrinkage-based models, we find the work of Fan _et al._[18], where they compared variants of the so-called soft autoencoder and found that the shrinkage-based model outperformed the ReLU variant.

#### V-A2 **Single/multi-scale designs**

Single-scale models have the advantage that they avoid aliasing because no down/up-sampling layers are used. Still this comes at the cost of more computations and memory. Furthermore, this approach may lead to models with larger filters and/or deeper networks to achieve the same receptive-field than multi-scale models, which may further increase the computation costs of single-scale models.

In the case of multi-scale models, the main consideration should be that the down/up-sampling structure should allow perfect signal reconstruction to avoid introducing aliasing and/or distortion to the image estimates (e.g. the discrete wavelet transform in the tight-frame U-Net and in the learned wavelet-frame shrinkage network).

#### V-A3 **(Non-) residual models**

Residual noise-reduction CNNs often perform better than their non-residual counterparts (e.g. the U-Net vs FBPConvNet and the LWFSN vs the rLWFSN). This may be caused because the trained models have more freedom to learn the filters, because the design does not need to learn to reconstruct the noiseless signal, but only to estimate the noise [12]. Also, it can be observed that non-residual models potentially need more parameters than residual networks, because the propagation/reconstruction of the noiseless signal is also dependent on the number of channels of the network.

### **State-of-the art**

Defining the state-of-the art in image denoising with CNNs is challenging for diverse reasons. First, there is a wide variety of available CNNs, which often are not compared to each other. Second, the suitability of a CNN for a given task may depend on the image and noise characteristics, such as noise distribution and (non-) stationarity. Third, the large amount of variables, in terms of e.g. optimization, data and data augmentation, adds reproducibility issues which further complicate making a fair comparison between all available models [11]. In addition, it should be noted that for many of
the existing models, the performance gap between state-of-the art models and other CNNs is often small.

Despite of the mentioned challenges, we have found some models that could be regarded as the state-of-the art. The first model to be addressed is the so-called DRU-Net [41], which is a bias-free model [14] that incorporates a U-Net architecture with residual blocks. In addition, the DRU-Net uses an additional input to indicate to the network the noise intensity, which increases its generalization to different noise levels. An additional state-of-the art model is DnCNN [7]. This network is residual, single-scale, while also using ReLU activations. Another state-of-the art model is the multi-level-wavelet CNN [42], which has a design very similar to the tight-frame U-Net [6]. Both of these models are based on the original U-Net design [34], but are deployed in residual configuration and the down/up-sampling structure is based on the discrete wavelet transform. Furthermore, in addition to use encoding-decoding CNNs stand-alone, have used CNNs as proximal operator within model-based methods [43, 41], which further improves the denoising power of non-model-based encoding-decoding CNNs.

## IX **Conclusions and Future Outlook**

In this paper, the widely used encoding-decoding CNN architecture has been analyzed from several signal processing principles. This analysis has revealed the following conclusions. (1) Multiple signal processing concepts converge in the mathematical formulation encoding-decoding CNNs models. For example, the convolution and down/up-sampling structure of the encoder-decoder structure is akin to the framelet decomposition, the activation functions are rooted on classical signal estimators. In addition, linear filtering also may happen within the model. (2) The activations implicitly assume noise and signal characteristics of the feature maps. (3) There are still many signal processing developments that can be integrated with current CNNs, further improving their performance in terms of accuracy, efficiency or robustness.

Despite of the signal processing nature of the encoding-decoding CNNs, at the moment of this publication, the integration of CNNs and existing signal processing algorithms is at an early stage. A clear example of the signal modeling limitations of current CNN denoisers are the activation functions, where the estimators provided by current activation layers neglect spatial correlation of the feature maps. Possible alternatives to solve this limitation could be to perform an activation function inspired on denoisers working on principles such as Markov random fields [44], the locally-spatial indicators [45] and multi-scale shrinkage [24]. Further ideas are provided by the extensive survey on denoising algorithms by Pizurika and Philips [46]. Additional approaches that can be further explored are non-local [47] and collaborative filtering [48]. Both techniques exploit the redundancy in natural images and only a few models are exploring these properties [49, 50].

To finalize this manuscript we encourage the reader to actively consider the properties of the signals processed, the design requirements and the existing signal processing algorithms when designing new CNNs. By doing so, we expect that next the generation of CNNs denoisers not only will be better performing, but also more interpretable, and reliable.

## Acknowledgement

We thank Dr. Ulugbek Kamilov and the anonymous reviewers for their valuable feedback and suggestions for this article.
