Unsupervised Image Restoration Using Partially Linear Denoisers

Natural images are non-Gaussian signals with fine structures such as sharp edges and textures. One of the main challenges in the image denoising task is to preserve these structures while removing the noise. Traditional linear denoisers such as Gaussian filtering and Tikhonov regularization usually do not achieve satisfactory denoising quality, as they tend to smooth the edges. Most of the existing image denoisers in the literature are therefore nonlinear.

Total variation (TV) denoising [30] is one of the most fundamental image denoising algorithms. The TV denoising solution is a minimizer of the optimization problem

in which $y\in\mathbb{R}^{m}$ is the noisy image and $\|x\|_{\rm TV}:=\|\nabla x\|_{1}$ is the total variation of $x\in\mathbb{R}^{m}$. TV denoising is known to preserve sharp edges thanks to the non-linearity introduced by the TV norm $\|\cdot\|_{\rm TV}$.

Patch based methods have gained popularity for image denoising tasks because of their capacity of capturing the self-similarity of images. The non-local means (NL-means) algorithm proposed by Buades et al. [3] is among the most successful methods in this category. The NL-means algorithm removes noise by calculating a weighted average of the pixel intensities, with weights defined based on a patch similarity measure which emphasizes the connection among pixels in similar patches. It is a nonlinear denoiser, different from the local mean filtering, as the weights are image-dependent. As another nonlinear patch based denoising method, the Block-matching and 3D filtering (BM3D) algorithm [6] divides image patches into groups based on a similarity criterion, and collaborative filtering on each of the groups is then performed to clean the patches. While patch based denoisers achieve promising denoising performance, they are often time-consuming due to the high computational complexity in calculating the weights or matching with similar patches.

In recent years, convolutional neural networks (CNN) based denoisers became the state of the art for image denoising [41, 42], thanks to the rapid development of deep learning techniques. In particular, CNN are known to be efficient in modeling image priors [35], which are crucial for the quality of the denoiser. A typical CNN denoiser can be formulated as a composition of mappings called layers, and a basic type of layers $y^{(k)}\rightarrow y^{(k+1)}$ has the form

where $*$ denotes the convolution operation, $(W,b)$ are parameters of the network that can be learned from the data, and $\sigma{\left(\cdot\right)}$ is an activation function which is often nonlinear. As such, CNN denoisers are in general nonlinear. There has been growing interest in developing CNN denoisers with new network architectures and building blocks being proposed, such as batch-normalization [41], residual connection [10], and residual dense block [42].

While deep CNN based models have great advantages in image denoising, the most standard learning techniques are limited to the availability of sufficiently many noise-free images. Recently, CNN-based learning algorithms for image denoising that do not require clean images as training data have been developed. Soltanayev et al. [32] propose to use Stein’s unbiased risk estimator (SURE) based loss function, which is computed without knowing ground truth images, as a replacement of the MSE loss function. In the setting of i.i.d. Gaussian noise, SURE provides an unbiased estimate of the MSE, and hence the minimization problems with respect to these two losses are equivalent. Consequently, CNN denoisers can be trained in an unsupervised manner, based on only one noisy realization of each training image. Nevertheless, SURE may not be identical to the MSE in the case of non-Gaussian noise, e.g., shot noise. The SURE training scheme has been extended to the setting of multiple noise realizations per image as well as imperfect ground truths [43].

The Noise2Noise approach developed by Lehtinen et al. [20] offers a different ground-truth free learning strategy for denoising, based on the assumption that for each image at least two different noisy observations are available. It is found that replacing the clean targets by their noisy observations in the MSE loss function leads to the same minimizers of the original supervised loss if the noise has zero mean and an infinite amount of training data are provided. Specifically, the Noise2Noise loss function is

where $y$ and $\hat{y}$ represent two independent noisy observations of the same image sample, and $f_{\theta}$ is the denoiser parameterized by $\theta$. In contrast to the MSE loss $\sum_{{\left(y,x\right)}\in{\left(\mathcal{Y},\mathcal{X}\right)}}\|f_{\theta}\mathopen{}\left(y\right)\mathclose{}-x\|^{2}$ where $x$ denotes the ground truth image, the loss function (1) is defined on a noisy pair ${\left(y,\hat{y}\right)}$ only. If the noise in $y$ and the noise in $\hat{y}$ are independent and have zero means, it is shown that the minimization of this loss function (1) with respect to $\theta$ is equivalent to minimizing the MSE [20]. This implies that the parameter $\theta$ can be computed from a training set of noisy pairs. Besides, the Noise2Noise approach does not rely on an explicit image prior or on structural knowledge about the noise models.

In certain denoising tasks, however, the acquisition of two or more noisy copies per image can be very expensive or impractical, in particular in medical imaging where patients are moving during the acquisition, or in videos with moving cars, etc. Several authors have investigated learning techniques that overcome this restriction [8, 36, 2, 15, 16]. The Frame2Frame [8] developed by Ehret et al. fine-tunes a denoiser for blind video denoising. The idea is to use optical flow to warp one video frame to match its neighboring frame, and then minimize the Noise2Noise loss (1) with $(y,\hat{y})$ being the pair of matched frames. The work [27] generalizes the Noise2Noise into the setting of a single noisy realization for each image. A synthetic noise, that is drawn from the same distribution as the underlying noise, is added to the noisy image $y$, and the new noisy image is then used to replace the second observation $\hat{y}$ in the training loss (1). At test time, the raw output of the network is post-processed to obtain the denoising results, by computing a linear combination of the output and the input. In this approach, the synthetic noise has to be of the same noise type as the underlying noise, and the post-processing can magnify the errors of the output of the network. Batson et al. [2] show that a denoiser $f_{\theta}{\left(\cdot\right)}$, satisfying a so called $\mathcal{J}$-invariant property for a partition $\mathcal{J}$ of the image pixels, can be trained without accessing a second observation of the image if the noise is independent across different regions in $\mathcal{J}$. A function $f_{\theta}\mathopen{}\left(\cdot\right)\mathclose{}$ is said to be $\mathcal{J}$-invariant if for each subset of pixels $J\in\mathcal{J}$, the pixel values of $f_{\theta}{\left(y\right)}$ at $J$ can be calculated without knowing the values of $y$ at $J$. By additionally assuming that the noises of $y$ at different elements of $\mathcal{J}$ are independent and have zero mean, one minimizes the loss

where the subscript $J$ in $y_{J}$ denotes the restriction of the image $y$ to the pixel collection $J$. It should be noted that, for a $\mathcal{J}$-invariant function $f_{\theta}{\left(\cdot\right)}$ and any $J\in\mathcal{J}$, the loss $\sum_{y\in\mathcal{Y}}{\left\|[f_{\theta}{\left(y\right)}]_{J}-y_{J}\right\|}^{2}$ can be interpreted as a variant of (1), given the fact that the noise contained in $y_{J}$ is independent of $[f_{\theta}{\left(y\right)}]_{J}$ and has zero mean. This enables the training of a model using a set of noisy images $y$ only. However, as the denoiser is a $\mathcal{J}$-invariant function, the images can not be perfectly reconstructed as the level of noise decreases to zero, i.e., it can not learn the identity mapping which is clearly not $\mathcal{J}$-invariant. The unused information from $y$ can be further leveraged if the noise model is known or can be estimated. Krull et al. [16] propose to combine $y$ with the network predictions, based on a probabilistic model for each pixel and a Bayesian formulation, to obtain the minimal MSE estimate. The integration of statistical inference effectively removes the noise remained in the predictions. A downside of this method is that it requires an explicit expression of the posterior distribution of the noisy images beforehand, and it is not end-to-end because the network outputs intermediate results that are improved in the post-processing step.

The motivation for this paper comes from the fact that, in practice the distribution of $\bm{x}$ is often unknown, and samples of $\bm{x}$ (i.e., the ground truth images) or of the noise $\bm{n}$ are not readily available. What can be easily accessed are noisy observations $\bm{y}$. With these alone, however, a direct evaluation of the MSE is not possible. In this work, we circumvent the need for $\bm{x}$ by introducing an auxiliary random vector and replacing the MSE $\mathcal{J}_{0}$ by an approximation. First, let $\bm{z}$ be a random vector satisfying assumptions

• (A2)

  the conditional mean $\mathbb{E}{\left(\bm{z}\mid\bm{x}\right)}=0$,

• (A3)

  $\bm{z}$ is independent of $\bm{n}$,

• (A4)

  the conditional covariance:

  $${\rm Cov}{\left(\bm{z},\bm{z}\mid\bm{x}\right)}={\rm Cov}{\left(\bm{n},\bm{n}\mid\bm{x}\right)}.$$

The auxiliary vector $\bm{z}$ does not necessarily need to follow the same distribution as $\bm{n}$. Samples of $\bm{z}$, therefore, can be easily generated from e.g., Gaussian distributions, once the variances of $\bm{z}$ are known. Then, associated with $\bm{z}$, we define

where $\alpha$ is a real constant. The new random vector $\bm{\hat{y}}=\bm{x}+\bm{\hat{n}}$ can be regarded as a noisy version of image $\bm{x}$ with noise vector $\bm{\hat{n}}$. In the following, the discussion will focus on denoising $\bm{\hat{y}}$ rather than denoising $\bm{y}$, but the objective remains unchanged, i.e., getting the same clean image $\bm{x}$. Specifically, if one can obtain a high-quality solution $\bm{x}$ from $\bm{\hat{y}}$, then there exists an algorithm taking $\bm{y}$ to the clean image because $\bm{\hat{y}}$ can be computed from $\bm{y}$. Such an algorithm can be achieved, for instance, by $R{\left(V{\left(\bm{y}\right)}\right)}$ where $R$ is a denoiser for the $\bm{\hat{y}}$ problem and $V(\bm{y}):=\bm{\hat{y}}$. It should be noted that the difficulty of the denoising problem is raised because of the additional uncertainty from the auxiliary random vector $\bm{z}$ encoded in the observed data $\bm{\hat{y}}$. However, one of the benefits of considering $\bm{\hat{y}}$ is that the quantity $\bm{z}$ is known and can be leveraged for constructing approximations to the MSE without knowing $\bm{x}$ or $\bm{n}$ as we will show in the following.

For the $\bm{\hat{y}}$ denoising problem, the MSE associated with the denoiser $R$ is defined as

where the expectation $\mathbb{E}$ is taken over random variables $\bm{x}$, $\bm{n}$ and $\bm{z}$. This is connected to $\mathcal{J}_{0}$ via $\mathcal{J}_{\rm mse}(R)=\mathcal{J}_{0}(R{\left(V\right)})$. Since $V$ is known, it can be shown that $\min_{R}\mathcal{J}_{0}(R)\leq\min_{R}\mathcal{J}_{\rm mse}(R)$. Nevertheless, if $\alpha$ is close to zero, the noise distributions of $\bm{y}$ and $V(\bm{y})$ are close, so the minimum of $\mathcal{J}_{0}$ can be approximated by the minimum of $\mathcal{J}_{\rm mse}$ which promises similar denoising quality.

Now, using the auxiliary vector $\bm{z}$, we define the following objective function for our proposed partially linear denoiser

for $\alpha\neq 0$. Indeed, as we will see later in Subsections III-B and III-C, if we consider the partially linear denoising model, then $\mathcal{J}$ provides a good estimate of $\mathcal{J}_{\rm mse}$. More precisely, according to the definition (4), we consider a set of denoisers for $\bm{\hat{y}}$ that have the form

where $\bm{\hat{n}}$ is defined in (5), and $L$ and $\bm{e}$ are depending on $\bm{x}$ and ${\left(\bm{x},\bm{\hat{n}}\right)}$, respectively, and the residual $\bm{e}$ is of small variance.

In the rest of this section, we show that for $\alpha\!\neq\!0$ and partially linear denoiser $R$, the term $\mathcal{J}{\left(R\right)}$ in (7) approximates the MSE $\mathcal{J}_{\rm mse}{\left(R\right)}$ up to an additive constant. This implies that an optimal denoiser of this class can be computed even if the distribution of ground truth images $\bm{x}$ and the distribution of $\bm{n}$ are not known. We then discuss unsupervised learning methods based on our proposed objective $\mathcal{J}{\left(R\right)}$ for image restoration tasks like denoising and deblurring.

To start with, we focus the discussion on the simplest case, in which the denoiser $R$ is linear. This is a subset of the set of partially linear denoisers (8) with $g=L$ and $\bm{e}=0$. The best estimator in this setting is the linear minimum mean square error (LMMSE) estimator, i.e., the minimizer of $\mathcal{J}_{\rm mse}{\left(R\right)}$ in (6) over all linear denoisers $R$. The following proposition establishes the equivalence between the MSE loss and $\mathcal{J}{\left(R\right)}$ in (7) over linear denoisers.

According to Proposition 1, for all linear denoisers, the term $\mathcal{J}{\left(R\right)}$ in (7) differs from the MSE by a constant $c$ not depending on $R$. Therefore a minimization of $\mathcal{J}{\left(R\right)}$ over all linear $R$ also leads to the LMMSE estimator. We remark that the objective function (7) is not defined based on the ground truth image $\bm{x}$, and the distribution of the random vector $\bm{n}$ is not necessarily known or equal to that of $\bm{z}$. In fact, the denoising problem can be reformulated as follows. Given $R$, the quantity $R\bm{\hat{y}}$ is known, and we want to decouple the term $R\bm{\hat{y}}-\bm{x}$ from

In the loss $\mathcal{J}{\left(R\right)}$, this is implemented by adding the scaled auxiliary vector $\bm{z}/\alpha$ to both sides of (14). The vector $\bm{z}/\alpha$ compensates the noise $\bm{n}$ in the sense that, when taking the expectation of ${\left\|R\bm{\hat{y}}-\bm{y}+\bm{z}/\alpha\right\|}^{2}$ over $\bm{n}$ and $\bm{z}$, the $\bm{n}$-related term $\mathbb{E}{\left(\langle R{\bm{n}},\bm{n}\rangle\mid\bm{x}\right)}$ is canceled by the $\bm{z}$-related term $\mathbb{E}{\left(\langle R{\bm{z}},\bm{z}\rangle\mid\bm{x}\right)}$ (i.e., the last equality in (12)), which leads to (9).

Though linear denoisers have nice properties that link (7) to the MSE, they may not be the best denoisers for imaging data. Nonlinearity is unavoidable in order to achieve good denoising quality and preserve fine structures such as edges in the image. To this end, we consider a more general class of denoisers that are partially linear, as expressed in (8), where $g$ is nonlinear and $\bm{e}$ is not necessarily zero.

In a special case where the denoiser outputs exactly the clean image, i.e. $R{\left(\bm{\hat{y}}\right)}:=\bm{x}$, $R$ is also partially linear with $L$ and $\bm{e}$ being both zero. However, such denoisers do not exist in most settings. Nevertheless, for good denoisers one could still expect the residual term $\bm{e}$ to be small.

The equivalence of the two objective functions stated in Proposition 1 does not hold for nonzero $\bm{e}$. It is therefore of interest to know how well $\eqref{eq:newmse}$ approximates the MSE in this general case. The next proposition quantifies the approximation error in the presence of a small nonzero residual term $\bm{e}$.

From the bounds given in Proposition 2, if the variance of $\bm{e}$ is small, then the error (7) provides a good estimate to the MSE. Note that making a Lipschitz continuity assumption on $\bm{e}$, the error bound in (16) can be further improved to (17) which is independent of $\alpha$ (in contrast to (16) where the bound has a $1/\alpha$ factor). In practice, this Lipschitz continuity assumption is not restrictive as we expect that the denoiser $R$ is stable with respect to small perturbations in $\bm{\hat{y}}$ and therefore has a small Lipschitz constant $K$ for $\bm{e}$.

Next, we take a closer look at the errors of the proposed denoisers. We first introduce the set $\mathcal{R}_{\epsilon}$ of partially linear denoisers. Second, the distance between the minimizers of $\mathcal{J}$ and $\mathcal{J}_{\rm mse}$ over $\mathcal{R}_{\epsilon}$ is estimated. Then we verify the convergence of the proposed denoisers to the best denoisers for corrupted images generated from a single ground truth image. Finally, focusing on a special subset of images that consists of constant patches, we demonstrate the partial linearity of the optimal denoisers with an example.

For $\epsilon>0$, let $\mathcal{R}_{\epsilon}$ be the set of denoiser satisfying (8) with variance of $\bm{e}$ less than or equal to $\epsilon^{2}$, i.e.,

We are interested in the best denoisers in $\mathcal{R}_{\epsilon}$ with as small MSE as possible. More precisely, for small $\delta\geq 0$, we define the MMSE denoiser, denoted by $R^{\epsilon,\delta}$, as one of the denoisers in $\mathcal{R}_{\epsilon}$ that satisfies

where $\mathcal{J}_{\rm mse}$ are defined in (6). Note that $R^{\epsilon,\delta}$ exists for arbitrarily small positive $\delta$. Based on the approximation properties stated in Proposition 2, we propose a denoiser $\hat{R}^{\epsilon,\delta}\in\mathcal{R}_{\epsilon}$ satisfying

as an alternative version of $R^{\epsilon,\delta}$. In light of the theoretical bound (16), the distance between $R^{\epsilon,\delta}$ and $\hat{R}^{\epsilon,\delta}$ has an upper bound described in the following proposition.

propositionpropositiond Let the denoisers $R^{\epsilon,\delta}$ and $\hat{R}^{\epsilon,\delta}$ be given by (22) and (23) respectively. If $\bm{n}$ and $\bm{z}$ satisfy Assumption (A1) - (A4) and $0<\delta<1$, it holds that

where the expectation is taken over $\bm{x}$, $\bm{n}$ and $\bm{z}$. The proof is given in Appendix A. The bound in Proposition 3 converges to $2\epsilon\sqrt{\mathbb{E}({\left\|\bm{n}-\bm{z}/\alpha\right\|}^{2})}$ as $\delta$ tends to $0$, and it converges to $0$ as $\epsilon$ and $\delta$ converge to zero.

Building on Proposition 3, we further analyze the convergence of the proposed solutions for a special case where all realizations of the corrupted images $\bm{\hat{y}}$ are generated from the same clean image. The next corollary shows that, in this setting, the approximation error is at most of the same order as $\epsilon^{2}$, regardless of the noise distributions or the noise variance. The proof of this corollary is given in Appendix A. {restatable}corollarypropositiondb Assume that the conditions in Proposition 3 hold. If $\bm{x}$ follows a delta distribution, i.e., the probability $P{\left(\bm{x}=x\right)}=1$ and $P{\left(\bm{x}\neq x\right)}=0$ for some image $x$, then

Proposition 3 and Corollary 3 illustrate that the minimizers of $\mathcal{J}$ are good estimates of the best partially denoiser when $\epsilon$ is small. In principle, $\bm{e}$ should be much smaller than the noise itself, as this is a necessary condition for the denoiser being robust to different realizations of $\bm{\hat{n}}$. We can provide some more insight into the size of the residual term for high quality denoisers when applied to a special subset of images. Smooth regions are one of the most important components in natural images, and they often account for a large fraction of the pixels. Here we study the basic case with constant images and pixelwise independent noise. We give an example to demonstrate the scale of the residual term of the minimum MSE denoiser for constant image patches.

For learning the denoiser we consider parametrized denoising models $R_{\theta}$, e.g., deep CNNs, that are parameterized by $\theta$. Suppose that a set $\left\{y_{i}\right\}$ of realizations of $\bm{y}$, e.g., the set of noisy images that one wants to denoise, is given. The noise in these noisy images satisfies Assumption (A1). Associated with each $y_{i}$, we define $z_{i}$ as a realization of the auxiliary random vector $\bm{z}$ that satisfies (A2) - (A4). In (A2) and (A4) only the first two moments of $\bm{z}$, rather than its distribution, are specified. Therefore without loss of generality, one can randomly generate $z_{i}$ from normal distributions. With the auxiliary vector $z_{i}$, the samples of $\bm{\hat{y}}$, denoted by $\hat{y}_{i}$, are computed according to Equation (5). Based on the objective function (7), we minimize the empirical loss function

under the condition that $R_{\theta}\in\mathcal{R}_{\epsilon}$. This condition can be implemented with the partial linearity constraint described below.

Partial linearity constraint. To preserve the partial linearity during training, restriction on the variance of the residual term $\bm{e}$ in (8) is required. Unfortunately, a direct evaluation of $\bm{e}$ is not feasible due to 1) the fact that the formulation of $\bm{e}$ depends on $g{\left(\bm{x}\right)}$ and $L$ which are unknown, and 2) the condition that only one noisy observation for each image is given. However, since $g{\left(\bm{x}\right)}$ and $L$ depend only on the image $\bm{x}$ and $L$ is linear, we can remove these two terms by perturbing $\bm{\hat{y}}$. One simple way of implementing this is to find $\bm{q}_{1}$ and $\bm{q}_{2}$, two perturbed versions of $\bm{\hat{y}}$, satisfying $\bm{\hat{y}}=\tau\bm{q}_{1}+(1-\tau)\bm{q}_{2}$ where $\tau\in[0,1]$. Let the residual terms associated with $R{\left(\bm{q}_{1}\right)}$ and $R{\left(\bm{q}_{2}\right)}$ be denoted by $\bm{e}_{1}$ and $\bm{e}_{2}$ respectively. Then, from Equation (8) we have

By the assumption that $\bm{e},\bm{e}_{1}$ and $\bm{e}_{2}$ are small, based on the representation of the residual terms in (26), one can therefore penalize $\psi{\left(R{\left(\bm{\hat{y}}\right)}-\tau R{\left(\bm{q}_{1}\right)}-(1-\tau)R{\left(\bm{q}_{2}\right)}\right)}$ for some metric $\psi$.

The proposed partially linear denoisers can be generalized to learning image reconstructors for other linear inverse problems such as image deblurring, in particular when the ground truth images are missing. We repeat the observation model (3) here

where $A$ is a linear forward operator (such as the Radon transform for CT reconstruction, a convolution operator for image deblurring, etc), and $\bm{n}$ models the zero-mean noise. For many inverse imaging problems, the naive reconstruction $A^{\dagger}\bm{y}$, is based on a direct inversion of the linear operator $A$ given by the pseudo-inverse operator $A^{\dagger}$. Due to ill-conditioning of $A$ such inversion is typically unstable and therefore not accurate in the presence of noise. If $A$, however, has full column rank, then the resulting artifacts in the reconstruction could be remedied by denoising $A^{\dagger}\bm{y}$. In fact, it is easy to see that $\bm{n}^{\dagger}:=A^{\dagger}\bm{y}-\bm{x}=A^{\dagger}\bm{n}$ has zero mean. The covariance of $\bm{n}^{\dagger}$ can be estimated as long as the covariance of $\bm{n}$ is known. Therefore a straightforward application of the partially linear denoisers to the noisy reconstruction $A^{\dagger}\bm{y}$ leads to an estimate of $\bm{x}$. Various approaches have been proposed for denoising $A^{\dagger}\bm{y}$ in a data-driven post-processing setup (see e.g., [13, 11]). These learned post-processing techniques require noisy-clean image pairs as training data. Another class of methods, in contrast, leverage denoisers to construct regularization for (27). For instance, the Regularization by Denoising (RED) methods [28, 12] minimize a variational objective function with a regularization term derived from the denoisers. The Plug-and-Play Prior framework [37, 33] is based on variable splitting algorithms for the Maximum a Posteriori (MAP) optimization problem with the proximal mapping associated with the regularization term being replaced by the denoisers. These denoiser-based regularization approaches require predefined denoisers. However, our framework does not require any noisy-clean pair or any predefined denoiser. In particular, our proposed partially linear denoiser allows training an estimator from noisy data alone.

Image deblurring is a special case of (27) where $A$ represents a convolutional kernel. In practice, $A$ is governed by different imaging factors including motion and camera focus. In the context of the proposed partial linear denoiser, one way to solve the deblurring problem is to train a single model $R_{\theta}$ that maps the measurement $\bm{y}$ directly to $\bm{x}$. Here, knowledge of $A$ is indirectly encoded in $R_{\theta}$. Assuming that the training data contains noisy samples of $\bm{\hat{y}}$ and the blurring kernel $A$ but no ground truth images, then similar to (7) one can minimize the loss function

and in this case, $A{R_{\theta}}$ acts as a partially linear denoiser. At test time, the deblurred image can be directly computed as ${R_{\theta}}{\left(\bm{\hat{y}}\right)}$ without knowing the operator $A$. These considerations also draw connections of the partial linear denoiser to the deep image prior approach [35]. There, an implicit regularizer is introduced, based on a convolutional neural network $R_{\theta}$ and a sole data-fitting loss function is minimized in the training. Early stopping is applied to prevent the estimator from over-fitting to the noise. With the partial linearity structure, however, we establish a connection between the cost (28) and the MSE of the noise free measurement $A\bm{x}$ which aims to get $A{R_{\theta}}{\left(\bm{\hat{y}}\right)}$ close to $A\bm{x}$ rather than its noisy versions.

In this test, we investigate the partial linearity of some existing standard denoising approaches, including the TV approach [30], BM3D [6] and DnCNN [41]. For the convenience of the readers, the formula of the partially linear denoiser (7) is repeated here

where $R$ is the underlying denoiser. In particular, the DnCNN is trained by minimizing the standard empirical MSE loss (6) with a training set of $400$ images with ground truth [41]. For a given image $\bm{x}$, we compute the decomposition (29) using Remark 1. Specifically, $g(\bm{x})\!=\!\mathbb{E}{\left(R{\left(\bm{\hat{y}}\right)}\mid\bm{x}\right)}$ and $L\!=\!\min_{L_{0}\in\mathcal{L}}\mathbb{E}{\left({\left\|R{\left(\hat{\bm{y}}\right)}-g{\left(\bm{x}\right)}-L_{0}\hat{\bm{n}}\right\|}^{2}\right)}$ where the expectation is taken over $\hat{\bm{n}}$, and $\mathcal{L}$ denotes the set of linear operators. In this example, we use a standard test image called parrot (cf. top middle of Fig. 1) as the ground truth image, i.e., a realization of $\bm{x}$. The noise $\bm{\hat{n}}$ is i.i.d. Gaussian noise with zero mean and standard deviation $25$ (corresponding to 256 gray levels). The expectations are evaluated using $20000$ random realizations of the pair ${\left(\bm{\hat{y}},R{\left(\bm{\hat{y}}\right)}\right)}$.

For the three denoisers, we report the variance of $\bm{e}$ averaged over the image pixels of the parrot image (i.e., $\epsilon^{2}/m$ where $m$ denotes the number of pixels in $\bm{e}$) in Table I. A comparison of the accuracy of the methods, measured in PSNR, is also given in the table. All figures reported in the table are averaged over $8000$ independent runs with different realizations of $\bm{\hat{y}}$. Based on the table, the DnCNN achieves the best denoising quality, and it outperforms BM3D by around $0.5$ dB and the TV approach by around $1.8$ dB. All three methods have $\epsilon^{2}/m$ less than $10^{-4}$, and the value for the CNN denoiser is about half of that of the TV method.

Given a denoiser $R$, in order to study the partial linearity we consider the pair

where $i$ denotes a fixed pixel (localized by the red dot on the first row of Fig. 3). To visualize the linearity for the three denoisers, the pair (30) is plotted on the second row of Fig. 3 under $8000$ realizations of the noise $\bm{\hat{n}}$. Each point in the plot corresponds to one realization. As shown in the figure, all the denoisers have certain degrees of partial linearity at the pixel $i$. For the TV denoiser, the residual $\bm{e}$ has a relatively larger variance (also illustrated by Table I) and there are some outliers for big $|[L\bm{\hat{n}}]_{i}|$.

We consider two types of synthetic noises (Gaussian noise and Poisson noise), and the performance of our method is compared with recent denoising methods that are not trained on ground truth images, such as SURE [32], Noise2Self [2] and Noisier2Noise [27], along with the classic BM3D approach [6]. Throughout the denoising experiments, we use the same network architecture DnCNN [41] for the fully supervised baseline (we will call it DnCNN in the subsequent), SURE, Noise2Self, Noisier2Noise and our approaches. The experimental setups for the five methods are the same, except that the clean images are used to train the DnCNN while they are unseen by the latter four methods in the learning phase. All models are trained on a benchmark denoising dataset [41] consisting of $400$ training images of size $180\!\times\!180$. In the training phase, we feed the CNNs with image patches of size $40\!\times\!40$ and set a batch size of $128$. Augmentations such as random flipping and random cropping are applied to the patches. In the inference phase, the inputs to the networks are the whole noisy images. In particular, in our approach we do not include the auxiliary vector $\bm{z}$ at this stage, and the denoised images are the outputs of the network without any post-processing. We evaluate the denoising quality on two different test image sets, namely the BSD68 (containing $68$ images) [26] and the $12$ wildly used images (Set12) [6].

Training. To train our denoising model, we minimize the loss function (25) using the Adam optimizer [14]. The minimization process consists of two stages. In the first stage, we fix $\alpha\!=\!1$ and minimize the loss function (25) without constraints. In the second stage, we randomly choose $\alpha\!\in\![0.1,0.5]$ for each of the noisy samples, and additionally, in order to control the variance of $\bm{e}$ defined in (8), we impose a partially linear constraint (26). The implementation details of the latter will be given in the next paragraph. Both stages contain $2\!\times\!10^{5}$ optimization steps with an initial learning rate $0.001$, which then drops to $0.0001$ and $0.00005$ at the $6\!\times\!10^{4}$th step and $1.2\!\times\!10^{5}$th step, respectively. Throughout the experiments, $z_{i}$, i.e., the samples of the auxiliary random vector $\bm{z}$ in (25), are generated from Gaussian distributions.

Partially linear constraint. We implement the partially linear constraint by perturbing the noisy images $\hat{y}_{i}$ (i.e., samples of the noisy image $\bm{\hat{y}}$) and by following the formula (26). Specifically, for each $\hat{y}_{i}$ let $\beta_{i}^{\left(1\right)}$ and $\beta_{i}^{\left(2\right)}$ be random numbers uniformly distributed in $[1,1.5]$, and let $q_{i}$ be a perturbation vector randomly generated from the same distribution as $z_{i}$. With $q_{i}$ we construct two perturbed versions of $\hat{y}_{i}$ as $q^{(1)}_{i}\!:=\!\hat{y}_{i}-\beta_{i}^{\left(1\right)}q_{i}$ and $q^{(2)}_{i}\!:=\!\hat{y}_{i}+\beta_{i}^{\left(2\right)}q_{i}$ respectively. Then we have the linear relationship

where $\tau_{1}:=\beta_{i}^{\left(2\right)}/{\left(\beta_{i}^{\left(1\right)}+\beta_{i}^{\left(2\right)}\right)}$ and $\tau_{2}:=1-\tau_{1}$. In practice, we let $q_{i}$ be independent of $z_{i}$. Since the noise levels of $q^{(1)}_{i}$ and $q^{(2)}_{i}$ depend on the pixel values of $q_{i}$, we modify $q_{i}$ to be sparse to avoid raising the noise level too much. To do this, we randomly select $1/25$ of the pixels, the pairwise distances of which are at least $4$ pixels. The remaining pixels of $q_{i}$ are set to zero, i.e., no perturbations are applied to these pixels. To avoid having outliers at some individual pixels caused by the perturbations, we also clip $q_{i}$ such that the pixel values $q_{i}^{\left(1\right)}$ and $q_{i}^{\left(2\right)}$ fit in the range $\left[1.2a-0.2b,1.2b-0.2a\right]$ where $a\!:=\!\min_{j}[\hat{y}_{i}]_{j}$ and $b\!:=\!\max_{j}[\hat{y}_{i}]_{j}$. Having $q_{i}^{\left(1\right)}$ and $q_{i}^{\left(2\right)}$, based on Property (26) we add the following penalty term to the loss function

where $\theta$ are the network parameters, $M$ is a diagonal matrix. The diagonal entries of $M$ are set as

if $[q_{i}]_{j}\!\neq\!0$, otherwise $M_{jj}\!:=\!0$, where $\sigma\!>\!0$ is the square root of the largest pixel-wise variance of the noise $\bm{n}$. Note that having $M$ in the loss (31) means that we penalize the nonlinearity at the perturbed pixels only. The term $0.1\sigma$ prevents division by very small numbers. In summary, the loss function is $\mathcal{L}+\gamma\mathcal{L}_{c}$ where $\mathcal{L}$ is defined in (25), and $\gamma$ is a hyperparameter which can be tuned based on the quality of the denoised images.