Analysis Operator Learning and Its Application to Image Reconstruction

I Introduction

I-B Synthesis Model and Dictionary Learning

For $d>n$, the dictionary is said to be overcomplete or redundant.

Now, using the knowledge that (2) allows a sparse solution, an estimation of the original signal in (1) can be obtained from the measurements $\mathbf{y}$ by first solving

with $0<p\leq 1$ and differentiable approximations of (4). As the signal is synthesized from the sparse coefficients, the reconstruction model (3) is called the synthesis reconstruction model .

To find the minimizer of Problem (3), various algorithms based on convex or non-convex optimization, greedy pursuit methods, or Bayesian frameworks exist that may employ different choices of $g$. For a broad overview of such algorithms, we refer the interested reader to . What all these algorithms have in common, is that their performance regarding the reconstruction quality severely depends on an appropriately chosen dictionary $\mathcal{D}$. Ideally, one is seeking for a dictionary where $\mathbf{s}$ can be represented most accurately with a coefficient vector $\mathbf{x}$ that is as sparse as possible. Basically, dictionaries can be assigned to two major classes: analytic dictionaries and learned dictionaries.

Analytic dictionaries are built on mathematical models of a general type of signal, e.g. natural images, they should represent. Popular examples include Wavelets , Bandlets, and Curvlets among several others, or a concatenation of various such bases/dictionaries. They offer the advantages of low computational complexity and of being universally applicable to a wide set of signals. However, this universality comes at the cost of not giving the optimally sparse representation for more specific classes of signals, e.g. face images.

It is now well known that signals belonging to a specific class can be represented with fewer coefficients over a dictionary that has been learned using a representative training set, than over analytic dictionaries. This is desirable for various image reconstruction applications as it readily improves their performance and accuracy . Basically, the goal is to find a dictionary over which a training set admits a maximally sparse representation. In contrast to analytic dictionaries, which can be applied globally to an entire image, learned dictionaries are small dense matrices that have to be applied locally to small image patches. Hence, the training set consists of small patches extracted from some example images. This restriction to patches mainly arises from limited memory, and limited computational resources.

Roughly speaking, starting from some initial dictionary the learning algorithms iteratively update the atoms of the dictionary, such that the sparsity of the training set is increased. This procedure is often performed via block-coordinate relaxation, which alternates between finding the sparsest representation of the training set while fixing the atoms, and optimizing the atoms that most accurately reproduce the training set using the previously determined sparse representation. Three conceptually different approaches for learning a dictionary became well established, which are probabilistic ones like , clustering based ones such as K-SVD , and recent approaches which aim at learning dictionaries with specific matrix structures that allow fast computations like . For a comprehensive overview of dictionary learning techniques see .

I-C Analysis Model

An alternative to the synthesis model (3) for reconstructing a signal, is to solve

Different analysis operators proposed in the literature include the fused Lasso , the translation invariant wavelet transform , and probably best known the finite difference operator which is closely related to the total-variation . They all have shown very good performance when used within the analysis model for solving diverse inverse problems in imaging. The question is: Can the performance of analysis based signal reconstruction be improved when a learned analysis operator is applied instead of a predefined one, as it is the case for the synthesis model where learned dictionaries outperform analytic dictionaries? In , it has been discussed that the two models differ significantly, and the naïve way of learning a dictionary and simply employing its transposed or its pseudo-inverse as the learned analysis operator fails. Hence, different algorithms are required for analysis operator learning.

I-D Contributions

In this work, we introduce a new algorithm based on geometric optimization for learning a patch based analysis operator from a set of training samples, which we name GOAL (GeOmetric Analysis operator Learning). The method relies on a minimization problem, which is carefully motivated in Section II-B. Therein, we also discuss the question of what is a suitable analysis operator for image reconstruction, and how to antagonize overfitting the operator to a subset of the training samples. An efficient geometric conjugate gradient method on the so-called oblique manifold is proposed in Section III for learning the analysis operator. Furthermore, in Section IV we explain how to apply the local patch based analysis operator to achieve global reconstruction results. Section V sheds some light on the influence of the parameters required by GOAL and how to select them, and compares our method to other analysis operator learning techniques. The quality of the operator learned by GOAL on natural image patches is further investigated in terms of image denoising, inpainting, and single image super-resolution. The numerical results show the broad and effective applicability of our general approach.

I-E Notations

Matrices are written as capital calligraphic letters like $\mathcal{X}$, column vectors are denoted by boldfaced small letters e.g. $\mathbf{x}$, whereas scalars are either capital or small letters like $n,N$. By $v_{i}$ we denote the $i^{\textit{th}}$ element of the vector $\mathbf{v}$, $v_{ij}$ denotes the $i^{\textit{th}}$ element in the $j^{\textit{th}}$ column of a matrix $\mathcal{V}$. The vector $\mathbf{v}_{:,i}$ denotes the $i^{\textit{th}}$ column of $\mathcal{V}$ whereas $\mathbf{v}_{i,:}$ denotes the transposed of the $i^{\textit{th}}$ row of $\mathcal{V}$. By $\mathcal{E}_{ij}$, we denote a matrix whose $i^{\textit{th}}$ entry in the $j^{\textit{th}}$ column is equal to one, and all others are zero. $\mathcal{I}_{k}$ denotes the identity matrix of dimension $(k\times k)$, ${\bm{0}}$ denotes the zero-matrix of appropriate dimension, and $\operatorname{ddiag}(\mathcal{V})$ is the diagonal matrix whose entries on the diagonal are those of $\mathcal{V}$. By $\|\mathcal{V}\|_{F}^{2}=\sum_{i,j}v_{ij}^{2}$ we denote the squared Frobenius norm of a matrix $\mathcal{V}$, $\operatorname{tr}(\mathcal{V})$ is the trace of $\mathcal{V}$, and $\operatorname{rk}(\mathcal{V})$ denotes the rank.

II Analysis Operator Learning

The topic of analysis operator learning has only recently started to be investigated, and only few prior work exists. In the sequel, we shortly review analysis operator learning methods that are applicable for image processing tasks.

An adaption of the widely known K-SVD dictionary learning algorithm to the problem of analysis operator learning is presented in . As in the original K-SVD algorithm, $G(\mathbf{\Omega}\mathcal{S})=\sum_{i}\|\mathbf{\Omega}\mathbf{s}_{i}\|_{0}$ is employed as the sparsifying function and the target cosparsity is required as an input to the algorithm. The arising optimization problem is solved by alternating between a sparse coding stage over each training sample while fixing $\mathbf{\Omega}$ using an ordinary analysis pursuit method, and updating the analysis operator using the optimized training set. Then, each row of $\mathbf{\Omega}$ is updated in a similar way as described in the previous paragraph for the method of . Interestingly, the operator learned on piecewise constant image patches by and closely mimics the finite difference operator.

In , the authors use $G(\mathbf{\Omega}\mathcal{S})=\sum_{i}\|\mathbf{\Omega}\mathbf{s}_{i}\|_{1}$ as the sparsity promoting function and suggest a constrained optimization technique that utilizes a projected subgradient method for iteratively solving (6). To exclude the trivial solution, the set of possible analysis operators is restricted to the set of Uniform Normalized Tight Frames, i.e. matrices with uniform row norm and orthonormal columns. The authors state that this algorithm has the limitation of requiring noiseless training samples whose analyzed vectors $\{\mathbf{\Omega}\mathbf{s}_{i}\}_{i=1}^{M}$ are exactly cosparse.

To overcome this restriction, the same authors propose an extension of this algorithm that simultaneously learns the analysis operator and denoises the training samples, cf. . This is achieved by alternating between updating the analysis operator via the projected subgradient algorithm and denoising the samples using an Augmented Lagrangian method. Therein, the authors state that their results for image denoising using the learned operator are only slightly worse compared to employing the commonly used finite difference operator.

An interesting idea related to the analysis model, called Fields-of-Experts (FoE) has been proposed in . The method relies on learning high-order Markov Random Field image priors with potential functions extending over large pixel neighborhoods, i.e. overlapping image patches. Motivated by a probabilistic model, they use the student-t distribution of several linear filter responses as the potential function, where the filters, which correspond to atoms from an analysis operator point of view, have been learned from training patches. Compared to our work and the methods explained above, their learned operator used in the experiments is underdetermined, i.e. $k<n$, the algorithms only works for small patches due to computational reasons, and the atoms are learned independently, while in contrast GOAL updates the analysis operator as a whole.

II-B Motivation of Our Approach

The algorithm presented here aims at minimizing the empirical expectation of a sparsifying function $g(\mathbf{\Omega}\mathbf{s}_{i})$ for all training samples $\mathbf{s}_{i}$, while additionally keeping the empirical variance moderate. In other words, we try to avoid that the analyzed vectors of many similar training samples become very sparse and consequently prevent $\mathbf{\Omega}$ from being adapted to the remaining ones. For image processing, this is of particular interest if the training patches are chosen randomly from natural images, because there is a high probability of collecting a large subset of very similar patches, e.g. homogeneous regions, that bias the learning process.

Concretely, we want to find an $\mathbf{\Omega}$ that minimizes both the squared empirical mean

of the sparsity of the analyzed vectors. We achieve this by minimizing the sum of both, which is readily given by

Using $g(\cdot)=\|\cdot\|_{p}^{p}$, and introducing the factor $\tfrac{1}{2}$ the function we employ reads as

with $0\leq p\leq 1$ and $\mathcal{V}=\mathbf{\Omega}\mathcal{S}$.

Certainly, without additional prior assumptions on $\mathbf{\Omega}$, the useless solution $\mathbf{\Omega}=\mathbf{0}$ is the global minimizer of Problem (6). To avoid the trivial solution and for other reasons explained later in this section, we regularize the problem by imposing the following three constraints on $\mathbf{\Omega}$.

The rows of $\mathbf{\Omega}$ have unit Euclidean norm, i.e. $\|{\bm{\omega}}_{i,:}\|_{2}=1$ for $i=1,\ldots,k$.

The analysis operator $\mathbf{\Omega}$ has full rank, i.e. $\operatorname{rk}(\mathbf{\Omega})=n$.

The analysis operator $\mathbf{\Omega}$ does not have linear dependent rows, i.e. ${\bm{\omega}}_{i,:}\neq\pm{\bm{\omega}}_{j,:}$ for $i\neq j$.

These constraints motivate the consideration of the set of full rank matrices with normalized columns, which admits a manifold structure known as the oblique manifold

Due to the full rank condition on $\mathcal{X}$, the product $\mathcal{X}\mathcal{X}^{\top}$ is positive definite, consequently the strict inequality $0<\det(\tfrac{1}{k}\mathcal{X}\mathcal{X}^{\top})$ applies. To see the second inequality of Lemma 1, observe that

which implies $\operatorname{tr}(\tfrac{1}{k}\mathcal{X}\mathcal{X}^{\top})=1$. Since the trace of a matrix is equal to the sum of its eigenvalues, which are strictly positive in our case, it follows that the strict inequality $0<\lambda_{i}<1$ holds true for all eigenvalues $\lambda_{i}$ of $\tfrac{1}{k}\mathcal{X}\mathcal{X}^{\top}$. From the well known relation between the arithmetic and the geometric mean we see

Now, since the determinant of a matrix is equal to the product of its eigenvalues, and with $\sum\lambda_{i}=\operatorname{tr}(\tfrac{1}{k}\mathcal{X}\mathcal{X}^{\top})=1$, we have

Regarding Condition (iii), the following result proves useful.

Thus, Condition (iii) can be enforced via the logarithmic barrier function of the scalar products between all distinctive rows of $\mathbf{\Omega}$, i.e.

Finally, combining all the introduced constraints, our optimization problem for learning the transposed analysis operator reads as

Let $\mathbf{\Omega}$ be a minimum of $h$ in the set of transposed oblique matrices, i.e.

then the condition number of $\mathbf{\Omega}$ is equal to one.

With other words, the minima of $h$ are uniformly normalized tight frames, cf. . From Lemma 3 we can conclude that with larger $\kappa$ the condition number of $\mathbf{\Omega}$ approaches one. Now, recall the inequality

with $\sigma_{\min}$ being the smallest and $\sigma_{\max}$ being the largest singular value of $\mathbf{\Omega}$. From this it follows that an analysis operator found with a large $\kappa$, i.e. obeying $\sigma_{\min}\approx\sigma_{\max}$, carries over distinctness of different signals to their analyzed versions. The parameter $\mu$ regulates the redundancy between the rows of the analysis operator and consequently avoids redundant coefficients in the analyzed vector $\mathbf{\Omega}\mathbf{s}$.

The difference between any two entries of the analyzed vector $\mathbf{\Omega}\mathbf{s}$ is bounded by

From the Cauchy-Schwarz inequality we get

Since by definition $\|\bm{\omega}_{i,:}\|_{2}=\|\bm{\omega}_{j,:}\|_{2}=1$, it follows that $\|\bm{\omega}_{i,:}-\bm{\omega}_{j,:}\|_{2}=\sqrt{2(1-\bm{\omega}_{i,:}^{\top}\bm{\omega}_{j,:})}$.∎

The above lemma implies, that if the $i^{\textit{th}}$ entry of the analyzed vector is significantly larger than then a large absolute value of $\bm{\omega}_{i,:}^{\top}\bm{\omega}_{j,:}$ prevents the $j^{\textit{th}}$ entry to be small. To achieve large cosparsity, this is an unwanted effect that our approach avoids via the log-barrier function $r$ in (16). It is worth mentioning that the same effect is achieved by minimizing the analysis operator’s mutual coherence $\max\limits_{i\neq j}|{\bm{\omega}}_{i,:}^{\top}{\bm{\omega}}_{j,:}|$ and that our experiments suggest that enlarging $\mu$ leads to minimizing the mutual coherence.

III Analysis Operator Learning Algorithm

Knowing that the feasible set of solutions to Problem (17) is restricted to a smooth manifold allows us to formulate a geometric conjugate gradient (CG-) method to learn the analysis operator. Geometric CG-methods have been proven efficient in various applications, due to the combination of moderate computational complexity and good convergence properties, see e.g. for a CG-type method on the oblique manifold.

To make this work self contained, we start by shortly reviewing the general concepts of optimization on matrix manifolds. After that we present the concrete formulas and implementation details for our optimization problem on the oblique manifold. For an in-depth introduction on optimization on matrix manifolds, we refer the interested reader to .

The concepts presented in this subsection are visualized in Figure 2 to alleviate the understanding.

Once $\alpha^{(i)}$ has been determined, the new iterate is computed by

Now, one straightforward approach to minimize $f$ is to alternate Equations (23), (24), and (25) using $\mathcal{H}^{(i)}=-\mathcal{G}^{(i)}$, with the short hand notation $\mathcal{G}^{(i)}:=\mathcal{G}(\mathcal{X}^{(i)})$, which corresponds to the steepest descent on a Riemannian manifold. However, as in standard optimization, steepest descent only has a linear rate of convergence. Therefore, we employ a conjugate gradient method on a manifold, as it offers a superlinear rate of convergence, while still being applicable to large scale optimization problems with low computational complexity.

III-B Geometric Conjugate Gradient for Analysis Operator Learning

Regarding geodesics, note that in general a geodesic is the solution of a second order ordinary differential equation, meaning that for arbitrary manifolds, its computation as well as computing the parallel transport is not feasible. Fortunately, as the oblique manifold is a Riemannian submanifold of a product of $k$ unit spheres $S^{n-1}$, the formulas for parallel transport and the exponential mapping allow an efficient implementation.

Let $\mathbf{x}\in S^{n-1}$ be a point on a sphere and $\mathbf{h}\in T_{\mathbf{x}}S^{n-1}$ be a tangent vector at $\mathbf{x}$, then the geodesic in the direction of $\mathbf{h}$ is a great circle

The associated parallel transport of a tangent vector $\bm{\xi}\in T_{\mathbf{x}}S^{n-1}$ along the great circle $\gamma(\mathbf{x},\mathbf{h},t)$ reads as

Now, to use the geometric CG-method for learning the analysis operator, we require a differentiable cost function $f$. Since, the cost function presented in Problem (17) is not differentiable due to the non-smoothness of the $(p,q)$-pseudo-norm (10), we exchange Function (10) with a smooth approximation, which is given by

The gradient of the rank penalty term (15) is

and the gradient of the logarithmic barrier function (16) is

Combining Equations (III-B), (41), and (42), the gradient of the cost function

which is used for learning the analysis operator reads as

Regarding the CG-update parameter $\beta^{(i)}$, we employ a hybridization of the Hestenes-Stiefel Formula (29) and the Dai Yuan formula (30)

which has been suggested in . As explained therein, formula (45) combines the good numerical performance of HS with the desirable global convergence properties of DY.

Finally, to compute the step size $\alpha^{(i)}$, we use an adaption of the well-known backtracking line search to the geodesic $\Gamma(\mathcal{X}^{(i)},\mathcal{H}^{(i)},t)$. In that, an initial step size $t^{(i)}_{0}$ is iteratively decreased by a constant factor $c_{1}<1$ until the Armijo condition is met, see Algorithm 1 for the entire procedure.

In our implementation we empirically chose $c_{1}=0.9$ and $c_{2}=10^{-2}$. As an initial guess for the step size at the first CG-iteration $i=0$, we choose

as proposed in . In the subsequent iterations, the backtracking line search is initialized by the previous step size divided by the line search parameter, i.e. $t^{(i)}_{0}=\frac{\alpha^{(i-1)}}{c_{1}}$. Our complete approach for learning the analysis operator is summarized in Algorithm 2. Note, that under the conditions that the Fletcher-Reeves update formula is used and some mild conditions on the step-size selection, the convergence of Algorithm 2 to a critical point, i.e. $\liminf_{i\to\infty}\|\mathcal{G}^{(i)}\|=0$, is guaranteed by a result provided in .

IV Analysis Operator based Image Reconstruction

Remember, that the size of $\mathbf{\Omega}^{\star}$ is very small compared to the size of the image, and it has to be applied locally to small image patches rather than globally to the entire image. Artifacts that arise from naïve patch-wise reconstruction are commonly reduced by considering overlapping patches. Thereby, each patch is reconstructed individually and the entire image is formed by averaging over the overlapping regions in a final step. However, this method misses global support during the reconstruction process, hence, it leads to poor inpainting results and is not applicable for e.g. Compressive Sensing tasks. To overcome these drawbacks, we use a method related to the patch based synthesis approach from and the method used in , which provides global support from local information. Instead of optimizing over each patch individually and combining them in a final step, we optimize over the entire image demanding that a pixel is reconstructed such that the average sparsity of all patches it belongs to is minimized. When all possible patch positions are taken into account, this procedure is entirely partitioning-invariant. For legibility, we assume square patches i.e. of size $(\sqrt{n}\times\sqrt{n}$) with $\sqrt{n}$ being a positive integer.

Formally, let $\mathbf{r}\subseteq\{1,\dots,h\}$ and $\mathbf{c}\subseteq\{1,\dots,w\}$ denote sets of indices with $r_{i+1}-r_{i}=d_{v}$, $c_{i+1}-c_{i}=d_{h}$ and $1\leq d_{v},d_{h}\leq\sqrt{n}$. Therein, $d_{v},d_{h}$ determine the degree of overlap between two adjacent patches in vertical, and horizontal direction, respectively. We consider all image patches whose center is an element of the cartesian product set $\mathbf{r}\times\mathbf{c}$. Hence, with $|\cdot|$ denoting the cardinality of a set, the total number of patches being considered is equal to $|\mathbf{r}||\mathbf{c}|$. Now, let $\mathcal{P}_{rc}$ be a binary $(n\times N)$ matrix that extracts the patch centered at position $(r,c)$. With this notation, we formulate the (global) sparsity promoting function as

being the global analysis operator that expands the patch based one to the entire image. We treat image boundary effects by employing constant padding, i.e. replicating the values at the image boundaries $\lfloor\frac{\sqrt{n}}{2}\rfloor$ times, where $\lfloor\cdot\rfloor$ denotes rounding to the smaller integer. Certainly, for image processing applications $\mathbf{\Omega}^{F}$ is too large for being applied in terms of matrix vector multiplication. Fortunately, applying $\mathbf{\Omega}^{F}$ and its transposed can be implemented efficiently using sliding window techniques, and the matrix vector notation is solely used for legibility.

According to , we exploit the fact that the range of pixel intensities is limited by a lower bound $b_{l}$ and an upper bound $b_{u}$. We enforce this bounding constraint by minimizing the differentiable function $\bm{b}(\mathbf{s}):=\sum\limits_{i=1}^{N}b(s_{i})$, where $b$ is a penalty term given as

Finally, combining the two constraints (48) and (57) with the data fidelity term, the analysis based image reconstruction problem is to solve

V Evaluation and Experiments

The first part of this section aims at answering the question of what is a good analysis operator for solving image reconstruction problems and relates the quality of an analysis operator with its mutual coherence and its condition number. This, in turn allows to select the optimal weighting parameters $\kappa$ and $\mu$ for GOAL. Using this parameters, we learn one general analysis operator $\mathbf{\Omega}^{\star}$ by GOAL, and compare its image denoising performance with other analysis approaches. In the second part, we employ this $\mathbf{\Omega}^{\star}$ unaltered for solving two classical image reconstruction tasks of image inpainting and single image super-resolution, and compare our results with the currently best analysis approach FoE , and state-of-the-art methods specifically designed for each respective application.

To quantify the reconstruction quality, as usual, we use the peak signal-to-noise ratio $\textit{PSNR}=10\log(255^{2}N/\sum_{i=1}^{N}(s_{i}-s_{i}^{\star})^{2})$. Moreover, we measure the quality using the Mean Structural SIMilarity Index (MSSIM) , with the same set of parameters as originally suggested in . Compared to PSNR, the MSSIM better reflects a human observer’s visual impression of quality. It ranges between zero and one, with one meaning perfect image reconstruction.

Throughout all experiments, we fixed the size of the image patches to $(8\times 8)$, i.e. $n=64$. This is in accordance to the patch-sizes mostly used in the literature, and yields a good trade-off between reconstruction quality and numerical burden. Images are reconstructed by solving the minimization problem (58) via the conjugate gradient method proposed in . Considering the pixel intensity bounds, we used $b_{l}=0$ and $b_{u}=255$, which is the common intensity range in $8$-bit grayscale image formats. The sparsity promoting function (48) with $p=0.4$ and $\nu=10^{-6}$ is used for both learning the analysis operator by GOAL, and reconstructing the images. Our patch based reconstruction algorithm as explained in Section IV achieves the best results for the maximum possible overlap $d_{h}=d_{v}=1$. The Lagrange multiplier $\lambda$ and the measurements matrix $\mathcal{A}$ depend on the application, and are briefly discussed in the respective subsections.

V-B Analysis Operator Evaluation and Parameter Selection

For evaluating the quality of an analysis operator and for selecting appropriate parameters for GOAL, we choose image denoising as the baseline experiment. The images to be reconstructed have artificially been corrupted by additive white Gaussian noise (AWGN) of varying standard deviation $\sigma_{\textit{noise}}$. This baseline experiment is further used to compare GOAL with other analysis operator learning methods. We like to emphasize that the choice of image denoising as a baseline experiment is not crucial neither for selecting the learning parameters, nor for ranking the learning approaches. In fact, any other reconstruction task as discussed below leads to the same parameters and the same ranking of the different learning algorithms.

For image denoising, the measurement matrix $\mathcal{A}$ in Equation (58) is the identity matrix. As it is common in the denoising literature, we assume the noise level $\sigma_{\textit{noise}}$ to be known and adjust $\lambda$ accordingly. From our experiments, we found that $\lambda=\frac{\sigma_{\textit{noise}}}{16}$ is a good choice. We terminate our algorithm after $6-30$ iterations depending on the noise level, i.e. the higher the noise level is the more iterations are required. To find an optimal analysis operator, we learned several operators with varying values for $\mu,\kappa$, and $k$ and fixed all other parameters according to Subsection V-A. Then, we evaluated their performance for the baseline task, which consists of denoising the five test images, each corrupted with the five noise levels as given in Table I. As the final performance measure we use the average PSNR of the 25 achieved results. The training set consisted of $M=200\;000$ image patches, each normalized to unit Euclidean norm, that have randomly been extracted from the five training images shown in Figure 3. Certainly, these images are not considered within any of the performance evaluations. Each time, we initialized GOAL with a random matrix having normalized rows. Tests with other initializations like an overcomplete DCT did not influence the final operator.

V-C Comparison with Related Approaches

The purpose of this subsection is to rank our approach among other analysis operator learning methods, and to compare its performance with state-of-the-art denoising algorithms. Concretely, we compare the denoising performance using $\mathbf{\Omega}^{\star}$ learned by GOAL with total-variation (TV) which is the currently best known analysis operator, with the recently proposed method AOL , and with the currently best performing analysis operator FoE . Note that we used the same training set and dimensions for learning the operator by AOL as for GOAL. For FoE we used the same setup as originally suggested by the authors. Concerning the required computation time for learning an analysis operator, for this setting GOAL needs about $10$-minutes on an Intel Core i7 3.2 GHz quad-core with 8GB RAM. In contrast, AOL is approximately ten times slower, and FoE is the computationally most expensive method requiring several hours. All three methods are implemented in unoptimized Matlab code.

The achieved results for the five test images and the five noise levels are given in Table I. Our approach achieves the best results among the analysis methods both regarding PSNR, and MSSIM. For a visual assessment, Figure 6 exemplarily shows some denoising results achieved by the four analysis operators.

To judge the analysis methods’ denoising performance globally, we additionally give the results achieved by current state-of-the-art methods BM3D and K-SVD Denoising , which are specifically designed for the purpose of image denoising. In most of the cases our method performs slightly better than the K-SVD approach, especially for higher noise levels, and besides of the "barabara" image it is at most $\approx 0.5$dB worse than BM3D. This effect is due to the very special structure of the "barbara" image that rarely occurs in natural images, which are smoothed by the learned operator.

V-D Image Inpainting

In image inpainting as originally proposed in , the goal is to fill up a set of damaged or disturbing pixels such that the resulting image is visually appealing. This is necessary for the restoration of damaged photographs, for removing disturbances caused by e.g. defective hardware, or for deleting unwanted objects. Typically, the positions of the pixels to be filled up are given a priori. In our formulation, when $N-m$ pixels must be inpainted, this leads to a binary $m\times N$ dimensional measurements matrix $\mathcal{A}$, where each row contains exactly one entry equal to one. Its position corresponds to a pixel with known intensity. Hence, $\mathcal{A}$ reflects the available image information. Regarding $\lambda$, it can be used in a way that our method simultaneously inpaints missing pixels and denoises the remaining ones.

As an example for image inpainting, we disturbed some ground-truth images artificially by removing $N-m$ pixels randomly distributed over the entire image as exemplary shown in Figure 7LABEL:sub@fig:plena10. In that way, the reconstruction quality can be judged both visually and quantitatively. We assumed the data to be free of noise, and empirically selected $\lambda=10^{-2}$. In Figure 7, we show exemplary results for reconstructing the "lena" image from $10\%$ of all pixels using GOAL, FoE, and the recently proposed synthesis based method . Table II gives a comparison of further images and further number of missing pixels. It can be seen that our methods performs best independent of the configuration.

V-E Single Image Super-Resolution

For our experiments, we artificially created a low resolution image by downsampling a ground-truth image by a factor of $d$ using bicubic interpolation. Then, we employed bicubic interpolation, FoE, the method from , and GOAL to magnify this low resolution image by the same factor $d$. This upsampled version is then compared with the original image in terms of PSNR and MSSIM. In Table III, we present the results for upsampling the respective images by $d=3$. The presented results show that our method outperforms the current state-of-the-art. We want to emphasize that the blur kernel used for downsampling is different from the blur kernel used in our upsampling procedure.

Note that many single image super-resolution algorithms rely on clean noise free input data, whereas the general analysis approach as formulated in Equation (58) naturally handles noisy data, and is able to perform simultaneous upsampling and denoising. In Figure 8 we present the result for simultaneously denoising and upsampling a low resolution version of the image "august" by a factor of $d=3$, which has been corrupted by AWGN with $\sigma_{\textit{noise}}=8$. As it can be seen, our method produces the best results both visually and quantitatively, especially regarding the MSSIM. Due to high texture this image is hard to upscale even when no noise is present, see the second column of Table III. Results obtained for other images confirm this good performance of GOAL but are not presented here due to space limitation.

VI Conclusion

References