---
abstract: 'In this paper we study the topic of signal restoration using complexity regularization, quantifying the compression bit-cost of the signal estimate. While complexity-regularized restoration is an established concept, solid practical methods were suggested only for the Gaussian denoising task, leaving more complicated restoration problems without a generally constructive approach. Here we present practical methods for complexity-regularized restoration of signals, accommodating deteriorations caused by a known linear degradation operator of an arbitrary form. Our iterative procedure, obtained using the alternating direction method of multipliers (ADMM) approach, addresses the restoration task as a sequence of simpler problems involving $ \ell _2$-regularized estimations and rate-distortion optimizations (considering the squared-error criterion). Further, we replace the rate-distortion optimizations with an arbitrary standardized compression technique and thereby restore the signal by leveraging underlying models designed for compression. Additionally, we propose a shift-invariant complexity regularizer, measuring the bit-cost of all the shifted forms of the estimate, extending our method to use averaging of decompressed outputs gathered from compression of shifted signals. On the theoretical side, we present an analysis of complexity-regularized restoration of a cyclo-stationary Gaussian signal from deterioration by a linear shift-invariant operator and an additive white Gaussian noise. The theory shows that optimal complexity-regularized restoration relies on an elementary restoration filter and compression spreading reconstruction quality unevenly based on the energy distribution of the degradation filter. Nicely, these ideas are realized also in the proposed practical methods. Finally, we present experiments showing good results for image deblurring and inpainting using the JPEG2000 and HEVC compression standards.'
author:
- |
    Yehuda Dar, Michael Elad, and Alfred M. Bruckstein\
    [^1]
bibliography:
- 'IEEEabrv.bib'
- 'complexity\_regularized\_inverse\_problems\_\_refs.bib'
title: Restoration by Compression
---

[ ]{}

 \
 

Complexity regularization, rate-distortion optimization, signal restoration, image deblurring, alternating direction method of multipliers (ADMM).

Introduction
============

restoration methods are often posed as inverse problems using regularization terms. While many solutions can explain a given degraded signal, using regularization will provide signal estimates based on prior assumptions on signals. One interesting regularization type measures the complexity of the candidate solution in terms of its compression bit-cost. Indeed, encoders (that yield the bit cost) rely on signal models and allocate shorter representations to more likely signal instances. This approach of complexity-regularized restoration is an attractive meeting point of signal restoration and compression, two fundamental signal-processing problems.

Numerous works [@saito1994simultaneous; @natarajan1995filtering; @chang1997image; @mihcak1999low; @rissanen2000mdl; @chang2000adaptive; @liu2001complexity] considered the task of denoising a signal corrupted by an additive white Gaussian noise using complexity regularization. In [@natarajan1995filtering; @liu2001complexity], this idea is translated to practically estimating the clean signal by employing a standard lossy compression of its noisy version. However, more complex restoration problems (e.g., deblurring, super resolution, inpainting), involving non-trivial degradation operators, do not lend themselves to a straightforward treatment by compression techniques designed for the squared-error distortion measure. Moulin and Liu [@moulin2000statistical] studied the complexity regularization idea for general restoration problems, presenting a thorough theoretical treatment together with a limited practical demonstration of Poisson denoising based on a suitably designed compression method. Indeed, a general method for complexity-regularized restoration remained as an open question for a long while until our recent preliminary publication [@dar2016image], where we presented a generic and practical approach flexible in both the degradation model addressed and the compression technique utilized.

Our strategy for complexity-regularized signal restoration relies on the alternating direction method of multipliers (ADMM) approach [@boyd2011distributed], decomposing the difficult optimization problem into a sequence of easier tasks including $ \ell_2 $-regularized inverse problems and standard rate-distortion optimizations (with respect to a squared-error distortion metric). A main part of our methodology is to replace the rate-distortion optimization with standardized compression techniques enabling an indirect utilization of signal models used for efficient compression designs. Moreover, our method relates to various contemporary concepts in signal and image processing. The recent frameworks of Plug-and-Play Priors [@venkatakrishnan2013plug; @sreehari2016plug] and Regularization-by-Denoising [@romano2017little] suggest leveraging a Gaussian denoiser for more complicated restoration tasks, achieving impressive results (see, e.g., [@venkatakrishnan2013plug; @sreehari2016plug; @romano2017little; @dar2016postprocessing; @dar2016reducing; @rond2016poisson]). Essentially, our approach is the compression-based counterpart for denoising-based restoration concepts from [@venkatakrishnan2013plug; @sreehari2016plug; @romano2017little].

Commonly, compression methods process the given signal based on its decomposition into non-overlapping blocks, yielding block-level rate-distortion optimizations based on block bit-costs. The corresponding complexity measure sums the bit-costs of all the non-overlapping blocks, however, note that this evaluation is shift sensitive. This fact motivates us to propose a shift-invariant complexity regularizer by quantifying the bit-costs of all the overlapping blocks of the signal estimate. This improved regularizer calls for our restoration procedure to use averaging of decompressed signals obtained from compressions of shifted signals. Our shift-invariant approach conforms with the Expected Patch Log-Likelihood (EPLL) idea [@zoran2011learning], where a full-signal regularizer is formed based on a block-level prior in a way leading to averaging MAP estimates of shifted signal versions. Our extended method also recalls the cycle spinning concept, presented in [@coifman1995translation] for wavelet-based denoising. Additional resemblance is to the compression postprocessing techniques in [@nosratinia2001enhancement; @nosratinia2003postprocessing] enhancing a given decompressed image by averaging supplementary compression-decompression results of shifted versions of the given image, thus, our method generalizes this approach to any restoration problem with an appropriate consideration of the degradation operator. Very recent works [@beygi2017compressed; @beygi2017efficient] suggested the use of compression techniques for compressive sensing of signals and images, but our approach examines other perspectives and settings referring to restoration problems as will be explained below.

In this paper we extend our previous conference publication [@dar2016image] with improved algorithms and new theoretical and experimental results. In [@dar2016image] we implemented our concepts in procedures relying on the half quadratic splitting optimization technique, in contrast, here we present improved algorithms designed based on the ADMM approach. The new ADMM-based methods introduce the following benefits (with respect to using half quadratic splitting as in [@dar2016image]): significant gains in the restoration quality, reduction in the required amount of iterations, and an easier parameter setting. In addition, in this paper we provide an extensive experimental section. While in [@dar2016image] we experimentally examined only the inpainting problem, in this paper we present new results demonstrating the practical complexity-regularized restoration approach for image deblurring. While deblurring is a challenging restoration task, we present compelling results obtained using the JPEG2000 method and the image compression profile of the HEVC standard [@RefWorks:112]. An objective comparison to other deblurring techniques showed that the proposed HEVC-based implementation provides good deblurring results. Moreover, we also extend our evaluation given in [@dar2016image] for image inpainting, where here we use the JPEG2000 and HEVC compression standards in our ADMM-based approach to restore images from a severe degradation of 80% missing pixels. Interestingly, our compression-based image inpainting approach can be perceived as the dual concept of inpainting-based compression of images and videos suggested in, e.g., [@galic2008image; @schmaltz2009beating; @andris2016proof] and discussed also in [@adam2017denoising].

Another prominent contribution of this paper is the new theoretical study of the problem of complexity-regularized restoration, considering the estimation of a cyclo-stationary Gaussian signal from a degradation procedure consisting of a linear shift-invariant operator and additive white Gaussian noise. We gradually establish few equivalent optimization forms, emphasizing two main concepts for complexity-regularized restoration: the degraded signal should go through a simple inverse filtering procedure, and then should be compressed so that the decompression components will have a varying quality distribution determined by the degradation-filter energy-distribution. We explain how these ideas materialize in the practical approach we propose, thus, establishing a theoretical reasoning for the feasible complexity-regularized restoration.

This paper is organized as follows. In section \[sec:Complexity-Regularized Restoration\] we overview the settings of the complexity-regularized restoration problem. In section \[sec:Proposed Methods\] we present the proposed practical methods for complexity-regularized restoration. In section \[sec:Rate-Distortion Theoretic Analysis for the Gaussian Case\] we theoretically analyze particular problem settings where the signal is a cyclo-stationary Gaussian process. In section \[sec:Experimental Results\] we provide experimental results for image deblurring and inpainting. Section \[sec:Conclusion\] concludes this paper.

Complexity-Regularized Restoration: Problem Settings {#sec:Complexity-Regularized Restoration}
====================================================

Regularized Restoration of Signals {#subsec:Regularized-Restoration Optimization}
----------------------------------

In this paper we address the task of restoring a signal $ {\mathbf{x}}_0 \in\nolinebreak \mathbb{R}^N $ from a degraded version, $ {\mathbf{y}} \in \mathbb{R}^M $, obeying the prevalent deterioration model:

[rCl]{} \[eq:corruption model\] = \_0 +

where $ {\mathbf{H}} $ is a $M\times N$ matrix being a linear degradation operator (e.g., blur, pixel omission, decimation) and ${\mathbf{n}}\in \mathbb{R}^M$ is a white Gaussian noise vector having zero mean and variance $\sigma _n ^2$.

Maximum A-Posteriori (MAP) estimation is a widely-known statistical approach forming the restored signal, $\hat{ {\mathbf{x}}}$, via

[rCl]{} \[eq:MAP estimation\] = \_ p( [|]{} )

where $  p\left( {{{\mathbf{x}}}|{{\mathbf{y}}}} \right) $ is the posterior probability. For the above defined degradation model (\[eq:corruption model\]), incorporating additive white Gaussian noise, the MAP estimate reduces to the form of

[rCl]{} \[eq:log MAP estimate - AWGN model\] = \_ \_2\^2 - p ( )

where $p({\mathbf{x}})$ is the prior probability that, here, evaluates the probability of the candidate solution.

Another prevalent restoration approach, embodied in many contemporary techniques, forms the estimate via the optimization

[rCl]{} \[eq:restoration using a general prior\] = \_ \_2\^2 + s()

where $s({\mathbf{x}})$ is a general regularization function returning a lower value for a more likely candidate solution, and $\mu\ge 0$ is a parameter weighting the regularization effect. This strategy for restoration based on arbitrary regularizers can be interpreted as a generalization of the MAP approach in (\[eq:log MAP estimate - AWGN model\]). Specifically, comparing the formulations (\[eq:restoration using a general prior\]) and (\[eq:log MAP estimate - AWGN model\]) exhibits the regularization function $s({\mathbf{x}})$ and the parameter $ \mu $ as extensions of $( - \log p \left( {{{\mathbf{x}}}} \right)) $ and the factor $2\sigma ^2_n$, respectively.

Among the various regularization functions that can be associated with the general restoration approach in (\[eq:restoration using a general prior\]), we explore here the class of complexity regularizers measuring the required number of bits for the compressed representation of the candidate solution. The practical methods presented in this section focus on utilizing existing (independent) compression techniques, implicitly employing their underlying signal models for the restoration task.

Operational Rate-Distortion Optimization {#subsec:Operational Rate-Distortion Optimization}
----------------------------------------

The practical complexity-regularized restoration methods in this section are developed with respect to a compression technique obeying the following conceptual design. The signal is segmented to equally-sized non-overlapping blocks (each is consisted of $ N_b $ samples) that are independently compressed. The block compression procedure is modeled as a general variable-rate vector quantizer relying on the following mappings. The compression is done by the mapping $ Q : \mathbb{R}^{N_b}  \rightarrow \mathcal{W} $ from the $ N_b $-dimensional signal-block domain to a discrete set $ \mathcal{W} $ of binary compressed representations (that may have different lengths). The decompression procedure is associated with the mapping $ F :  \mathcal{W} \rightarrow \mathcal{C} $, where $ \mathcal{C} \subset \mathbb{R}^{N_b} $ is a finite discrete set (a codebook) of block reconstruction candidates. For example, consider the block ${\mathbf{x}}_{block} \in \mathbb{R}^{N_b}$ that its binary compressed representation in $ \mathcal{W} $ is given via $ b = Q\left( {\mathbf{x}}_{block} \right) $ and the corresponding reconstructed block in $ \mathcal{C} $ is $ \hat{{\mathbf{x}}}_{block} = F\left(  b \right) $. Importantly, it is assumed that shorter codewords are coupled with block reconstructions that are, in general, more likely.

The signal $ {\mathbf{x}} $ is compressed based on its segmentation into a set of blocks $\{{\mathbf{x}}_i\}_{i\in\mathcal{B}}$ (where $ \mathcal{B}$ denotes the index set of blocks in the non-overlapping partitioning of the signal). In addition we introduce the function $ r({\mathbf{z}}) $ that evaluates the bit-cost (i.e., the length of the binary codeword) for the block reconstruction $ {\mathbf{z}}\in\mathcal{C} $. Then, the operational rate-distortion optimization corresponding to the described architecture and a squared-error distortion metric is

[rCl]{} \[eq:rate-distortion optimization - blocks\] {\_i}\_[i]{} = \_[{\_i}\_[i]{} ]{} \_[i]{}[ \_2\^2]{} + \_[i]{}[r(\_i)]{} ,

where $ \lambda \ge 0 $ is a Lagrange multiplier corresponding to some total compression bit-cost. Importantly, the independent representation of non-overlapping blocks allows solving (\[eq:rate-distortion optimization - blocks\]) separately for each block [@shoham1988efficient; @ortega1998rate].

Our mathematical developments require the following algebraic tools for block handling. The matrix $ {\mathbf{P}}_i $ is defined to provide the $ i^{th} $ block from the complete signal via the standard multiplication ${\mathbf{P}}_i{\mathbf{x}} = {\mathbf{x}}_i$. Note that $ {\mathbf{P}}_i $ can extract any block of the signal, even one that is not in the non-overlapping grid $ \mathcal{B} $. Accordingly, the matrix $ {\mathbf{P}}_i^T $ locates a block in the $ i^{th} $ block-position in a construction of a full-sized signal and, therefore, lets to express the a complete signal as ${\mathbf{x}} = \mathop \sum\limits_{i\in\mathcal{B}} {\mathbf{P}}_i^T{\mathbf{x}}_i$.

Now we can use the block handling operator $ {\mathbf{P}}_i $ for expressing the block-based rate-distortion optimization in its corresponding full-signal formulation:

[rCl]{} \[eq:rate-distortion optimization - full signal\] = \_[ \_ ]{} [ \_2\^2]{} + r\_[tot]{}() .

where $ \mathcal{C}_{\mathcal{B}} $ is the full-signal codebook, being the discrete set of candidate reconstructions for the full signal, defined using the block-level codebook $ \mathcal{C} $ as

[rCl]{} \[eq:group of solutions based on nonoverlapping blocks\] \_ = .

Moreover, the regularization function in (\[eq:rate-distortion optimization - full signal\]) is the total bit cost of the reconstructed signal defined for $ {\mathbf{v}} \in \mathcal{C}_{\mathcal{B}} $ as $ r_{tot}({\mathbf{v}}) \triangleq \nolinebreak \sum\limits_{i\in \mathcal{B}} r({\mathbf{P}}_i {\mathbf{v}})$.

Complexity-Regularized Restoration: Basic Optimization Formulation {#subsec:Complexity-Regularized Restoration}
------------------------------------------------------------------

While the regularized-restoration optimization in (\[eq:restoration using a general prior\]) is over a continuous domain, the operational rate-distortion optimization in (\[eq:rate-distortion optimization - full signal\]) is a discrete problem with solutions limited to the set $ \mathcal{C}_{\mathcal{B}} $. Therefore, we extend the definition of the block bit-cost evaluation function such that it is defined for any $ {\mathbf{z}} \in \mathbb{R}^{N_b} $ via

[rCl]{} \[eq:extended block bit-cost\] |r () = {

[\*[20]{}[c]{}]{} [r( )]{}&[, ]{}\
&[, ]{}

. ,

and the corresponding extension of the total bit-cost $\bar{r}_{tot}({\mathbf{x}})\triangleq\nolinebreak \sum\limits_{i\in \mathcal{B}} \bar r({\mathbf{P}}_i {\mathbf{x}})$ is defined for any $ {\mathbf{x}} \in \mathbb{R}^{N} $.

Now we define the complexity regularization function as

[rCl]{} \[eq:complexity regularization function\] s()=|[r]{}\_[tot]{}()

and the corresponding restoration optimization is

[rCl]{} \[eq:complexity-regularized restoration - full signal\] = \_ \_2\^2 + |[r]{}\_[tot]{}().

Due to the definition of the extended bit-cost evaluation function, $\bar{r}_{tot}({\mathbf{x}})$, the solution candidates of (\[eq:complexity-regularized restoration - full signal\]) are limited to the discrete set $ \mathcal{C}_{\mathcal{B}}  $ as defined in (\[eq:group of solutions based on nonoverlapping blocks\]).

Examining the complexity-regularized restoration in (\[eq:complexity-regularized restoration - full signal\]) for the Gaussian denoising task, where $ {\mathbf{H}} = {\mathbf{I}} $, shows that the optimization reduces to the regular rate-distortion optimization in (\[eq:rate-distortion optimization - full signal\]), namely, the compression of the noisy signal $ {\mathbf{y}} $. However, for more complicated restoration problems, where $ {\mathbf{H}} $ has an arbitrary structure, the optimization in (\[eq:complexity-regularized restoration - full signal\]) is not easy to solve and, in particular, it does not correspond to standard compression designs that are optimized for the regular squared-error distortion metric.

Proposed Methods {#sec:Proposed Methods}
================

In this section we present three restoration methods leveraging a given compression technique. The proposed algorithms result from two different definitions for the complexity regularization function. While the first approach regularizes the total bit-cost of the non-overlapping blocks of the restored signal, the other two refer to the total bit-cost of all the overlapping blocks of the estimate.

Regularize Total Complexity of Non-Overlapping Blocks {#subsec:Total Complexity of Non-Overlapping Blocks}
-----------------------------------------------------

Here we establish a practical method addressing the optimization problem in (\[eq:complexity-regularized restoration - full signal\]) based on the alternating direction method of multipliers (ADMM) approach [@boyd2011distributed] (for additional uses see, e.g., [@venkatakrishnan2013plug; @sreehari2016plug; @dar2016postprocessing; @rond2016poisson; @afonso2010fast]). The optimization (\[eq:complexity-regularized restoration - full signal\]) can be expressed also as

[rCl]{} \[eq:complexity-regularized restoration - expressing blocks\] = \_ \_2\^2 + \_[i]{} |r(\_i),

where the degradation matrix $ {\mathbf{H}} $, having a general structure, renders a block-based treatment infeasible.

Addressing this structural difficulty using the ADMM strategy [@boyd2011distributed] begins with introducing the auxiliary variables $ \left\lbrace {\mathbf{z}}_i \right\rbrace_{i\in\mathcal{B}} $, where $ {\mathbf{z}}_i $ is coupled with the $ i^{th} $ non-overlapping block. Specifically, we reformulate the problem (\[eq:complexity-regularized restoration - expressing blocks\]) into

[rCl]{} \[eq:complexity-regularized restoration - with splitting - constrained form\] &&( , \_i \_[i]{} ) = \_[,\_i \_[i]{}]{} \_2\^2 + \_[i]{} [|[r]{}(\_i)]{}\
&&       \_i = \_i   i .

Then, reformulating the constrained optimization (\[eq:complexity-regularized restoration - with splitting - constrained form\]) using the augmented Lagrangian (in its scaled form) and the method of multipliers (see [@boyd2011distributed Ch. 2]) leads to the following iterative procedure

[rCl]{} \[eq:complexity-regularized restoration - with splitting\] ( \^[(t)]{}, \_i\^[(t)]{} \_[i]{} ) &=& \_[,\_i \_[i]{}]{} \_2\^2 + \_[i]{} |r(\_i)       \
&& + \_[i]{} [ \_2\^2]{}\
\_i\^[(t+1) ]{} & = & \_i\^[(t) ]{} + ( \_i \^[(t)]{} - \_i\^[(t)]{} )    i ,

where $ t $ is the iteration number, $ \beta $ is a parameter originating in the augmented Lagrangian, and $ {\mathbf{u}}_i^{(t) } \in \mathbb{R}^{N_b} $ is the scaled dual variable corresponding to the $ i^{th} $ block (where $ i\in \mathcal{B} $).

Each of the optimization variables in (\[eq:complexity-regularized restoration - with splitting\]) participates only in part of the terms of the cost function and, therefore, employing one iteration of alternating minimization (see [@boyd2011distributed Ch. 2]) leads to the ADMM form of the problem, where the included optimizations are relatively simple. Accordingly, the $ t^{th} $ iteration of the proposed iterative solution is

[rCl]{} \[eq:complexity-regularized restoration - iterative solution - inversion\] && \^[(t)]{} = \_ \_2\^2 + \_[i]{} [ \_2\^2]{}    \
\[eq:complexity-regularized restoration - iterative solution - compression\] && \_i\^[(t)]{} = \_[\_i]{} + |r(\_i),    i\
\[eq:complexity-regularized restoration - iterative solution - beta update\] && \_i\^[(t+1)]{} = \_i\^[(t)]{} + ( \_i \^[(t)]{} - \_i\^[(t)]{} ) ,    i ,

where $ {\tilde{{\mathbf{z}}}}_{i}^{(t)} \triangleq {\hat{{\mathbf{z}}}}_{i}^{(t-1)} - {\mathbf{u}}_i^{(t) } $ and $ {\tilde{{\mathbf{x}}}}_{i}^{(t)} \triangleq {\mathbf{P}}_i \hat{{\mathbf{x}}}^{(t)} + {\mathbf{u}}_i^{(t)} $ for $  i\in\mathcal{B} $.

The analytic solution of the first stage optimization in (\[eq:complexity-regularized restoration - iterative solution - inversion\]) is

[rCl]{} \[eq:complexity-regularized restoration - iterative solution - inversion - analytic form\] \^[(t)]{} = ( \^[T]{} + )\^[-1]{} ( \^T + \_[i]{}[ \_i\^T \_[i]{}\^[(t)]{} ]{} )\

rendering this stage as a weighted averaging of the deteriorated signal with the block estimates obtained in the second stage of the previous iteration. While the analytic solution (\[eq:complexity-regularized restoration - iterative solution - inversion - analytic form\]) explains the underlying meaning of the $ \ell_2 $-constrained deconvolution stage (\[eq:complexity-regularized restoration - iterative solution - inversion\]), it includes matrix inversion that, in general, may lead to numerical instabilities. Accordingly, in the implementation of the proposed method we suggest to address (\[eq:complexity-regularized restoration - iterative solution - inversion\]) via numerical optimization techniques (for example, we used the biconjugate gradients method).

The optimizations in the second stage of each iteration (\[eq:complexity-regularized restoration - iterative solution - compression\]) are rate-distortion optimizations corresponding to each of the non-overlapping blocks of the signal estimate $ \hat{ {\mathbf{x}}}^{(t)} $ obtained in the first stage. Accordingly, the set of block-level optimizations in (\[eq:complexity-regularized restoration - iterative solution - compression\]) can be interpreted as a single full-signal rate-distortion optimization with respect to a Lagrange multiplier value of $ {\lambda} = \frac{2\mu}{\beta} $. We denote the compression-decompression procedure that replaces (\[eq:complexity-regularized restoration - iterative solution - compression\]) as

[rCl]{} \[eq:compression-decompression function\] \^[(t)]{}=CompressDecompress\_ ( \^[(t)]{} ),

where $  \tilde{{\mathbf{x}}}^{(t)} \triangleq  \mathop \sum\limits_{i\in\mathcal{B}} {\mathbf{P}}_i^T \tilde{{\mathbf{x}}}_i^{(t)} $ is the signal to compress, assembled from all the non-overlapping blocks, and $ \hat{{\mathbf{z}}}^{(t)} $ is the corresponding decompressed full signal. Moreover, by defining a full-sized scaled dual variable ${{\mathbf{u}}}^{(t)} \triangleq  \mathop \sum\limits_{i\in\mathcal{B}} {\mathbf{P}}_i^T{{\mathbf{u}}}_i^{(t)} $ we get that $  \tilde{{\mathbf{x}}}^{(t)} =  \hat{{\mathbf{x}}}^{(t)} + {{\mathbf{u}}}^{(t)}$. Then, using the definitions established here we can translate the block-level computations (\[eq:complexity-regularized restoration - iterative solution - inversion\])-(\[eq:complexity-regularized restoration - iterative solution - beta update\]) into the full-signal formulations described in Algorithm \[Algorithm:Proposed Method Non-overlapping\].

We further suggest using a standardized compression method as the compression-decompression operator (\[eq:compression-decompression function\]). While many compression methods do not follow the exact rate-distortion optimizations we have in our mathematical development, we still encourage utilizing such techniques as an approximation for (\[eq:complexity-regularized restoration - iterative solution - compression\]). Additionally, since many compression methods do not rely on Lagrangian optimization, their operating parameters may have different definitions such as quality parameters, compression ratios, or output bit-rates. Accordingly, we present the suggested algorithm with respect to a general compression-decompression procedure with output bit-cost directly or indirectly affected by a parameter denoted as $ \theta $. These generalizations are also implemented in the proposed Algorithm \[Algorithm:Proposed Method Non-overlapping\]. In Section \[sec:Experimental Results\] we elaborate on particular settings of $ \theta $ that were empirically found appropriate for utilization of the HEVC and the JPEG2000 standard. In cases where the compression method significantly deviates from a Lagrangian optimization form, it can be useful to appropriately update the compression parameter in each iteration (this is the case for JPEG2000 as explained in Section \[sec:Experimental Results\]).

Importantly, Algorithm \[Algorithm:Proposed Method Non-overlapping\] does not only restore the deteriorated input image, but also provides the signal estimate in a compressed form by employing the output of the compression stage of the last iteration.

Inputs: ${\mathbf{y}}$, $ \beta $, $ \theta $. Initialize $ {\hat{{\mathbf{z}}}}^{(0)} $ (depending on the deterioration type). $t = 1$ and $ {{\mathbf{u}}}^{(1)} = {\mathbf{0}} $

$\tilde{{\mathbf{z}}}^{(t)} =  {\hat{{\mathbf{z}}}}^{(t-1)} - {\mathbf{u}}^{(t)}$ Solve the $ \ell_2 $-constrained deconvolution:$~~~~~~~~\hat{ {\mathbf{x}}}^{(t)} = \mathop {\text{argmin}}\limits_{{\mathbf{x}}} \left\| {  {\mathbf{H}} {\mathbf{x}} - {\mathbf{y}} } \right\|_2^2 + \frac{\beta}{2}  {\left\| {   {\mathbf{x}} - \tilde{{\mathbf{z}}}^{(t)} } \right\|_2^2} $

$\tilde{{\mathbf{x}}}^{(t)} =  {\hat{{\mathbf{x}}}}^{(t)} + {\mathbf{u}}^{(t)}$ $ \hat{{\mathbf{z}}}^{(t)} = {CompressDecompress}_{\theta}\left( \tilde{{\mathbf{x}}}^{(t)} \right) $

${\mathbf{u}}^{(t+1)} = {\mathbf{u}}^{(t)} + \left( \hat{{\mathbf{x}}}^{(t)} - \hat{{\mathbf{z}}}^{(t)} \right)$ $ t \gets t + 1$

Regularize Total Complexity of All Overlapping Blocks {#subsec:Complexity of Overlapping Blocks}
-----------------------------------------------------

Algorithm \[Algorithm:Proposed Method Non-overlapping\] emerged from complexity regularization measuring the total bit-cost of the estimate based on its decomposition into non-overlapping blocks (see Eq. (\[eq:complexity-regularized restoration - expressing blocks\])), resulting in a restored signal available in a compressed form compatible with the compression technique in use. Obviously, the above approach provides estimates limited to the discrete set of signals supported by the compression architecture, thus, having a somewhat reduced restoration ability with respect to methods providing estimates from an unrestricted domain of solutions. This observation motivates us to develop a complexity-regularized restoration procedure that provides good estimates from the continuous unrestricted domain of signals while still utilizing a standardized compression technique as its main component.

As before, our developments refer to a general block-based compression method relying on a codebook $ \mathcal{C} $ as a discrete set of block reconstruction candidates. We consider here the segmentation of the signal-block space, $ \mathbb{R}^{N_b} $, given by the voronoi cells corresponding to the compression reconstruction candidates, namely, for each $ {\mathbf{c}} \in \mathcal{C} $ there is a region

[rCl]{} \[eq:voronoi cell of codebook element\] V\_ \^[N\_b]{}    |[   = \_ - \_2\^2 ]{}

defining all the vectors in $ \mathbb{R}^{N_b} $ that $ {\mathbf{c}} $ is their nearest member of $ \mathcal{C} $. We use the voronoi cells in (\[eq:voronoi cell of codebook element\]) for defining an alternative extension to the bit-cost evaluation of a signal block (i.e., the new definition, $\bar{r}_{v}({\mathbf{z}})$, will replace $\bar r({\mathbf{z}})$ given in (\[eq:extended block bit-cost\]) that was used for the development of Algorithm \[Algorithm:Proposed Method Non-overlapping\]). Specifically, we associate a finite bit-cost to any $ {\mathbf{z}} \in \mathbb{R}^{N_b} $ based on the voronoi cell it resides in, i.e.,

[rCl]{} \[eq:voronoi extended block bit-cost\] |[r]{}\_[v]{}() = r() V\_

where $ r({\mathbf{c}}) $ is the regular bit-cost evaluation defined in Section \[subsec:Operational Rate-Distortion Optimization\] only for blocks in $\mathcal{C}$.

The method proposed here emerges from a new complexity regularization function that quantifies the total complexity of all the overlapping blocks of the estimate. Using the extended bit-cost measure $ \bar{r}_{v}(\cdot) $, defined in (\[eq:voronoi extended block bit-cost\]), we introduce the full-signal regularizer as

[rCl]{} \[eq:complexity regularizer of all overlapping blocks\] s\^[\*]{}()=\_[i\^\*]{} |[r]{}\_v (\_i)

where $ {\mathbf{x}} \in \mathbb{R}^N $, and $ \mathcal{B}^* $ is a set containing the indices of all the overlapping blocks of the signal. The associated restoration optimization is

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks\] = \_ \_2\^2 + \_[i\^\*]{} |[r]{}\_v (\_i).

Importantly, in contrast to the previous subsection, the function $ s^{*}({\mathbf{x}}) $ evaluates the complexity of any $ {\mathbf{x}} \in \mathbb{R}^N $ with a finite value and, thus, does not restrict the restoration to the discrete set of codebook-based constructions, $ \mathcal{C}_{\mathcal{B}} $, defined in (\[eq:group of solutions based on nonoverlapping blocks\]).

While the new regularizer in (\[eq:complexity-regularized restoration - overlapping blocks\]) is not separable into complexity evaluation of non-overlapping blocks, the ADMM approach can accommodate it as well. This is explained next. We define the auxiliary variables $ \left\lbrace {\mathbf{z}}_i \right\rbrace_{i\in\mathcal{B}^*} $, where each $ {\mathbf{z}}_i $ is coupled with the $ i^{th} $ overlapping block. Then, the optimization (\[eq:complexity-regularized restoration - overlapping blocks\]) is expressed as

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - with splitting - constrained form\] &&( , \_i \_[i\^\*]{} ) = \_[,\_i \_[i\^\*]{}]{} \_2\^2 + \_[i\^\*]{} [|[r]{}\_v (\_i)]{}\
&&       \_i = \_i   i\^\* .

As in Section \[subsec:Total Complexity of Non-Overlapping Blocks\], employing the augmented Lagrangian (in its scaled form) and the method of multipliers results in an iterative solution provided by the following three steps in each iteration (as before $ t $ denotes the iteration number):

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - iterative solution - inversion\] && \^[(t)]{} = \_ \_2\^2 + \_[i\^[\*]{}]{} [ \_2\^2]{}    \
\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression\] && \_i\^[(t)]{} = \_[\_i]{} + |r(\_i),    i\^[\*]{}\
\[complexity-regularized restoration - overlapping blocks - iterative solution - beta update\] && \_i\^[(t+1)]{} = \_i\^[(t)]{} + ( \_i \^[(t)]{} - \_i\^[(t)]{} ) ,    i\^[\*]{} ,

where $ {\mathbf{u}}_i^{(t)} $ is the scaled dual variable for the $ i^{th} $ block, $ {\tilde{{\mathbf{z}}}}_{i}^{(t)} \nolinebreak \triangleq \nolinebreak  {\hat{{\mathbf{z}}}}_{i}^{(t-1)} - {\mathbf{u}}_i^{(t) } $ and $ {\tilde{{\mathbf{x}}}}_{i}^{(t)} \triangleq {\mathbf{P}}_i \hat{{\mathbf{x}}}^{(t)} + {\mathbf{u}}_i^{(t)} $ for $  i\in\mathcal{B}^{*} $.

While the procedure above resembles the one from the former subsection, the treatment of overlapping blocks has different interpretations to the optimizations in (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - inversion\]) and (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression\]). Indeed, note that the block-level rate-distortion optimizations in (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression\]) are not discrete due to the extended bit-cost evaluation $ \bar{r}_v(\cdot) $ defined in (\[eq:voronoi extended block bit-cost\]). Due to the definition of $ \bar{r}_v(\cdot) $, the rate-distortion optimizations (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression\]) can be considered as continuous relaxations of the discrete optimizations done by the practical compression technique. Since we intend using a given compression method without explicit knowledge of its underlying codebook, we cannot construct the voronoi cells defining $ \bar{r}_v(\cdot) $ and, thus, it is impractical to accurately solve (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression\]). Consequently, we suggest to approximate the optimizations (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression\]) by the discrete forms of

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression - discrete\] && \_i\^[(t)]{} = \_[\_i]{} + |[r]{}(\_i),   i\^\*

where $ \bar{r}(\cdot) $ is the discrete evaluation of the block bit-cost, defined in (\[eq:extended block bit-cost\]), letting to identify the problems as operational rate-distortion optimizations of the regular discrete form.

Each block-level rate-distortion optimization in (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - compression - discrete\]) is associated with one of the overlapping blocks of the signal. Accordingly, we interpret this group of optimizations as multiple applications of a full-signal compression-decompression procedure, each associates to a shifted version of the signal (corresponding to different sets of non-overlapping blocks). Specifically, for a signal $ {\mathbf{x}} $ and a compression block-size of $ N_b $ samples, there are $ N_b $ shifted grids of non-overlapping blocks. For mathematical convenience, we consider here cyclic shifts such that the $ j^{th} $ shift ($ j=1,...,N_b $) corresponds to a signal of $ N $ samples taken cyclically starting at the $ j^{th} $ sample of $ {\mathbf{x}} $ (in practice other definitions of shifts may be used, e.g., see Section \[sec:Experimental Results\] for a suggested treatment of two-dimensional signals). We denote the $ j^{th} $ shifted signal as $ shift_{j}\left\lbrace {\mathbf{x}} \right\rbrace $. Moreover, we denote the index set of blocks included in the $ j^{th} $ shifted signal as $ \mathcal{B}^{j} $ (noting that $ \mathcal{B}^{1} =  \mathcal{B} $), hence, $ \mathcal{B}^{*} = \cup_{j=1}^{N_{b}}\mathcal{B}^{j} $. Therefore, the decompressed blocks $\left\lbrace{\hat{{\mathbf{z}}}_i^{(t)}}\right\rbrace_{i\in\mathcal{B}^*}$ can be decomposed into $ N_{b} $ subsets, $\left\lbrace{\hat{{\mathbf{z}}}_i^{(t)}}\right\rbrace_{i\in\mathcal{B}^j}$ for $ j=1,...,N_b $, each contains non-overlapping blocks corresponding to a different shifted grid. Moreover, the $ j^{th} $ set of blocks, $\left\lbrace{\hat{{\mathbf{z}}}_i^{(t)}}\right\rbrace_{i\in\mathcal{B}^j}$, is associated with the full signal $ \hat{{\mathbf{z}}}^{j,(t)} \triangleq \sum\limits_{i\in \mathcal{B}^j}{ {\mathbf{P}}_i^T  \hat{{\mathbf{z}}}_{i}^{(t)}}  $. Then, the set of full signals $\left\lbrace{\hat{{\mathbf{z}}}^{j,(t)}}\right\rbrace_{j=1}^{N_b}$ can be obtained by multiple full-signal compression-decompression applications, namely, for $ j\nolinebreak=\nolinebreak 1,...,N_b $: $ \hat{{\mathbf{z}}}_{shifted}^{j,(t)}\nolinebreak=\nolinebreak CompressDecompress_{ \lambda} \left( \tilde{{\mathbf{x}}}^{j,(t)}_{shifted} \right) $, where the Lagrangian multiplier value is $ {\lambda} = \frac{2\mu}{\beta} $ and

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - shifted compression input - definition\] && \^[j,(t)]{}\_[shifted]{} shift\_[j]{}\^[(t)]{} + \^[j,(t)]{}

is the compression input formed as the $ j^{th} $ shift of $ \hat{{\mathbf{x}}}^{(t)} $ combined with the full-sized dual variable defined via

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - shifted dual variable - definition\] && \^[j,(t)]{} \_[i\^j]{}[ \_i\^T \_[i]{}\^[(t)]{}]{}

assembled from the block-level dual variables corresponding to the $ j^{th} $ grid of non-overlapping blocks. Notice that inverse shifts are required for obtaining the desired signals, i.e.,

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - shifted decompression output - definition\] && \^[j,(t)]{} = shift\_[j]{}\^[-1]{}\_[shifted]{}\^[j,(t)]{}

where $ shift^{-1}_{j} \left\lbrace \cdot \right\rbrace$ is the inverse shift operator that (cyclically) shifts back the given full-size signal by $ j $ samples.

The deconvolution stage (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - inversion\]) of the iterative process can be rewritten as

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - iterative solution - inversion - full signals\] && \^[(t)]{} = \_ \_2\^2 + \_[j=1]{}\^[N\_b]{} [ \_2\^2]{}    

where the regularization part (the second term) considers the distance of the estimate from the $ N_b $ full signals defined via

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - iterative solution - inversion - full signals - auxiliary definition\] && \^[j,(t)]{} \^[j,(t)]{} - \^[j,(t)]{}

for $ j=1,...,N_b $, where $ {\hat{{\mathbf{z}}}}^{j,(t)} $ and $ {{\mathbf{u}}}^{j,(t)} $ were defined above. The analytic solution of the optimization (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - inversion - full signals\]) is

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - inversion - analytic form - full signals\] &&\^[(t)]{} = ( \^[T]{} + N\_b )\^[-1]{} ( \^T + \_[j=1]{}\^[N\_b]{}[ \^[j,(t)]{} ]{} ) ,   

showing that the first stage of each iteration is a weighted averaging of the given deteriorated signal with all the decompressed signals (and the dual variables) obtained in the former iteration. It should be noted that the analytic solution (\[eq:complexity-regularized restoration - overlapping blocks - inversion - analytic form - full signals\]) is developed here for showing the essence of the $ \ell_2 $-constrained deconvolution part of the method. Nevertheless, the possible numerical instabilities due to the matrix inversion appearing in (\[eq:complexity-regularized restoration - overlapping blocks - inversion - analytic form - full signals\]) motivate the practical direct treatment of (\[eq:complexity-regularized restoration - overlapping blocks - iterative solution - inversion\]) via numerical optimization techniques.

Algorithm \[Algorithm:Proposed Method Overlapping Blocks\] summarizes the practical restoration method for a compression technique operated by the general parameter $ {\theta} $ for determining the bit-cost (see details in Section \[subsec:Total Complexity of Non-Overlapping Blocks\]).

The computational cost of Algorithm \[Algorithm:Proposed Method Overlapping Blocks\] stems from its reliance on repeated applications of compressions, decompressions, and $ \ell_2 $ - constrained deconvolution procedures. While the actual run-time depends on the computational complexity of the utilized compression technique, we can generally state that the total run-time will be of at least the run-time of compression and decompression processes for a total number of applications equal to the product of the number of iterations and the number of shifts considered.

Inputs: $ {\mathbf{y}} $, $\beta$, $ \theta $. Initialize $ \left\lbrace {\hat{{\mathbf{z}}}}^{j,(0)} \right\rbrace_{j=1}^{N_b} $ (depending on the deterioration type). $t = 1$ and $ {{{\mathbf{u}}}}^{j,(1)} = {\mathbf{0}} $ for $ j = 1,...,{N_b} $.

${\tilde{{\mathbf{z}}}}^{j,(t)} = {\hat{{\mathbf{z}}}}^{j,(t-1)} - {\mathbf{u}}^{j,(t) }  \text{~~for~}j=1,...,N_b  $ Solve the $ \ell_2 $-constrained deconvolution:$~~~~\hat{ {\mathbf{x}}}^{(t)} = \mathop {\text{argmin}}\limits_{{\mathbf{x}}} \left\| {  {\mathbf{H}} {\mathbf{x}} - {\mathbf{y}} } \right\|_2^2 + \frac{\beta}{2} \sum\limits_{j=1}^{N_b} {\left\| {  {\mathbf{x}} - {\tilde{{\mathbf{z}}}}^{j,(t)}} \right\|_2^2} $

$\tilde{{\mathbf{x}}}_{shifted}^{j,(t)} = shift_{j}\left\lbrace \hat{{\mathbf{x}}}^{(t)}  + {{\mathbf{u}}}^{j,(t)} \right\rbrace$

$ \hat{{\mathbf{z}}}_{shifted}^{j,(t)} = {CompressDecompress}_{\theta}\left( \tilde{{\mathbf{x}}}_{shifted}^{j,(t)}  \right) $

$\hat{{\mathbf{z}}}^{j,(t)} = shift^{-1}_{j}\left\lbrace { \hat{{\mathbf{z}}}_{shifted}^{j,(t)} } \right\rbrace$

${\mathbf{u}}^{j,(t+1)} = {\mathbf{u}}^{j,(t)} + \left( \hat{{\mathbf{x}}}^{(t)} - \hat{{\mathbf{z}}}^{j,(t)} \right) $

$ t \gets t + 1$

The ADMM is known for promoting distributed optimization structures [@boyd2011distributed]. In Algorithm \[Algorithm:Proposed Method Overlapping Blocks\] the distributed nature of the ADMM is expressed in the separate optimization of each of the shifted block-grids (see stages 8-11). In particular, the dual variables $ \left\lbrace {\mathbf{u}}^{j,(t)} \right\rbrace_{j=1}^{N_b} $, associated with the various grids (see stages 5,8, and 11 in Algorithm \[Algorithm:Proposed Method Overlapping Blocks\]), are updated independently in stage 11 such that each considers only its respective $ \hat{{\mathbf{z}}}^{j,(t)} $. However, the dual variables $ \left\lbrace {\mathbf{u}}^{j,(t)} \right\rbrace_{j=1}^{N_b} $ essentially refer to the same data based on different block-grids. Accordingly, we suggest to merge the independent dual variables to form a single, more robust, dual variable defined as

[rCl]{} \[eq:complexity-regularized restoration - overlapping blocks - robust dual variable\] && \_[total]{}\^[(t)]{} = \_[j=1]{}\^[N\_b]{}[ \^[j,(t)]{} ]{}

where the averaging tends to reduce particular artifacts that may appear due to specific block-grids. We utilize the averaged dual variable (\[eq:complexity-regularized restoration - overlapping blocks - robust dual variable\]) to extend Algorithm \[Algorithm:Proposed Method Overlapping Blocks\] into Algorithm \[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\]. Notice stages 5,8, and 13 of Algorithm \[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\], where the averaged dual variable is used instead of the independent ones.

Inputs: $ {\mathbf{y}} $, $\beta$, $ \theta $. Initialize $ \left\lbrace {\hat{{\mathbf{z}}}}^{j,(0)} \right\rbrace_{j=1}^{N_b} $ (depending on the deterioration type). $t = 1$ and $ {\mathbf{u}}_{total}^{(1)} = {\mathbf{0}} $.

${\tilde{{\mathbf{z}}}}^{j,(t)} = {\hat{{\mathbf{z}}}}^{j,(t-1)} - {\mathbf{u}}_{total}^{(t)}  \text{~~for~}j=1,...,N_b  $ Solve the $ \ell_2 $-constrained deconvolution:$~~~~\hat{ {\mathbf{x}}}^{(t)} = \mathop {\text{argmin}}\limits_{{\mathbf{x}}} \left\| {  {\mathbf{H}} {\mathbf{x}} - {\mathbf{y}} } \right\|_2^2 + \frac{\beta}{2} \sum\limits_{j=1}^{N_b} {\left\| {  {\mathbf{x}} - {\tilde{{\mathbf{z}}}}^{j,(t)}} \right\|_2^2} $

$\tilde{{\mathbf{x}}}_{shifted}^{j,(t)} = shift_{j}\left\lbrace \hat{{\mathbf{x}}}^{(t)}  + {\mathbf{u}}_{total}^{(t)} \right\rbrace$

$ \hat{{\mathbf{z}}}_{shifted}^{j,(t)} = {CompressDecompress}_{\theta}\left( \tilde{{\mathbf{x}}}_{shifted}^{j,(t)}  \right) $

$\hat{{\mathbf{z}}}^{j,(t)} = shift^{-1}_{j}\left\lbrace { \hat{{\mathbf{z}}}_{shifted}^{j,(t)} } \right\rbrace$

${\mathbf{u}}^{j,(t+1)} = {\mathbf{u}}^{j,(t)} + \left( \hat{{\mathbf{x}}}^{(t)} - \hat{{\mathbf{z}}}^{j,(t)} \right) $

${\mathbf{u}}_{total}^{(t+1)} = \frac{1}{N_b} \sum\limits_{j=1}^{N_b}{ {\mathbf{u}}^{j,(t+1)}  }$

$ t \gets t + 1$

In Section \[sec:Experimental Results\] we further discuss practical aspects of the proposed Algorithms \[Algorithm:Proposed Method Non-overlapping\]-\[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\] and evaluate their performance for deblurring and inpainting of images.

Rate-Distortion Theoretic Analysis for the Gaussian Case {#sec:Rate-Distortion Theoretic Analysis for the Gaussian Case}
========================================================

In this section we theoretically study the complexity-regularized restoration problem from the perspective of rate-distortion theory. While our analysis is focused on the particular settings of a cyclo-stationary Gaussian signal and deterioration caused by a linear shift-invariant operator and additive white Gaussian noise, the results clearly explain the main principles of complexity-regularized restoration.

In general, theoretical studies of rate-distortion problems for the Gaussian case provide to the signal processing practice optimistic beliefs about which design concepts perform well for the real-world non-Gaussian instances of the problems (see the excellent discussion in [@donoho1998data Sec. 3]). Moreover, theoretical and practical solutions may embody in a different way the same general concepts. Therefore, one should look for connections between theory and practice in the form of high-level analogies.

The optimal solution presented in this section considers the classical framework of rate-distortion theory and a particular, however, important case of a Gaussian signal and a linear shift-invariant degradation operator. Our rate-distortion analysis below will show that the optimal complexity-regularized restoration consists of the following two main ideas: pseudoinverse filtering of the degraded input, and compression with respect to a squared-error metric that is weighted based on the degradation-filter squared-magnitude (considering a processing in the Discrete Fourier Transform (DFT) domain). In Subsection \[subsec:Rate-Distortion Theoretic - Conceptual Relation to The Proposed Approach\] we explain how these two concepts connect to more general themes having different realizations in the practical approach proposed in Section \[sec:Proposed Methods\].

In this section, consider the signal $ {\mathbf{x}} \in \mathbb{R}^N $ modeled as a zero-mean cyclo-stationary Gaussian random vector with a circulant autocorrelation matrix $ {\mathbf{R}}_{{\mathbf{x}}} $, i.e., ${\mathbf{x}} \sim  \mathcal{N}\left(0,{\mathbf{R}}_{{\mathbf{x}}}\right) $. The degradation model studied remains

[rCl]{} \[eq:theoretic Gaussian analysis - degradation model\] = + ,

where here $ {\mathbf{H}} $ is a real-valued $N\times N$ circulant matrix representing a linear shift-invariant deteriorating operation and ${\mathbf{n}} \sim  \mathcal{N}\left(0,\sigma_n^2{\mathbf{I}}\right) $ is a length $N$ vector of white Gaussian noise. Clearly, the degraded observation $ {\mathbf{y}} $ is also a zero-mean cyclo-stationary Gaussian random vector with a circulant autocorrelation matrix $ {\mathbf{R}}_{{\mathbf{y}}}\nolinebreak=\nolinebreak{\mathbf{H}}  {\mathbf{R}}_{{\mathbf{x}}} {\mathbf{H}}^* + \sigma_n^2{\mathbf{I}} $.

Prevalent Restoration Strategies {#subsec:Alternative Restoration Strategies}
--------------------------------

We precede the analysis of the complexity-regularized restoration with mentioning three well-known estimation methods. The restoration procedure is a function

[rCl]{} \[eq:theoretic Gaussian analysis - restoration function\] = f ( ),

where $ f $ maps the degraded signal $ {\mathbf{y}} $ to an estimate of $ {\mathbf{x}} $ denoted as $ \hat{{\mathbf{x}}} $. In practice, one gets a realization of $ {\mathbf{y}} $ denoted here as $ {\mathbf{y}}_r $ and forms the corresponding estimate as $ {\hat{{\mathbf{x}}}}_r \nolinebreak = \nolinebreak f \left( {\mathbf{y}}_r \right)  $.

### Minimum Mean Squared Error (MMSE) Estimate

This restoration minimizes the expected MSE of the estimate, i.e.,

[rCl]{} \[eq:theoretic Gaussian analysis - MMSE estimate optimization\] f\_[MMSE]{} =   E \_2\^2 ,

yielding that the corresponding estimate is the conditional expectation of $ {\mathbf{x}} $ given $ {\mathbf{y}} $

[rCl]{} \[eq:theoretic Gaussian analysis - MMSE estimate\] \_[MMSE]{} = f\_[MMSE]{} ( ) = E .

Nicely, for the Gaussian case considered in this section, the MMSE estimate (\[eq:theoretic Gaussian analysis - MMSE estimate\]) reduces to a linear operator, presented below as the Wiener filter.

### Wiener Filtering

The Wiener filter is also known as the Linear Minimum Mean Squared Error (LMMSE) estimate, corresponding to a restoration function of the form

[rCl]{} \[eq:theoretic Gaussian analysis - Wiener restoration function\] = f\_[Wiener]{} ( ) = + ,

optimized via

[rCl]{} \[eq:theoretic Gaussian analysis - Wiener estimate optimization\] =   E \_2\^2 .

In our case, where $ {\mathbf{x}} $ and $ {\mathbf{y}} $ are zero mean, $ \hat{{\mathbf{b}}} = {\mathbf{0}} $ and

[rCl]{} \[eq:theoretic Gaussian analysis - Wiener filter matrix\] = \_ \^\* ( \_ \^\* + \_n\^2 )\^[-1]{} .

If the distributions are Gaussian, this linear operator coincides with the optimal MMSE estimator.

### Constrained Deconvolution Filtering

This approach considers a given degraded signal $ {\mathbf{y}}_r = {\mathbf{H}} {\mathbf{x}}_0 + {\mathbf{n}}_r $, with the noise vector a realization of a random process while the signal ${\mathbf{x}}_0$ is considered as a deterministic vector, with perhaps some known properties. Then, the restoration is carried out by minimizing a carefully-designed penalty function, $ g $, that assumes lower values for ${\mathbf{x}}$ vectors that fit the prior knowledge on ${\mathbf{x}}_0$. Note that for a sufficiently large signal dimension we get that $ \left\| { {\mathbf{y}}_r  - {\mathbf{H}} {\mathbf{x}}_0 } \right\|_2^2  = \left\| {\mathbf{n}}_r \right\|_2^2 \approx N \sigma_n^2 $. The last result motivates to constrain the estimate, $ \hat{{\mathbf{x}}} $, to conform with the known degradation model (\[eq:theoretic Gaussian analysis - degradation model\]), by demanding the similarity of $  { {\mathbf{y}}_r  - {\mathbf{H}} \hat{{\mathbf{x}}} } $ to the additive noise term via the equality relation $ \left\| { {\mathbf{y}}_r  - {\mathbf{H}} \hat{{\mathbf{x}}} } \right\|_2^2  = N \sigma_n^2 $. The above idea is implemented in an optimization of the form

[rCl]{} \[eq:theoretic Gaussian analysis - constrained deconvolution filtering\]

& & & [ g( ) ]{}\
& & & \_2\^2 = N \_n\^2 .

Our practical methods presented in Section \[sec:Proposed Methods\] emerge from an instance of the constrained deconvolution optimization (\[eq:theoretic Gaussian analysis - constrained deconvolution filtering\]), in its Lagrangian version, where the penalty function $ g $ is the cost in bits measuring the complexity in describing the estimate $ \hat{{\mathbf{x}}} $. In the remainder of this section, we study the complexity-regularized restoration problem from a statistical perspective.

The Complexity-Regularized Restoration Problem and its Equivalent Forms {#subsec:Reformulations of the Problem}
-----------------------------------------------------------------------

Based on rate-distortion theory (e.g., see [@cover2012elements]), we consider the estimate of $ {{\mathbf{x}}} $ as a random vector $ \hat{{\mathbf{x}}} \in \mathbb{R}^N $ with the probability density function (PDF) $ p_{\hat{{\mathbf{x}}}}\left(  {{\hat{x}}} \right) $. The estimate characterization, $ p_{\hat{{\mathbf{x}}}}\left(  {{\hat{x}}} \right) $, is determined by optimizing the conditional PDF $ p_{\hat{{\mathbf{x}}} | {\mathbf{y}} }\left( {{\hat{x}}} | {{y}}  \right) $, statistically representing the mapping between the given data $ {\mathbf{y}} $ and the decompression result $ \hat{{\mathbf{x}}} $. Moreover, the rate is measured as the mutual information between $\hat{{\mathbf{x}}}$ and $ {\mathbf{y}} $, defined via

[rCl]{} \[eq:theoretic Gaussian analysis - mutual information definition\] I( ; ) = .

Then, the basic form of the complexity-regularized restoration optimization is expressed as

\[problem: basic form\]

[rCl]{}

& & & [ I( ; ) ]{}\
& & & E \_2\^2 = N \_n\^2 .

Here the estimate rate is minimized while maintaining suitability to the degradation model (\[eq:theoretic Gaussian analysis - degradation model\]) using a distortion constraint set to achieve an a-priori known total noise energy. In general, Problem \[problem: basic form\] is complicated to solve since the distortion constraint considers $ \hat{{\mathbf{x}}} $ through the degradation operator $ {\mathbf{H}} $, while the rate is directly evaluated for $ \hat{{\mathbf{x}}} $.

The shift invariant operator $ {\mathbf{H}} $ is a circulant $ N \times N $ matrix, thus, diagonalized by the $ N\times N $ Discrete Fourier Transform (DFT) matrix $ {\mathbf{F}} $. The $ (k,l) $ component of the DFT matrix ($ k,l=0,...,N-1$) is ${\mathbf{F}}_{k,l}= W_N^{kl} $ where $ W_N \triangleq \frac{1}{\sqrt{N}} e^{-i2\pi/N} $. Then, the diagonalization of $ {\mathbf{H}} $ is expressed as

[rCl]{} \[eq:theoretic Gaussian analysis - diagonlization of H\] \^\* = \_H ,

where $ {\mathbf{\Lambda}}_H $ is a diagonal matrix formed by the components $ h^F_k $ for $ k=0,...,N-1 $. Using $ {\mathbf{\Lambda}}_H $ we define the pseudoinverse of $ {\mathbf{H}} $ as

[rCl]{} \[eq:theoretic Gaussian analysis - pseudoinverse of H\] \^[+]{} = \^\* \_H\^[+]{} ,

where $ {\mathbf{\Lambda}}_H^{+} $ is the pseudoinverse of $ {\mathbf{\Lambda}}_H $, an $ N \times N $ diagonal matrix with the $ k^{th} $ diagonal element:

[rCl]{} \[eq:theoretic Gaussian analysis - pseudoinverse of Lambda\_H - diagonal elements\] h\^[F,+]{}\_k = {

[\*[20]{}[c]{}]{} &[, h\^F\_k0 ]{}\
0&[, h\^F\_k = 0 . ]{}

.

We denote by $N_{H}$ the number of nonzero diagonal elements in $ {\mathbf{\Lambda}}_H $, the rank of $ {\mathbf{H}} $.

The first main result of our analysis states that Problem \[problem: basic form\], being the straightforward formulation for complexity-regularized restoration, is equivalent to the next problem.

\[problem: pseudoinverse filtered input\]

[rCl]{}

& & & [ I( ; ) ]{}\
& & & E ( [ - ]{} ) \_2\^2 = N\_[H]{} \_n\^2 ,

where

[rCl]{} \[eq:theoretic Gaussian analysis - pseudoinverse filtered input definition\] = \^[+]{}

is the pseudoinverse filtered version of the given degraded signal $ {\mathbf{y}} $. One should note that Problem \[problem: pseudoinverse filtered input\] has a more convenient form than Problem \[problem: basic form\] since the distortion is an expected weighted squared error between the two random variables determining the rate. The equivalence of Problems \[problem: basic form\] and \[problem: pseudoinverse filtered input\] is proved in Appendix \[appendix:Equivalence of Problems 1 and 2\].

In this section, $ {\mathbf{x}} $ is a cyclo-stationary Gaussian signal, hence, having a circulant autocorrelation matrix $ {\mathbf{R}}_{{\mathbf{x}}} $. Consequently, and also because $ {\mathbf{H}} $ is circulant, the deteriorated signal $ {\mathbf{y}} $ is also a cyclo-stationary Gaussian signal. Moreover, $ {\mathbf{H}}^{+} $ is also a circulant matrix, thus, by (\[eq:theoretic Gaussian analysis - pseudoinverse filtered input definition\]) the pseudoinverse filtering result, $ {\tilde{{\mathbf{y}}}} $, is also cyclo-stationary and zero-mean Gaussian. Specifically, the autocorrelation matrix of $ {\tilde{{\mathbf{y}}}} $ is

[rCl]{} \[eq:theoretic Gaussian analysis - pseudoinverse filtered input - autocorrelation matrix\] \_ & = & \^[+]{} \_ [\^[+]{}]{}\^\*\
& = & \^[+]{} \_ \^\* [\^[+]{}]{}\^\* + \_n\^2 \^[+]{} [\^[+]{}]{}\^\* ,

and, as a circulant matrix, it is diagonalized by the DFT matrix yielding the eigenvalues

[rCl]{} \[eq:theoretic Gaussian analysis - pseudoinverse filtered input - autocorrelation eigenvalues\] \^[()]{}\_k = {

[\*[20]{}[c]{}]{} [\^[()]{}\_k + ]{}&[, h\^F\_k0 ]{}\
0&[, h\^F\_k = 0 . ]{}

.

The DFT-domain representation of $ {\tilde{{\mathbf{y}}}} $ is

[rCl]{} \[eq:theoretic Gaussian analysis - pseudoinverse filtered input - in DFT domain\] \^F = ,

consisted of the coefficients $ \left\lbrace {\tilde{y}}^F_k \right\rbrace_{k=0}^{N-1} $, being independent zero-mean Gaussian variables with variances corresponding to the eigenvalues in (\[eq:theoretic Gaussian analysis - pseudoinverse filtered input - autocorrelation eigenvalues\]).

Transforming Problem \[problem: pseudoinverse filtered input\] to the DFT domain, where $ {\tilde{{\mathbf{y}}}} $ becomes a set of independent Gaussian variables to be coded under a joint distortion constraint, simplifies the optimization structure to the following separable form (see proof sketch in Appendix \[appendix:Equivalence of Problems 2 and 3\]).

\[problem: Separable Form in DFT domain\]

[rCl]{}

& & & [ \_[k=0 ]{}\^[N-1]{} I( [\^F\_k]{}; \^F\_k ) ]{}\
& & & \_[k=0]{}\^[N-1]{} [| h\^F\_k | \^2]{} E = N\_[H]{} \_n\^2 ,

where $ \left\lbrace {\hat{x}}^F_k \right\rbrace_{k=0}^{N-1} $ are the elements of $ {\hat{{\mathbf{x}}}}^F = {\mathbf{F}}{\hat{{\mathbf{x}}}} $. Nicely, the separable distortion in Problem \[problem: Separable Form in DFT domain\] considers each variable using a squared error that is weighted by the squared magnitude of the corresponding degradation-filter coefficient.

The rate-distortion function of a single Gaussian variable with variance $ \sigma^2 $ has the known formulation [@cover2012elements]:

[rCl]{} \[eq:theoretic Gaussian analysis - rate-distortion function of a Gaussian variable\] R(D) = \_+

evaluating the minimal rate for a squared-error allowed reaching up to $ D \ge 0 $. In addition, the operator $ \left[\cdot\right]_+ $ is defined for real scalars as $ \left[\alpha\right]_+ \triangleq \max\left\{\alpha, 0\right\} $, hence, $ R \left( D \right) = 0 $ for $ D \ge \sigma^2 $. Accordingly, the rate-distortion function of the Gaussian variable ${\tilde{y}^F_k}$ is

[rCl]{} \[eq:theoretic Gaussian analysis - rate-distortion function of DFT component\] R\_k(D\_k) = \_+

where $ D_k $ denotes the maximal squared-error allowed for this component. Now, similar to the famous case of jointly coding independent Gaussian variables with respect to a regular (non-weighted) squared-error distortion [@cover2012elements], we explicitly express Problem \[problem: Separable Form in DFT domain\] as the following distortion-allocation optimization.

\[problem: Gaussian distortion allocation\]

[rCl]{} \[eq:theoretic Gaussian analysis - Gaussian distortion allocation\]

& & & \_[k=0]{}\^[N-1]{} [ \_+ ]{}\
& & & \_[k=0]{}\^[N-1]{} [ | h\^F\_k |\^2 D\_k ]{} = N\_[H]{} \_n\^2\
& & & D\_k 0    , k=0,...,N-1.

The optimal distortion-allocation satisfying the last optimization is

[rCl]{} \[eq:theoretic Gaussian analysis - optimal distortion allocation\] D\_k\^[opt]{} = {

[\*[20]{}[c]{}]{} &[, h\^F\_k0 ]{}\
0&[, h\^F\_k = 0 ]{}

.

and the associated optimal rates are

[rCl]{} \[eq:theoretic Gaussian analysis - optimal rate allocation\] R\_k\^[opt]{} = {

[\*[20]{}[c]{}]{} [ ( | h\^F\_k |\^2 + 1 )]{}&[, h\^F\_k0 ]{}\
0&[, h\^F\_k = 0 . ]{}

.

Results (\[eq:theoretic Gaussian analysis - optimal distortion allocation\]) and (\[eq:theoretic Gaussian analysis - optimal rate allocation\]) are proved in Appendix \[appendix:Solution of Problem 4\].

Demonstration of The Explicit Results {#subsec:Rate-Distortion Theoretic - Demonstration of The Explicit Results}
-------------------------------------

Let us exemplify the optimal rate-distortion results (\[eq:theoretic Gaussian analysis - optimal distortion allocation\])-(\[eq:theoretic Gaussian analysis - optimal rate allocation\]) for a cyclo-stationary Gaussian signal, $ {\mathbf{x}} $, having the circulant autocorrelation matrix presented in Fig. \[fig:signal\_autocorrelation\_matrix\], corresponding to the eigenvalues $ \{ {\lambda}^{({\mathbf{x}})}_k \} _{k=0}^{N-1} $ (Fig. \[fig:signal\_autocorrelation\_matrix\_DFT\_decomposition\]) obtained by a DFT-based decomposition. We first examine the denoising problem, where the signal-domain degradation matrix is $ {\mathbf{H}} = {\mathbf{I}} $ (Fig. \[fig:gaussian\_denoising\_\_degradation\_filter\_signal\_domain\]) and its respective DFT-domain spectral representation consists of $ h^F_k = 1 $ for any $ k $ (see Fig. \[fig:gaussian\_denoising\_\_degradation\_filter\_DFT\_coefficient\_magnitudes\]). The additive white Gaussian noise has a sample variance of $ \sigma_n^2 = 5 $. Fig. \[fig:gaussian\_denoising\_\_effective\_waterfilling\_demonstration\] exhibits the optimal distortion allocation using a reverse-waterfilling diagram, where the signal-energy distribution $ \{ {\lambda}^{({\mathbf{x}})}_k \} _{k=0}^{N-1} $ (black solid line) and the additive noise energy (the light-red region) defining together the noisy-signal energy level (purple solid line) corresponding to $\lambda^{\left(\tilde{{\mathbf{y}}}\right)}_k = \lambda^{({\mathbf{x}})}_k + \sigma_n^2$. The blue dashed line in Fig. \[fig:gaussian\_denoising\_\_effective\_waterfilling\_demonstration\] shows the water level associated with the uniform distortion allocation. The optimal rate-allocation, corresponding to Fig. \[fig:gaussian\_denoising\_\_effective\_waterfilling\_demonstration\] and Eq. (\[eq:theoretic Gaussian analysis - optimal rate allocation\]), is presented in Fig. \[fig:gaussian\_denoising\_\_optimal\_rate\_allocation\] showing that more bits are spent on components with higher signal-to-noise ratios.

Another example considers the same Gaussian signal described in Fig. \[Fig:Gaussian signal for demonstration\] and the noise level of $ \sigma_n^2 = 5 $, but here the degradation operator is the circulant matrix shown in Fig. \[fig:gaussian\_restoration\_\_degradation\_filter\_signal\_domain\] having a DFT-domain representation given in magnitude-levels in Fig. \[fig:gaussian\_restoration\_\_degradation\_filter\_DFT\_coefficient\_magnitudes\] exhibiting its frequency attenuation and amplification effects. The waterfilling diagram in Fig. \[fig:gaussian\_restoration\_\_effective\_waterfilling\_demonstration\] includes the same level of signal energy (black solid line) as in the denoising experiment, but the effective additive noise levels and the allocated distortions are clearly modulated in an inversely proportional manner by the squared magnitude of the degradation operator. For instance, frequencies corresponding to degradation-filter magnitudes lower than 1 lead to increase in the effective noise-energy addition and in the allocated distortion. The optimal rate allocation (Fig. \[fig:gaussian\_restoration\_\_optimal\_rate\_allocation\]) is affected by the signal-to-noise ratio and by the squared-magnitude of the degradation filter (see also Eq. (\[eq:theoretic Gaussian analysis - optimal rate allocation\])), e.g., components that are attenuated by the degradation operator get less bits in the rate allocation.

Conceptual Relation to The Proposed Approach {#subsec:Rate-Distortion Theoretic - Conceptual Relation to The Proposed Approach}
--------------------------------------------

As explained at the beginning of this section, theoretical and practical solutions may include different implementations of the same general ideas. Accordingly, connections between theory and practice should be established by pointing on high-level analogies. Our rate-distortion analysis (for a Gaussian signal and a LSI degradation operator) showed that the optimal complexity-regularized restoration relies on two prominent ideas: pseudoinverse filtering of the degraded input, and compression with respect to a squared-error metric that is weighted based on the degradation-filter squared-magnitude (considering the DFT-domain procedure). We will now turn to explain how these two concepts connect to more general themes having different realizations in the practical approach proposed in Section \[sec:Proposed Methods\] [^2].

$\bullet$ ***Design Concept \#1:** Apply simple restoration filtering*. The general idea of using an elementary restoration filter is implemented in the Gaussian case as pseudoinverse filtering. Correspondingly, our practical approach relies on a simple filtering mechanism, extending the pseudoinverse filter as explained next. Stage 6 of Algorithm \[Algorithm:Proposed Method Non-overlapping\] is an $ \ell_2 $-constrained deconvolution filtering that its analytic solution can be rewritten, using the relation $ {\mathbf{H}}^*\left( {\mathbf{I}} - {\mathbf{H}}{\mathbf{H^*}} \right) \nolinebreak=\nolinebreak {\mathbf{0}} $, as (see proof in Appendix \[appendix:Equivalent Forms of Stage 4 of Algorithm 1\])

[rCl]{} \[eq:theory-practice relation - interpreting first stage\] \^[(t)]{} = ( \^\* + )\^[-1]{} ( \^\* + \^[(t)]{} ) . 

As before, $ \tilde{{\mathbf{y}}} = {\mathbf{H}}^+ {\mathbf{y}} $, i.e., the pseudoinverse-filtered version of $ {\mathbf{y}} $. The expression (\[eq:theory-practice relation - interpreting first stage\]) can be interpreted as an initial pseudoinverse filtering of the degraded input, followed by a simple weighted averaging with $\tilde{{\mathbf{z}}}^{(t)}$ (that includes the decompressed signal obtained in the last iteration). Evidently, the filtering in (\[eq:theory-practice relation - interpreting first stage\]) is determined by the $\beta$ value, specifically, for $ \beta=0 $ the estimate coincides with the pseudoinverse filtering solution and for a larger $\beta$ it is closer to $\tilde{{\mathbf{z}}}^{(t)}$.

$\bullet$ ***Design Concept \#2:** Compress by promoting higher quality for signal-components matching to higher $ h $-operator magnitudes*. This principle is realized in the theoretic Gaussian case as weights attached to the squared-errors of DFT-domain components (see Problems \[problem: Separable Form in DFT domain\] and \[problem: Gaussian distortion allocation\]). Since the weights, ($ \left\lvert h_0^F \right\rvert ^2 , ..., \left\lvert h_{N-1}^F \right\rvert ^2 $), are the squared magnitudes of the corresponding degradation-filter coefficients, in the compression of the pseudoinverse-filtered input the distortion is spread unevenly being larger where the degradation filter-magnitude is lower. Remarkably, this concept is implemented differently in the proposed procedure (Algorithm \[Algorithm:Proposed Method Non-overlapping\]) where regular compression techniques, optimized for the squared-error distortion measure, are applied on the filtering result of the preceding stage. We will consider the essence of the effective compression corresponding to these two stages together. Let us revisit (\[eq:theory-practice relation - interpreting first stage\]), expressing stage 6 of Algorithm \[Algorithm:Proposed Method Non-overlapping\]. Assuming $ {\mathbf{H}} $ is a circulant matrix, we can transform (\[eq:theory-practice relation - interpreting first stage\]) into its Fourier domain representation

[rCl]{} \[eq:theory-practice relation - interpreting first stage - circulant H\] \^[F,(t)]{} = ( \_H\^\* \_H + )\^[-1]{} ( \_H\^\* \_H \^[F]{} + \^[F,(t)]{} )\

where $ \hat{ {\mathbf{x}}}^{F,(t)} $ and $ \tilde{{\mathbf{z}}}^{F,(t)} $ are the Fourier representations of $ \hat{ {\mathbf{x}}}^{(t)} $ and $ \tilde{{\mathbf{z}}}^{(t)} $, respectively. Furthermore, (\[eq:theory-practice relation - interpreting first stage - circulant H\]) reduces to the componentwise formulation

[rCl]{} \[eq:theory-practice relation - interpreting first stage - circulant H - component\] \^[F,(t)]{}\_k =

where $ \hat{ {x}}^{F,(t)}_k $ and $ \tilde{z}^{F,(t)}_k $ are the $ k^{th} $ Fourier coefficients of $ \hat{ {\mathbf{x}}}^{F,(t)} $ and $ \tilde{{\mathbf{z}}}^{F,(t)} $, respectively. Equation (\[eq:theory-practice relation - interpreting first stage - circulant H - component\]) shows that signal elements (of the pseudoinverse-filtered input) corresponding to degradation-filter components of weaker energies will be retracted more closely to the respective components of $ \tilde{z}^{F,(t)}_k $ – thus, will be farther from $\tilde{y}^{F}_k$, yielding that the corresponding components in the standard compression applied in the next stage of this iteration will be of a relatively lower quality with respect their matching components of $\tilde{y}^{F}_k$ (as they were already retracted relatively far from them in the preceding deconvolution stage).

To conclude this section, we showed that the main architectural ideas expressed in theory (for the Gaussian case) appear also in our practical procedure. The iterative nature of our methods (Algorithms \[Algorithm:Proposed Method Non-overlapping\]-\[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\]) as well as the desired shift-invariance property provided by Algorithms \[Algorithm:Proposed Method Overlapping Blocks\]-\[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\] are outcomes of treating real-world scenarios such as non-Gaussian signals, general linear degradation operators, and computational limitations leading to block-based treatments – these all relate to practical aspects, hence, do not affect the fundamental treatment given in this section.

Experimental Results {#sec:Experimental Results}
====================

In this section we present experimental results for image restoration. Our main study cases include deblurring and inpainting using the image-compression profile of the HEVC standard (in its BPG implementation [@hevc_software_bpg]). We also provide evaluation of our method in conjunction with the JPEG2000 technique for the task of image deblurring.

We empirically found it sufficient to consider only a part of all the shifts, i.e., a portion of $ \mathcal{B}^* $. When using HEVC, the limited amount of shifts is compensated by the compression architecture that employs inter-block spatial predictions, thus, improves upon methods relying on independent block treatment. The shifts are defined by the rectangular images having their upper-left corner pixel relatively close to the upper-left corner of the full image, and their bottom-right corner pixel coincides with that of the full image. This extends the mathematical developments in Section \[sec:Proposed Methods\] as practical compression handles arbitrarily sized rectangular images.

Many image regularizers have visual interpretation, for example, the classical image-smoothness evaluation. In our framework, the regularization part in (\[eq:complexity-regularized restoration - full signal\]) measures the complexity in terms of the compression bit-cost with respect to a specific compression architecture, designed based on some image model. Our complexity regularization also has a general visual meaning since, commonly, low bit-cost compressed images tend to be overly smooth or piecewise-smooth.

Image Deblurring {#subsec:Deblurring}
----------------

Here we consider two deterioration settings taken from [@danielyan2012bm3d]. The first setting, denoted here also as ‘Set. 1’, considers a noise variance $ \sigma_n^2 = 2 $ coupled with a blur operator defined by the two-dimensional point-spread-function (PSF) $ h(x_1,x_2)=1/{(1+x_1^2+x_2^2)} $ for $ x_1,x_2=-7,...,7 $, and zero-valued otherwise. The second setting, denoted here also as ‘Set. 2’ (named in [@danielyan2012bm3d] as ‘Scenario 3’), considers a noise variance $ \sigma_n^2 \approx 0.3 $ joint with a blur operator defined by the two-dimensional uniform blur PSF of size $ 9\times 9 $.

We precede the deblurring experiments with empirical evaluations of four important aspects of the proposed method.

### Iterative Reduction of the Fundamental Restoration Cost

In Section \[sec:Proposed Methods\] we established the basic optimization problems for restoration by regularizing the bit-costs of the non-overlapping and the overlapping blocks of the estimate (see (\[eq:complexity-regularized restoration - expressing blocks\]) and (\[eq:complexity-regularized restoration - overlapping blocks\]), respectively). As explained above, these two fundamental optimization tasks cannot be directly addressed and, therefore, we developed the ADMM-based Algorithms \[Algorithm:Proposed Method Non-overlapping\]-\[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\], that iteratively employ simpler optimization problems. Figures \[fig:empirical\_analysis\_cost\_evolution\_deblurring\_shift1\_cameraman\] and \[fig:empirical\_analysis\_cost\_evolution\_deblurring\_shift3\_cameraman\] demonstrate that, for appropriate parameter settings, the fundamental optimization cost reduces in each iteration. The provided figures also show the fidelity term, $ \left\| {  {\mathbf{H}} {\mathbf{x}} - {\mathbf{y}} } \right\|_2^2 $, and the regularizing bit-cost (multiplied by $ \mu $) of each iteration.

### Iterative Improvement of the Restored Image

The fundamental optimization costs in (\[eq:complexity-regularized restoration - expressing blocks\]) and (\[eq:complexity-regularized restoration - overlapping blocks\]) include the fidelity term $ \left\| {  {\mathbf{H}} {\mathbf{x}} - {\mathbf{y}} } \right\|_2^2 $ that considers the candidate estimate ${\mathbf{x}}$ and the given degraded signal $ {\mathbf{y}} $. However, the ultimate goal of the restoration process is to produce an estimate ${\mathbf{x}}$ that will be close to the original (unknown!) signal ${\mathbf{x}}_0$. It is common to evaluate proposed methods in experiment settings where ${\mathbf{x}}_0$ is known and used only for the evaluation of the squared error $ \left\| {  {\mathbf{x}} - {\mathbf{x}}_0 } \right\|_2^2 $ or its corresponding PSNR. Accordingly, it is a desired property that our iterative methods will provide increment in the PSNR along the iterations and, indeed, Figures \[fig:empirical\_analysis\_psnr\_evolution\_deblurring\_shift1\_cameraman\] and \[fig:empirical\_analysis\_psnr\_evolution\_deblurring\_shift3\_cameraman\] show that this is achievable for appropriate parameter settings (the use of improper parameters may lead to unwanted decrease of the PSNR starting at some unknown iteration that, however, can be detected in many cases by heuristic divergence rules based on the dual variables used in the ADMM process). Interestingly, for some parameter settings, the PSNR may increase with the iterations, whereas the fundamental restoration cost will not necessarily consistently decrease. The last behavior may result from the fact the our optimization problem (with respect to a standard compression technique) is discrete, non-linear, and usually not convex and, therefore, the convergence guarantees of the ADMM [@boyd2011distributed] do not hold here for the fundamental restoration cost.

### The Optimal Compression Parameter

Another question of practical importance is the value of the parameter $ \theta $, determining the compression level of the standard technique utilized in the proposed Algorithms \[Algorithm:Proposed Method Non-overlapping\]-\[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\]. Recall that the ADMM-based developments in Section \[sec:Proposed Methods\] led to an iterative procedure including a stage of Lagrangian rate-distortion optimization operated for a Lagrange multiplier $ \lambda \triangleq \frac{2\mu}{\beta} $ and, then, we replaced this optimization with application of a standard compression-decompression technique operated based on a parameter $ \theta $. It is clear that $ \theta $ is a function of $ \lambda $. In the particular case where the standard compression has the Lagrangian form from our developments, then $ \theta = \lambda $, however, this is not the general case. For an arbitrary compression technique, we assume that its parameter $ \theta $ has $ K $ possible values $ \theta_1,...,\theta_K $ (for example, the HEVC standard supports 52 values for its quantization parameter), then, for a given $ \lambda \triangleq \frac{2\mu}{\beta} $ the required $ \theta $ value in stage 8 of Algorithm \[Algorithm:Proposed Method Non-overlapping\] can be determined via

[rCl]{} \[eq:experiments - general setting of theta - rate-distortion optimization\] \^[(t)]{}\_[,opt]{} = \_[\_1,...,\_K ]{} [ \_2\^2]{} + r\_[tot,]{} .

where $  \hat{{\mathbf{z}}}^{(t)}_{\theta } = {CompressDecompress}_{\theta}\left( \tilde{{\mathbf{x}}}^{(t)} \right) $ is the decompressed signal and $ r_{tot,\theta} $ is the associated compression bit-cost. We present here experiments (see Figs. \[Fig:empirical\_analysis\_for\_deblurring\_shift1\_cameraman\] and \[Fig:empirical\_analysis\_for\_deblurring\_shift3\_cameraman\]) for Algorithm \[Algorithm:Proposed Method Non-overlapping\] and \[Algorithm:Proposed Method Overlapping Blocks\] that in each iteration optimize the $ \theta $ value corresponding to $ \lambda \triangleq \frac{2\mu}{\beta} $ based on procedures similar to (\[eq:experiments - general setting of theta - rate-distortion optimization\]). Nicely, it is shown in Figures \[fig:empirical\_analysis\_compression\_factor\_evolution\_deblurring\_shift1\_cameraman\] and \[fig:empirical\_analysis\_compression\_factor\_evolution\_deblurring\_shift3\_cameraman\] that, for the HEVC compression used here, the best $ \theta $ values along the iterations are nearly the same (for a specific restoration task). This important property may be a result of the fact that HEVC extensively relies on Lagrangian rate-distortion optimizations (although in much more complex forms than those presented in Section \[sec:Proposed Methods\]). Accordingly, in order to reduce the computational load, in the experiments shown below we will use a constant compression parameter given as an input to our methods. Interestingly, when examining the optimal compression parameters (compression ratios in this case) for the JPEG2000 method that applies wavelet-based transform coding, there is a decrease in the optimal compression ratio along the iterations (see Fig. \[fig:empirical\_analysis\_deblurring\_ADMM\_JPEG2000\_compression\_ratios\_cameraman\]). Accordingly, in order to reduce the computational load in the experiments below, we employed a predefined rule for reducing the JPEG2000 compression ratio along the iterations. Importantly, when we use the sub-optimal predefined rules for setting the compression parameter $ \theta $ values (see Table \[table:Experimental Results - Parameter Setting\] for the settings used in our evaluation comparisons in Table \[table:Deblurring - PSNR Comparison\]) we do not longer need to set a value for $ \mu $ (since, in these cases, $ \mu $ is practically unused).

\[\]

\[table:Experimental Results - Parameter Setting\]

\[\]

  --------------------------------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- -----------
                                                                                                                        
                                      Set. 1      Set. 2      Set. 1      Set. 2      Set. 1      Set. 2      Set. 1      Set. 2
             Input PSNR                22.23       20.76       25.61       24.11       27.25       25.84       23.34       22.49
   ForWaRD [@neelamani2004forward]     28.99       28.10       32.96       33.67       33.30       32.81       27.03       26.51
     SV-GSM [@guerrero2008image]       29.68       28.09       34.25       33.15         -           -         30.19       27.56
      BM3DDEB [@dabov2008image]      **30.42**     29.10     **34.93**     34.96     **35.20**     33.81     **31.14**   **28.35**
    TVMM [@oliveira2009adaptive]       29.64       29.30       33.59       34.50       33.61       33.31       26.44       25.98
   CGMK [@chantas2010variational]      30.03       29.91       33.92       34.86       34.01       33.70       25.79       26.04
    IDD-BM3D [@danielyan2012bm3d]    **31.08**   **31.21**   **35.56**   **37.00**   **35.22**   **34.75**   **30.98**   **28.54**
      EPLL [@zoran2011learning]        29.40       29.54       33.88       35.87     **34.64**     34.39       27.62       26.78
                                       28.99       27.48       33.55       33.89       34.11       32.56       26.42       25.14
                                       29.28       28.10       33.07       34.01       34.38       32.60       26.31       25.18
                                       30.19     **30.20**     33.17     **36.47**     33.32     **34.40**     29.83     **27.74**
                                     **30.35**   **30.14**   **34.37**   **36.57**     34.55     **34.41**   **30.20**     27.72
  --------------------------------- ----------- ----------- ----------- ----------- ----------- ----------- ----------- -----------

### Restoration Improvement for Increased Number of Image Shifts

In our experiments we noticed that the restoration quality improves as more shifts are used, however, at some point the added gain due to the added shifts becomes marginal. As an example to the benefits due to shifts see Figs. \[fig:empirical\_analysis\_psnr\_evolution\_deblurring\_shift1\_cameraman\] and \[fig:empirical\_analysis\_psnr\_evolution\_deblurring\_shift3\_cameraman\], where the PSNR obtained for deblurring the Cameraman image using Algorithm 2 and 9 shifts is about 2dB higher than the PSNR obtained using Algorithm 1 (i.e., without additional shifts).

In our main evaluation, we examined the proposed Algorithms 2 and 3 for image deblurring in conjunction with the JPEG2000 and HEVC compression techniques (see the parameter settings in Table \[table:Experimental Results - Parameter Setting\]). Table \[table:Deblurring - PSNR Comparison\] shows a comparison between various deblurring methods tested in the above two settings for four grayscale images[^3]. In 7 out of the 8 cases, the proposed Algorithms 2 or 3 utilizing the HEVC standard provided one of the best three results. Visual results are presented in Figures \[Fig:Deblurring results - cameraman\] and \[Fig:Deblurring results - barbara\].

Image Inpainting {#subsec:Inpainting}
----------------

We presented in [@dar2016image] experimental results for the inpainting problem, in its noisy and noiseless settings. Here we focus on the noiseless inpainting problem, where only pixel erasure occurs without an additive noise. The degradation is represented by a diagonal matrix $ {\mathbf{H}} $ of $ N\times N $ size with main diagonal values of zeros and ones, indicating positions of missing and available pixels, respectively. Then, the product $ {\mathbf{H}}{\mathbf{x}} $ equals to an $ N $-length vector where its $ k^{th} $ sample is determined by $ {\mathbf{H}} $: if $ {\mathbf{H}}[k,k] = 0$ then it is zero, and for $ {\mathbf{H}}[k,k] = 1$ it equals to the corresponding sample of $ {\mathbf{x}} $. The structure of the pixel erasure operator let us to simplify the optimization in step 6 of Algorithm \[Algorithm:Proposed Method Overlapping Blocks\]. We note that $ {\mathbf{H}} $ is a square diagonal matrix and, therefore, $ {\mathbf{H}}^T = {\mathbf{H}} $ and $ {\mathbf{H}}^T {\mathbf{y}} $ is equivalent to a vector $ {\mathbf{y}} $ with zeroed components according to $ {\mathbf{H}} $’s structure. Additional useful relation is $ {\mathbf{H}}^T {\mathbf{H}} = {\mathbf{H}} $. Consequently, step 6 of Algorithm \[Algorithm:Proposed Method Overlapping Blocks\] facilitates a componentwise computation that is interpreted to form the $ k^{th} $ sample of $  \hat{{\mathbf{x}}}^{(t)}  $ as

[rCl]{} \[eq:noiseless inpainting - k-th component of x\] \^[(t)]{} \[k\] = {

[ll]{}        \[k,k\]=1\
\_[j=1]{}\^[N\_b]{} \^[j,(t)]{}\[k\]            \[k,k\]=0

   .

We initialize the shifted images $ \left\lbrace {\hat{{\mathbf{z}}}}^{j,(0)} \right\rbrace _{j=1}^{N_b} $ as the given image with the missing pixels set as the corresponding local averages of the available pixels in the respective $ 7\times 7 $ neighborhoods. When the iterative processing ends, we use the fact that the available pixels are noiseless and set them in the reconstructed image. The rest of the procedure remains as before.

We present here implementations of Algorithms \[Algorithm:Proposed Method Overlapping Blocks\] and \[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\] utilizing the JPEG2000 and HEVC image compression (the parameter settings are described in Table \[table:Experimental Results - Parameter Setting\]). We consider the experimental settings from [@ram2013image], where 80% of the pixels are missing (see Fig. \[fig:inpainting - barbara - degraded image\] and \[fig:inpainting - house - degraded image\]). Five competing inpainting methods are considered: cubic interpolation of missing pixels via Delaunay triangulation (using Matlab’s ’griddata’ function); inpainting using sparse representations of patches of $ 16\times 16 $ pixels based on an overcomplete DCT (ODCT) dictionary (see method description in [@elad2010sparse Ch. 15]); using patch-group transformation [@li2008patch]; based on patch clustering [@yu2012solving]; and via patch reordering [@ram2013image]. The PSNR values of images restored using the above methods (taken from [@ram2013image]) are provided in Table \[table:Experimental Results - Inpainting\] together with our results. For two images our HEVC-based implementation of Algorithm 3 provides the highest PSNR values. Visually, Figures \[fig:inpainting - barbara - proposed restoration\] and \[fig:inpainting - house - proposed restoration\] exhibit the effectiveness of our method in repairing the vast amount of absent pixels.

\[\]

\[table:Experimental Results - Inpainting\]

   Image           ODCT                                                       
  ------- ------- ------- ------- ----------- ------- ------- ------- ------- -----------
           30.25   29.97   31.62     32.22     31.96   30.31   30.96   32.14   **32.55**
           22.88   27.15   25.40   **30.94**   29.71   24.25   24.83   26.06     28.80
           29.21   29.69   32.87     33.05     32.71   29.49   30.50   32.42   **33.10**

Conclusion {#sec:Conclusion}
==========

In this paper we explored the topic of complexity-regularized restoration, where the likelihood of candidate estimates are determined by their compression bit-costs. Using the alternating direction method of multipliers (ADMM) approach we developed three practical methods for restoration using standard compression techniques. Two of the proposed methods rely on a new shift-invariant complexity regularizer, evaluating the total bit-cost of the signal shifted versions. We explained few of the main ideas of our approach using an insightful theoretical-analysis of complexity-regularized restoration of a cyclo-stationary Gaussian signal from deterioration of a linear shift-invariant operator and additive white Gaussian noise. Experiments for deblurring and inpainting of images using the JPEG2000 and HEVC technique showed good results.

Proofs for the theory section
=============================

Equivalence of Problems \[problem: basic form\] and \[problem: pseudoinverse filtered input\] {#appendix:Equivalence of Problems 1 and 2}
---------------------------------------------------------------------------------------------

We start by showing the equality between the distortion constraints of Problems \[problem: basic form\] and \[problem: pseudoinverse filtered input\]. We develop the distortion of Problem \[problem: basic form\] as follows:

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - distortions\] \_2\^2 & = & \_2\^2\
& = & \_2\^2\
& = & \_2\^2 + \_2\^2\
&& + ( \^[+]{} - )\^[\*]{} \^[\*]{} ( - \^[+]{} )\
&& + \^[\*]{} ( - \^[+]{} )\^[\*]{} ( \^[+]{} - )\
& = & \_2\^2 + \_2\^2     

where the last equality follows from

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - distortions - auxiliary result\] \^[\*]{} ( - \^[+]{} ) =

that can be easily proved, e.g., by using the DFT-based diagonalization of $ {\mathbf{H}} $ and $ {\mathbf{H}}^{+} $.

The first term in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - distortions\]) can be further developed:

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - distortions - first part - deterministic\] \_2\^2 & = & \_2\^2\
& = & \_2\^2\
& = & \_2\^2\
& = & \_[k: h\_k\^F = 0]{} y\_k\^F \^2 =\
& = & \_[k: h\_k\^F = 0]{} n\_k\^F \^2

where $ {\mathbf{y}}^{F} \triangleq {\mathbf{F}} {\mathbf{y}} $ is the DFT-domain representation of $ {\mathbf{y}} $ (correspondingly, we use these notations to any vector), and the last equality is implied from the DFT-component relation $ y_k^F = h_k^F x_k^F + n_k^F $ that reduces to $ y_k^F = n_k^F $ for components with $ h_k^F = 0 $. Consequently,

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - distortions - first part\] E \_2\^2 & = & E\_[k: h\_k\^F = 0]{} n\_k\^F \^2\
& = & (N - N\_[H]{}) \_n\^2

where $ N_{H} $ was defined in Section \[subsec:Reformulations of the Problem\] as the rank of $ {\mathbf{H}} $. Accordingly, and also using (\[eq:appendix - proofs - equivalence of problem 1 and 2 - distortions\]), the distortion constraint of Problem \[problem: basic form\], i.e.,

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - distortion of problem 1\] E \_2\^2 = N \_n\^2

equals to (recall that $ \tilde{{\mathbf{y}}} = {\mathbf{H}}^{+} {\mathbf{y}} $)

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - distortion of problem 2\] E ( [ - ]{} ) \_2\^2 = N\_[H]{} \_n\^2 ,

that is, the distortion constraint of Problem \[problem: pseudoinverse filtered input\].

We now turn to prove the equivalence of Problems \[problem: basic form\] and \[problem: pseudoinverse filtered input\]. Our proof sketch conforms with common arguments in rate-distortion function proofs (see [@cover2012elements]): first, we lower bound the mutual information $ I\left( {\mathbf{y}}; \hat{{\mathbf{x}}} \right) $, which is the cost function of Problem \[problem: basic form\]; then, we provide a statistical construction achieving the lower bound while obeying the distortion constraint.

The proposed lower bound for $ I\left( {\mathbf{y}}; \hat{{\mathbf{x}}} \right) $ is established by noting that $ \tilde{{\mathbf{y}}} = {\mathbf{H}}^{+} {{\mathbf{y}}} $ and, therefore, the data processing inequality [@cover2012elements] implies here that

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - data processing inequality\] I( ; ) I( ; ) ,

where $ I\left( \tilde{{\mathbf{y}}} ; \hat{{\mathbf{x}}} \right) $ is the cost function of Problem \[problem: pseudoinverse filtered input\]. The relation in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - data processing inequality\]) is known to be attained with equality when $ {\mathbf{y}} $ and $\hat{{\mathbf{x}}} $ are independent given $ \tilde{{\mathbf{y}}} $. The next construction shows that this is indeed the case.

We will now show the achievability of the lower bound in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - data processing inequality\]) by describing a two-stage setting that statistically represents $ \tilde{{\mathbf{y}}} $ as an outcome of $ \hat{{\mathbf{x}}} $, and $ {{\mathbf{y}}} $ as a consequence of $ \tilde{{\mathbf{y}}} $. This layout is an instance of the construction concept known as the (backward) test channel [@cover2012elements]. The first stage of our construction is based on

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - x \] & \~& ( , \^[+]{} \_ \^[+\*]{} - \_n\^2 \^[+]{}\^[+\*]{} )\
\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - z \] & \~& ( , \_n\^2 \^[+]{}\^[+\*]{} )

where $ \hat{{\mathbf{x}}} $ and $ {\mathbf{z}} $ are independent. Consequently, we define

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - y\_tilde as a sum \] = + ,

implying $ \tilde{{\mathbf{y}}}  \sim \mathcal{N} \left( {\mathbf{0}}, {\mathbf{H}}^{+} {\mathbf{R}}_{{\mathbf{y}}} {\mathbf{H}}^{+*} \right) $ that, indeed, conforms with $ \tilde{{\mathbf{y}}} = {\mathbf{H}}^{+} {\mathbf{y}} $ where ${\mathbf{y}} \sim \mathcal{N} \left( {\mathbf{0}}, {\mathbf{R}}_{{\mathbf{y}}} \right)$. Moreover, the construction (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - x \])-(\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - y\_tilde as a sum \]) yields

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - distortion constraint satisfies\] E ( [ - ]{} ) \_2\^2 & = & E \_2\^2\
& = & E\^[\*]{} \^[\*]{}\
& = & ETrace\^[\*]{} \^[\*]{}\
& = & ETrace \^[\*]{} \^[\*]{}\
& = & Trace \_ \^[\*]{}\
& = & \_n\^2 Trace \^[+]{}\^[+\*]{} \^[\*]{}\
& = & \_n\^2 Trace\_[[H]{}]{} \_[[H]{}]{}\^[\*]{}\
& = & \_n\^2 N\_H

where $ {\mathbf{P}}_{{H}} \triangleq {\mathbf{H}}  {\mathbf{H}}^{+}  $ is the matrix projecting onto the range of $ {\mathbf{H}} $, note it is also a circulant matrix diagonalized by the DFT matrix to the diagonal matrix $ {\mathbf{\Lambda}}_{{\mathbf{P}}_{{H}}} \triangleq {\mathbf{\Lambda}}_{H} {\mathbf{\Lambda}}_{H}^{+}  $. The last computation of the trace is due to the structure of ${\mathbf{\Lambda}}_{{\mathbf{P}}_{H}}$, having ones in main-diagonal entries corresponding to the DFT-domain indices of the range of $ {\mathbf{H}} $, and zeros elsewhere. The result in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - distortion constraint satisfies\]) shows that the distortion constraint (\[eq:appendix - proofs - equivalence of problem 1 and 2 - distortion of problem 2\]) is satisfied.

Let us consider the second stage of the construction, awaiting to prove that $ {\mathbf{y}} $ and $\hat{{\mathbf{x}}} $ are independent given $ \tilde{{\mathbf{y}}} $. We precede the construction with examining the following decomposition of $ {\mathbf{y}} $

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - decomposition of y\] & = & \_[H]{} + ( - \_[H]{} )\
& = & \^[+]{} + ( - \_[H]{} ) ( + )\
& = & + ( - \_[H]{} )

where the second equality uses the degradation model, and the third equality is due to $ \left( {\mathbf{I}} - {\mathbf{P}}_{H} \right) {\mathbf{H}} = {\mathbf{0}} $. Importantly, Eq. (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - decomposition of y\]) describes $ {\mathbf{y}} $ as a linear combination of two independent random vectors: $ \tilde{{\mathbf{y}}} $ and $ \left( {\mathbf{I}} - {\mathbf{P}}_{H} \right) {\mathbf{n}} $. Since $ \tilde{{\mathbf{y}}} $ and $ \left( {\mathbf{I}} - {\mathbf{P}}_{H} \right) {\mathbf{n}} $ are Gaussian random vectors, their independence is proved by showing they are uncorrelated via

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - uncorrelated backward variables - new \] && E( - \_[H]{} ) \^[\*]{}\
&& = E( - \_[H]{} ) \^[\*]{} \^[+\*]{}\
&& = E( - \^[\*]{} \_[\_[H]{}]{} ) \^[\*]{} \^[\*]{} \_[H]{}\^[+\*]{}\
&& = E\^[\*]{} ( - \_[\_[H]{}]{} ) ( )\^[\*]{} \_[H]{}\^[+\*]{}\
\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - uncorrelated backward variables - new - explain\] && = \^[\*]{} E( - \_[\_[H]{}]{} ) \^[F]{} ( \_[H]{}\^[+]{} \^[F]{} )\^[\*]{}\
&& = \^[\*]{}\
&& =

where in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - uncorrelated backward variables - new - explain\]) we used the facts that $\left( {\mathbf{I}} -  {\mathbf{\Lambda}}_{{\mathbf{P}}_{H}}  \right) {\mathbf{n}}^{F}$ is a DFT-domain vector with zeros in components corresponding to the range of $ {\mathbf{H}} $, and $ {\mathbf{\Lambda}}_{H}^{+} {\mathbf{y}}^{F} $ is a DFT-domain vector with zeros in entries corresponding to the nullspace of $ {\mathbf{H}} $, hence, these zero patterns yield the outer-product matrix which is all zeros.

The decomposition in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - decomposition of y\]) motivates us to consider $ {\mathbf{y}} $ to emerge from $ \tilde{{\mathbf{y}}} $ via the statistical relation

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - y from y\_tilde \] = +

where $ {\mathbf{w}} \sim \mathcal{N} \left( {\mathbf{0}}, \sigma_n^2  \left({\mathbf{I}} - {\mathbf{P}}_{H}  \right) \left({\mathbf{I}} - {\mathbf{P}}_{H}  \right)^{*} \right) $ and is independent of $ \tilde{{\mathbf{y}}} $ (and also of $ {\mathbf{z}} $). Note that $ {\mathbf{w}} $ takes the role of $ \left( {\mathbf{I}} - {\mathbf{P}}_{H} \right) {\mathbf{n}} $ appearing in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - decomposition of y\]), e.g., they have the same distribution. One can further examine the suitability of the construction (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - y from y\_tilde \]) to the considered problem by noting it satisfies the distortion constraint of Problem \[problem: basic form\], namely,

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - distortion constraint satisfies - Problem 2\] E \_2\^2 & = & E \_2\^2\
& = & E \_2\^2\
& = & E \_2\^2\
& = & E \_2\^2 + E \_2\^2\
& = & \_n\^2 Trace \^[+]{}\^[+\*]{} \^[\*]{}\
&& + \_n\^2 Trace( - \_[H]{} ) ( - \_[H]{} )\^[\*]{}\
& = & \_n\^2 Trace\_[H]{} \_[H]{}\^[\*]{}\
&& + \_n\^2 Trace( - \_[H]{} ) ( - \_[H]{} )\^[\*]{}\
& = & \_n\^2 N\_H + \_n\^2 ( N - N\_H )\
& = & N \_n\^2

as required. Furthermore, the constructed $ {\mathbf{y}} $ indeed obeys ${\mathbf{y}} \nolinebreak \sim \nolinebreak \mathcal{N} \left( {\mathbf{0}}, {\mathbf{R}}_{{\mathbf{y}}} \right)$. Specifically, its autocorrelation matrix stems from the calculation

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - Problem 2 - y autocorrelation\] E \^[\*]{} & = & E( [ + ]{} ) ( [ + ]{} )\^[\*]{}\
& = & E( [ ]{} ) ( [ ]{} )\^[\*]{} + E\^[\*]{}\
& = & \^[+]{} \_ \^[+\*]{} \^[\*]{} + \_\
& = & \_[H]{} ( \_ \^[\*]{} + \_n\^2 ) \_[H]{}\^[\*]{} + \_\
& = & \_[H]{} \_ \^[\*]{} \_[H]{}\^[\*]{} + \_n\^2 \_[H]{} \_[H]{}\^[\*]{} + \_\
& = & \_ \^[\*]{} + \_n\^2 \_[H]{} \_[H]{}\^[\*]{} + \_n\^2 ( - \_[H]{} ) ( - \_[H]{} )\^[\*]{}\
& = & \_ \^[\*]{} + \_n\^2\
& = & \_

where we used the auxiliary result

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - Problem 2 - y autocorrelation - auxiliary result\] && \_[H]{} \_[H]{}\^[\*]{} + ( - \_[H]{} ) ( - \_[H]{} )\^[\*]{}\
&& = \^[\*]{} \_[\_[H]{}]{} \_[\_[H]{}]{}\^[\*]{} + \^[\*]{} ( - \_[\_[H]{}]{} ) ( - \_[\_[H]{}]{} )\^[\*]{}\
&& = \^[\*]{} \_[\_[H]{}]{} + \^[\*]{} ( - \_[\_[H]{}]{} )\
&& = .

Joining the two stages of the construction, presented in (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - y\_tilde as a sum \]) and (\[eq:appendix - proofs - equivalence of problem 1 and 2 - backward channel construction - y from y\_tilde \]), exhibits $ \hat{{\mathbf{x}}} \rightarrow \tilde{{\mathbf{y}}} \rightarrow {\mathbf{y}} $ as a Markov chain and, therefore, $ {\mathbf{y}} $ and $\hat{{\mathbf{x}}} $ are independent given $ \tilde{{\mathbf{y}}} $. This evident construction turns (\[eq:appendix - proofs - equivalence of problem 1 and 2 - data processing inequality\]) into

[rCl]{} \[eq:appendix - proofs - equivalence of problem 1 and 2 - mutual information equality\] I( ; ) = I( ; )

that completes proving the equivalence of Problems \[problem: basic form\] and \[problem: pseudoinverse filtered input\].

Equivalence of Problems \[problem: pseudoinverse filtered input\] and \[problem: Separable Form in DFT domain\] {#appendix:Equivalence of Problems 2 and 3}
---------------------------------------------------------------------------------------------------------------

The rate-distortion function for a Gaussian source with memory (i.e., correlated components) is usually derived in the Principle Component Analysis (PCA) domain where the components are independent Gaussian variables (see, e.g., [@berger1971rate]). In our case, where the signal is cyclo-stationary, the PCA is obtained using the DFT matrix. As in the usual case,

[rCl]{} \[eq:appendix - proofs - equivalence of problem 2 and 3 - mutual information equality under DFT\] I( ; ) & = & I( \^[F]{} ; \^[F]{} )\
\[eq:appendix - proofs - equivalence of problem 2 and 3 - mutual information reduction to sum\] & = & \_[k=0 ]{}\^[N-1]{} I( [\^F\_k]{} ; \^F\_k )

where (\[eq:appendix - proofs - equivalence of problem 2 and 3 - mutual information equality under DFT\]) emerges from the reversibility of the transformation, and (\[eq:appendix - proofs - equivalence of problem 2 and 3 - mutual information reduction to sum\]) is due to the independence of the $ \left\lbrace {\tilde{y}^F_k} \right\rbrace_{k=0}^{N-1} $ components [@cover2012elements] and the backward-channel construction (see Appendix \[appendix:Equivalence of Problems 1 and 2\]) that can be translated into a form of independent DFT-domain component-level channels.

The main difference from the well-known rate-distortion analysis is that here, in Problem \[problem: pseudoinverse filtered input\], the distortion constraint is not a regular squared error – but a weighted one, that will be developed next. Since DFT is a unitary transformation, its energy preservation property yields

[rCl]{} \[eq:appendix - proofs - equivalence of problem 2 and 3 - distortion energy preservation\] E ( [ - ]{} ) \_2\^2 & = & E \_[H]{} ( [ \^F - \^F ]{} ) \_2\^2\
& = & \_[k=0]{}\^[N-1]{} [ | h\^F\_k |\^2 E| [ \^F\_k - \^F\_k ]{} |\^2 ]{}

where the last equality is due to the diagonal structure of $ {\mathbf{\Lambda}}_{H} $. Hence, we got that the two expected-distortion expressions in Problems \[problem: pseudoinverse filtered input\] and \[problem: Separable Form in DFT domain\] are equal.

Solution of Problem \[problem: Gaussian distortion allocation\] {#appendix:Solution of Problem 4}
---------------------------------------------------------------

For a start, the transition between Problem \[problem: Separable Form in DFT domain\] and Problem \[problem: Gaussian distortion allocation\] is analogous to the familiar case of jointly coding a set of independent Gaussian variables [@cover2012elements]. Accordingly, and also due to lack of space, we do not elaborate here on this problem-equivalence proof.

Problem \[problem: Gaussian distortion allocation\] is compelling as it is a distortion-allocation optimization, where the distortion levels $ \left\lbrace D_k \right\rbrace_{k=0}^{N-1} $ are allocated under the joint distortion constraint. We address Problem \[problem: Gaussian distortion allocation\] via its Lagrangian form (temporarily ignoring the constraints of non-negative distortions)

[rCl]{} \[eq:appendix - proofs - solution of problem 4 - Lagrangian form\] \_[k=0]{}\^[N-1]{} [ \_[+]{} ]{} + \_[k=0]{}\^[N-1]{} [ | h\^F\_k |\^2 D\_k ]{}\

where $ \mu \ge 0 $ is the Lagrange multiplier. Recalling that some components may correspond to $ h^F_k = 0 $ and, by (\[eq:theoretic Gaussian analysis - pseudoinverse filtered input - autocorrelation eigenvalues\]), also $\lambda^{\left(\tilde{{\mathbf{y}}}\right)}_k = 0$ – meaning they are deterministic variables. These deterministic components do not need to be coded (i.e., $ R_k \nolinebreak = \nolinebreak 0 $) while still attaining $ D_k = 0 $. Accordingly, the Lagrangian optimization (\[eq:appendix - proofs - solution of problem 4 - Lagrangian form\]) is updated into

[rCl]{} \[eq:appendix - proofs - solution of problem 4 - Lagrangian form - only components with nonzero h\_k\^F\] \_[k: h\^F\_k 0]{} [ ( ) ]{} + \_[k: h\^F\_k 0]{} [ | h\^F\_k |\^2 D\_k ]{}\

where we used the expression from (\[eq:theoretic Gaussian analysis - pseudoinverse filtered input - autocorrelation eigenvalues\]), and assumed that the distortions are small enough such that the operator $ \left[\cdot\right]_+ $ can be omitted (a correct assumption as will be later shown). Now, the optimal $ D_k $ value can be determined by equating the respective derivative of the Lagrangian cost to zero, leading to allocated distortion (still as a function of $ \mu $)

[rCl]{} \[eq:appendix - proofs - solution of problem 4 - optimal distortion allocation - function of mu\] D\_k\^[opt]{} = k: h\^F\_k 0

and by setting the $ \mu $ satisfying the total distortion constraint from Problem \[problem: Gaussian distortion allocation\] we get

[rCl]{} \[eq:appendix - proofs - solution of problem 4 - optimal distortion allocation\] D\_k\^[opt]{} = k: h\^F\_k 0.

Expressing a nonuniform distortion-allocation (for components with nonzero $ h^F_k $), being inversely proportional to the weights $ \{ \left| h^F_k \right|^2 \} _{k=0}^{N-1} $. One should note that the assumption on small-enough distortions is satisfied as $ D_k^{opt} \le \lambda^{\left(\tilde{{\mathbf{y}}}\right) }_k $ for any $ k $ obeying $ h^F_k \ne 0 $, and that all the distortions are non-negative as required. The optimal distortions established here (for $ k $ where $ h^F_k \ne 0 $) are set in the rate formula (\[eq:theoretic Gaussian analysis - rate-distortion function of a Gaussian variable\]), providing the optimal rate allocation

[rCl]{} \[eq:appendix - proofs - solution of problem 4 - optimal rate allocation\] R\_k\^[opt]{} = {

[\*[20]{}[c]{}]{} [ ( | h\^F\_k |\^2 + 1 )]{}&[, h\^F\_k0 ]{}\
0&[, h\^F\_k = 0 . ]{}

.

Equivalent Form of Stage 6 of Algorithm \[Algorithm:Proposed Method Non-overlapping\] {#appendix:Equivalent Forms of Stage 4 of Algorithm 1}
-------------------------------------------------------------------------------------

The analytic solution of Stage 6 of Algorithm \[Algorithm:Proposed Method Non-overlapping\] is considered here with the conjugate-transpose operator, $ ^* $, extending the regular transpose:

[rCl]{} \[eq:appendix - proofs - Equivalent Forms of Stage 4 of Algorithm 1 - signal domain\] \^[(t)]{} & = & ( \^[\*]{} + )\^[-1]{} ( \^\* + \^[(t)]{} ) .

We note that

[rCl]{} \[eq:appendix - proofs - Equivalent Forms of Stage 4 of Algorithm 1 - H\^\*y\] \^\* & = & \^\* ( \^[+]{} + ( - \^[+]{} ) )\
& = & \^[\*]{} + \^[\*]{} ( - \^[+]{} )\
& = & \^[\*]{}

where the last equality results from the relation $ {\mathbf{H}}^{*} \left({\mathbf{I}} - {\mathbf{H}} {\mathbf{H}}^{+} \right) = {\mathbf{0}} $. Consequently, (\[eq:appendix - proofs - Equivalent Forms of Stage 4 of Algorithm 1 - signal domain\]) becomes

[rCl]{} \[eq:appendix - proofs - Equivalent Forms of Stage 4 of Algorithm 1 - signal domain - 2\] \^[(t)]{} & = & ( \^[\*]{} + )\^[-1]{} ( \^[\*]{} + \^[(t)]{} ) ,

which is the form presented in (\[eq:theory-practice relation - interpreting first stage\]).

[^1]: The authors are with the Department of Computer Science, Technion, Israel. E-mail addresses: {ydar, elad, freddy}@cs.technion.ac.il.

[^2]: Since the differences between Algorithm \[Algorithm:Proposed Method Non-overlapping\] and Algorithms \[Algorithm:Proposed Method Overlapping Blocks\] and \[Algorithm:Proposed Method Overlapping Blocks with Robust Dual Variables\] are for a shift-invariance purpose, an issue that we do not concern in this section, we compare our theoretic results only to Algorithm \[Algorithm:Proposed Method Non-overlapping\].

[^3]: The results in Table \[table:Deblurring - PSNR Comparison\] for the methods from [@neelamani2004forward; @guerrero2008image; @dabov2008image; @oliveira2009adaptive; @chantas2010variational; @danielyan2012bm3d] were taken as is from [@danielyan2012bm3d].
