Learning Sparse and Low-Rank Priors for Image Recovery via Iterative Reweighted Least Squares Minimization
Authors: Stamatios Lefkimmiatis
Iaroslav Sergeevich Koshelev
Venue: ICLR 2023 poster (2023)
Keywords: sparsity
low-rank
IRLS
inverse problems
optimization
recurrent networks

Abstract: In this work we introduce a novel optimization algorithm for image recovery under learned sparse and low-rank constraints, which are parameterized with weighted extensions of the $\ell_p^p$-vector and $\mathcal{S}_p^p$ Schatten-matrix quasi-norms for $0\!<p\!\le1$, respectively. Our proposed algorithm generalizes the Iteratively Reweighted Least Squares (IRLS) method, used for signal recovery under $\ell_1$ and nuclear-norm constrained minimization. Further, we interpret our overall minimization approach as a recurrent network that we then employ to deal with inverse low-level computer vision problems. Thanks to the convergence guarantees that our IRLS strategy offers, we are able to train the derived reconstruction networks using a memory-efficient implicit back-propagation scheme, which does not pose any restrictions on their effective depth. To assess our networks' performance, we compare them against other existing reconstruction methods on several inverse problems, namely image deblurring, super-resolution, demosaicking and sparse recovery. Our reconstruction results are shown to be very competitive and in many cases outperform those of existing unrolled networks, whose number of parameters is orders of magnitude higher than that of our learned models. 

Published as a conference paper at ICLR 2023
LEARNING SPARSE AND LOW-RANK PRIORS FOR IMAGE RECOVERY VIA ITERATIVE REWEIGHTED LEAST SQUARES MINIMIZATION
Stamatios Lefkimmiatis, laroslav Koshelev Huawei Noah's Ark Lab {stamatios lefkimmiatis, koshelev.iaroolay com
ABSTRACT
We introduce a novel optimization algorithm for image recovery under learned sparse and low-rank constraints, which we parameterize as weighted extensions of the ve-vector and Sp. Schatten-mate matrix quasi- quasi-norms for 0< p 1, respectively. Our proposed algorithm generalizes the Iteratively Reweighted Least Squares (IRLS) method, used for signal recovery under l and nuclear-nore constrained mini- mization. Further, we interpret our overall minimization approach as a recurrent network that we then employ to deal with inverse low-level computer vision prob- lems. Thanks to the convergence guarantees that our IRLS strategy offers, we are able to train the derived reconstruction networks using memory-efficiee implicit back- back-propagation scheme, which does not pose any restrictions on their effective depth. To assess networks' performance, other existing reconstruction methods on several inverse problems, namely image deblurring, super-resolution, demosaicking and sparse recovery. Our reconstruc- tion results are shown to be very competitive and in many cases outperform those of existing unrolled networks, whose number of parameters orders of magnitude higher than that of our learned models. The code is available at this link
compare them against
INTRODUCTION
With the advent of modern imaging techniques, we are witnessing significant rise of interest in inverse problems, which appear increasingly in host of applications ranging from microscopy and medical imaging to digital photography, 2D&3D computer vision, and astronomy (Bertero & Boccacci] 1998). An inverse imaging problem amounts to estimating a latent image from set of possibly incomplete and distorted indirect measurements. In practice, such problems are typical ill- posed (Tikhonov 1963; Vogel 2002), which implies that the equations relating the underlying image with the measurements (image formation model) are not adequate by themselves to uniquely char- acterize the solution. Therefore, in order recover approximate solutions, which a statistical or physical sense, from the set of solutions that are consistent with the image formation model, it is imperative exploit prior knowledge about certain properties of the underlying image.
meaningful in
Among the key approaches for solving ipl-posed inverse problems are variational methods (Benning & Burger 2018), which entail the minimization of objective function. crucial part of such an objective function is the regularization term, whose role is to promote those solutions that fit best our prior knowledge about the latent image. Variational methods have also direct links Bayesian methods and can be interpreted as seeking the penalized maximum likelihood or the maximum posteriori (MAP) estimator (Figueiredo et al.] 2007), with the regularizer matching the negative log-prior. Due to the great impact of the regularizer in the reconsturction quality, significant re- search effort has been put in the design of suitable priors. Among the overwhelming number of existing priors in the literature, sparsity and low-rank (spectral-domaal have received considerable attention. This is mainly due to their solid mathematical foundation and the competitive results they achieve (Bruckstein et al. 2009 Mairal al. 2014).
sparsity) promoting priors
Nowdays, thanks to the advancements of deep learning there is a plethora of networks dedicated to image reconstruction problems, which significantly outperform conventional approaches. Nev- ertheless, they are mostly specialized and applicable to a single task. Further, they are difficult to analyze and interpret since they do not explicitly model any of the well-studied image properties,

---

Published as a conference paper at ICLR 2023
successfully utilized in the past (Monga et al. 2021). In this work, we aim to harness the power of supervised leanniig but at the same time rely on the rich body of modeling and algorithmic ideas that have been developed in the past for dealing with inverse problems. To this end our contributions are: (1) We introduce a generalization of the Iterative Reweighted Least Squares (IRLS) method based on novel tight upper-bounds that we derive. (2) We design recurrent network architecture that plicitly models
image priors and is applicable to a wide range of reconstruction problems. (3) We propose a memory efficient training strategy based on implicit back- pack-propagttion that does not restrict in any way our network's ' effective depth.
2
IMAGE RECONSTRUCTION
Let us first focus on how one typically deals with inverse imaging problems of the form:
y Ax n.
(1)
where x E RR- is the multichannel latent image of e channels, that we seek to recover, A TR - Rm. is a linear operator that models the impulse response of the sensing device, y E mmm is the observation vector, and E IRm is noise vector that models all approximation errors of the forward model and mcasurement noise. Hereafter, we will assume that n consists of i.i.d samples drawn from a Gaussian distribution of zero mean and variance o Note that despite of the seeming simplicity of this observation model, it is widely used in the literature, since it can accurately enough describe plethora of practical problems. Specifically, by varying the form of A, Eq. ( can cover many different inverse imaging problems such as denoising, deblurring, demosaicking, inpainting, super-resolution, MRI reconstruction, etc. If we further define the objective function:
I (x)
_. ly Aa | R (x²);
where R: + R+ E RR 2 0} is the regularizer (image prior), we can recover an estimate of the latent image x² as the minimizer of the optimization problem: x' arg ming J (z). Since the type of the regularizer R (x) can significantly affect the reconstruction quality, it is of the utmost importance to employ a proper regularizer for the reconstruction task at hand.
2.1
SPARSE AND LOW-RANK IMAGE PRIORS
Most of the existing image regularizers in the literature can be written in the generic form:
R (())
Σ Φ (G
(3)
where - is the regularization operator that transforms x, G. = M,G, M,=L, with & denoting the Kronecker product, and e; is the unit vector of the standard RR basis. Further, : Rd + R2 is a potential function that penalizes the response of the d-dimensional transform- domain feature vector, 2; = G. E Rd. Among such regularizers, widely used are those that promote sparse and low-rank responses by utilizing as their potential functions the l and nuclear norms (Rudin et al., 1992 Figueiredo et al.) 2007; Lefkimmiatis al. 2013 2015). Enforcing spar- sity of the solution in some Transform-doman domain has been studied in-depth and is supported both by solid mathematical theory (Donoho 2006 Candes & Wakin 2008 Elad 2010)
well as strong empirical results, which indicate that distorted images do not typically exhibit sparse or low-rank representations, as opposed their clean counterparts. More recently it has also been advocated that non-convex penalties such as the Eg vector and St Schatten-matrix quasi-norms enforce spar- sity better and lead improved image reconstruction results (Chartrand) [007]] Lai al. 2013 Candes et al. 2008; Gu et al. 2014: Liu et al. (0044) Xie et al.] 2016; Kümmerle & Verdun 2021
Based on the above, we consider two expressive parametric forms for the potential function o (-), which correspond to weighted and smooth extensions of the Ey and the Schatten matrix Sy quasi- norms with < p S 1, respectively. The first one is a sparsity- sparsity-promoting promoting penalty, defined as:
Ssp (z; w. w, wy (z) + +7) &
2, w E Rd
(4)
while the second one a low-rank (spectral- domain sparsity) promoting penalty, defined as:
dir (Z: w. w.p) p) = Σ w (o (z) + +7) Z E I mxn w E RC, with r = min (m, n).
(5)
1=1
2

---

Published as a conference paper at ICLR 2023
In both definitions vy is small fixed constant that ensures the smoothness of the penalty functions. Moreover, in Eq. ( o ( denotes the vector with the singular values of Z sorted in decreasing order, while the weights better promote low-rank solutions we need to penalize more the smaller singular values of the matrix than its larger ones. Next, we define our proposed sparse and low-rank promoting image priors as:
sorted in increasing order. The reason for the latter order is that to
Rop (a)) = Σ sp (z); W.,
(6a)
Rir (x)
Øir (ZN; P),
(6b)
where in Eq. 6ap 2; = G E Rd, while in Eq. (6b) Z E Rex0 is matrix whose dependence on is expressed as: vec (Zi) = G. E F., with d= c. q. In words, the j-th row of Z, is formed by the g-dimensional feature vector z (..)) extracted from the image channel x'. = 1....... The motivation for enforcing the matrices Z, to have low-rank is that the channels of natural images are typically highly correlated. Thus, it is reasonable to expect that the features stored in the rows of Z; would be dependent to each other. We note that a similar low-rank enforeing regularization strategy is typically followed by regularizers whose goal is to model the non-local similarity property of natural images Gu et al. (2014); ; Xie et al. (2016). In these cases the matrices Z; are formed such way that each of their rows holds the elements of patch extracted from the image, while the entire matrix consists of groups of structurally similar patches.
2.2 MAJORIZATION MAJORIAATIONNINMMIATATOOO MINIMIZATION STRATEGY
One of the challenges in the minimization of the overall objective function in Eq. (2) is that the image priors introduced Eq. (6) are generally non-convex w.r.t æc This precludes any guarantees of reaching a global minimum, and we can only opt for a stationary point. One way to handle the minimization task would be to employ the gradient descent method. Potential problems in such case are the slow convergence as well as the need to adjust the step-size in every iteration of the algorithm, so as to avoid divergence from the solution. Other possible minimization strategies are variable splitting (VS) techniques such as the Alternating Method of Multipliers (Boyd et al., 2011) and Half- Quadratic splitting (Nikolova & Ng 2005), or the Fast Iterative Shrinkage Algorithm (FISTA) (Beck) & Teboulle 2009). The underlying idea of such methods is to transform the original problem in easier to solve sub-problems. However, VS techniques require finetuning of additional algorithmic parameters, to ensure that satisfactory convergence rate achieved, while FISTA works-well under the assumption that the proximal map of the regularizer (Boyd & Vandenberghe (2004) is not hard to compute. Unfortunately, this not the case for the regularizers considered in this work.
For all the above reasons, here we pursue majorization-miiiiizzaoo (MM) approach (Hunter & Lange (004)), which does not pose such strict requirements. Under the MM approach, instead of trying to minimize the objective function J (x) directly, we follow an iterative procedure where each iteration consists of two steps: (a) selection of a surrogate function that serves as a tight upper- bound of the original objective (mujorization-stepp and (b) computation of current estimate of the solution by minimizing the surrogate function (minimization-sttep) Specifically, the iterative algorithm for solving the minimization problem x' arg min, J (x) takes the form: =+1
arg min, Q z*), where Q (x; is the majorizer of the objective function I (x)) at some point æ², satisfying the two conditions:
Q (xx 2 J (z), x* and Q (a*;=*) = (z*).
Given these two properties of the majorizer, it can be easily shown that iteratively minimizing Q (x:=*) also monotonically decreases the objective function J (x) (Hunter & Lange 2004). In fact, to ensure this, it only suffices to find x++1 that decreases the value of the majorizer, i.e., Q (x++ S Q (ar*; a+). Moreover, given that both @ (a; a²*) and I (ar) are bounded from below,
safely state that upon convergence we will reach stationary point.
The success of the described iterative strategy solely relies on our ability to efficiently minimize the chosen majorizer. Noting that the data fidelity term of the objective is quadratic, we proceed by seeking a quadratic majorizer for the image prior R (z). This way the overall majorizer will be of quadratic form, which is amenable to efficient minimization by solving the corresponding normal equations. Below we provide two results that allow us to design such tight quadratic upper bounds both for the
(6a) and the low-rank promoting (6b) regularizers. Their proofs are
3

---

Published as a conference paper at ICLR 2023
provided in Appendix A.1 We note that, to the best of our knowledge, the derived upper-bound
presented in Lemma 2 is novel and it can find use in a wider range of applications, which also utilize low-rank inducing penalties (Hu et al.] 202T), than the ones focus here.
Lemma 1. Let æ, y E IR" and w E RT. The weighted- Eg, function Osp w.;w)) defined in Eq. (4) can be upper-bounded as:
Øap (2; wp S P W , + p2 tr (W,))
2-P 2 sp (y: w. p). Ya,y
2
where W, = diag (w) I (yy) (1]
- and o denotes the Hadamard product. The equality
in is attained when x =
Lemma 2. Let X, Y E Rmx= and o (Y), E R with r= min (m, n). The vector o (Y) holds the singular values of Y in decreasing order while the elements of w are sorted in increasing order. The weighhed-Scchaten--nee Schatten-matri function our (X; w. defined in Eq. (5) can be upper-bounded as:
Φtr (X; p) S 2 XXXX 2 tr (Wy) + 2 oi (Y; Y.W,Pl Y
where WY U diag (w) UT (YYT +11) and Y = U diag (a (Y)) VT, with U E RRxr
and V E RRxx The equality in 2 is altained when X =Y.
Next, with the help of the derived tight upper-bounds we obtain the quadratic majorizers for both of our regularizers as:
Cep (z: (r="")= = 2 2 (10a)
2ir xx) = tr Z W z 2.)
(10b)
where and Z; are defined as in Eq. 6 W x W 2 are defined in Lemmas and respectively, and we have ignored all the constant terms that do not depend on æ. In both cases, by adding the majorizer of the regularizer, 2reg (z:; x²), ) to the quadratic data fidelity term we obtain the overall majorizer Q (z: x²). of the objective function J (æ)). Since this majorizer is quadratic, it is now possible to obtain the (k+1)-th update of the MM iteration by solving the normal equations:
k+1 A + o 2 Σ
G. G. G.
+ +
a1)
(A y + ta*t) = (s* + +al) b²,
(1I)
where a= 5o W E for Rep (z), and - Wzt 2 E Roaxx for Rir (z). We note that, to ensure that the system matrix in Eq. (IT) is invertible, we use an augmented majorizer that includes the additional term & || to Appendix I. for a justification of the validity of this strategy). This leads to the presence of the extra term af additionally to the system matrix S and of ax in b² Based on the above, the minimization of J (x), incorporating any of the two regularizers of Eq. (6). boils down to solving sequence of re-weighted least squares problems, where the weights W of the current iteration are updated using the solution of the previous one. Given that our regularizers in Eq. (6) include the weights un f 1, our proposed algorithm generalizes the IRLS methods introduced by (Dabbechies et al. (2010) and Mohan & Fazel (2012), which only consider the case where w, 1 and have been successfully applied in the past on sparse and low-rank recovery problems, respectively.
fir-x* x* [|. with > 0 being fixed small constant (we refer
3 LEARNING PRIORS USING IRLS RECONSTRUCTION NETWORKS
To deploy the proposed IRLS algorithm, we have to specify both the regularization operator G and the parameters w [w,) / i= PP of the potential functions in Eqs. 4 and (5). which constitute the image regularizers of Eq. (6). Manually selecting their values, with the goal of achieving satisfactory
reconstructions, can be cumbersome task. Thus, instead of explicitly defining them, we pursue the idea of implementing IRLS as a recurrent network, and learn their values in a supervised way. Under this approach, our learned IRLS-based reconstruction network (LIRLS) can be described as: x k+1 = fe (x=; 1)), with 0= (G. w, p} denoting its parameters. The network itself consists of three main components: (a) A feature extraction layer that accepts input the current reconstruction estimate, x²,, and outputs the feature maps {2 = G
} i= (b) The weight module that acts on = and the parameters w; to construct the weight matrix W which is part of the system matrix S in Eq. (IT). (c) The least- least-squares (LS) layer whose role is to produce refined reconstruction
estimate, a as the solution of Eq. II] The overall architecture of LIRLS is shown in Fig.
4

---

Published as a conference paper at ICLR 2023
2 or a 1
On WPN w
{7.p)
{A. n o.p)
3.1 NETWORK TRAINING
common training strategy for re- current networks is to unroll the net- work using a fixed number of iter- ations and update its parameters ei- ther by means of back-propagation through time (BPTT) or by its trun- cated version (TBPTT) (Robinson & Kokkinos & Letkim-
Feature extraction
IRLS weight matrix
W
layer
x+1
Eq. (1)) (PCG)
(Lemmas 1 or 2)
G (shared weights)
feedbackkloop
Figure 1: The proposed LIRLS recurrent architecture.
Fallside 1987 7
miatis| 2019). However, this strategy cannot be efficiently applied to LIRLS. The reason is that the system matrix S in Eq. [T] typically is of huge dimensions and its direct inversion is generally in- feasible. Thus, to compute x k+1 via Eq. (11) we need to rely on matrix-free iterative linear solver such as the conjugate gradient method (CG) (Shewchuk) 1994). This means that apart from the IRLS iterations we need also to take into account the internal iterations required by the linear solver. Therefore, unrolling both types of iterations would result in very deep network, whose efficient training would be extremely challenging for two main reasons. The first is the required amount of memory, which would be prohibitively large. The second one is the problem of vanishing/exploding gradients that appears in recurrent architectures (Pascanu et al., 2013), which would be even more pronounced the case of LIRLS due to its very deep structure.
To overcome these problems, we rely on the fact that our IRLS strategy guarantees the convergence to a fixed point x'. Practically, this means that upon convergence of IRLS, it will hold x' = fe (x Considering the form of our IRLS iterations, as shown in Eq. [II, this translates to:
9 (x 0) x² fe (x S (x", 000x'' A y
(12)
where we explicitly indicate the dependency of S on xx and O. To update the network's parameters during training, we have compute the gradients VoL (ax) Vox" V.. f (ax), where L is loss function. Now, if we differentiate both sides of Eq. 12) w.r.t 0 then we get:
0g (x" 0) 00
0g (x:", 0) Ox
Vox Veg (x* 0) (V w g (a" 0))
(13)
Ox
00
Thus, we can now compute the gradient of the loss function w.r.t the network's parameters as:
VeL (x*) = - - (x" 0)v.
(14)
where is obtained as the solution of the linear problem Va-g (a* V-. (x²) and all the necessary auxiliary gradients can be computed via automatic differentiation. Based on the above, train the LIRLS network without having restrict its overall effective depth save any intermediate results that would significantly increase the memory requirements. The implementation details of our network for both its forward and backward passes, are described in Algorithm | in Sec. A.4 Finally, while our training strategy is similar in spirit with the one used for training Deep Equilibrium Models (DEQ) (Bai et al. 2019). an important difference is that our recurrent networks are guaranteed to converge to a fixed point, while in general this is not true for DEQ models.
3.2 LIRLS NETWORK IMPLEMENTATION
In this section we discuss implementation details for the LIRLS variants whose performance we will assess on different grayscale/color image reconstruction tasks. As mentioned earlier, among the parameters that we aim to learn is the regularization operator G. For all the different networks we parameterize G with valid convolution layer that consists of 24 filters of size of color images these filters are shared across channels, Further, depending on the particular LIRLS instance that we utilize, we either fix the values of the parameters w, p or we learn them during training. Hereafter, we will use the notations Ep and Sp.w to refer to the networks that employ learned sparse-pronee prommotin prior and learned low-rank promoting prior, respectively. The different LIRLS instances that we consider in this work are listed below:
In the case
1. I (nuclear): fixed p=1, fixed w w=1, where 1 is a vector ones.
Weighted 1S" (weighted nuclear): fixed p= 1, weights w are computed by a weight prediction neural network (WPN). WPN accepts as inputs either the features 2n = Gx for the grayscale case their singluar values for the color case, well the noise standard deviation on. The
5

---

Published as a conference paper at ICLR 2023
Table 1: Comparisons on grayscale image deblurring.
Datnset
TV-2
TV
41
'1"
2
RED
IRCNN
FDN
Sun et al, 2013
PSNR
30.81
30.96
31,98
32.36
32.70
32.67
31.68
32.55
32.63
SSIM
0.8272
0.8351
0.8759
0.8849
0.8965
0.8960
0.8573
0.8860
0.8894
Levin et al
2009
PSNR
34.74
34.92
35.19
35.33
35.75
35.60
35.19
32,86
35,08
SSIM
0,9552
0.9568
0.9600
0.9617
0.9647
0,9628
0.9584
0.9118
0.9609
Table 2: Comparisons on color image deblurring.
TVN
VTV
Si
Sp
S1"
Sp
RED
IRCNN
FDN
DWDN
PSNR
31.71
31.38
33.09
34.17
33,95
34.24
31.38
33.19
32.51
34.09
SSIM
0,8506
0.8440
0.8995
0.9227
0,9179
0,9234
0.8233
0.9123
0.8857
0.9197
vector Do is the estimate obtained after IRLS steps of the pretrained / networks, while G is their learned regularization operator. For all studied problems except of MRI reconstruction, we use a lightweight RFDN architecture proposed by Liu et al. (2020) to predict the weights, while for MRI reconstruction we use lightweight UNet from [Roneeergger al. (2015), which we have found to be more suitable for this task. For the S network, in order to enforce the predicted weights to be sorted in increasing order, we apply across channels of the ouput of WPN a cumulative sum. In all cases, the number of parameters of WPN does not exceed 0.3M.
3. e:spp learned p [0.4, 0.9), fixed w=1.
4. p ep."/Sp.. learned p E [0.4, 0.9), weights w are computed as described in item 2 with the only difference being that both 2g and G are now obtained from the pretarained E/SS networks.
We note that the output of the LS layer, which corresponds to the solution of the normal equations in Eq. (T). is computed by utilizing preconditioned version of CG (PCG) Hestenes & Stiefel 1952). This allows for an improved convergence of the linear solver. Details about our adopted preconditioning strategy are provided in Appendix/A.2 Finally, in all the reported cases, the constant ó related to the parameter in Eq. (II) is set to Se -. while as initial solution for the linear solver we use the output of an independently trained fast FFT-based Wiener filter.
3.3 TRAINING DETAILS
The convergence of the forward pass of LIRLS to a fixed point x" is determined according to the criterion: ||S* (x* 0) x - A' y y/l2/11 A ylll2 < le- which needs to be satisfied for 3 consecutive IRLS steps. If this criterion is not satisfied, the forward pass is terminated after 400 IRLS steps during training and 15 steps during inference. In the LS layers we at most 150 CG iterations during training and perform an early exit if the relative tolerance of the residual falls below le while during inference the maximum amount of CG iterations is reduced to 50. When training LIRLS, in the backward pass, as shown in Eq. (14) we are required to solve a linear problem whose system matrix corresponds the Hessian of the objective function (2). This symmetric matrix is positive definite when J (x) is convex and indefinite otherwise. In the former case we utilize the CG algorithm to solve the linear problem, while in the latter one we use the Minimal Residual Method (MINRES) (Paige & Saunders 1975). We perform early exit if the relative tolerance of the residual below le - and limit the maximum amount of iterations to 2000. should be noted that we have experimentally found these values to be adequate for achieving stable training.
All our models are trained using as loss function the negative peak-to-signal-nne -oise--aioos (PSNR) between the ground truth and the network's output. We use the Adam optimizer (Kingma & Ba 2015) with a learning rate 5e -3 for all the models that do not involve a WPN and Le otherwise. The learning rate is decreased by factor of 0.98 after each epoch. On average we set the batch size to 8 and train our models for 100 epochs, where each epoch consists of 500 batch passes.
4 EXPERIMENTS
In this section assess the performance of all the LIRLS instances described in Sec. 3.2 four different image reconstruction tasks, namely image deblurring, super-resolution, demosaicking and MRI reconstruction. In all these cases the only difference in the objective function (a)) is the form of the degradation operator A. Specifically, the operator A has one of the following forms: (a) low-pass valid convolutional operator (deblurring). (b) composition of low-pass valid convo- lutional operator and a decimation operator (super-resolution), (c) color filter array (CFA) operator (demosaicking). and (d) sub-sampled Fourier operator (MRI reconstruction). The first two recovery tasks are related to either grayscale or color images, demosaicking is related to color images, and MRI reconstruction related to single-channel images.

---

Published as a conference paper at ICLR 2023
Table 3: Comparisons on grayscale image super-resolution.
Scale
Noise
TV-2
TV
f
2
41"
Bicubic
RED
IRCNN
USRNet
0%
PSNR
26.23
28.04
28.20
27.98
28.76
28.31
24.74
28.34
28.58
SSIM
0.6937
0.8028
0.8064
0,7953
0,8186
0,8067
0.6699
0,8213
0.8324
0.8188
x2
FSNR
17.37
26.64
26.81
26.64
27.32
26,94
24.63
26.88
26.22
27.26
1%
SSIM
0,2739
0.7320
0.7379
0,7292
0,7610
0,7470
0.6521
0,7275
0.6799
0.7607
PSNR
24.50
25,55
25.63
25.44
26.02
25.67
23.34
25.70
25.91
26.06
0%
SSIM
0,5875
0.6825
0.6857
0.6741
0.6975
0,6832
0,5872
0,7006
0,7144
0.7101
x3
1%
FSNK
14.36
24.02
24.67
24.48
25.03
24.79
23,26
24.81
24.70
25.07
SSIM
0,1812
0.6276
0.6306
0.6195
0.6472
0.6353
0.5739
0.6417
0.6219
0.6586
PSNR
22.51
24.07
24.09
23.92
24,16
24.10
22.27
24.15
23.72
23.30
0%
SSIM
0.4773
0.6031
0.6052
0.5960
0.6061
0.6010
0,5342
0.6164
0.6168
0.6203
x4
1%
PSNR
11.81
23.31
23.32
23.13
23.35
23.38
22.21
23.35
21.41
SSIM
0.1228
0,5627
0.5643
0.5551
0,5661
0,5648
0.5247
0,5731
0,5713
0.4945
Table 4: Comparisons on color image super-resolution.
Scale
Noise
TVN
VTV
S1
S
Sq
S
Bicubic
RED
IRCNN
USRNet
0%
PSNR
28.01
28.13
28.21
28.62
28.42
28.78
24.73
28.68
28.33
28,86
SSIM
0.8035
0.8042
0.8088
0.8203
0.2122
0.8217
0,6661
08210
0,8331
0.8359
x2
1%
FSNK
26.87
26.86
27.06
27.43
27.46
17.48
24.62
25.74
26.37
27.97
SSIM
0.7404
0.7392
0.7469
0.7646
0.7640
0.7687
0.6486
0.6215
0.6843
0.7924
FSNR
25.46
25.02
25.50
25.60
25.35
25.00
23.35
25,90
25.74
26.17
0%
SSIM
0,6812
0,6823
0.6828
0.6904
0,6877
0,6928
0,5833
0,7030
0,7092
0.7203
1%
PSNR
24.76
24.76
24.83
25.09
25.08
25.11
23.25
24.72
24.77
23.65
SSIM
0.6322
0,6314
0.6348
0,6468
0,6452
0.6497
0,5702
0.6078
0,6203
0.6857
0%
PSNR
24.05
24.13
24.15
24.28
24.22
24.27
22.25
24.33
23.46
23,30
SSIM
0.6018
0.6024
0.6045
0.6111
0.6079
0.6101
0,5305
0.6191
0.6073
0.6238
x4
1%
PSNK
23.43
23.42
23.46
23,66
23,53
23.66
22.19
23.64
23.40
23.57
SSIM
0.5655
0.5645
0.5666
0.5752
0.5698
0.5774
0,5210
0.5681
0.5643
0.6028
4.1
TRAIN AND TEST DATA
To train our models for the first three recovery tasks, we use random crops of size 64 x 64 taken from the train and val subsets of the BSD500 dataset provided by Martin al. (200)) while for MRI reconstruction we use 320 x 320 single-coil 3T images from the NYU fastMRI knee dataset (nolll et al. (2020). In order to train the deblurring networks we use synthetically created blur kernels varying in size from 13 to 35 pixels according to the procedure described by [Boracchi & Foi 2012). For the super-resolution task we use scale factors 2, 3 and 4 with 25 25 kernels randomly synthesized using the algorithm provided by |Bell-Kligler et al. (2019). For accelerated mnnjugate-symmette
MRI we consider x4 and x8 undersampling factors in k-space using so that training and inference is performed in the real domain. For the demosaicking problem we use the RGGB CFA pattern. All our models are trained and evaluated on a range of noise levels that are typical in the literature for each considered problem. Based on this, we consider on to be up to 1% of the maximum image intensity for deblurring, super-resolution and MRI reconstruction tasks, while for demosaicing we use wider range of noise levels, with on being up to 3%.
masks
For our evaluation purposes we use common benchmarks proposed for each recovery task. In par- ticular, for deblurring we use the benchmarks proposed by (Snn et al. (2013) and Levin et al. (2009). For super-resolution we use the BSD100RK dataset, which consists of 100 test images from the BSD500 dataset, and the degradation model proposed by Bell- et al. 2019). For demosaick- ing we use the images from the McMaster (Zhang et al. 2011)) dataset. For MRI reconstruction we use the 3T scans from the val subset of the NYU fastMRI knee dataset, where from each scan we take the central slice and two more slices located 8 slices apart from the central one. None of our training data intersect with the data we use for evaluation purposes. We should further note, that all benchmarks we use for evaluation contain diverse sets of images with relatively large resolutions. This strategy is not common for other general purpose methods, which typically report results on a limited set of small images, due to their need to manually fine-tune certain parameters. In our case all network parameters are learned via training and remain fixed during inference. In order to provide fair comparison, for the methods that do require manual tuning of parameters, perform a grid search on a smaller subset of images to find the values that lead to the best results. Then we use these values fixed during the evaluation for the entire set of images per each benchmark. Finally, note that for all ba-based regularizers we compute their solutions using an IRLS strategy.
4.2 RESULTS
Deblurring. In Table I we report the average results in terms of PSNR and structure-ste similarit index measure (SSIM) for several of our LIRLS instances, utilizing different pparsity-promoot priors, and few competing methods, namely anisotropic (TV-(; (Zach et al. 2007) and isotropic
7

---

Published as a conference paper at ICLR 2023
Table 5: Comparisons on image demosaicking.
Dataset
Noise
TVN
VTV
S1
5
S1
Sp
Bilinear
RED
IRONN
0%
33.33
32.64
37.66
32.27
37.52
SSIM
0.9377
0.9285
0.9605
0.9649
0.9606
0.9632
0,9259
0.9338
0.9612
1% SSIM
PSNR
33.32
32.16
35.14
36.43
35.49
56.37
31,76
34,40
36.12
McMaster
0.9198
0,9076
0,9416
0,9498
0.9451
0.9490
0,9014
0.9208
0,9403
Zhang et al. 2011
2011
2% SSIM
PSNR
32.21
31.16
33.37
34.35
33.95
34.80
30.64
33.08
34.42
0.8843
0,8638
0.9151
0.9224
0.9226
0,9281
0,8443
0.8821
0.9178
PSNR
31.19
30.13
31.39
32,44
32,60
33.34
29.34
31.39
33:23
3%
0.8539
088164
0.8808
0.7745
0.8308
0.8931
Table 6: Comparisons on MRI reconstruction.
Acceleration
TV-
4
2
PSNR
26.31
x4
SSIM
0.5883
0.5921
0.6585
0.6701
0,7081
0.6958
0.6014
PSNR
24.43
24.51
25.86
26.19
28.02
27.30
24.34
X8 SSIM
0,5092
0.5120
0.5518
0.5620
0.6168
0.5973
0.5237
Total Variation (Rudin et al.
1992). RED (Romano et al. 2017), IRCNN (Zhang et al. 2017).
and FDN (Kruse et al.J 2017). From these results we observe that our LIRLS models lead very competitive performance, despite the relatively small number of learned parameters. In fact, the e and fgww variants obtain the best results among all methods, including FDN, which is the current grayscale sota method, and IRCNN that involves 4.7M parameters in total (for its 25 denoising networks). Similar comparisons for color images are provided in Table 2 where we report the reconstruction results of our LIRLS variants, which utilize different learned low-rank promoting priors. For these comparisons we further consider the vector-valued TV (VTV) (Blomgren & Chan 1998), Nuclear TV (TVN) Lefkimmiatis et al. (2015) and DWDN (Dong et al. 2020). From these results we again observe that the LIRLS models perform very well. The most interesting observation though, is that our Sy model with total of 601 learned parameters manages to outperform the DWD) network which involves 7M learned parameters. It is also better than St" and very close to Spe" which indicates that the norm order p plays more crucial role in this task than the weights w.
Super-resolution.
Similarly to the image deblurring problem, we provide comparisons among com- peting methods both for grayscale and color images. For these comparisons we still consider the RED and IRCNN methods as before and we further include results from bicubic upsampling and the USRNet network Zhang et al. (2020). The latter network, unlike RED and IRCNN, specifically designed to deal with the problem of super-resolution and involves 17M of learned parameters. From the reported results we can see that for both the grayscale and color cases the results we obtain with our LIRLS instances are very competitive and they are only slightly inferior than the specialized USRNet network. For a visual assessment of the reconstructions quality we refer to Fig.
Demosaicking, For this recovery task we compare our low-rank promoting LIRLS models against the same general-purpose classical and deep learning approaches as color deblurring, plus bicubic upsamplping. The average results of all the competing methods are reported in Table 5 From these results we come to a similar conclusion as in the debluring case. Specifically, the best performance is achieved by the S and Sp.w LIRLS models. The first one performs best for lower noise levels while the latter achieves the results for the highest ones. Visual comparisons among the different methods are provided in Fig. 4
MRI Reconstruction. In Table we compare the performance of our LIRLS models with anisotropic and isotropic TV reconstruction and the Fourier back-projection (FBP) method. While
Input
TVN
SY
IRCNN
DWDN
Figure 2: Visual comparisons among several methods on a real blurred color image. Image and blur kernel were taken from Pan et al. (2016). For more visual examples please refer to Appendix[A.7
8

---

Published as a conference paper at ICLR 2023
Input (Ricuhie)
TVN
$
S1
RED
USRNer
Figure 3: Visual comparisons among several methods real low-resolution image enlarged by scale factor Image from Cai et al. (2019), downscaling kernel was obtained using the method by Bell- Kligler et al. (2019). For more visual examples please refer to Appendix AA.7
these three methods are conventional ones, they can be quite competitive on this task and are still relevant due to their interpretability, which very important in medical applications. Based the reported results and the visual comparisons provided in Fig. | we can conclude that the LIRLS also worth noting that similarly
models do very good job in terms of reconstruction quality. to the other two gray-scale recovery tasks, we observe that the best performance among the LIRLS models achieved by Z This can be attributed the fact that the learned prior of this model while being adaptive, thanks to the presence of the weights, is still convex, unlike Ey and fp.w. There- fore, its output is not significantly affected by the initial solution. Additional discussions about our networks reconstruction performance related to the choice of and the weights is provided in Appendix A.3 along with some possible extensions that we plan to explore the future.
5 CONCLUSIONS
In this work we have demonstrated that our proposed IRLS method, which covers rich family of sparsity and low-rank promoting priors, can be successfully applied to deal with wide range of practical inverse problems. In addition, thanks its convergence guarantees have managed to use in the context of supervised learning and efficiently train recurrent networks that involve only a small number of parameters, but can still lead very competitive reconstructions. Given that most of the studied image priors are non-convex, an interesting open research topic is to analyze how the initial solution affects the output the LIRLS models and how we can exploit this knowledge to further improve their reconstruction performance. Finally, it would be interesting to explore whether learned priors with spatially adaptive norm orders can lead to additional improvements.
Bilincar: 32.05
TVN: 32.00
S 35.07
S1º: 32.57
RED: 32.12
IRCNN: 33.93
Target
Figure 4: Visual comparisons among several methods on mosaicked image with 1 1% noise. For each image its PSNR value is provided in dB. For more visual examples please refer to Appendix A.7
FBP: 22.95
25.25
Ex: 27.77
28.59
EMM: 30.30
Target
Figure 5: Visual comparisons among several methods simulated MRI with x4 acceleration. For each image its PSNR value is provided in dB. For more visual examples please refer to Appendix A.7

---

Published as a conference paper at ICLR 2023
REFERENCES
Shaojie Bai, J. Zico Kolter, and Vladlen Koltun. Deep equilibrium models. In Advances in Neural Information Processing Systems, volume 32, 2019.
A. Beck and M. Teboulle. A fast iterative threesooiggr algorithm for linear inverse prob- lems. SIAM J. Imaging Sci., 2:183-202, 2009.
Sefi Bell-Kligler, Assaf Shocher, and Michal Irani. Blind super-resolution kernel estimation using an internal-gan. Advances Neural Information Processing Systems, 32, 2019.
Martin Benning and Martin Burger. Modern regularization methods for inverse problems. Acta Numerica, 27:1-111, 2018
Mario Bertero and Patrizia Boccacci. Introduction to Inverse Problems in Imaging. IOP Publishing, 1998.
P. Blomgren and T.F. Chan. Color TV: Total Variation methods for restoration of vector-valued images. IEEE Trans. Image Process., 7:304-309, 1998.
Giacomo Boracchi and Alessandro Foi. Modeling the performance of image restoration from motion blur. IEEE Transactions on Image Processing, 21(8):3502-3511 2012.
S. Boyd and L. Vandenberghe. Convex Optimization. Kluwer Academic Publishers, 2004.
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Now Publishers, 2011.
Alfred M Bruckstein, David L Donoho, and Michael Elad. From sparse solutions of systems of equations to sparse modeling of signals and images. SIAM review, 51(1):44-81 2009.
Jianrui Cai, Hui Zeng, Hongwei Yong. Zisheng Cao, and Lei Zhang. Toward real-worke single image super-resolution: A new benchmark and a new model. In Proceedings of the IEEE International Conference Computer Vision, 2019.
Emmanuel J. Candes and Michael B. Wakin. An introduction to compressive sampling. IEEE Signal Processing Magazine, 25(2):21-30, 2008
Emmanuel J Candes, Michael B Wakin, and Stephen P Boyd. Enhancing sparsity by reweighted & minimization. Journal of Fourier analysis and applications, 14(5):877-905, 2008
Rick Chartrand. Exact reconstruction of sparse signals via nonconvex minimization. IEEE Signal Processing Letters, 14(10):707-710, 2007.
Ingrid Daubechies, Ronald Verore, Massimo Fornasier, and Sinan Güntürk. Iteratively reweighted least squares minimization for sparse recovery. Communications on Pure and Ap- plied Mathematics, 63(1):1-38, 2010.
Jiangxin Dong, Stefan Roth, and Bernt Schiele. Deep wiener deconvolution: Wiener meets deep learning for image deblurring. Advances in Neural Information Processing Systems, 33:1044 10:4 1059, 2020.
David I. Donoho. Compressed sensing. IEEE Transactions on information theory, 52(4):1289 1306, 2006.
lain S Duff and Jacko Koster. On algorithms for permuting large entries to the diagonal of a sparse matrix. SIAM Journal Matrix Analysis and Applications, 22(4):973-996, 2001.
Michael Elad. Sparse and redundant representations: from theory to applications in signal and image processing, volume 2. Springer, 2010.
M.A.T. Figueiredo, J.M. Bioucas-Dee Dias, and R.D. Nowak. Majorization Majirizationnmmiiiizztaae minimization algorithms for wavelet-based image restoration. IEEE Trans. Image Process., 16:2980-2999 2007.
10

---

Published as a conference paper at ICLR 2023
Shuhang Gu, Lei Zhang, Wangmeng Zuo, and Xiangchu Feng. Weighted nuclear norm minimization with application to image denoising. Proc. IEEE Int. Conf. Computer Vision and Pattern Recognition, pp. 2862-2869, 2014.
Magnus R Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving. Journal of research of the National Bureau of Standards, 49(6):409, 1952.
Zhanxuan Hu, Feiping Nie, Rong Wang, and Xuelong Li. Low rank regularization: A review. Neural Networks, 136:218-233 2222, 2021.
D. Hunter and K. Lange. A tutorial on MM algorithms. The American Statistician, 58:30-37, 2004.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR (Poster), 2015.
Florian Knoll, Jure Zbontar, Anuroop Sriram, Matthew J Muckley, Mary Bruno, Aaron Defazio, Mare Parente, Krzysztof J Geras, Joe Katsnelson, Hersh Chandarana, et al. fastmri: A publicly available raw k-space and dicom dataset of knee images for accelerated mr image reconstruction using machine learning. Radiology. Artificial intelligence, 2(1), 2020.
Filippos Kokkinos and Stamatios Lefkimmiatis. Iterative joint image demosaicking and denoising using a residual denoising network. IEEE Transactions on Image Processing, 28())41177-4188, 2019.
Jakob Kruse, Carsten Rother, and Uwe Schmidt. Learning to push the limits of efficient fft-based image deconvolution. In Proceedings of the IEEE International Conference Computer Vision, pp. 4586-4594, 2017.
Christian Kümmerle and Claudio M Verdun. A scalable second order method for ill-conditioned matrix completion from few samples. In International Conference on Machine Learning, pp. 5872-5883, PMLR, 2021.
Ming-Jun Lai, Yangyang Xu, and Wotao Yin. Improved iteratively reweighted least squares for unconstrained smoothed l, minimization. SIAM Journal on Numerical Analysis, 51(2):927-957, 2013.
S. Lefkimmiatis, A. Bourquard, and M. Unser. Hessian-based norm regularization for image restora- tion with biomedical applications. IEEE Trans. Image Process., 21(3):883-99 995, 2012.
S. Lefkimmiatis, J. Ward, and M. Unser. Hessian Schatten-nowe regularization for linear inverse problems. IEEE Trans. Image Process., 1265):1873-1888, 2013.
Lefkimmiatis, A. Roussos, P. Maragos, and M. Unser. Structure tensor total variation. SIAM J. Imaging Sci., 8:1090-1122 2015.
A. Levin, Y. Weiss, F. Durand, and W. T. Freeman. Understanding and evaluating blind decon- volution algorithms. In IEEE Conf. Comput. Vision and Patt. Recogn. (CVPR), pp. 1964-1971, 2009.
Jie Liu, Jie Tang. and Gangshan Wu. Residual feature distillation network for lightweight image super-resolution. In European Conference on Computer Vision, pp. 41-55. Springer, 2020.
Lu Liu, Wei Huang, and Di-Rong Chen. Exact minimum rank approximation via schatten p-norm minimization. Journal of Computational and Applied Mathematics, 267:218-22 2222118227711 2014.
Julien Mairal, Francis Bach, and Jean Ponce. Sparse modeling for image and vision processing.
Foundations and Trends@ in Computer Graphics and Vision, 8(2-3):85-283, 2014.
A. W Marshall, I. Olkin, and C Arnold, Inequalities: theory of majorization and applications, volume 143. Springer, 1979.
D. Martin, C. Fowlkes, D. Tal, and J. Malik. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. In Proc. IEEE Int. Conf. Computer Vision, pp. 416-423, 2001.

---

Published as a conference paper at ICLR 2023
Karthik Mohan and Maryam Fazel. Iterative reweighted algorithms for matrix rank minimization. The Journal of Machine Learning Research, 13(1):3441-33 2012.
Vishal Monga, Yuelong Li, and Yonina Eldar. Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing, IEEE Signal Processing Magazine, 38(2):18-44, 2021.
Mila Nikolova and Michael K Ng. Analysis of half-quadratti minimization methods for signal and image recovery. SIAM Journal on Scientific computiing 27(3):937-966) 2005.
Christopher C Paige and Michael A Saunders. Solution of sparse indefinite systems of linear equa- tions. SIAM journal on numerical analysis, 12(4):617-629, 1975.
Jinshan Pan, Zhouchen Lin, Zhixun Su, and Ming- Musng- Yang. Robust kernel estimation with outliers handling for image deblurring. In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1800-2808, 2016.
Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In International conference on machine learning, pp. 1310-1318. PMLR, 2013.
AJ Robinson and Frank Fallside. The utility driven dynamic error propagation network. University of Cambridge Department of Engineering Cambridge, 1987.
Yaniv Romano, Michael Elad, and Peyman Milanfar. The little engine that could: Regularization by denoising (red). SIAM Journal on Imaging Sciences, 10(4): 000:-180448444 2017.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedi- cal image segmentation. In International Conference Medical image computing and computer- assisted intervention, pp. 234-241. Springer, 2015.
L. Rudin, Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. Physica D, 60:259 268, 1992.
R. Shewchuk. An introduction to the conjugate gradient method without the agonizing pain, 1994.
URL BtRhhttp:///et
CS cmu. wdw. arseddd////aaeaeeeeeeeeeeeaeeeee jrs wwap:rssrirssapprr. http://sppppp.p.eeeooo html
Libin Sun. Sunghyun Cho, Jue Wang, and James Hays. Edge-based blur kernel estimation using patch priors. In Proc. IEEE International Conference on Computational Photography, 2013.
Y. Sun, P. Babu, and D. P. Palomar. communications, and machine learning. IEEE Transactions on Signal Processing, 65(3):794 816, 2017.
algorithms in signal processing,
A. Tikhonov. Solution of incorrectly formulated problems and the regularization method. In Sovier Math. Dokl., volume 5, pp. 1035-1038, 1963.
Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky. Instance normalization: The missing in- gredient for fast stylization. arXiv preprint arXiv: 1rXiv:1607.08022 2016.
C. Vogel. Computational Methods for Inverse Problems. SIAM, 2002.
Yuan Xie, Shuhang Gu, Yan Liu, Wangmeng Zuo, Wensheng Zhang, and Lei Zhang. Weighted schatten p-norm minimization for image denoising and background subtraction. IEEE transac- fions image processing, 25(10):4842 4857, 2016.
Christopher Zach, Thomas Pock, and Horst Bischof. A dduliit based approach for realtime TV-L1 optical flow. In Pattern Recognition, pp. 214-223. Springer, 2007.
Kai Zhang, Wangmeng Zuo, Shuhang Gu, and Lei Zhang. Learning deep CNN denoiser prior for image restoration. In Proc. IEEE Int. Conf. Computer Vision and Pattern Recognition, July 2017.
Kai Zhang, Luc Van Gool, and Radu Timofte. Deep unfolding network for image super-resolution. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp-
3217-3226, 2020.
Lei Zhang, Xiaolin Wu, Antoni Buades, and Xin Li. Color demosaicking by local directional in- terpolation and nonlocal adaptive thresholding. Journal of Electronic imaging, 20(2):023016, 2011.
12

---

Published as a conference paper at ICLR 2023
A APPENDIX
A.1
PROOFS
In this section we derive the proofs for the tight upper-bounds presented in Lemmal aand Lemma 2 For our proofs use the following inequality for function |rr|", 0 < p S 2, which is given by Sun et al. (017)):
1 1x < 2
P
2 2-P lyl", Va E IR and E \ (0) x² +
(15)
and the Ruhe's trace inequality (Marshall et al. 1979). which is stated in the following theorem:
Theorem 1. Ler A and B be n x positive semidefinite Hermitrian matrices. Then it holds that:
tr (r 2 Σ os (A) (B):
(16)
where oi (A) denotes the i-th singular value of A and the singular values are sorted in a decreasing
order, that is a; (A) > 0x+1 (A).
The proof is straightforward and relies on the inequality of Eq. (5) Specifically, let us consider the positive scalars v xt + 7.
Proof of Lemmal 7
v + 7, with y > If we plug them in 15 we get:
(x.
P ( + 7) 8 =
E +
p (v. + +7)²
(17)
2
2
Multiplying both sides by non-negative scalar wi leads to:
w ( +)) & < J/2 (n? +7) - wi (x) + 7)
2 p w;
(v? + +7)8
(18)
2
The above inequality is closed under summation and, thus, it further holds that:
w. (x² + 7) &
Dep (x: w. p)
(19)
wi (y +7) 222
(x? +7) +
2 2- p
Σ
w. (y² + +7)
2
2
1=1
i=1
P Σ
wi (y) +7)
a
+
P7 2
wi
(y² +7)
n=2
+
p
wi (y² + 7)
i=1
i=1
1=1
x
T W y +
Py
tr (W₂) v)
2
dap (y: w. p), Vx,
(20)
2
+ +7) 4
2
with W, diag
wi (vi
wn (v + + = diag (w) [I (yy + tz)]
By substitution and carrying over the algebraic operations on the r.h.s of Eq. 20 we can show that when æ y, the inequality reduces to equality.
We note that it is possible to derive the IRLS algorithm that minimizes J (x) of Eq. 2 under the weighted & regularizers, without relying on the MM framework. In particular, we can redefine the regularizer &ep (.) as:
n
wi (x) + +7) 2-² 2
3ap (x: w.
(21)
from where the weights Wy y of Lemma | can be inferred. Then, the convergence of the IRLS strategy stationary point can be proven according (Daubechies al. 2010).
13

---

Published as a conference paper at ICLR 2023
Unfortunately, this strategy doesn't seem to apply for the weights Wy of Lemma 2 that the weighted Sp. regularizers don't 't apply directly on the matrix X but instead on its singular values. Specifically, it holds that:
The reason is
(rrww 61r (X; = Σ wi (x) tr (W (X) (XXX +11)) 8).
(22)
Unlike the previous case, here the weights w are included in the matrix W (X), which directly depends on X as: W (X) = (X) diag (w) U (X), with U (X) being the left singular vectors of X. Therefore it is unclear how the approach by Mohan & Fazel (012)) would apply in this case
U
and how the convergence of IRLS can be established.
Proof of Lenuna 2
Let consider the positive scalars (o, + Vo, +
If we plug
them in 15 we get:
? (x)+))
P (o (Y) 2
(o (x)
2 2 (o? (Y) +7)
(23)
Multiplying both sides by a non-negative scalar w0 leads to:
w; (o (X) + +7)² 8 ك P' 2 w. (o (Y) +7) ² (o: +2) +
2 P w. (X) (24)
2
The above inequality is closed under summation and, thus, it further holds that:
6ir 1: p)
w. (o, (x)+
wi (o? (Y) 7) (o (X) +7)
2 p
wi (0, (Y) + +7)$
(Y,(
Σ wi (o? (Y) +7) of (x) - m Σ w; (o, (Y) + +7) p=3
2 2-p
Sir (Y; w.P)
2
2
2
i=1
i=1
r- (wy)o (xx") + py 2
+ 2-p
tr (Wy)
2 Øir (Y: w.w)p p),
(25)
2
where Wy = U diag (w) UT (YYT + 71 and Y admits the singular value decomposition
Y = U diag (o (Y) VT with U E Rmxr, V E Raxr, and r= min (m, n). Further, we show that
holds:
Wy
U diag (w) U (YYT + +イI)
wi (oi (Y) +7) E
= U
UT
w, (0 (Y) +7) 2
w, (o (1) _
(26)
wi (o (Y) +7) 2
= Ô diag (a (Wr)) ÜT -
where Ü-UJJ with J denoting the exchange matrix (row-reversed identity matrix). We note that the vector a (Wy) R similarly to (X) and (Y), holds the singular values of Wy in decreasing order, given that
Wia1 (o (Y) + 7) 2 wi (o? (Y) +7) m Vi.
(27)
14

---

Published as a conference paper at ICLR 2023
This is true because according to the definition of w, it holds w++1 w; Vi, while it also holds (o (V) (o? 7)
mn Vi, since 0+41 (Y) o; and m=? 0.
Finally, given that both WY and X x are positive semidefinite symmetric matrices, we can invoke Ruhe's trace inequality from Theorem
| and combine it with Eq. (25) to get:
Sr (X; u. p) - P ir (Wr + PY tr (Wy)
2-P oer (Y; w. VX Y.
(28)
By substitution and carrying over the algebraic operations the r.h. of Eq. when X Y the inequality reduces to equality.
28 we can show that
A.1.1 THEORETICAL JUSTIHICATION FOR USING AN AUGMENTED MAJORIZER
In Section 2.2 where we consider the solution of the normal equations Eq. (T]. instead of the majorizers that we derived in Eqs. (TO). 2reg. we consider their augmented counterparts which are of the form:
(x:; x²)) = Oreg (2; x*)
+
a |
xx
(29)
2
The reason that, under this choice the system matrix of Eq. (1I) is guaranteed to be non-singular and, thus, a unique solution of the linear system always exists. To verify that this choice doesn't 't compromise the convergence guarantees of our IRLS approach, we note that Ô (æ; x*) is still valid majorizer and satisfies both properties of Eq. (7). required by the MM framework. Specifically, it is straightforward to show that:
Q (z: x) Qreg (x) and 2 (z*; 2 Qreg Vx, x*
(30)
Finally, we note that the use of the augmented majorizer serves an additional purpose. In particular, due to the term & | x 2 the majorizer ② enforces the IRLS estimates between two successive IRLS iterations, x² and =+1 not to differ significantly. Both the unique solution the linear system and the closeness of the the successive IRLS estimates play important role for the stability of the training stage of our LIRLS networks.
2.
A.2 MATRIX EQUILIBRATION PRECONDITIONING
During the training and inference of LIRLS, both the network parameters as well the samples in the input batches vary significantly. This results in a convergence behavior that is not consistent, which mostly attributed to the varying convergence rate of the linear solver at each IRLS step. Indeed, it turns out that the main term S of system matrix, defined in (I. in certain cases can be poorly conditioned. To deal with this issue and improve the overall convergence of LIRLS we apply preconditioning strategy. In particular, we employ matrix equilibration (Duff & Koster, (001) such that the resulting preconditioned matrix has a unit diagonal, while its off-diagonal entries are not greater than 1 in magnitude. In our case all the components that form the system matrix are given in operator form, and thus we do not have access to the individual matrix elements of S S* aI. For this reason, describe below the practical technique of forming diagonal matrix preconditioner that equilibrates the matrix S E c
We start by noting that such matrix can be decomposed as:
A Sp-onwi/20 Val
S =
[AT
GT vall
B
(31)
where B E [ with n1-c We further note simple fact, that any diagonal matrix D diag (d) multiplied with B from the right (BD) equally scales all the elements of each column B.j with corresponding diagonal element dj. and the same diagonal matrix multiplied with TT from the left (DBT) scales the same way each matrix row of BT. Let diagonal elements d, of matrix D to be the form d, = B.jlBl,ll i.c the inverse of the & norm of the corresponding column of the matrix B. In this case, the product B D results in matrix
select the
15

---

Published as a conference paper at ICLR 2023
with normalized columns, while the product DB results in matrix with normalized rows. The product of the two, f.e., the preconditioned matrix D B BD, becomes equilibrated since it has a unit diagonal and all of its non-diagonal elements are smaller or equal to 1. The task then becomes to develop a simple way of calculating the vector that holds the norms of the matrix B columns.
From the definition of matrix B in Eq (1 the squared norm of its j-th column be computed as:
By 12 Σ
11
12 W
2 G 1
||
B -
Σ 2 + p
1/2
+a.
(32)
All restoration problems under study are large-scale, meaning the matrices A and G and W are structured, but only available in an operator form. Depending on the problem at hand, A is either valid convolution matrix (deblurring), strided valid convolution matrix (super-resolution) diagonal binary matrix (demosaicing) or orthonormal subsampled FFT matrix (MRI reconstruction), G is block matrix with valid convolution matrices as blocks, while the matrix W is either a diagonal one, for the case of sparse promoting priors, or a block diagonal matrix with each block being matrix of dimensions e x c, for the case of low-rank promoting priors.
2 and [wI/G]] appear- We utilize the following trick in order to calculate both terms Σ
ing in Eq. 32 In particular, considering any arbitrary matrix C' and a vector of ones 1 we note that Ci ([co2]
T 1),- where C C C' is the Hadamard square operation and denotes
the Hadamard product. The computation of the Haramard square for operators A is straightforward for possibly strided valid convolution matrices (the following also holds for convolutions with pe- riodic and zero boundaries) used in deblurring and super-resolution problems the Hadamard square operator can be obtained by squaring element-wise all the elements of the convolution kernel, for the demosaicing problem we have A = A, and for the MRI reconstruction with subsampled orthonormal FFT matrix, we compute the squared norms of the columns directly, as it is equal to the corresponding acceleration (sampling) rate. Since the matrix W has a the specific structure de- scribed above, and the convolution filter bank operator G is applied independently for each color channel, it is straightforward to show, that for both sparse and low-rank cases the following holds:
w G LJj
(((mirG))
(T)
02 T
1)
02 - (G=((1
(33)
As already discussed above, Go2 can be obtained by squaring element-wise all the elements of the convolution filter bank. The computation of (W1/2) is trivial for the sparse case where W is a diagonal matrix for which (w1/2)0² c2 = W. For the low-rank case we construct W from its cigendecomppsstionn so all of its eigenvalues and eigenvectors are at hand, meaning that we can easily obtain 1/22 by computing the square root of eigenvalues of W and composing them with is computed easily by squaring element-wise all of the elements
its eigenvectors. Then, (W1/2) o2 of 1/2
A.3 DISCUSSION ON THE PERFORMANCE OF THE LIRLS MODELS
In Section we reported the reconstruction results that different LIRLS models have achieved for the studied reconstruction tasks both for grayscale/single-hhan results we observe that for different recovery problems. the best performance is not always achieved by the same model. In particular, for the grayscale reconstruction tasks we can see that the e and fo.w LIRLS models perform better on average than the & and Pg models. This is a strong indication that the presence of the learned weights w lead to more powerful sparsity-promotin and have a positive impact on the reconstruction quality. Moreover, we also observe that the best performance among all the models is accomplished by the eT, which can be somehow counter- intuitive since in theory the choice of p < 1 should promote sparse solutions better. A possible explanation for this is that the fl" regularizer is a convex one and, thus, is amenable to efficient minimization and I LIRLS will converge to the global minimum of the objective function J (ar) in
and color images. From these
regularizers
16

---

Published as a conference paper at ICLR 2023
Eq. 2 On the other hand, the e regularizer is non-convex, which means that the stationary point reached by LIRLS can be sensitive to the initialization and might far from the global minimum.
Regarding the color reconstruction tasks and the learned low-rank promoting regularizers that have considered, we observe that the best performance on average is achieved by the S and Sp.u LIRLS models with p 1. To better interpret this results, need have in mind that the only convex regularizer out of this family is the S (nuclear norm). Given that the nuclear norm is not as expressive as the rest of the low-rank promoting regularizers, it is expected that this LIRLS model will be the least performing, Also note that the weighted nuclear norm, ST unlike to the e. is non- convex and thus does not benefit from the convex optimization guarantees. Based on the above, and having in mind the results we report in Section 4.2 it turns out that in this case the choice of p < plays a more important role the reconstruction quality than the learned weights wn.
Another important issue that is worth of discussing, is the fact that in this work we have applied searsityppromotig regularization on grayscale reconstruction tasks and low-rank promoting reg- ularization on color recovery tasks. We note that by no means this is strict requirement and it possible to seek for low-rank solutions some transform domain when dealing with grayscale images as in Lefkimmiatis et al. (2012) 2013 2015). or sparse solutions when dealing with color images. Due to space limitations we have not explored this cases in this work, but we plan to include related results in an extended version of this paper.
Algorithm 1: Forward and backward passes of LIRLS networks.
Inputs:
initial solution, y:: degraded image, A: degradation operator Input parameters: 0 = [G, w. p}: network parameters, o o,
Forward Pass
Initialize: k = repeat
1. Compute the feature maps {2*
G, (Z vvec for the low-rank case).
2. Compute the updated W weight matrices based on the current estimate "
Sparse case: W diag (wv,) To (**x** + +11)
2 where
Low-rank case: W. = Ig U diag (w.) U (zi Zt +71)
z U diag (z()) v.at
Find the updated solution t+1 by solving linear system:
1
A p.on E c.w.a.+ar) "A'' + caat).
4. k = until the convergence criterion is satisfied; Return x" " Backward Pass
1. Use xt to compute W = W. (G, x²) following steps and 2 in the Forward Pass. Then use both to define the following auxiliary network with parameters 0:
g (a"
AA E G. (G,z') x' AT
2. Compute (V..9)-* p by solving the linear system V -.000000 P. where p = a. L and C the training loss function.
3. Obtain the gradient V by computing the product Veg
4. Use VoL to update the network's parameters 0 or backpropagate further into their parent leafs.
17

---

Published as a conference paper at ICLR 2023
A.4 ALGORITHMIC IMPLEMENTATION OF LIRLS NETWORKS
In Algorithm | we provide the pseudo-code for the forward and backward passes of the LIRLS models, where we distinguish between the learned low-rank and sparsity promoting scenarios. The gradients in the backward pass can be easily computed using any of the existing autograd libraries.
A.5 EMPIRICAL CONVERGENCE OF LIRLS TO A FIXED POINT
As we have explained the manuscript, relying on Lemmas] and We have managed to find valid quadratic majorizers for both the sparsity- and low-rank promoting regularizers defined in Eq. O Given that these majorizers satisfy all the necessary conditions required by the MM framework, we can safely conclude that the proposed IRLS strategy will converge to fixed point. In this section we provide further empirical evidence which support our theoretical justification. In particular, we have conducted several evaluations of the trained LIRLS models. In the first scenario models for 30 steps for color deblurring and simulated MRI with x4 acceleration on the correspond- ing datasets described in Subsection 4.1 After that, we calculate the mean PSNR and SSIM scores individually for each step across all images in each dataset. We provide the resulting plots in Fig, 6 As we can notice, after approximately 25 iterations, both PSNR and SSIM curves for all different learned regularizers have reached an equillibrium and their values do not change. For comparison reasons we also plot the evolution of PSNR and SSIM for standard TV-based regularizers that exist in the literature. These results provide a strong indication that our LIRLS models indeed reach fixed point, which well aligned with the theory.
LIRLS
Additionally to the previous averaged convergence results, provide some representative examples of convergence to a fixed point per individual images and models. For this reason we have selected images from grayscale and color reseutresolttion benchmarks and provide the inference results of fp." and sp.", respectively in Figs
8 The plots in the top row of each figure depict the evolution
of the relative tolerance rtol by a constant value, and the PSNR score across the number of the performed IRLS iterations. The middle rows show the estimated solutions at specific steps, while the bottom rows include the corre- sponding relative error, i.e. the difference between the current latent estimate and the one from the previous step. The corresponding PSNR values and relative errors are provided for each image. For visualization purposes the images with the difference are normalized by the maximum relative error. From these figures it is clearly evident that the relative error between the current estimate and the one from the previous IRLS iteration gradually decreases and approaches zero. At the same time the value of the objective function approaches stationary point, and the PSNR value starts saturating at the later iterations. Please note that in Fig. S we include results obtained by employing more than the 15 IRLS iterations that used for the comparisons reported in Sec. 4.2 The reason for this, is that our main purpose here is to experimentally demonstrate that our Sp.w LIRLS model indeed converges fixed point and not to increase the computational efficiency of the model.
the value of the objective function J (x) const shifted
18

---

Published as a conference paper at ICLR 2023
Color deblurring
SF
S!"
5.
TVN
VTV
0.925
34
0.900
33
PSNR, dB
0.875
32
SSIM
0.850
31
0.825
30
0.800
29
0.775
0
5
10
15
20
25
10
15
20
25
iterations
iterations
MRI reconstruction
TV-4
TV
30
0.70
29
0.68
PSNR, dB
SSIM
0.60
28
0.64
27
62
26
0.60
10
# iterations
#
iterations
Figure 6: Convergence of LIRLS models fixed point. Top: color deblurring Sun al. (2013)
dataset with 1% noise, bottom: simulated MRI with x4 acceleration and 1% noise benchmark based on Knoll et al. 2020) and described in Subsection 4.]
19

---

Published as a conference paper at ICLR 2023
1720
1500
24.6
1250
PSNR. dB
2000
J(
700
15.1
500
10-1
10
30
100
150
200
30
100
150
200
50
100
150
200
iterations
# iterations
iteracions
Input
Step
Step 2
Step
Step 4
Step
Step
Step 15
PSNR-22.83
PSNR:24.11
PSNR:24.36
PSNR:24.48
PSNR-24.55
PSNR-24.58
PSNR:24.61
PSNR:24.66
rlol:1.55 10- 1o0.8.9-100 Hol:5.40
1ol3376-100
10.22.6 10-
nol:2.05 10-3
not:3.50-10
Figure 7: Demonstration of convergence to fixed point of the the f: LIRLS model for the task of grayscale super-resolution. The input corresponds to a synthetically downscaled image by scale factor of 3 with 1 % noise. The image is taken from the BSDI0ORK dataset.
20

---

Published as a conference paper at ICLR 2023
37.6
1050
(() ronnt
27.4
2000
PSNR. JB
750
15al
3000
200
06.8
100
200
200
:00
200
200
410
500
150
200
300
600
500
# 1te alaozus
# iterrations
itezationss
Input
Step
Step 3
Step
Step 7
Step 9
Step 11
Step 50
PSNR:22.59
PSNR:26.67
PSNR:26.80
PSNR:26.94
PSNR-27.04
PSNR:27.11
PSNR:27.17
PSNR:27.55
nll:4425 1144255-10-
10 nol:3.766-10 10 nol:2.77-100 10 mol:2.38 10
101:2.12-10-3
nol:1.87-10" 10- nol:7.46-10 10-1
Figure 8: Demonstration of convergence fixed point of the the Sp.w LIRLS model for the task of color super-resolution. The input corresponds to a synthetically downscaled image by a scale factor of 4 without noise. The image is taken from the BSDI00RK dataset.
21

---

Published as a conference paper at ICLR 2023
A.6 WEIGHT PREDICTION NETWORKS ARCHITECTURES
In the main manuscript we have specified a neural network, whose role is to predict the weights w from some initial solution æp. that will then be used in the weighted Eg and Sy norms during the optimization stage (see Fig. D. This weight prediction network is chosen to be either a lightweight RFDN architecture proposed by Liu al. (2020) for the deblurring. super-resolution and demosaic- ing problems, or a lightweight UNet from Ronneberger et al. (2013) for MRI reconstruction. Below in Table 7 and Table we present detailed per-layer structure of both networks that we have used in our reported experiments.
Table 7: Detailed per-layer structure of the RFDN Liu et al. (020)) weight prediction network (WPN) that was used predict the weights ur of the weighted Ey (eT, (g.w) and weighted S (S1". Spi") norms for the deblurring, super-resolution and demosaicing problems. Here "conv" denotes a convolution layer, "relu" denotes a rectified linear unit function (ReLU), "Irelu" denotes leaky ReLU function with negative slope of 0.05, "sigmoid" denotes sigmoid function, "sc" denotes a skip-connection, "cat" denotes a concatenation along channels dimension, "maxpool" denotes max-pooling operation, "interp" denotes bilinear interpolation, "mul" denotes a point- wise multiplication. For sparsity promoting priors the number of input and output channels is 25 and 24 respectively, while for low-rank promoting priors repeated structures (but not shared weights) are denoted with
and 3, respectively. Blocks with For a more detailed architecture
description refer the codes released by the authors of Liu al. (2220). which we have used
without any modifications: https: github. com uull 1/ R/FDDN
Kernel Size
Input Channels
Output Channels
Block
Layer
Stride
Padding
Bias
conv
x 3
1 1
1 x
25/4
40
True
conv+Irelu
3
40
20
True
conv+sc+hrelu
3 x 3
1 x 1
1 1
40
40
True
conv+Irelu
3
40
20
True
conv+sc+Irelu
3 x 3
1 x 1
1 1
40
40
True
conv+Irelu
x 3
40
True
conv+sc+trelu
3 3
1 X 1 1 X 1 1 1X1
1 1 1x1
40
40
True
Residual
conv+Irelu
x 3
40
20
True
Feature
cat+conv conv conv maxpool
3 x 3
80
40
True
Distillation
X 1
1
40
10
True
Block
x 3
2 2
10
10
True
7 X 7
3 3
conv+relu
3 3
1X1
1x1
10
10
True
conv+relu
3 3
1
1 1x1
10
10
True
conv+interp
3 3
1x1
10
10
True
conv+sc
X 1
1 X
10
10
True
con++sigmoid+mul
10
40
True
Residual
Feature
Distillation
Block
Residual
Feature
Distillation
Block
Residual
Feature
Distillation
Block
cat+conv++reiu conv+sc
1x
1x]
160
40
True
3x3 3
1
40
40
True
conv
3x3
1x1
1X
40
24/3
True
22

---

Published as a conference paper at ICLR 2023
Table Detailed per-layer structure of the U-Net Ronnebbrerrr et al. (2015) weight prediction net- work (WPN) that used to predict the weights of the weighted E norm (ep.") for the MRI reconstruction problem. Here "conv" denotes a convolution layer, "up-conv" denotes transpose convolution layer, "relu" denotes rectified linear unit function (ReLU), "norm" denotes an in- stance normalization layer [Ulannov et al. (2016). "sc" denotes a skip-connection, "cat" denotes a concatenation along channels dimension, "maxpool" denotes max-pooling operation.
Kernel Size
Input
Output
Block
Layer
Stride
Bias
Channels
Channels
con-+norm--elu
3 x 3
1 x
25
25
True
con-+norm+relu
x 3
1 1
25
25
True
maxpool
x 2
2 x 2
Down
con-+norm+-elu
x 3
1 x 1
25
32
True True
conv+norm+relu
3x3 3
1 X 1
32
32
maxpool
2x2 2
2 x 2
Down
conv+norm+relu
3
32
64
True
conv+norm+relu
3x3
1 X 1
64
64
True
maxpool
x 2
2 x 2
Down
conv+norm+relu
3x3 3
1
64
64
True
con-+norm+relu
3 3x3
1 X 1
64
64
True
up-conv+cat
2x2
2 2
64
64
True
Up
conv+norm+relu
3 3
1 X 1
128
32
True
conv+norm+relu
3 3
1
32
32
True
up-conv+cat
2 2
2 2
32
32
True
Up
conv+norm+relu
3 3
1 X 1
64
25
True
conv+norm+relu
3x3
X 1
25
25
True
up-conv+cat
2 2
2 x 2
25
25
True
Up
conv+norm+relu
3 3
1
50
25
True
conv+norm+relu
3 3
1 X
25
25
True
conv
1x1
1 X
25
24
True
23

---

Published as a conference paper at ICLR 2023
A.7 VISUAL RESULTS
In this section we provide additional visual comparisons among the competing methods for all the inverse problems under study.
Input: 20.43
TV- 33.75
TV: 33.89
Zy 34.13
Ey 34,28
c 34.50
MI 34,38
RED: 33.19
IRCNN: 26.06
FDN: 34.04
Target
Figure 9: Visual comparisons among several methods on an optically blurred image from the Levin et al. 2009) dataset. For each reconstructed image its PSNR value is provided in dB.
24

---

Published as a conference paper at ICLR 2023
Input: 24.62
TV-4 30.32
TV: 30.39
Ex: 31.23
22: 31.82
:1.96 31.96
EP 1P
IRCNN: 31.22
FDN: 31.92
Target
Figure 10: Visual comparisons among several methods on a synthetically blurred image from the Sun al. (203)) dataset. For each reconstructed image its PSNR value provided in dB.
25

---

Published as a conference paper at ICLR 2023
Input: 23.54
29.2
TVN: 29.46
Sy 29.90
S 30.99
S 30.95
30 31.9
FDN: 29.33
DWDN: 30.50
Target
Figure 11: Visual comparisons among several methods on a synthetically blurred image from the Sun al. (2033) dataset. For each reconstructed image its PSNR value is provided in dB.
26

---

Published as a conference paper at ICLR 2023
Input (Bicubic):
TV-41
Eg 36.33
C. 36.45
38.55
8P.10: 37.62
RED: 35.73
IRCNN: 34.71
A
A
USRNet:: 38.77
Target
Figure 12: Visual comparisons among several methods on x3 synthetically downscaled image with 1% noise from the BSDI0ORK dataset. For each reconstructed image its PSNR value is provided in dB.
27

---

Published as a conference paper at ICLR 2023
Input (Bicubic): 23.79
VTV: 24.09
TVN: 24.10
24:14.13
S4.24.17
24.11
ST.2 24.20
RED: 24.19
IRCNN: 23.98
USRNet: 23.09
Target
Figure 13: Visual comparisons among several methods on x4 synthetically downscaled image with 1% noise from the BSD100RK dataset. For each reconstructed image its PSNR value is provided in dB.
28

---

Published as a conference paper at ICLR 2023
Bilincar: 29.51
3LV5
TVN: 32.90
S1: 34,82
S 36.11
STc: 35.22
SP 36.02
RED: 34.45
IRCNN: 36.09
Target
Figure 14: Visual comparisons among several methods on mosaicked image without noise. For each demosaicked image its PSNR value is provided in dB.
Bilincar: 28.88
29.5
TVN: 30.08
S:: 29.54
SP. 30.75
Sth: 30.83
S 31.69
RED: 29.53
IRCNN: 31.06
Target
Figure 15: Visual comparisons among several methods on mosaicked image with 3% noise. For each demosaicked image its PSNR value is provided in dB.
29

---

Published as a conference paper at ICLR 2023
Sampling mask (k-sc)
FBP: 22.95
TV. 25.52
TV: 25.99
61 27.77
Eg 28.59
1:0 30.82
g0.50 g 30.30
Target
Figure 16: Visual comparisons among several methods on a simulated MRI with x4 acceleration. For each reconstructed image its PSNR value is provided in dB.
30

---

Published as a conference paper at ICLR 2023
Sampling mask (k (poce)
FBP: 23.69
TV.4;: 25.29
26.85
EG: 27.36
FW: 30.43
29.31
Target
Figure 17: Visual comparisons among several methods on a simulated MRI with x8 acceleration. For each reconstructed image its PSNR value is provided in dB.
31


Reviews:
[{"review_id": "nWi9IqXTXL4", "reviewer": "Reviewer_vB8U", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a new image processing method based on combination of several technics. First the authors introduce a generalization of the iterative reweigthed least square (IRLS) optimization scheme. Based on the majorization-minimization scheme, it solved a classical optimization problem build from a least square loss function with regularizers. Second, they proposed to learn the filters used in the regularization using deep learning method, for that they proposed to use the previous method as a layer of a recurrent network. The results show the effectiveness of the framework on different inverse problems in image processing.  ", "summary_of_the_review": "To sum up the paper, it is good paper that introduce a  image processing method. The two main contribution are the generalization of the IRLS  with a extended set of regularizers and a new neural network layer based on the optimization framework. While some part still need to be clarified, the full method is exposed. The results are among the best, but shows minor improvement compare to state-of-art. ", "strength_and_weaknesses": "# Strength\n\nThe proposed method is a complete framework to solve inverse problem for image processing. It relies on regularizers that are well known to be effective for such problems and leads to pretty good results. The biggest asset of the framework is the construction of a layer which relies on the IRLS with learned filters for regularizations. Using in some way the best of two worlds it leads to a fairly strong neural network that is able to solve inverse problems.\n\n# Weakness\n\nThey are several weakness in the proposed framework. \n1. The image formation model (with a Gaussian additive noise) is very general, but fails on some non-additive noise like Poisson noise in low-luminosity setting. As changing the noise change the loss function, the construction of a surrogate may become delicate.\n2. The IRLS layer asks for solving a optimization problem. Even if conjugate gradient is a very powerful method it takes times to apply, thus with deep network the learning process may become untrackable...\n3. While the results are among the best, both quantitative and qualitative results show little improvements.\n4. The results are computed for the specific level of noise. It would have been interesting to try a blind version, where the noise level is random during the training.\n5. The reproduction of the methods could be difficult (see below).", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "kLGJ5yKlPt9", "reviewer": "Reviewer_gm4f", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "# Summary of The Paper\nThe paper considers the problem of image recovery with sparse or\nlow-rank regularizers. The method considered to solve this problem is\niteratively reweighted least-squares (IRLS), as well as its learned\nvariant. The methods are implemented and applied to several image\nprocessing problems, e.g., super-resolution, demosaicking, and MRI reconstruction.", "summary_of_the_review": "The paper has an interesting idea of unrolling IRLS with several applications to image processing. \n\nHowever, it also has several problems when it comes to related works and technical clarity. In its current form, the paper appears suboptimal for an ICLR publication and a significant revision is needed in my opinion. This is why I intend to give a score of reject.\n", "strength_and_weaknesses": "# Weakness\nThe present manuscript has several drawbacks.\n\n## 1. Insufficient Review of Related Works\nThere are many papers on IRLS with $\\ell_p$\nregularization/minimization, which are overlooked by the paper (except\nDaubechies et al. 2010, Mohan & Fazel 2012). This leads to several\ninaccurate claims and developments, e.g.,\n1. the paper claims that the proposed algorithm is a generalization\n   of IRLS.\n2. the paper claims that the bounds in Lemma 1 and Lemma 2 are novel.\n3. The paper attempts to learn the weights for IRLS in the experiments.\n\nPlease allow me to elaborate on why the two claims are inaccurate.\n1. First it should be noted that the proposed algorithm with $\\ell_p$\n   regularization has been studied in\n   - [1] \"Improved Iteratively Reweighted Least Squares for Unconstrained\n     Smoothed $\\ell_q$ Minimization\" (SIAM Journal on Numerical\n     Analysis, 2013)\n\n   There, Lemma 1 has been proved. In fact, Lemmas 1 and 2 are already\n   known by experts who have read Daubechies et al. (2010), and Mohan &\n   Fazel (2012). And the reweighting strategy of the paper follows\n   directly from Mohan & Fazel (2012).\n2. The reweighting strategy derived by Mohan & Fazel (2012) could be\n   sub-optimal, but the provably optimal weights have been found\n   in a recent paper:\n   - [2] \"A Scalable Second Order Method for Ill-Conditioned Matrix\n     Completion from Few Samples\" (ICML 2021)\n\n   Thus it is a bit unclear to me why the paper proposes to learn\n   the weights.\n3. The paper believes the proposed IRLS variant is a generalization, perhaps\n   because Eq. (11) has an additional augmented Lagrangian term, which\n   is to make a matrix in (11) invertible. However, this modification\n   is problematic for several reasons. First, note that (11) requires\n   knowing the variance of the noise, which is typically unavailable.\n   Anyway, let us suppose the variance is known. The second issue is\n   that, the matrix would be invertible in general. Finally, even if\n   it is not invertible, the least-squares problem has infinitely many\n   solutions, a natural way to handle this is to take the minimum-norm\n   solution, or to project $x^k$ onto the space defined by normal\n   equations. With this mysterious regularization term\n   $\\|x-x^k\\|_2^2$, it is theoretically unclear whether the algorithm\n   will be convergent (even to stationary points), as the\n   least-squares problem is not solved exactly in each iteration.\n   Later we will see this creates a lot of issues.\n\nI would recommend reading more on the recent developments of IRLS\nsince Daubechies et al. (2010), appreciate prior works more, and write\nmore humbly regarding actual contributions.\n\n## 2. Technical Clarity\nThe paper has great potential of confusing the readers. Details are\ngiven below:\n1. In Section 2.1, the reader would get confused at the very\n   beginning, as there is no motivation of having the regularization\n   operator $G$ in (3). Is $G$ given? Where does it come from? Why\n   should we have it there? Can there be any example of $G$?\n\n2. Section 2.2 should undergo a complete revision, and should be\n   turned into a review of IRLS. Note that the\n   majorization-minimization interpretation of IRLS has been known in\n   the literature, and the two lemmas are not new in my opinion. There\n   seems to be no need to emphasize them so much. Instead, this\n   section should have an algorithmic listing for IRLS (which is very\n   easy to understand but now missing) with an explanation of each\n   step. Then write down clearly the objective\n   function that IRLS is minimizing; only after this it makes sense to talk\n   about stationary points of that objective, and cite some papers\n   which prove stationary convergence. (Note that the paper\n   keeps talking about stationary points and local minima, but it is not\n   clear which optimization problem is referred to.)\n3. To my knowledge, Section 3 is the main novelty of the paper. But it\n   is poorly written.\n   - The first question to which the paper should answer intuitively\n     is **why would it be a good idea of learning/unrolling IRLS**? Is it\n     because others are doing this kind of thing for different\n     optimization algorithms and getting cited? Note that IRLS is a well-defined\n     algorithm that converges to some stationary point of the objective\n     that it minimizes (which is not the case for the proposed IRLS\n     variant, as explained above). It could be the case that learning IRLS\n     destroys this guarantee, as the convolution matrix $G$ is\n     changing over time (it remains the same only if the gradient is\n     zero), and as the weight prediction module does not compute the\n     weights in an IRLS way.\n   - In Section 3.1, the way that the paper addresses the challenge of\n     training is not well justified. It has not been proven that the\n     proposed IRLS variant converges to stationary points or some\n     fixed points (which are two different concepts but the paper\n     confused them at Eq 12), hence (12) is mainly a heuristic choice.\n     And it has not been proved that $x^*=f_{\\theta}(x^*;y)$.\n## 3. Experiments\nFirst of all, my impression is that there are too many manually chosen\nfactors that impact the performance, and different experiments have\ndifferent settings. Hence I will not take this section into serious account.\n\nHowever, I do want to point out the importance of ablation study,\nwhich seems missing in the paper.\n\n1. The first thing that the paper should convince me is that, **the learned IRLS\n   performs significantly much better than IRLS**. Without such an experiment, I would\n   not be able to agree that the proposed method is a good idea.\n   Please note that the learned IRLS is significantly more complicated than\n   IRLS when it comes to tuning parameters. And also note that\n   nowadays many of our researchers are using a lot of electricity,\n   money, and trees to train deep networks. We have to avoid training\n   deep networks for simply this reason, unless absolutely necessary.\n\n2. The second experiment the paper should conduct is to **show that, at\n   least empirically, the learned IRLS is actually convergent**, e.g.,\n   the optimality condition holds approximately. This is easy to\n   detect and plot. Otherwise, no such claim should be made in the paper.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "HoVWunIw_wm", "reviewer": "Reviewer_FciL", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This work introduces an algorithm for incorporating a learned sparse and low-rank penalty for image recovery tasks. A majorization bound of the penality functions (as quadratic functions) are derived for efficient optimization by Iteratively Reweighted Least\nSquares (IRLS). Due to the recurrent nature of the algorithm, the authors propose to learn the prior function $G_i$, the weights and $p$ by \na recurrent network, while the other parts follow the IRLS framework. The proposed algorithm is shown to achieve competitive performance on several image recovery tasks.", "summary_of_the_review": "Overall, this paper has made a solid contribution to the area of deep image recovery. It uses a novel combination of unrolling networks and the majorization minimization framework to efficiently incorporate sparsity and low-rank prior. In the meantime, I encourage the authors to include more discussions on the choices of regularizers and $p$. ", "strength_and_weaknesses": "Strength:\n\n1. The idea of the paper is interesting and very natural. Many existing unrolling networks require a large number of parameters and data, whereas the proposed network only learns a few number of functions as components of IRLS framework. This is a nice and natural way to combine neural networks with a regularized least squares, where the regularizers lie in some restricted spaces that can be efficiently represented by neural networks. Moreover, the majorizors that the authors developed enables incorporating sparse or low rank prior for deep image recovery by combining IRLS with an \"unrolling\" strategy. \n\n2. The paper is overall well-written. The idea is clearly motivated and easy to follow. I really enjoyed reading this paper.\n\n3. The experiments look impressive. The proposed method achieves top performance on a variety of datasets for different image recovery tasks.\n\nWeaknesses and questions:\n\n1. There seems to be no regularizer that consistently outperforms others, which makes sense. However, it raises another question that how to choose the \"searching space\" of regularizers for different datasets. Could author(s) add some discussions on whether fixing $p$ or not, and whether promote sparsity or low-rankness? The authors should also explain the reason why sometimes fixing $p=1$ (instead of learning $p$) can be more advantageous, which is somewhat counterintuitive (is that a case of overfitting?).\n\n2. In some image recovery and matrix recovery problems, the ground truth could be both low rank and sparse. Is there anyway to jointly emphasize the sparsity and low-rankness?\n\n3. In the literature review, the authors missed some related works in unrolling networks, and I name a few:\n\n[1] Algorithm Unrolling: Interpretable, Efficient Deep Learning for Signal and Image Processing\n\n[2] Deep Algorithm Unrolling for Blind Image Deblurring\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster