Insights into analysis operator learning: From patch-based sparse models to higher-order MRFs

I Introduction

I-B Patch-based analysis operator learning

In the case of the synthesis model, the learning of an optimized dictionary has become ubiquitous. However, in analysis-based models, fixed operators inspired from variational methods such as the discrete total variation have been used for a long time. It is only recently that people started to develop customized algorithms to learn in some sense optimal analysis operators.

For the noise aware case, the objective is to learn an optimal analysis operator $A$, which enforces the coefficient vector $Ax_{i}$ to be sparse, while $\|{x_{i}-y_{i}}\|_{2}^{2}\leq\varepsilon$ for each training sample ($\varepsilon$ is an error tolerance, which is derived from the noise level). This requires solving a problem of the form

Using a Lagrange multiplier $\lambda>0$ this can be equivalently expressed as

where $\phi$ is again a sparsity promoting function and $\|{\cdot}\|_{F}$ denotes the Frobenius norm.

Unfortunately, the above optimization problems suffer from the problem of trivial solutions. Indeed, if no constraints are imposed on $A$, it is easy to see that the trivial solution $A\equiv 0$ is the global minimizer of (I.3), (I-B) and (I.5). A possible solution to exclude the trivial solution is to impose additional assumptions on $A$, i.e., restricting the solution set to an admissible set $\mathcal{C}$. The following constraints have been investigated in and :

row norm constraints. All the rows of $A$ have the same norm, i.e., $\|{A_{i}}\|_{2}=c$ for the $i^{th}$ row of operator $A$.

row norm + full rank constraints. The analysis operator $A$ has full rank, i.e., $rk(A)=m$.

As pointed out in , each individual constraint presented above does not lead to satisfactoy results. Therefore, in a constraint called the Uniform Normalized Tight Frame (UNTF) was proposed, which is a combination of the unit row norm and the tight frame constraint. The authors of employed a constraint combining the unit row norm and the full rank constraint with an additional consideration that the analysis operator $A$ doesn’t have trivially linear dependent rows, i.e., $A_{i}\neq A_{j}$ for $i\neq j$.

In , Hawe et al. exploited the above constraints - full rank matrices with normalized rows, and a non-convex sparsity measurement function called the mixed $(p,q)$-pseudo-norm to minimize problem (I.3). They employed a conjugate gradient method on manifolds to solve this optimization problem. Their experimental results for classical image restoration problems show competitive performance compared to state-of-the-art techniques.

Ophir et al. proposed a simple analysis operator learning algorithm, where analysis “atoms” are learned sequentially by identifying directions that are orthogonal to a subset of the training data.

Apart from the above analysis operator learning algorithms, Peyré and Fadili proposed an attractive learning approach in . They considered the analysis operator from a particular viewpoint. They interpreted the behavior of the analysis operator as a convolution with some finite impulse response filters. Keeping this idea in mind, they formulated the analysis operator learning as a bi-level programming problem which was solved using a gradient descent algorithm. However, their work only considered a simple case - one filter and 1D signals. Following this direction, a preliminary attempt to apply this idea to 2D image processing was done in .

I-C Motivation and contributions

Among the existing algorithms for analysis operator learning, only few prior works have been evaluated based on natural images . Moreover, most of these algorithms have to impose a non-convex constraint on the analysis operator $A$, making the corresponding optimization problems hard to solve. Thus a question arises: Is it possible to introduce a more principled technique to learn optimized analysis operators without the need to impose additional constraints on the operators?

In this paper, we give an answer to this question. First, we extend the patch-based analysis model to a global image regularization term, which allows to consider also more general inverse problems such as image deconvolution and image inpainting. Then, we show that this model is equivalent to higher-order filter-based MRF models such as the FoE model . Motivated by this observation, we apply a loss-function based training scheme and show that this approach excludes the trivial solution of the analysis operator learning problem without imposing any additional constraints. Furthermore, we carefully investigate the effect of different aspects of the analysis based model. We show that the choice of the sparsity promoting function is the most important aspect. We present various experimental results for standard image restoration problems to demonstrate the effectiveness of our training model. Numerical results show that our trained model significantly outperforms existing analysis operator based models and is on par with specialized image denoising algorithms while being computationally very efficient. Therefore, our training procedure provides an attractive alternative to existing approaches.

A shorter version of this paper was presented in GCPR .

I-D Notation

II Insights into analysis based models

In this section, we first show the equivalence between the patch-based analysis model and filter-based probabilistic image patch modeling - Product of Experts (PoE) . Then we extend the patch-based analysis model to the image-based model and show connections to higher order MRFs .

The patch-based analysis model in (I.2) focuses on modeling small image patches, which is formulated as a matrix-vector multiplication ($Ax$). This procedure can be interpreted as projecting a signal $x$ (an image patch) onto a set of linear components $\{A_{i}\}_{i=1}^{n}$, where each component $A_{i}$ is a row of the matrix $A$. Note that projecting an image patch onto a linear component ($A_{i}x$) is equivalent to filtering the patch with a linear filter given by $A_{i}$.

The PoE model provides a prior distribution on small image patches by taking the product of several expert distributions, where each expert works on a linear filter and the expert function. The PoE model is formally written as $p(x)=\frac{1}{Z(\Theta)}\text{exp}(-E_{PoE}(x,\Theta))$ with

where $\rho_{i}$ is the potential function, $Z(\Theta)$ is the normalization and $\Theta$ are the parameters of this model.

Comparing the analysis prior given in (I.2), ${\phi(Ax)=\sum\nolimits_{i=1}^{n}\phi_{i}(A_{i}x)}$ with the above PoE model, we can see they are actually the same model if we choose the penalty function as $\phi_{i}=-\text{log}\rho_{i}$. In this case, if we consider the analysis operator learning problem based on the strategy which focuses on the modeling of small image patches rather than defining a prior model over an entire image, the learning problem is tantamount to learning filters in the PoE model.

II-B From patch-based to image-based model

Patch-based models are only valid for the reconstruction of a single patch. When they are applied to full image recovery, a common strategy is patch averaging . All the patches in the entire image are treated independently, reconstructed individually and then integrated to form the final reconstruction result by averaging the overlapping regions. While this method is simple and intuitive, it clearly ignores the coherence between over-lapping patches, and thus misses global support during image reconstruction. To overcome these drawbacks, an extension to the whole image is necessary where patches are not treated independently but each of them is a part of the image.

A promising direction to formulate an image-based model is to make use of the formalism of higher-order MRFs which enforce coherence across patches . The basic idea is to modify the patch-based analysis model in (I.2) such that all possible patches in the entire image and the corresponding coefficient vectors $Ax$ are considered at once. This leads to an image-based prior model of the form:

A key characteristic of the model (II.2) is that it explicitly models the overlapping of image patches, which are highly correlated. Intuitively, it is a better strategy for image modeling compared to the patch averaging approach. In Subsection V-A we will provide experimental results to support this claim.

II-C Equivalence between the image-based analysis model and the FoE model

If we consider in (II.2) each row of $A$ ($A_{i}$) as a 2-D filter ($\sqrt{m}\times\sqrt{m}$), we can rewrite this term as

where $(A_{i}*u_{p})$ denotes the result of convolving the patch at pixel $p$ with filter $A_{i}$. After having a closer look at this prior term, interestingly we find that it is the same as the FoE model proposed by Roth and Black . The FoE models the prior probability of an image by using a set of linear filters and a potential function. The probability density function for the entire image is written as $p(u)=\frac{1}{Z(\Theta)}\text{exp}(-E_{FoE}(u,\Theta))$ with

where $\rho_{i}$ is the potential function, $Z(\Theta)$ is the normalization and $\Theta$ is a vector holding the parameters of this model. Based on the observation that responses of linear filters applied to natural images typically exhibit heavy tailed distribution, two types of heavy tailed potential functions, the Student-t distribution (ST) and generalized Laplace distribution (GLP), are commonly considered:

Comparing (II.3) and (II.4), we can see that they are exactly the same if we choose the penalty function $\phi_{i}=-\text{log}\rho_{i}$,

Note that these choices lead to commonly used non-convex sparsity promoting functions .

In conclusion, the FoE model can be seen as an extension of the co-sparse analysis model from a patch-based formulation to an image-based formulation. It comes along with the advantage of inherently capturing the coherence between overlapping patches which has to be enforced explicitly in patch-based models. As we will see in the next section, the image-based model also allows to learn optimized analysis operators without the need for additional constraints.

III Learning

In this section, we first present a loss-based training procedure based on bi-level optimization. Our algorithm is closely related to the algorithm proposed in but we propose to solve the lower level problems with high accuracy which leads to improved gradient directions for minimizing the loss function with respect to the model parameters. As a result of this seemingly minor modification, we achieve significantly better results compared to previous work. For more details about the refined training algorithm we refer to .

Existing approaches to learn the parameters in the FoE model fall into two main types: (1) probabilistic learning using sampling-based algorithms, e.g., ; (2) bi-level training based on MAP estimation, e.g., . Reviewing all algorithms is beyond the scope of this paper. Here we focus on the the bi-level training scheme and refer the interested reader to for a survey.

Bi-level optimization is a popular and effective technique for selecting hyper parameters, cf. . The problem of learning the analysis operator can be written as the following bi-level optimization problem

Given a noisy observation $f$ and the ground truth $g$, our goal is to find optimal hyper parameters $\vartheta$ such that the minimizer $u^{*}(\vartheta)$ of the lower level problem $E(u,f,\vartheta)$ minimize the higher level problem $L(u^{*}(\vartheta),g)$. In our case, the hyper parameters will be used to parametrize the analysis operators as well as the potential functions, the lower level problem is given by the energy function and the higher level problem is given by a certain loss function that compares the solution of the lower level problem with the ground truth solution.

where $\phi(\mathcal{A}_{i}u)=\sum\nolimits_{p=1}^{N_{p}}\phi((\mathcal{A}_{i}u)_{p})$, $\mathcal{A}_{i}$ is an $N_{p}\times N_{p}$ highly sparse matrix, which makes the convolution of the filter $A_{i}$ with a two-dimensional image $u$ equivalent to the product of the matrix $\mathcal{A}_{i}$ with the vectorization of $u$, i.e., $A_{i}*u\Leftrightarrow\mathcal{A}_{i}u$, and $\alpha_{i}\geq 0$ is the weight parameter associated to the filter $\mathcal{A}_{i}$. In our training model, we express the filter $\mathcal{A}_{i}$ as a linear combination of a set of basis filters $\{B_{1},\cdots,B_{N_{B}}\}$, i.e.,

The loss function is defined to penalize the difference (loss) between the optimal solution of the energy minimization problem and the ground-truth. In this paper, we make use of the following differentiable function as in :

where $g$ is the ground-truth image and $u^{*}$ is the minimizer of energy function (III-A). This loss function has an interpretation of pursuing as high PSNR as possible.

Given the training samples $\{f_{s},g_{s}\}_{s=1}^{S}$, where $g_{s}$ and $f_{s}$ are the $s^{th}$ clean image and the associated noisy version respectively, our aim is to learn an optimal analysis operator or a set of filters which are defined by parameters $\vartheta=(\alpha,\beta)$ (we group the coefficients $\beta_{ij}$ and weights $\alpha_{i}$ into a single vector $\vartheta$), such that the overall loss function for all samples is as small as possible. Therefore, our learning model is formulated as the following bi-level optimization problem:

We eliminate $\lambda$ for simplicity since it can be incorporated into the weights $\alpha$. Our analysis operator training model has two advantages over existing analysis operator learning algorithms.

It is completely unconstrained with respect to the analysis operator $A$. Normally, existing approaches such as have to impose some non-convex constraints over the analysis operator. On the one hand, this makes the corresponding optimization problem difficult to solve, and on the other hand it decreases the probability of learning a meaningful analysis operator, because as indicated in , there is no evidence to prove that the introduced constraints are the most suitable choices. The reason why constraints are indispensable for these approaches lies in the need to exclude the trivial solution $A=0$. However, looking back at our training model, this trivial solution can be avoided naturally. If $A=0$, the optimal solution of the lower-level problem in (III.5) is certainly $u_{s}^{*}=f_{s}$, which makes the loss function still large; thus this trivial solution is not acceptable in that the goal of our model is to minimize the loss function. Therefore, the optimal operator $A$ must comprise some meaningful filters such that the minimizer of the lower-level problem is close to the ground-truth.

The learned analysis operator inherently captures the properties of overlapping patches. In , their approaches present a patch-based prior, and thus for global reconstruction of an entire image, the common strategy consists of two stages: (i) extract overlapping patches, reconstruct them individually by synthesis-prior or analysis-prior based model, and (ii) form the entire image by averaging the final reconstruction results in the overlapping regions. This strategy clearly misses global support during the reconstruction process; however our approach can overcome these drawbacks.

In the work of , the authors employ the patch-based model to train the analysis operator, but use it in the manner of an image based model. Clearly, if the final intent is to use the analysis operator in an image-based model, a better strategy is to train it also in the same framework.

III-B Solving the bi-level problem

In this subsection, we consider the bi-level optimization problem from a general point of view. For convenience, we only consider the case of a single training sample and we show how to extend the framework to multiple training samples in the end.

According to the optimality condition, the solution of the lower-level problem in (III.5) is given by $u^{*}$, such that $\nabla_{u}E(u^{*})=0$. Therefore, we can rewrite problem (III.5) as following constrained optimization problem

In principle, we can continue to calculate the second derivatives of (III.7), i.e., the Jacobian matrix of G, with which we can then employ a Newton’s method to solve the necessary optimality system (III.8) as in . However, for this problem, calculating the Jacobian of G is computationally expensive; thus in this paper we do not consider the second derivatives and only make use of the first derivatives. An efficient Newton’s method is subject to the future work.

In our training model, what we are interested in is the parameters $\vartheta=\{\alpha,\beta\}$. We can reduce unnecessary variables in (III.8) by solving for $p$ and $u$ in (III.8), and substituting them into the second and the third equation. We arrive at the following gradients of the loss function with respect to the parameters $\vartheta$:

$H_{E}(u)$ denotes the Hessian matrix of $E(u)$,

Note that in (III.9) we also eliminated the Lagrange multiplier $\mu$ associated to the inequality constraint $\alpha\geq 0$ since we utilize a quasi-Newton’s method for optimization, which can easily handle this type of box constraints; therefore we do not need to consider the inequality constraint in the derivatives. The derivatives in (III.9) are equivalent to the results presented in , which used implicit differentiation for the derivation.

Considering the case of $S$ training samples, in fact it turns out that the derivatives of the overall loss function in (III.5) with respect to the parameters $\vartheta$ are just the sum of gradients given in (III.9) with respect to all the training samples.

III-C Bi-level learning algorithm

In (III.9), we have collected all the necessary information to compute the gradients of the loss function with respect to the parameters $\vartheta$, so we can now employ gradient descent based algorithms, e.g., the steepest descent method, for optimization. Although this type of algorithm is very easy to implement, it is not efficient. In this paper, we turn to a more efficient non-linear optimization method - the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) quasi-Newton’s method . We summarize our bi-level learning scheme in Algorithm 1.

In our work, step (ii) in Algorithm 1 is completed using the L-BFGS algorithm, as this problem is smooth, to which L-BFGS is perfectly applicable. We solve this minimization problem to a very high accuracy with $\|\nabla_{u}E(u^{*})\|_{2}\leq 10^{-3}$ (gray-values in the range ), i.e., we use a more conservative convergence criterion in this inner loop than previous work . The training algorithm is terminated when the relative change of the loss is less than a tolerance, e.g., $tol=10^{-5}$, a maximum number of iterations e.g., $maxiter=500$ is reached or L-BFGS can not find a feasible step to decrease the loss.

IV Training experiments

We conducted our training experiments using the training images from the BSDS300 image segmentation database . We used the whole 200 training images, and randomly sampled one $64\times 64$ patch from each training image, resulting in a total of 200 training samples. We then generated the noisy versions by adding Gaussian noise with standard deviation ${\sigma=25}$. Figure 1 shows an exemplary subset of the training data together with the noisy version.

In order to evaluate the performance of the learned analysis operators, we applied them to the image denoising experiment over a validation dataset consisting of 68 images from Berkeley database . This is a common denoising test dataset for natural images, which was selected by Roth and Black . The performance of an image denoising algorithm varies greatly for different image contents. We therefore consider the average performance over the whole test dataset as performance measure.

IV-B Training experiments

For the preliminary experiment, we initialized the analysis operator using 48 random filters having unified norms and weights, which are 0.01 and 1, respectively. Finally training result shows that all the coefficients with respect to the first atom of DCT-7 (an atom with constant entries) are approximately equal to zero, implying that the first atom isn’t necessary to construct the filters. Therefore the learned filters are undoubtedly zero-mean because all the remaining atoms are zero-mean; this makes the analysis prior based model (III-A) illumination invariant. This result is coherent with the findings in the work that meaningful filters should be zero-mean. Then we explicitly exclude the first atom in DCT-7 to speed up the training process for the following experiments.

We then conducted training experiments based on three different penalty functions. In this paper, the regularization parameter $\varepsilon$ in (IV.1) was set to $\varepsilon=10^{-2}$. Smaller $\varepsilon$ implies a better fitting to the absolute function, but it makes the lower-level problem harder to solve and the training algorithm fail. Just like in the preliminary experiment, we also learned 48 filters. We initialized the filters using the modified DCT-7 basis with unified norms and weights. The optimal analysis operator learned by using the penalty function $\text{log}(1+z^{2})$ is shown in Figure 2. The final loss function values (normalized by the number of training images) of these three experiments are presented in Table I (first three columns), together with the average denoising PSNR results based on 68 test images with $\sigma=25$ Gaussian noise.

As shown in Figure 2, the learned filters present some special structures. We can find high-frequency filters as well as derivative filters including the first derivatives along different directions, the second and the third derivatives. These filters make the analysis prior based model (III-A) a higher-order model which is able to capture the structures in natural images that cannot be captured by using only the first derivatives as in the total variation based methods.

In our training model, the size and the number of filters are free parameters; thus we can train filters of various sizes and numbers. Our current implementation is unoptimized Matlab code. The training time for 48 filters of size $7\times 7$ was approximately 24 hours on a server (Intel X5675, 3.07GHz), 98 filters of size $7\times 7$ took about 80 hours. However, the training time for larger filter size $9\times 9$ was much longer; it took about 20 days. Fortunately, the training procedure is off-line; thus the training time does not matter too much in practice. Demo training code can be downloaded from our homepage www.gpu4vision.org.

IV-C The influence of the penalty function

In order to further investigate how important the non-convex penalty function is for the analysis prior based model, we considered an analysis model consisting of 48 fixed and predefined filters (DCT-7 filters excluding the filter with uniform entries) and making use of the $\text{log}(1+z^{2})$ penalty function. We only optimized the norm and weight of each filter by using our bi-level training algorithm. The training loss value and the denoising test result of this model are shown in Table I (the sixth column entitled “direct DCT-7”). The image denoising test result is surprisingly good, even though this analysis model only utilizes a predefined analysis operator DCT-7. We will see in Table II of Section V that the performance of this model is already on par with the currently best analysis operator learning model - GOAL , which involves much more carefully trained filters. The success of this model lies in the non-convex penalty function $\text{log}(1+z^{2})$.

IV-D The influence of the number of filters

From Table I, one can see that the improvement achieved by over-complete analysis operator is marginal. Therefore for the analysis model, under-complete operators already work sufficiently well. An increase of the number of filters can not bring large improvements.

IV-E The influence of filter size

Intuitively the size of filters should be an important factor for the analysis model. In order to investigate the influence of filter size, we conducted training experiments for several different analysis models, where the filter size varies from $3\times 3$ to $9\times 9$. The training and evaluation results of these models are presented in Table I and Figure 4.

One can see that increasing the filter size yields some improvements. However, the performance is close to saturation when the filter size is increasing to $7\times 7$. The improvement brought by increasing the filter size to $9\times 9$ is negligible. This implies that we can not expect large improvements by increasing the filter size to $11\times 11$ or larger.

IV-F The robustness of our training scheme

As our training model (III.5) is a non-convex optimization problem, we can only find local minima. Thus a natural question about the initialization arises. We did have experiments for different initializations, such as random initialization. The final learned analysis operators are surely different, but all of them have almost the same training loss, which is the goal of our optimization problem. In addition, these operators perform similarly in evaluation experiments.

Another issue about the robustness of our training scheme is the influence of the training dataset. Since the training patches were randomly selected, we could run the training experiment multiple times by using a different training dataset. Finally, we found that the deviation of test PSNR values based on 68 test images is within 0.02dB, which is negligible.

V Application results using learned operators

An important question for a learned prior model is how well it generalizes. To evaluate this, we directly applied our learned analysis operators, which where trained for the image denoising task, to various image restoration problems such as image deconvolution, inpainting and super-resolution, as well as denoising. To start with, we first express the image restoration model by using our learned analysis operator, which is formulated as:

where $K$ is a linear operator which depends on the respective application to be handled.

where $A_{i}*u$ denotes the convolution of image $u$ with filter $A_{i}$, and $A_{i}^{-}$ denotes the filter obtained by mirroring $A_{i}$ around its center pixel (in practice, we need to carefully handle the boundaries as we need to ensure that the filters $A_{i}^{-}$ correspond to $\mathcal{A}_{i}^{T}$).

We provide Matlab demo code for training and denoising with penalty function $\text{log}(1+z^{2})$ on our homepage www.gpu4vision.org.

We first apply the analysis model based on our learned operators to the image denoising problem. In the case of image denoising, $K$ is simply the identity matrix, i.e., $K=\mathcal{I}$. Since the image denoising performance of one method varies greatly for different image contents, in order to make a fair comparison, we conducted denoising experiments over a standard test dataset - 68 Berkeley test images identified by Roth and Black . We used exactly the same noisy version of each test image for different methods and different test images were added with distinct noise realizations. All results were computed per image and then averaged over the test dataset.

We considered image denoising for various noise levels $\sigma=\{15,25,50\}$. For noise levels other than $\sigma=25$, we need to tune the parameter $\lambda$ in (V.1). An empirical choice of $\lambda$ is: $\sigma=15$, $\lambda=25/\sigma\times 1.15$; $\sigma=50$, $\lambda=25/\sigma\times 0.8$.

Table II shows the summary of denoising results achieved by different penalty functions. One can clearly see that two non-convex penalty functions lead to similarly good results and they significantly outperform the results of the convex function $|z|$. In addition, we can also see that the over-complete operator can not improve the performance too much and larger filters ($9\times 9$) can only achieve slightly better performance. Both of these two models are more time consuming than the model with 48 filters for inference; therefore, the analysis model based on 48 filters of size $7\times 7$ offers the best trade-off between computational cost and performance. In the following experiments, we only consider the model of 48 filters and the penalty function $\text{log}(1+z^{2})$. We prefer the penalty function $\text{log}(1+z^{2})$, since it is completely smooth, making the corresponding minimization problem easier to solve. We present two denoising examples obtained by three different penalty functions in Figure 5.

V-A2 Comparison to other analysis models

An interesting result in Table II is that the performance of the direct DCT-7 model, which only utilizes a predefined analysis operator DCT-7 (48 filters of size $7\times 7$), is already on par with the GOAL model, which involves much more carefully trained filters (98 filters of size $7\times 7$). For this direct DCT-7 model, we used the $\text{log}(1+z^{2})$ penalty function, and only optimized the norms and weights of the filters using our bi-level training algorithm. This result demonstrates the importance of non-convex penalty functions and the effectiveness of our bi-level training scheme.

As the analysis operator of GOAL model is trained using a patch-based model, we can also use it in the manner of patch-averaging to conduct image denoising like K-SVD . We embedded the learned analysis operator $\Omega$ into the patch-based analysis model (I.2), and used it to denoise each patch extracted from an image. We also considered overlapped windows and averaged the results in the overlapping regions to form the final denoised image. As expected, we got inferior results (average PSNR 28.25 over 68 test images) to the model formulated under the FoE framework (average PSNR 28.45).

V-A3 Comparison to state-of-the-art methods

In order to evaluate how well our analysis models work for the denoising task, we compared their performance with leading image denoising methods, including three state-of-the-art methods: (1) BM3D ; (2) LSSC ; (3) GMM-EPLL along with three leading generic methods: (4) a MRF-based approach, FoE ; (5) a synthesis sparse representation based method, KSVD trained on natural image patches; and (6) the currently published best analysis operator learning method, GOAL . All implementations were downloaded from the corresponding authors’ homepages. We conducted denoising experiments over 68 Berkeley test images with various noise levels $\sigma=\{15,25,50\}$. All results were computed per image and then averaged over the number of images.

Table II shows the summary of results. One can see that our trained model based on the penalty function $\text{log}(1+z^{2})$ (48 learned filters, $7\times 7$) outperforms three leading generic methods and is on par with three state-of-the-art methods for any noise level. To the best of our knowledge, this is the first time that a MRF model based on generic priors of natural images has achieved such clear state-of-the-art performance. Figure 6 gives a detailed comparison between our learned analysis model and three state-of-the-art methods over 68 test images for $\sigma=25$. We can see that all the points surround the diagonal line “$y=x$” closely, i.e., all considered methods achieve very similar results. Therefore, it is clear that our learned analysis models based on non-convex penalty functions are state-of-the-art. We present an image denoising example of the considered methods in Figure 7.

Our model is well-suited to GPU parallel computation since it solely contains convolution of filters with an image. Our GPU implementation based on a NVIDIA Geforce GTX 580 accelerates the inference procedure significantly; for a denoising task with $\sigma=25$, typically it takes 0.87s for image size $512\times 512$, 0.60s for $481\times 321$ and 0.29s for $256\times 256$, i.e., using our GPU based implementation, image denoising can be conducted in near real-time at 3.4fps for an $256\times 256$ image sequence, with state-of-the-art performance. In Table IV, we show the average running time of the considered denoising methods on $481\times 321$ images.

Considering the speed and quality of our model, it is a perfect choice as a base method in the image restoration framework proposed in , which leverages advantages of existing methods.

V-B Single image super-resolution

For single image super-resolution, the linear operator $K$ is constructed by a decimation operator $\Phi$ and a blurring operator $B$, i.e., $K=\Phi B$. In order to perform a better comparison with the latest analysis model GOAL , we conducted the same single image super-resolution experiment. We artificially created a low resolution image by downsampling a ground-truth image by a factor of 3 using bicubic interpolation. Then the low resolution image was corrupted by Gaussian noise with $\sigma=8$. We magnified the noisy low resolution image by the same factor using (a) bicubic interpolation, (b) GOAL method , (c) our learned analysis model based on penalty function log$(1+z^{2})$ respectively. Figure 8 shows the results for different methods. One can see that two analysis models present similar results, which are visually and quantitatively better than the bicubic method.

V-C Non-blind image deconvolution

V-D Image inpainting

Image inpainting is the process of filling in lost image data such that the resulting image is visually appealing. Typically, the positions of the pixels to be filled up are given. In our formulation, the linear operator $K$ is simply a sampling matrix, where each row contains exactly one entry equal to one. Its position indicates a pixel with given value. The parameter $\lambda$ corresponds to joint inpainting and denoising, and the choice $\lambda\rightarrow+\infty$ means pure inpainting. In our experiment since we assumed the test images are noise free, we empirically selected $\lambda=10^{3}$. Due to space limitation, we only considered a classical image inpainting task here.

We destroyed the ground-truth “ Lena” image ($512\times 512$) artificially by masking 90% of the entire pixels randomly as shown in Figure 10(a). Then we reconstructed the incomplete image using our learned analysis model - log($1+z^{2}$)-based model. In order to present a comparison, we also give the inpainting result of the GOAL model and FoE model . From Figure 10, one can see that the result of our learned analysis model based on log($1+z^{2}$) penalty achieves equivalent results with respect to the GOAL model.

VI Conclusion and outlook

In this paper, we have expressed our insights into the co-sparse analysis model. We propose to go beyond existing patch-based models, and to exploit the framework of FoE model to define a image prior over the entire image, rather than image patches. We have pointed out that the image based analysis model is equivalent to the FoE model. Starting from this conclusion, we have introduced a bi-level training approach for analysis operator learning, which is solved effectively with L-BFGS algorithm. By using our training framework, we have carefully investigated the effect of different aspects of the analysis prior model including the filter size, the number of filters and the penalty function.

Since our training scheme directly optimizes the MAP-inference based analysis model, the learned model is an optimal MAP inference for image restoration problems. Numerical results have confirmed the good performance of our learned analysis model. For classic tasks such as image denoising, image deconvolution, image super-resolution and image inpainting, our learned analysis model has achieved strongly competitive performance with current state-of-the-art methods, and clearly outperforms existing analysis learning approaches.

For future work, focusing on generic priors of natural images, we expect that our learned analysis model could be improved potentially in two aspects: (1) consider more flexible penalty function. In our current model, the penalty function is fixed to the same form for every filter. If we free the shape of the penalty function, our model will possess more freedom, which might increase the performance. A feasible way to consider alterable penalty function is to make use of the GSMs prior . (2) make use of larger training dataset. Our training is conducted based on 200 training samples, which is only a very small part of the natural images. Consequently, the learned filters may over-fit on the training dataset. However, our current training scheme is not available for large training dataset, e.g., $\sim 10^{6}$, because it needs to solve the lower-level problem for each training sample. Feasible methods may include making use of stochastic optimization.

VII Acknowledgments

The authors wish to thank Qi Gao (TU Darmstadt) for sharing details of the FoE model; Daniel Zoran (Hebrew University of Jerusalem) for sharing testing details for the GMM-EPLL algorithm; Simon Hawe (TU Munich) for useful discussion about analysis operator learning.

References