
==== Front
Numer Funct Anal Optim
Numer Funct Anal Optim
Numerical Functional Analysis and Optimization
0163-0563
1532-2467
Taylor & Francis

10.1080/01630563.2024.2384849
2384849
Version of Record
Research Article
Research Article
Iteratively Refined Image Reconstruction with Learned Attentive Regularizers
Numerical Functional Analysis and Optimization
https://orcid.org/0000-0002-4791-6388
Pourya Mehrsa a
https://orcid.org/0000-0002-9041-7373
Neumayer Sebastian b
https://orcid.org/0000-0003-1248-2513
Unser Michael a
a Biomedical Imaging Group, EPFL, Lausanne, Switzerland
b Professorship of Inverse Problems, TU Chemnitz, Chemnitz, Germany
CONTACT Sebastian Neumayer sebastian.neumayer@mathematik.tu-chemnitz.de Professorship of Inverse Problems, TU Chemnitz, Chemnitz, Germany.
11 8 2024
2024
11 8 2024
45 7-9 411440
13 2 2024
9 7 2024
22 7 2024
KnowledgeWorks Global Ltd.9 8 2024
© 2024 The Author(s). Published with license by Taylor & Francis Group, LLC
2024
The Author(s)
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.

Abstract

We propose a regularization scheme for image reconstruction that leverages the power of deep learning while hinging on classic sparsity-promoting models. Many deep-learning-based models are hard to interpret and cumbersome to analyze theoretically. In contrast, our scheme is interpretable because it corresponds to the minimization of a series of convex problems. For each problem in the series, a mask is generated based on the previous solution to refine the regularization strength spatially. In this way, the model becomes progressively attentive to the image structure. For the underlying update operator, we prove the existence of a fixed point. As a special case, we investigate a mask generator for which the fixed-point iterations converge to a critical point of an explicit energy functional. In our experiments, we match the performance of state-of-the-art learned variational models for the solution of inverse problems. Additionally, we offer a promising balance between interpretability, theoretical guarantees, reliability, and performance.

KEYWORDS

Convex regularization
data-driven priors
fixed-point equations
inverse problems
majorization minimization
solution-driven models
European Research Council 10.13039/100019180 Swiss National Science Foundation 10.13039/501100001711 200020_219356/1 The research leading to this publication was supported by the European Research Council (ERC) under European Union’s Horizon 2020 (H2020), Grant Agreement - Project No 101020573 FunLearn, by the Swiss National Science Foundation, Grant 200020_219356/1, and the German Research Foundation (DFG) within the SPP2298 under the project number 543939932.
==== Body
pmc1 Introduction

In biomedical imaging [1], including magnetic resonance imaging (MRI) and computed tomography, reconstructions are often achieved via the resolution of an inverse problem. Its task is to recover an unknown signal x∈RN from noisy measurements y=Hx+n∈RM, where H∈RM,N encodes the data-acquisition process and the noise n∈RM accounts for imperfections in this description. From a variational perspective [2], one defines the reconstruction as the solution to the minimization problem (1) arg minx∈RN(E(Hx,y)+λR(x)),

which involves a data-fidelity term E:RM×RM→R≥0 and a regularizer R:RN→R≥0. In (1), the data fidelity term ensures the consistency of the reconstruction with the measurements, while the regularization, whose strength is controlled by λ∈R>0, imposes some regularity constraints (prior information) on the solution.

For a large variety of data-acquisition and noise models, a well-studied zoo of data fidelities E can be found in the literature. While an instance-specific E is natural, it is desirable that the regularizer R is agnostic to H and n and solely depends on the properties of the underlying images. Hence, a regularizer that captures these inherent properties would be of great interest. Attempts can be traced as far back as to the Tikhonov regularization [3], where images are modeled as smooth signals. Later, this approach was outperformed by compressed sensing [4]. Such models either assume that the signal is sparse in some latent space (e.g., wavelet decomposition [5]) or involve a filter-based regularizer R such as the total variation (TV) [4, 6] and its generalizations [7]. These classic signal-processing approaches achieve a baseline performance with the advantage that they provide stability and robustness guarantees [8].

With the emergence of deep-learning techniques for the solution of inverse problems [9], the traditional approaches have been outperformed in many applications. The end-to-end training achieves state-of-the-art performance in terms of quantitative metrics such as the peak signal-to-noise ratio (PSNR). However, such models are often neither interpretable nor trustworthy for sensitive applications such as biomedical imaging [10, 11]. Therefore, a recent line of research [12–15] is focusing on the use of deep learning for the solution of inverse problems within the variational framework (1). There, instead of learning the whole reconstruction pipeline in an end-to-end manner, one only learns the regularizer R. Up to now, these models have relied mostly on deep architectures to parameterize R, which makes an interpretation difficult. To bypass this issue, the authors in [16] have proposed to parameterize the learnable R as (2) R:x↦∑c=1NC〈1N,ψc(Wcx)〉,

with channel-wise data-driven convolutional matrices Wc∈RN,N and ψc(x)=(ψc(xk))k=1N, where the convex and symmetric profiles ψc are members of C≥01,1(R), the space of nonnegative differentiable functions with Lipschitz-continuous derivatives. Based on the architecture (2), the authors of [16] obtain the best performance among known convex regularizers in their experiments. Moreover, (2) has a clear interpretation as a filter-based regularizer. To further improve the reconstruction performance, we need to look beyond convexity. As an extension of the model (2), the authors of [17] have proposed to learn symmetric potentials ψc∈C≥01,1(R) with ψc″≥−ρ a.e., namely ρ-weakly convex ones. This relaxation significantly improves over the convex setting. In particular, it gets close to the performance of the DRUNet-based model [18], which is among the best-performing methods with a (loose) energy interpretation.

1.1 Outline and contribution

First, we introduce the theoretical concepts in Section 2. Then, we establish in Section 3.1 a link between the use of a ρ-weakly convex ψc within (2) and spatially-adaptive regularization [19–22]. To this end, we investigate the regularizer (3) RMMR:x↦∑c=1NC〈1N,ψc(Bc|Wcx|)〉,

where the convolutional matrices Wc∈RN,N and Bc∈R≥0N,N, and the ψc(x)=(ψc(xk))k=1N with concave potentials ψc∈C≥01,1(R≥0) are data-driven. In (3), |·| is applied component-wise to the vector Wcx and Bc is constrained to have normalized rows. As shorthand, we introduce the notations W=(Wc)c=1NC, B=(Bc)c=1NC, and ψ=(ψc)c=1NC. For the regularizer (3), we show in Theorem 3 that the variational problem (1) is guaranteed to have at least one minimizer. To reach the latter, we propose to use the iterative majorization-minimization regularization (MMR) characterized by (4) xk+1∈arg minx∈RN(E(Hx,y)+λRMMR,k(x)),

with initialization x1∈RN and (5) RMMR,k:x↦∑c=1NC〈Λc(xk),|Wcx|〉,

where Λ=(Λc)c=1NC:RN→(R≥0N)NC with Λc(x)=Bc⊤ψc′(Bc|Wcx|). If E is strictly convex and differentiable, then (4), namely the majorized problem at the k-th step, is strictly convex. Its unique minimum can be computed using the forward-backward splitting (FBS) algorithm [23]. Hence, we can rewrite (4) using the associated update operator TΛ,W,y:RN→RN as (6) xk+1=TΛ,W,y(xk).

In Theorem 4, we prove that the iterations (6) converge to a critical point of the underlying problem (1).

In (5), we can interpret Λc as a channel-wise spatial adaption of the regularization strength, which is attentive to image structures. This viewpoint of solution-driven spatial adaptivity [24, 25] serves as a starting point for the generalization of the MMR model in Section 3.2. More precisely, we propose to replace Λ in (5) with a more expressive convolutional neural network Λ˜. For this, we relax the constraints on the activation functions and linear operators in the mask generator Λ associated with (3), see Figure 1. This leads to the solution-adaptive fixed-point iterations (SAFI) as reconstruction scheme, which involves the regularizers (7) RSAFI,k:x↦∑c=1NC〈Λ˜c(xk),|Wcx|〉.

Fig. 1 Mask-generation architecture of the majorization-minimization (top) and the solution-driven (below) setting. Above each arrow, we denote the signal dimension at the corresponding stage.

In (7), Λ˜:RN→([0,1]N)NC is a 3-layer network with Λ˜c(x)=ϕ3,c(B^cϕ2(B˜ϕ1(W˜x))). Its convolutional operators have dimensions W˜∈R(NC·N),N, B˜∈R(NC·N),(NC·N), and B^c∈RN,(NC·N). The activation functions ϕ1(x)=(ϕ1,⌈k/N⌉(xk))k=1NC·N and ϕ2(x)=(ϕ2,⌈k/N⌉(xk))k=1NC·N share linear splines ϕr,c∈C(R) on input blocks of size N. The final activation functions are ϕ3,c(x)=(ϕ3,c(xk))k=1N, where each ϕ3,c∈C(R) is composed of a linear spline and a Sigmoid function. The latter enforces that the entries of each Λ˜c remain in [0,1]. For the regularizers (7), the minimization problem (4) is still strictly convex. Therefore, each update (6) in the pipeline is numerically tractable and gives rise to an update operator TΛ˜,W,y:RN→RN. In Theorem 5, we prove that TΛ˜,W,y:RN→RN admits at least one fixed point. In this relaxed setting, the convergence of the SAFI scheme to a fixed point is encouraged by the use of regularization techniques during training [26]. The parameterization details for the architectures (5) and (7) are given in Section 3.3. The learning of the associated parameters on denoising problems is discussed in Section 4.

Our numerical evaluation for both denoising and MRI reconstruction in Section 5 indicates that the learning of the parameters W, B, and ψ in (3) leads to a reconstruction performance similar to that of the weakly convex model from [17] with ρ = 1. By setting Bc=Id in (5), we obtain a weakly convex regularizer without a bound on ρ as special case. Hence, our theoretical analysis, corroborated by the numerical results, leads to yet another reasonable explanation for the performance gain of the weakly convex model [17] over the convex one [16]. With the more general regularizer (7) associated to SAFI, the performance gets similar to that of [27], despite the much simpler mask generator Λ˜. As in previous works, we observe that the learned regularizers generalize well to the previously unseen inverse problems. Finally, conclusions are drawn in Section 6.

1.2 Relation to previous work

Our regularizer RMMR relies on the architecture (2) from [17], where we add the inner activation |·| and the nonnegative convolutional matrices Bc. The decomposition ψc=μψc,cvx+ψc,ccv is proposed in [17] with a convex ψc,cvx∈ C≥01,1(R), a concave ψc,ccv∈ C≥01,1(R) with (−ρ)≤ψc,ccv″≤0 a.e., and μ∈R≥0. The convex part ψc,cvx of the learned ψc is necessary to maintain differentiability at 0. Since our ψc only takes positive inputs, we do not have this issue and can drop the term ψc,cvx. Further, we relax ρ from 1 to ∞ to fully explore the role of concavity. Given that the experimental results are similar, we do not expect that the inclusion of a convex part ψc,cvx in (3) leads to a significant gain in performance.

For the regularizer RSAFI,k, we use the absolute value |·| instead of non-convex potentials, as proposed in [27, 28]. Hence, the subproblem (4) for each SAFI update is convex and the deployed optimization algorithm converges to a minimizer. This is stronger than the mere convergence to stationary points of [27]. Since we learn W, our RSAFI,k generalizes the data-adaptive total-variation model in [20]. Moreover, in contrast to these approaches, we iteratively refine the mask in RSAFI,k based on xk+1=TΛ˜,W,y(xk). This leads to implicit depth, which is a possible explanation for why complex generators Λ˜ are not required in our framework.

The majorization-minimization (MM) perspective also shows up in [21], which deploys MM iterations to minimize a spatially adaptive model that is similar to [27, 28]. To ensure closed-form solutions for the minimization problems (4), the authors deploy |·|2 as potentials instead of |·| in (7). In contrast to the SAFI approach, their masks Λ˜(xk) for the MM iterations are induced completely by the underlying regularizer.

2 Preliminaries

Throughout this work, X⊆RN denotes a closed convex set.

2.1 Concave functions

A function f:X→R is said to be concave if it satisfies (8) f(αx1+(1−α)x2)≥αf(x1)+(1−α)f(x2),  ∀x1,x2∈X, ∀α∈[0,1].

If X is open and f∈C1(X), then f is concave if and only if its gradient ∇f satisfies (9) 〈∇f(x1)−∇f(x2),x1−x2〉≤0,  ∀x1,x2∈X.

In the special case N = 1, condition (9) simply states that the derivative f′ is non-increasing on X. Another useful property is that any differentiable concave function f is upper-bounded by its first-order Taylor expansion (10) f(x1)≤f(x2)+〈∇f(x2),x1−x2〉,  ∀x1,x2∈X.

A function f:X→R is convex if and only if (−f) is concave.

2.2 Majorization-minimization algorithm

For a deeper exposition to MM algorithms, we refer to [29–31]. Here, we only collect some basic definitions and the core results. For a continuous f:X→R, we investigate the problem (11) arg minx∈Xf(x).

The idea behind MM algorithms is to replace f by a sequence of (approximating) majorizations g(·,xk), xk∈X for which the computation of a (global) minimizer is tractable. A function g:X×X→R is said to be a majorization of f:X→R if it satisfies the upper-bound f(x)≤g(x,xk), ∀x,xk∈X;

and the local tight bound g(xk,xk)=f(xk), ∀xk∈X.

Next, we introduce the formal MM algorithm together with a convergence result [31, 32].

Theorem 1. For a continuous f:X→R with majorization g:X×X→R and a starting point x1∈X, the MM sequence is given by (12) xk+1∈arg minx∈Xg(x,xk),

and the function values f(xk) are non-increasing. If g is continuous, f and every g(·,xk) is continuously differentiable, and the sub-level set {x∈X:f(x)≤f(x1)} is compact, then all accumulation points of {xk}k∈N are critical points of f. Moreover, if the set X*={x:〈∇f(x),z−x〉≥0, ∀z∈X} is a singleton or if X* is discrete and limk→∞||xk+1−xk‖2→0, then the MM iterations {xk}k∈N converge to a critical point of f.

Remark 1. The condition that X* is a singleton is met if f is strongly convex, namely if (f−σ2||·||22) is convex for some σ∈R+. Hence, we get in this setting global convergence guarantees that are similar to those of convex-minimization algorithms.

2.3 Γ-Convergence

Here, we recall the basic concepts of Γ-convergence within our Euclidean framework and refer to [33] for a more detailed exposition. A family of functions {Jk}k∈N with Jk:X→[0,∞] is said to Γ-converge to J:X→[0,∞] if the following two conditions are fulfilled for every x∈X: for all xk→x, it holds that J(x)≤liminfk→∞Jk(xk);

for every x∈X, there is a sequence {xk}k∈N with xk→x and limsupk→∞Jk(xk)≤J(x).

The importance of Γ-convergence is captured by Theorem 2. Recall that a family of functions Jk:X→R is equi-coercive if it is bounded from below by a coercive function.

Theorem 2 (Theorem of Γ-convergence [33]). Let {Jk}k∈N be an equi-coercive family of functions Jk:X→R. If Jk Γ-converges to J, then it holds that the optimal function values converge limk→∞infx∈XJk(x)=infx∈XJ(x);

all accumulation points of the minimizers of Jk are minimizers of J.

In particular, if all the Jk and J have unique minimizers, then Theorem 2 directly implies convergence of the minimizers of the Jk to the one of J.

3 New perspectives on ridge-based regularization

First, we provide a novel perspective on weakly convex ridge regularizers [17] through the MMR model. Based on this perspective, we then derive our more general SAFI reconstruction scheme.

3.1 Majorization-minimization regularization

For the MMR model, we specify R in the generic problem (1) as (3) and choose E as the squared norm. Moreover, we allow for linear constraints by minimizing over a closed convex polytope X⊂RN. This leads to the problem (13) arg minx∈Xf(x):=(12‖Hx−y‖22+λ∑c=1NC〈1N,ψc(Bc|Wcx|)〉).

First, we establish the existence of minimizers for (13).

Theorem 3. Let ψc:R→R≥0, c=1,…,NC, be continuous and piecewise-polynomial functions with finitely many pieces, and let X⊂RN be a closed convex polytope. Then, problem (13) admits a minimizer.

Proof. Each ψc partitions R into finitely many closed1 intervals (Icm)m=1Lc on which it is a polynomial. Hence, if we denote the n-th row of Bc by Bc,n, each ψc(Bc,n|Wc·|) partitions X into Lc closed unions of polytopes Ωc,nm={x∈RN:Bc,n|Wcx|∈Icm}. Based on these, we can further partition X into finitely many closed polytopes, each of which is contained in one of the ∩c,n=1NC,NΩc,nmc,n, where mc,n∈{1,…,Lc}, and on which all the Bc,n|Wc·| are linear. The infimum in (13) is the infimum of f on (at least) one of these polytopes, say P.

Now, we pick a minimizing sequence (xk)k∈N⊂P. Due to the coercivity of ‖·‖22, we get that the sequence (Hxk)k∈N remains bounded. By construction, there exist diagonal matrices Dc∈RN,N such that Bc,n|Wcx|=Bc,nDcWcx for every n=1,…,N and x∈P. Let M be the matrix which is the vertical concatenation of H and all the Bc,nDcWc with c and n such that (Bc,nDcWcxk)k∈N remains bounded. Since the sequence (Mxk)k∈N is bounded, we can extract a convergent subsequence with limit u∈ran(M). The associated set (14) Q={x∈RN:Mx=u}={M†u}+ker(M)

is a closed polytope. It holds that (15) dist(xk,Q)=dist(M†Mxk+Pker(M)(xk),Q)≤dist(M†Mxk,M†u)→0

as k→+∞ and, thus, that dist(P,Q)=0. The distance of P and Q is 0 if and only if P∩Q≠∅ [34, Theorem 1]. For the Bc,nDcWc that were not added to M, it holds that Bc,n|Wcxk|→∞. Hence, the interval Icmc,n has to be unbounded. Since ψc is a nonnegative polynomial on it, ψc(Bc,n|Wc·|) has to be constant2 on P and ψc(Bc,n|Wcxk|)=ψc(Bc,n|Wcx|) for every x∈P∩Q. Hence, any x∈P∩Q is a minimizer for (13). □

Remark 2. A crucial ingredient for our proof is the architecture (3) with |·| as the inner nonlinearity. In general, it is much harder to guarantee the existence of minimizers for piecewise-polynomial functions [35].

The f in (13) is not necessarily convex. Hence, one should not attempt to solve (13) using conventional convex-optimization algorithms. Instead, one can use the majorization-minimization (MM) algorithm defined in (12) to search for stationary points. When f is convex, this algorithm converges to a minimizer. To apply the MM algorithm, we first show that the concavity of the ψc∈C≥01,1(R) implies the concavity of gc(x)=〈1N,ψc(Bcx)〉. Based on this property, we then construct a majorization of RMMR.

Lemma 1. If ψc∈ C≥01,1(R),  c=1,…,NC, then gc is differentiable with ∇gc(x)=Bc⊤ψc′(Bx). Moreover, if ψc is also concave, then gc is concave as well.

Proof. We have that gc(x)=h(ψc(Bcx)) with h(x)=〈1N,x〉. Hence, the Jacobian Jg is given through the chain rule as (16) Jgc(x)=Jh°ψc°Bc(x)=Jh(ψc(Bx))Jψc(Bx)B=1N⊤diag(ψc′(Bcx))Bc=ψc′(Bcx)⊤Bc,

where diag:RN→RN,N returns a diagonal matrix whose diagonal elements are the input vector. As ∇g(x)=Jg(x)⊤, the first claim readily follows. Further, it holds for x1,x2∈RN that (17) 〈∇gc(x1)−∇gc(x2),x1−x2〉=〈Bc⊤ψc′(Bcx1)−Bc⊤ψc′(Bcx2),x1−x2〉=〈Bc⊤(ψc′(Bcx1)−ψc′(Bcx2)),x1−x2〉=〈ψc′(Bcx1)−ψc′(Bcx2),Bc(x1−x2)〉=〈ψc′(Bcx1)−ψc′(Bcx2),Bcx1−Bcx2〉≤0,

where the inequality stems from the concavity of ψc. By (9), the gc are concave and the proof is complete. □

Now, we majorize gc using its first-order Taylor expansion, see (10), and get that (18) 〈1N,ψc(Bcx)〉≤〈1N,ψc(Bcxk)〉+〈Bc⊤ψc′(Bxk),x−xk〉,  ∀xk∈RN.

With the change of variables x↦|Wcx| and by summing over all c, we then get for any xk∈RN that (19) RMMR(x)=∑c=1NC〈1N,ψc(Bc|Wcx|)〉≤∑c=1NC〈1N,ψc(Bc|Wcxk|)〉+∑c=1NC〈Λc(xk),|Wcx|−|Wcxk|〉,

where Λc(x)=Bc⊤ψc′(Bc|Wcx|). Regarding the notation from Section 2, we choose (20) g(x,xk)=12‖Hx−y‖22+λ∑c=1NC〈1N,ψc(Bc|Wcxk|)〉+λ∑c=1NC〈Λc(xk),|Wcx|−|Wcxk|〉.

It is easy to verify that the chosen g(x,xk) is a valid majorization of f in (13). Therefore, to compute a stationary point of f based on (12), we have to compute the estimates (21) xk+1∈arg minx∈X(12‖Hx−y‖22+λ∑c=1NC〈Λc(xk),|Wcx|〉).

As ψc′≥0, these majorizations of the original problem can be interpreted as spatially reweighted l1-analysis regularization, where the strength of the convex summands ||Wcxk‖1 is reweighted by Λc(xk). Accordingly, we rewrite the convex problem of (21) in a more compact form as (22) xk+1∈arg minx∈X(12‖Hx−y‖22+λ‖Lkx‖1)  with Lk=[diag(Λc(xk))Wc]c=1NC.

In Algorithm 1, we provide an iterative procedure based on FBS [23, 36] to compute (22). To this end, we choose 12‖H·−y‖22 as the differentiable part of the objective and λ‖Lk·‖1 for the non-differentiable one. The most time-consuming part in Algorithm 1 for generic L is the evaluation of the proximal operator proxαλ‖L·‖1 defined as (23) proxαλ‖L·‖1(z)=arg minw∈X(12‖w−z‖22+αλ‖Lw‖1).

For computational purposes, it is better to consider the dual problem of (23), which we derive as in [37].

Proposition 1. Let ProjX:RN→RN denote the orthogonal projection onto X. If u^ solves the problem (24) arg minu∈RNC·N(12||L⊤u−z‖22−12||ProjX{L⊤u−z}−(L⊤u−z)‖22) subject to ‖u^‖∞≤αλ,

then ProjX{z−L⊤u^} equals (23).

Proof. By duality, we have that (25) αλ‖Lw‖1=maxu{u⊤(Lw):‖u‖∞≤αλ}.

Plugging this into (23) leads to (26) minw∈Xmaxu(12‖w‖22+12‖z‖22+w⊤(L⊤u−z)) subject to ‖u‖∞≤αλ=minw∈Xmaxu(12‖w−(z−L⊤u)‖22−12‖z−L⊤u‖22+12‖z‖22) subject to ‖u‖∞≤αλ.

Now, we can swap the min and max because the objective is convex in w and concave in u [38, Cor. 37.3.2]. Then, we directly get that w=ProjX{z−L⊤u} is optimal for the inner minimization. By removing the constant term 12‖z‖22 and a change of sign, we obtain (24). □

To solve (24), we apply once again FBS with the objective as the differentiable part and the constraints for the non-differentiable one, see Algorithm 2. Note that the subtrahend in (24) is the concatenation of a Moreau envelope with the affine map u↦(L⊤u−z). Hence, its gradient reads (L(L⊤u−z)−LProjX{L⊤u−z}), and the overall gradient of the objective in (24) with respect to u is (27) L(z−L⊤u)−L((z−L⊤u)−ProjX{z−L⊤u})=LProjX{z−L⊤u}.

Our last ingredient is the saturating function clip[κ1,κ2]:RN→RN, which is defined component-wise as (28) [clip[κ1,κ2](a)]k=clip[κ1,κ2](ak)={κ1,ak<κ1ak,κ1≤ak≤κ2κ2,ak>κ2.

Algorithm 1 FBS for solving (22)

1: Input: filter matrix L, previous minimizer x1, current iteration kout

2: Parameters: maximal iteration number KFBS, dynamic tolerance ϵFBS=fϵ,FBS(kout)>0

3: Initialize: t1=1, α=1/||H‖22, x˜1=x1

4: for k = 1 to KFBS do

5: xk+1=Proxαλ||L·‖1(x˜k−αH⊤(Hx˜k−y),kout,k)

6: tk+1=(k+5)/3

7: x˜k+1=xk+1+tk−1tk+1(xk+1−xk);

8: if ‖xk+1−xk‖2<ϵFBS‖xk‖2 then

9: break

10: end if

11: end for

12: return xk+1

Algorithm 2 Computation of proxγ||L·‖1 based on the dual (24) using FBS

1: Input: vector z∈RN, current iteration kout, current iteration kFBS

2: Parameters: maximal iteration number Kprox, dynamic tolerance ϵprox=fϵ,prox(kout,kFBS)>0

3: Initialize: u1=Lz, v1=Lz, x1=ProjX{z−L⊤u1}, t1=1, α=1/||L‖22

4: for k = 1 to Kprox do

5: uk+1=clip[−γ,γ](vk−αLProjX{L⊤vk−z})

6: tk+1=(k+5)/3

7: vk+1=uk+1+tk−1tk+1(uk+1−uk)

8: xk+1=ProjX{z−L⊤uk+1}

9: if ‖xk+1−xk‖2<ϵprox‖xk‖2 then

10: break

11: end if

12: end for

13: return xk+1

Algorithm 3 MMR scheme for (13)

1: Parameters: maximal iteration number Kout, tolerance ϵout>0

2: Initialize: x1=0,L1=[Wc]c=1NC

3: for k = 1 to Kout do

4: xk+1=FBS(Lk,xk,k)

5: Compute Lk+1=[diag(Λc(xk+1))Wc]c=1NC

6: if ‖xk+1−xk‖2<ϵout‖xk‖2 then

7: break

8: end if

9: end for

10: return xk+1

Our proposed MMR scheme is summarized in Algorithm 3. It deploys the FBS (Algorithm 1) to solve the majorization-minimization problems (22). If H=Id, Algorithm 1 can be terminated after one step. The involved operator Proxαλ||Lk·‖1 is computed using again the FBS (Algorithm 2). For both algorithms, our choice of {tk}k∈N ensures the convergence of the iterates [36]. Under the assumption of infinite precision in the sub-routines, Algorithm 3 finds indeed a critical point of f.

Theorem 4. Assume that the estimates (22) are obtained exactly within Algorithm 3. Then, {f(xk)}k∈N is non-increasing. If H is invertible, then f is coercive and all accumulation points of {xk}k∈N are in the set of critical points (29) X*={x1∈X:〈H⊤(Hx1−y),x2−x1〉   +λ∑c=1NC〈Λc(x1),|Wcx2|−|Wcx1|〉≥0, ∀x2∈X}.

Moreover, if X* is a singleton or if X* is discrete and limk→∞||xk+1−xk‖2→0, then the MM iterates {xk}k∈N converge to a critical point of f.

Proof. First, we introduce the auxiliary variable z∈R≥0NC·N with grouped components zc=|Wcx|∈RN in (13) and investigate the equivalent problem (30) arg minx∈Xminz∈R≥0NC·Nf˜(x,z):=12‖Hx−y‖22+λ∑c=1NC〈1N,ψc(Bczc)〉 subject to zc=|Wcx|.

Then, the majorizations will take the form (31) g˜((x,z),(xk,zk))=12‖Hx−y‖22+λ∑c=1NC〈1N,ψc(Bczk,c)〉+λ∑c=1NC〈Bc⊤ψc′(Bczk,c),z−zk,c〉,

and their minimization subject to zc=|Wcx| leads indeed to (4). Observe that g˜ are continuous. Further, both f˜ and the g˜(·,(xk,zk)) are differentiable. Hence, we can apply Theorem 1 and the claim follows. □

3.2 Solution-adaptive fixed-point iterations

For the MMR model with (4), the mask generator Λ:RN→ (R≥0N)NC allows for a successive spatial adaption of the regularization strength. So far, the architecture of each Λc:RN→R≥0N is motivated by the MMR perspective. One might wonder if a more generic Λ˜:RN→([0,1]N)NC leads to improvements. This leads to the SAFI scheme based on (7). As the masks generated by Λ˜c are nonnegative, the resulting SAFI updates (32) xk+1∈arg minx∈XJ(x,xk):=12‖Hx−y‖22+λ∑c=1NC〈Λ˜c(xk),|Wcx|〉

are minimizers of convex problems. Moreover, if H is invertible, then (32) is a singleton. Hence, this case gives rise to an update operator TΛ,W,y:X→X. For the MMR framework in Section 3.1 with Λc(x)=Bc⊤ψc′(Bc|Wcx|), we are guaranteed that the fixed-point iterates (33) xk+1=TΛ,W,y(xk)

are convergent. Moreover, the resulting fixed point is a critical point of problem (13).

If one is only interested in obtaining convergence of the iterates (33), this choice is overly constraining. For example, convergence can be guaranteed whenever TΛ,W,y is nonexpansive. Any fixed point x* of (33) is a critical point of (32) with xk=x*, and we cannot improve x* by updating the Λ˜c anymore. In contrast to [27], the spatial adaptivity of the SAFI is driven by every estimate (32), and not only by an initial reconstruction based on the data y. For the special case of total-variation regularization, this was therefore coined as solution-driven adaptivity instead of data-driven adaptivity [24, 25]. In general, the updates (33) are unrelated to the critical points of the non-convex minimization problem (34) arg minx∈X(12‖Hx−y‖22+λ∑c=1NC〈 Λ˜c(x),|Wcx|〉),

where one also minimizes over the input of Λ˜c.

Following the discussed ideas, we proceed as outlined in Figure 1 and Section 1, and replace Λc(x) by the richer architecture Λ˜c(x)=ϕ3,c(B^cϕ2(B˜ϕ1(W˜x))). As observed in [27], one expects that Λ˜c dampens the response of the Wc to structure and leaves it unchanged for noise or artifacts. Under some conditions, we can show that TΛ˜,W,y admits indeed at least one fixed point. Hence, the definition of reconstructions as fixed points of the operator TΛ˜,W,y makes sense.

Theorem 5. Let H be invertible and let σmin denote its smallest singular value. Then, TΛ˜,W,y:X→X maps X into a ball centered at 0 with radius 2‖y‖2/σmin. Further, TΛ˜,W,y admits a fixed point.

Proof. First, we investigate the range of TΛ˜,W,y. By definition of TΛ˜,W,y, it holds for any x∈X that (35) 12||HTΛ˜,W,y(x)−y‖22≤J(TΛ˜,W,y(x),x)≤J(0,x)=12‖y‖22.

From this, we conclude that (36) ||TΛ˜,W,y(x)‖2≤1σmin||HTΛ˜,W,y(x)‖2≤2‖y‖2σmin.

For the second part, we want to apply Brouwer’s fixed-point theorem. To this end, we additionally need to prove that TΛ˜,W,y is continuous. Due to Theorem 2, it suffices to check equi-coercivity and the conditions for Γ-convergence of the family J(·,x) parameterized by x∈X. First, note that it holds for any z∈X that (37) σmin24||z‖22≤14||Hz‖22≤12||Hz−y‖22+||y‖22≤J(z,x)+||y‖22,

which implies equi-coercivity. Now, let xk→x* and zk→z*. By the triangle inequality, we get that (38) |J(z*,x*)−J(zk,xk)|≤|J(z*,x*)−J(z*,xk)|+|J(z*,xk)−J(zk,xk)|.

To obtain the liminf inequality, it suffices to prove that the first two terms converge to 0 as k→∞. For the first one, we have that (39) |J(z*,x*)−J(z*,xk)|≤λ∑c=1NC〈|Λ˜c(x*)−Λ˜c(xk)|,|Wcz*|〉→0

because Λc is continuous. For the second one, we have that (40) |J(z*,xk)−J(zk,xk)|≤12|||Hz*−y‖22−||Hzk−y‖22|+λ∑c=1NC〈|Λ˜c(xk)|,||Wcz*|−|Wczk||〉≤12|||Hz*−y‖22−||Hzk−y‖22|+λ∑c=1NC〈1N,|Wcz*−Wczk|〉.

Again, we conclude that this quantity converges to zero. Hence, we have established the liminf inequality. For the limsup inequality, we use the constant recovery sequence xk=x*, for which the claim follows as in (40). In summary, this implies that TΛ˜,W,y is continuous and that a fixed point exists. □

Remark 3. Based on a quasi-variational inequality perspective, the authors of [25] prove the uniqueness of fixed points for certain problems of the form (33). Unfortunately, their assumptions are hard to verify in practice for Λ˜c. Hence, we do not pursue this direction further and only provide a proof of existence.

For finding a fixed point of the SAFI operator TΛ˜,W,y, we propose to use the fixed-point iterations (33) detailed in Algorithm 4. Unfortunately, a proof of convergence for these iterations is highly nontrivial. In practice, we encourage this property by using a random number of iterations for the training of the model, as detailed in Section 4.

Algorithm 4 SAFI scheme for (33)

1: Parameters: maximal iteration number Kout, tolerance ϵout>0

2: Initialize: x1=0,L1=[Wc]c=1NC

3: for k = 1 to Kout do

4: xk+1=FBS(Lk,xk,k)

5: Compute Lk+1=[diag(Λ˜c(xk+1))Wc]c=1NC

6: if ‖xk+1−xk‖2<ϵout‖xk‖2 then

7: break

8: end if

9: end for

10: return xk+1

Imposing Lipschitz constraints on the masks could potentially be helpful for proving the convergence of the fixed-point iterations. Note that, for our simple generator Λ˜, we can efficiently enforce such constraints, as detailed in [39]. For Theorems 4 and 5, we require the invertibility of the forward operator H to define a single-valued update operator TΛ˜,W,y. For its set-valued generalization, which naturally arises if we drop the invertibility assumption, a stability analysis of the defining problem (32) was recently established in [28]. Note that in the single-valued case, such results are often key to establish the existence of fixed points. Independent of any theoretical considerations, we observed a converging behavior of both MMR and SAFI for the compressed-sensing MRI experiment in Section 5, where H is not invertible.

3.3 Parameterization of the learnable parameters

We now provide details of the parameterization for our two solution-adaptive regularizers. The regularization strength λ in (13) and (32) is learnable for the corresponding reconstruction models. For the MMR model (13) from Section 3.1, the remaining parameters are the linear operators {W, B} and the concave potentials in Ψ. For the SAFI problem (32) from Section 3.2, the remaining parameters are the linear operators {W˜, B˜, B^} and the activation functions {ϕ1,ϕ2,ϕ3}. Taking a closer look at Algorithms 1–3, we observe that we actually only need access to Ψ′ and not to Ψ itself. Hence, we directly parameterize the derivatives Ψ′ instead.

3.3.1 Parameterization of linear operators

All linear operators are constructed with the Conv2d module from PyTorch. Here, we only detail the construction for the output dimension NC. More specifically, we decompose each operator into S stacked Conv2d modules; each with NC output channels, a kernel size (ks×ks), and a group size G. This was observed to be more effective than the direct use of a single Conv2d module with a larger kernel size [16, 17]. Here, the group size G controls the potential transfer of information across the different channels. In particular, if G = 1, then each kernel of the kth layer, k∈2,…,S, is convolved with all the ones of the (k−1)th layer. If G=NC, then each kernel is only convolved with the one of its channel.

3.3.2 Constrained linear operators

We impose constraints on some convolution kernels. All the kernels of W and W˜ should have zero mean. To ensure this, let w∈Rks2 contain the vectorized elements of the respective kernel. Then, we can use the parameterization w↦(w−(1⊤w)/ks2), and optimize over unconstrained variables. For B, we impose that the kernel elements are positive and sum to one. Let b∈Rks2 be the vectorized kernel elements. Here, the implementation of the constraint is nonnegative with the parameterization b↦(|b|(1⊤|b|)). Note that |·| is applied element-wise to b.

3.3.3 Learnable activation functions

For the {ϕ1,ϕ2, ϕ3} in Λ˜, we rely on the learnable linear-spline framework introduced in [40]. More precisely, we use a uniform grid centered at 0 with stepsize Δ and 2M+1 points, M∈N, and the B-spline of degree one defined as (41) β1(x)={1−|x|,x∈[−1,1]0,otherwise.

Then, we parameterize each ϕp,c based on the vector dp,c∈R2M+1 of function values at the grid points as (42) ϕp,c(x)={dp,c,1+dp,c,2−dp,c,1Δ(x+MΔ),x∈(−∞,−MΔ)∑k=−MMdp,c,k+M+1β1(x/Δ−k),x∈[−MΔ,MΔ]dp,c,2M+1+dp,c,2M+1−dp,c,2MΔ(x−MΔ),x∈(MΔ,∞).

In particular, ϕp,c is nonlinear on [−MΔ,MΔ] and extrapolated linearly outside of this interval.

3.3.4 Concave potentials

For the MMR model, we parameterize the ψc′, c=1,…,NC, as (43) ψc′(x)=clip[0,1](σc(rcx)),

where σc:R≥0→R are learnable linear splines and rc∈R>0 are learnable scaling constants that adapt the range. To parameterize {rc}c=1NC, we use the nn.Parameter module of PyTorch. To ensure their positivity, we use |rc| instead of rc in the implementation. As σc is only defined on R≥0, we parameterize it with its M + 1 values on the nonnegative part of the grid from (42) denoted by dc∈RM+1. As ψc must be concave, its derivative ψc′ is constrained to be non-increasing on R. This can be achieved by using a non-increasing σc with σc(0)=1. To enforce the condition σc(0)=1, it suffices to fix dc,0=1. Let D∈RM,M+1 be defined via (Ddc)m−1=(dc,m−dc,m−1), m=2,…,M+1. If all elements of Ddc are non-positive, then σc is non-increasing. To directly embed this constraint into the parameterization, we define (44) P↓(dc)=Sclip[−∞,0](Ddc)+1M+1,

where S∈RM+1,M with (Sdc)m=∑k=1m−1dc,k, m=1,…,M+1. By projecting the unconstrained coefficients dc∈RM+1 to P↓(dc), we ensure that the corresponding σc is non-increasing. With the proposed parameterization, the associated concave profiles ψc:R≥0→R≥0 satisfy the following properties. They are piecewise-quadratic, nonnegative, and increasing.

We have that 0≤ψc′(x)≤1 for all x∈R≥0 and ψc′(0)=1.

4 Architecture and training

For our regularizers (5) and (7), we now describe the learning of the parameters detailed in Section 3.3. For both architectures and their respective reconstruction routines Algorithm 3 and 4, we learn them by solving a denoising problem with additive white Gaussian noise of standard deviation σ∈{5/255,15/255,25/255}. Since the training procedure is exactly the same for both architectures, we restrict our discussion to (5).

Let {xm}m=1M with xm∈R40×40 be a set of clean patches from the grayscale BSD500 dataset [41], and let {ym}m=1M={xm+nσm}m=1M be some noisy versions, where nσm is a realization of the noise. For all experiments, we train with M = 238400 patches. In the following, we collect all the learnable parameters of (5) in the variable θ. Given a noisy patch ym, we obtain its denoised version Dθ,σn1,n2,n3(ym) by applying Algorithm 3 with Kout=n1 steps. As discussed in Section 3.1, fixing KFBS=n2=1 for Algorithm 1 suffices to guarantee convergence in the denoising case. For calculating the involved proxγ||L·‖1 based on Algorithm 2, we use KFBS=n3 steps. During the training phase, all the tolerances are set to (−1) to ensure that the maximum number of steps is used. Now, we propose to learn the optimal parameters θ^ in Dθ,σn1,1,n3 based on the empirical risk (45) θ^∈arg minθ∑n1=46∑n3=1012∑m=1M||Dθ,σn1,1,n3(xm+nσm)−xm‖22.

To solve (45), we use the ADAM optimizer [42] with a learning rate of 10−3 and a batch size of 128 patches that are reconstructed for a single pair (n1, n3). This pair is uniformly drawn at random from one of the possible values for each batch. As documented in [26], using a random numbers of iterations has a regularizing effect when unrolling fixed-point iterations. In particular, this prevents the models from getting overfitted to a specific number of iterations. We perform 40 training epochs and reduce the learning rate by a factor of 0.1 at the 5th and 10th epoch. After each epoch, we evaluate the performance of the model on the Set12 validation data, and choose the output model as the one with the best performance. As a result, we obtain the regularization strength in (13) as well as the linear layers and the potentials that appear in (5).

Remark 4. Instead of pursuing an unrolling approach for training, one can also aim to minimize (45) for n1=n3=∞ with implicit-differentiation techniques [43, 44]. However, as already observed in [16], it is usually unnecessary to fully compute the involved fixed points in Dθ,σn1,1,n3(ym) to learn good parameters θ for the regularizer. Moreover, as we have two nested fixed-point problems, namely the problems (22) and (23), this easily gets prohibitively expensive.

4.1 Architecture and initialization

We use the default nn.Conv2D initialization for every linear layer, and initialize λ as 10−4. Below, we discuss the remaining hyperparameters and initializations.

MMR model: For the operators {Wc}c=1NC, we proceed as described in Section 3.3 with NC = 64, S = 2, ks = 7, and G = 1. Further, we force the kernels to be zero-mean. For the modeling of {Bc}c=1NC, we use a linear layer with NC = 64, S = 2, ks = 7, and G = 64. Here, we enforce that the kernels are positive and normalized. For each concave potential ψc, we use 21 gridpoints, which corresponds to M = 20 and Δ=0.05. We initialize the expansion coefficients of the splines with zero, except dc,0=1. Every rc is initially set to one.

SAFI scheme: To model the linear layers {Wc}c=1NC, we choose NC = 64, S = 2, ks = 7, and G = 1. We enforce that the kernels are zero-mean. To model {Wc,1}c=1NC, {Wc,2}c=1NC, and {Wc,3}c=1NC, we use linear layers with NC = 64, S = 1, ks = 7 and G = 1. Further, we use linear splines with no constraints to parameterize {ϕ1,c}c=1NC, {ϕ2,c}c=1NC, and {ϕ3,c}c=1NC, as described in Section 3.3. For this, we use 21 knots, which correspond to M = 10 and Δ=0.1. We initialize all expansion coefficients of the splines with zero.

4.2 Fine tuning

The interpretability of the learned denoiser Dθ,σn1,1,n3 with the small n1 and n2 from the training stage is, however, limited. In particular, we only perform a partial minimization (unrolling) of (13). To remain within our theoretical setup, we need to iterate Algorithms 1–3 until convergence during the evaluation phase. Doing so without modifying the regularization strength λ in (13) has led to over-smoothing in our experiments. Moreover, the training of the model is for denoising only and not necessarily adapted to other inverse problems with H≠Id. To deal with these issues, we propose to deploy (5) with the previously learned parameters for (13) and solely fine-tune λ on a small set of task-specific validation data with a coarse-to-fine grid search, as described in [16]. Hence, we get two different denoisers for our numerical evaluation: first, the unrolled version Dθ,σn1,1,n3, which is exactly what we have trained for, but which is not necessarily a fixed point of (33); second, the exact fixed point Dθ,σ∞,1,∞ with an adapted λ, which uses more iterations and for which our theoretical analysis holds. For the inverse problems, we use the parameters of the denoising models that are trained with σ=15/255. Here, we only have the fixed-point-based reconstruction operator as we do not train for the task.

4.3 Algorithm hyperparameters for evaluation

We aim to iterate Algorithms 1–3 until convergence, namely, up to machine precision. Still, we enforce an upper bound on the number of iterations for all algorithms. This ensures that we always remain within a reasonable computational budget. Independent of H, we set Kprox=500 and Kout=10. For the denoising case with H=Id, we know that KFBS=1 suffices for convergence. There, we use the iteration-dependent tolerance (46) fϵ,prox(kout,kFBS)={10−3(0.01)kout5,kout≤510−5,kout>5.

For H≠Id, we set KFBS=1000 and use the iteration-dependent tolerances (47) fϵ,FBS(kout)={10−3(0.01)kout5,kout≤510−5,kout>5, andfϵ,prox(kout,kFBS)={3ϵFBS(19)kFBS50,kFBS≤50ϵFBS3,kFBS>50.

To summarize, for efficiency, the inner subproblems are solved with lower precision early on, while the precision for the later stages is higher to ensure convergence. This is a common technique to accelerate majorization minimization models [31] and the FBS algorithm [45, 46].

5 Numerical results

First, we present denoising results as this is our training problem. Then, we deploy the regularizers (5) and (7), which we learned for denoising, to a MRI problem without additional training. For this, we need to adapt the λ in (13) and (32) on some (small) validation set. With this task shift, we want to underline the universality of our approach. The code for our experiments is available on GitHub3. In this section, the images of each row in a figure are plotted with the same grayscale.

5.1 Denoising

Before investigating the qualitative behavior of the proposed regularizers (5) and (7), we first compare their quantitative performance with competing learned regularization methods. The achieved PSNR values on the BSD68 test set are given in Table 1. There, we compare our approach with BM3D, which is a popular baseline [47]. We also compare with the WCRR model of [17] and its spatially adaptive extension SARR [27], which both motivated our approach. Finally, we include Prox-DRUNet [18] as a regularizer with a deeper parameterization and some (loose) theoretical guarantees. For the SAFI, we report results for both the training (SAFI5) and the evaluation configuration (SAFI) with the λ adaption. The performance difference between them is negligible. Hence, from now on, we solely use the evaluation configuration. Additionally, we provide a visual denoising comparison for the castle image in Figure 2. Here, we observe that our SAFI scheme recovers the tip of the spire, which is in general hard to achieve for σ=25/255. The spatial adaptivity helps to preserve sharp edges in the image. Still, all but the Prox-DRUNet method tend to slightly smooth the image.

Fig. 2 Denoising of the castle image corrupted by additive white Gaussian noise with σ=25/255.

Table 1 Denoising performance (in terms of PSNR) on the BSD68 test set.

Method	BM3D [47]	WCRR [17]	MMR	SARR [27]	SAFI	SAFI5	Prox-DRUNet [18]	
σ=5/255	37.54	37.65	37.67	37.84	37.90	37.91	37.97	
σ=15/255	31.13	31.20	31.05	31.55	31.56	31.60	31.70	
σ=25/255	28.61	28.68	28.62	29.07	29.05	29.10	29.18	
The average standard deviation of the PSNR for each image (based on 5 reconstructions) is similar for all settings and is roughly 0.02.

Regarding the qualitative behavior, we provide a solution path for MMR and SAFI in Figures 3 and 4, respectively. Somewhat surprisingly, the algorithm outputs a blurred reconstruction after the first step, in which all the noise is removed at the onset. This initial reconstruction is then progressively sharpened throughout the remaining iterations. This behavior is particularly striking as SAFI still recovers the tip of the spire, see Figure 2, which only reemerges in the later iterations. This is only possible since we update the mask iteratively based on the previous reconstruction, which is not the case for the one step method SARR. As guaranteed by Theorem 4, the residuals along the path for MMR in Figure 3 become small. The same is the case for SAFI where the iterates seem to converge to a fixed point, which necessarily exists due to Theorem 5. We observe the same converging behavior for all images in the BSD68 test set. Also, the visual behavior along the path is very similar in terms of an initial strong smoothing followed by a later recovery of sharp features.

Fig. 3 Solution path of the MMR method for denoising with σ=25/255. Each image (k, ek, PSNRk) represents xk+1 at the kth step of Algorithm 3, with relative error ek=‖xk+1−xk‖2‖xk‖2.

Fig. 4 Solution path of the SAFI scheme for denoising with σ=25/255. Each image (k, ek, PSNRk) represents xk+1 at the kth step of Algorithm 4, with relative error ek=‖xk+1−xk‖2‖xk‖2.

Now, we provide some intuition for the superiority of our approach over its nonadaptive counterpart. If ‖W·‖1 is a well-performing regularizer, the Wc should not respond to the distinctive properties of an image. To investigate this for both MMR and SAFI, we display the response of the respective ∑c=1NC|Wc·| to the noisy image y in Figure 5. For both cases, the structure of the image is also triggered in addition to the noise. This leaves some room for improvements of the reconstruction results. In particular, we can dampen this undesirable response using the masks. Then, the effect of image structure on the regularization cost becomes less pronounced. In Figure 6, we see how the masks become progressively more attentive to the image structure. Overall, the richer parameterization of the mask generator Λ˜ for SAFI captures the image structure better. In particular, the masks for SAFI can still impose a high penalization in the vicinity of edges, whereas this is impossible for the masks from MMR. Overall, this results in a regularizer for which the image structures are less penalized. To conclude, the SAFI scheme leads to a better reconstruction performance than MMR model.

Fig. 5 Masks and responses for the learned regularization architectures (5) and (7). Black corresponds to lower values and white to higher ones. Note that {Wc}c=1NC is learned within the MMR and SAFI frameworks for the first and last three figures (from the left), respectively.

Fig. 6 Evaluation of the masks for MMR (top) and SAFI (bottom). Both models become successively attentive to image structure. Still, the extracted structure in the MMR masks is far less pronounced.

5.2 Magnetic resonance imaging

Now, we deploy the proposed regularizers (5) and (7) to solve MRI-reconstruction problems. We use the single- and 15-coil MRI setups detailed in [16]. For each setup, the ground-truth images consist of proton-density-weighted knee images from the fastMRI dataset [11], both with fat suppression (PDFS) and without fat suppression (PD). In total, this leads to four evaluation tasks. For each task, we use a validation set of ten images to fine-tune the regularization strength λ in (13) and (32), respectively. We then report the test performance of the calibrated models on the remaining fifty test images. To generate the ground-truth image, we use the fully sampled k-space measurements. For the single-coil setup, we generate the measurements through a direct masking of the Fourier measures. In the 15-coil setup, we subsample the Fourier transforms of the ground-truth images multiplied by the respective sensitivity maps. For this, we use the BART [48] implementation of the ESPIRiT algorithm [49]. The subsampling rate of each setup is determined by the acceleration factor Macc with the number of columns kept in the k-space being proportional to 1/Macc. Our single-coil setup is 4-fold (Macc=4) and our multi-coil setup is 8-fold (Macc=8). The measurements are then corrupted with additive white Gaussian noise of standard deviation σ=2·10−3. In Table 2, we provide both the PSNR and structural-similarity index measure (SSIM) values on centered (320×320) patches. Here, we compare against the popular TV regularization, the CRR as a state-of-the-art convex regularizer, its weakly convex extension WCRR, and the Prox-DRUNet as a popular PnP approach. Note that all of these methods are universal in the sense that they can be deployed without additional training. The full implementation details for the CRR and WCRR can be found in the respective papers. For Prox-DRUNet, we deploy the DRS-PnP algorithm proposed in [18], which was previously adapted to our experimental setups in [27].

Table 2 PSNR (first columns) and SSIM (second columns) values for the MRI experiment.

 	4-fold single coil	 	 	 	8-fold multi-coil	 	 	 	
 	PD	PDFS	PD	PDFS	PD	PDFS	PD	PDFS	
Zero-fill (H⊤y)	27.40	29.68	0.729	0.745	23.80	27.19	0.648	0.681	
TV [23]	32.44	32.67	0.833	0.781	32.77	33.38	0.850	0.824	
CRR [16]	33.99	33.75	0.880	0.831	34.29	34.50	0.881	0.852	
WCRR [17]	35.78	34.63	0.899	0.838	35.57	35.16	0.894	0.856	
SARR [27]	36.25	34.77	0.904	0.839	35.98	35.26	0.901	0.858	
Prox-DRUNet [18]	36.20	35.05	0.901	0.847	35.78	35.12	0.894	0.857	
MMR	35.63	34.49	0.896	0.833	35.33	34.97	0.891	0.849	
SAFI	36.43	34.92	0.908	0.844	36.06	35.36	0.901	0.860	

As we observe in Table 2, the MMR model achieves a performance close to that of the weakly convex model introduced in [17]. This underlines again the strong relationship between the two regularization architectures and the associated models. The proposed SAFI regularizer achieves the best performance in three out of the four tasks and is second-best in the other one. Overall, these results indicate that our regularizers (5) and (7) generalize well to inverse problems with the model parameters that were obtained by training on a denoising task. If enough data and compute resources are available, task-specific fine-tuning (second training stage) of the model parameters using the actual data or the forward operator H can help to further increase the performance.

In Figure 7, we provide multi-coil MRI reconstructions for a PD-type image. There, we observe that MMR results in a reconstruction that is sharper than the one with CRR, while the SAFI scheme yields even better results. Most importantly, these improvements do not come at the price of artifacts in the reconstruction. In terms of quantitative metrics, the Prox-DRUNet solution is comparable to SAFI. However, as we observe in the insets, this solution represents poorly the original texture of the images. In particular, it allows for sharp transitions but smooths out the textured parts of the image in this example. In Figure 8, we investigate the single-coil setup for a PDFS image. This is the only case where the Prox-DRUNet is best on average. Although the Prox-DRUNet solution achieves higher PSNR than the SAFI solution, it is hard to observe pronounced visual differences between them.

Fig. 7 Reconstructions for multi-coil MRI (PD). The reported metric is (PSNR, SSIM). The second row contains the zoomed-in insets. The last row shows the squared value of the residuals, which are cutoff at 0.003.

Fig. 8 Reconstructions for single-coil MRI (PDFS). The reported metric is (PSNR, SSIM). The second row contains the zoomed-in insets. The last row shows the squared value of the residuals, which are cutoff at 0.003.

5.3 Algorithmic aspects for MMR and SAFI

5.3.1 Initialization

In principle, the proposed MMR and SAFI reconstructions depend on the initialization of the schemes. The initializations are required to compute the first masks Λ and Λ˜ for MMR and SAFI, respectively. In Algorithms 3 and 4, we initialize with the solution of the nonadaptive convex problems (21) and (32) with Λc=1N and Λ˜c=1N, respectively. These plain experiments are denoted by CVX. To evaluate the sensitivity to this choice, we compare it against two alternatives. First, we perturb the proposed initialization by additive white Gaussian noise with σ=15/255. These noisy experiments are denoted by Perturbed CVX. As an even stronger deviation, we use a random initialization, where each entry is drawn from the standard normal distribution. These challenging experiments are denoted by Random. The PSNR values of the respective reconstructions for the MRI experiment from Figure 8 are given in Table 3. Most importantly, we observe that all variants eventually lead to the same PSNR. This indicates that the fixed point does not depend on the initialization. Unsurprisingly, the convergence to this fixed point occurs faster with a better initialization. We also observe the same behavior for other images. Finally, as indicated in Algorithms 1 and 2, the involved convex subproblems are always initialized with the minimizer of the previous one to accelerate the convergence.

Table 3 Robustness study: PSNR value after each MMR/SAFI update for the multi-coil MRI (PD) reconstruction experiment in Figure 7 depending on the initialization.

Update	0	1	2	3	4	5	6	Final	
MMR, CVX	34.24	35.93	36.30	36.38	36.41	36.43	36.43	36.44	
MMR, Perturbed CVX	24.16	35.10	36.17	36.35	36.41	36.42	36.43	36.43	
MMR, Random	0	34.17	36.08	36.33	36.40	36.42	36.43	36.43	
SAFI, CVX	33.70	35.55	35.84	35.91	35.94	35.95	35.96	35.97	
SAFI, Perturbed CVX	24.10	35.09	35.76	35.90	35.94	35.95	35.96	35.97	
SAFI, Random	0	21.66	34.70	35.59	35.82	35.89	35.93	35.96	

5.3.2 Computational complexity

For the discussed MRI setups, the iterative SAFI approach is on average five times slower than the Prox-DRUNet approach, which does not incorporate any refinement steps. Reconstruction methods with a similar regularization architecture that do not incorporate a mask refinement (such as WCRR) can be even 50 times faster than SAFI. Memory-wise, SAFI has almost 10 times fewer parameters than ProxDRUNet and about 100 times more than WCRR. Since our approach brings valuable insights regarding weakly convex and spatially adaptive regularization, future work should focus on the improvement of the computational effectiveness of the approach. Since the subproblems (21) and (32) are convex, we can choose from a rich pool of methods for this goal. Moreover, we can draw from the literature on accelerating MM iterations [30, 31].

6 Conclusion

We have proposed to use an iterative majorization-minimization regularization (MMR) along with solution-adaptive-fixed-point iterations (SAFI) as new families of data-driven regularizers. They give rise to a sequence of convex reconstruction problems. Numerically, the minimizers associated with this sequence converged to a fixed point in all of our experiments. Overall, this leads to a robust, universal, and interpretable regularization method for inverse problems. A benefit of our simple mask generator for SAFI is that it is well-suited to the enforcement of Lipschitz constraints, which are in turn important to obtain stability estimates. Such constraints might be the key to the proof of the convergence of the fixed point iterations. Finally, it could also be interesting to explore other architectural constraints for generating the masks to obtain theoretical guarantees.

Acknowledgments

The authors thank Alexis Goujon for fruitful discussions.

Disclosure statement

None of the authors have competing interests to disclose.

1 Such a partition with closed interval exists because ψc is continuous.

2 A non-constant polynomial cannot have a finite limit at ±∞.

3 https://github.com/mehrsapo/MMR_SAFI
==== Refs
References

1 McCann, M. T., Unser, M. (2019). Biomedical image reconstruction: From the foundations to deep neural networks. Found. Trends® Signal Process. 13 (3 ):283–359. DOI: 10.1561/2000000101.
2 Scherzer, O., Grasmair, M., Grossauer, H., Haltmeier, M., Lenzen, F. (2009). Variational Methods in Imaging, volume 167 of Applied Mathematical Sciences. New York: Springer.
3 Tikhonov, A. N. (1963). Solution of incorrectly formulated problems and the regularization method. Soviet Math. 4 :1035–1038.
4 Donoho, D. L. (2006). Compressed sensing. IEEE Trans. Inf. Theory 52 (4 ):1289–1306. DOI: 10.1109/TIT.2006.871582.
5 Mallat, S. G. (1989). A theory for multiresolution signal decomposition: The wavelet representation. IEEE Trans. Pattern Anal. Machine Intell. 11 (7 ):674–693. DOI: 10.1109/34.192463.
6 Rudin, L. I., Osher, S., Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms. Phys. D: Nonlinear Phenom. 60 (1–4 ):259–268.
7 Lefkimmiatis, S., Bourquard, A., Unser, M. (2011). Hessian-based norm regularization for image restoration with biomedical applications. IEEE Trans. Image Process. 21 (3 ):983–995. DOI: 10.1109/TIP.2011.2168232.21937351
8 del Aguila Pla, P., Neumayer, S., Unser, M. (2023). Stability of image-reconstruction algorithms. IEEE Trans. Comput. Imag. 9 :1–12. DOI: 10.1109/TCI.2023.3236161.
9 Arridge, S., Maass, P., Öktem, O., Schönlieb, C.-B. (2019). Solving inverse problems using data-driven models. Acta Numerica 28 :1–174. DOI: 10.1017/S0962492919000059.
10 Antun, V., Renna, F., Poon, C., Adcock, B., Hansen, A. C. (2020). On instabilities of deep learning in image reconstruction and the potential costs of AI. Proc. National Acad. Sci. 117 (48 ):30088–30095. DOI: 10.1073/pnas.1907377117.
11 Knoll, F., Zbontar, J., Sriram, A., Muckley, M. J., Bruno, M., Defazio, A., Parente, M., Geras, K. J., Katsnelson, J., Chandarana, H., Zhang, Z., Drozdzalv, M., Romero, A., Rabbat, M., Vincent, P., Pinkerton, J., Wang, D., Yakubova, N., Owens, E., Zitnick, C. L., Recht, M. P., Sodickson, D. K., and Lui, Y. W. (2020). fastMRI: A publicly available raw k-space and DICOM dataset of knee images for accelerated MR image reconstruction using machine learning. Radiol. Artif. Intell. 2 (1 ):e190007. DOI: 10.1148/ryai.2020190007.32076662
12 Duff, M. A. G., Campbell, N. D. F., Ehrhardt, M. J. (2024). Regularising inverse problems with generative machine learning models. J. Math. Imag. Vision 66 (1 ):37–56. DOI: 10.1007/s10851-023-01162-x.
13 Kobler, E., Effland, A., Kunisch, K., Pock, T. (2020). Total deep variation for linear inverse problems. In: 2020 Conference on Computer Vision and Pattern Recognition, Virtual, June 14–19 , pp. 7549–7558.
14 Li, H., Schwab, J., Antholzer, S., Haltmeier, M. (2020). NETT: Solving inverse problems with deep neural networks. Inverse Problems 36 (6 ):065005. DOI: 10.1088/1361-6420/ab6d57.
15 Lunz, S., Öktem, O., Schönlieb, C.-B. (2018). Adversarial regularizers in inverse problems. In: Advances in Neural Information Processing Systems, Vol. 31 .
16 Goujon, A., Neumayer, S., Bohra, P., Ducotterd, S., Unser, M. (2023). A neural-network-based convex regularizer for inverse problems. IEEE Trans. Comput. Imaging 9 :781–795. DOI: 10.1109/TCI.2023.3306100.
17 Goujon, A., Neumayer, S., Unser, M. (2024). Learning weakly convex regularizers for convergent image-reconstruction algorithms. SIAM J. Imaging Sci. 17 (1 ):91–115. DOI: 10.1137/23M1565243.
18 Hurault, S., Leclaire, A., Papadakis, N. (2022). Proximal denoiser for convergent Plug-and-Play optimization with nonconvex regularization. In: 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, Hawaii HI, USA, July 23–29, 2022, pp. 9483–9505.
19 Hintermüller, M., Papafitsoros, K., Rautenberg, C. N. (2017). Analytical aspects of spatially adapted total variation regularisation. J. Math. Anal. Appl. 454 (2 ):891–935. DOI: 10.1016/j.jmaa.2017.05.025.
20 Kofler, A., Altekrüger, F., Antarou Ba, F., Kolbitsch, C., Papoutsellis, E., Schote, D., Sirotenko, C., Zimmermann, F. F., Papafitsoros, K. (2023). Learning regularization parameter-maps for variational image reconstruction using deep neural networks and algorithm unrolling. SIAM J. Imaging Sci. 16 (4 ):2202–2246. DOI: 10.1137/23M1552486.
21 Lefkimmiatis, S., Koshelev, I. S. (2023). Learning sparse and low-rank priors for image recovery via iterative reweighted least squares minimization. In: ICLR.
22 Van Chung, C., De los Reyes, J. C., Schönlieb, C. (2017). Learning optimal spatially-dependent regularization parameters in total variation image denoising. Inverse Probl. 33 (7 ):074005. DOI: 10.1088/1361-6420/33/7/074005.
23 Beck, A., Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Imaging Sci. 2 (1 ):183–202. DOI: 10.1137/080716542.
24 Lenzen, F., Berger, J. (2015). Solution-driven adaptive total variation regularization. In: SSVM . Cham: Springer, pp. 203–215.
25 Lenzen, F., Lellmann, J., Becker, F., Schnörr, C. (2014). Solving quasi-variational inequalities for image restoration with adaptive constraint sets. SIAM J. Imaging Sci. 7 (4 ):2139–2174. DOI: 10.1137/130938347.
26 Anil, C., Pokle, A., Liang, K., Treutlein, J., Wu, Y., Bai, S., Kolter, J. Z., Grosse, R. B. (2022). Path independent equilibrium models can better exploit test-time computation. In: Adv. Neural Inf. Process. Syst., Vol. 35 , pp. 7796–7809.
27 Neumayer, S., Pourya, M., Goujon, A., Unser, M. (2023). Boosting weakly convex ridge regularizers with spatial adaptivity. In: NeurIPS Workshop on Deep Learning and Inverse Problems.
28 Neumayer, S., Altekrüger, F. (2024). Stability of data-dependent ridge-regularization for inverse problems. arXiv:2406.12289.
29 Figueiredo, M. A. T., Bioucas-Dias, J. M., Nowak, R. D. (2007). Majorization–minimization algorithms for wavelet-based image restoration. IEEE Trans. Image Process. 16 (12 ):2980–2991. DOI: 10.1109/tip.2007.909318.18092597
30 Hunter, D. R., Lange, K. (2004). A tutorial on MM algorithms. Amer. Stat. 58 (1 ):30–37. DOI: 10.1198/0003130042836.
31 Sun, Y., Babu, P., Palomar, D. P. (2017). Majorization-minimization algorithms in signal processing, communications, and machine learning. IEEE Trans. Signal Process. 65 (3 ):794–816. DOI: 10.1109/TSP.2016.2601299.
32 Jacobson, M. W., Fessler, J. A. (2007). An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms. IEEE Trans. Image Process. 16 (10 ):2411–2422. DOI: 10.1109/tip.2007.904387.17926925
33 Braides, A. (2002). Γ-Convergence for Beginners, volume 22 of Oxford Lecture Series in Mathematics and Its Applications. Oxford: Oxford University Press.
34 Willner, L. B. (1968). On the distance between polytopes. Q. Appl. Math. 26 (2 ):207–212. DOI: 10.1090/qam/231106.
35 Pham, T.-S. (2023). Tangencies and polynomial optimization. Math. Program. 199 (1–2 ):1239–1272. DOI: 10.1007/s10107-022-01869-6.
36 Chambolle, A., Dossal, C. H. (2015). On the convergence of the iterates of” fista”. J. Optim. Theory Appl. 166 (3 ):25.
37 Beck, A., Teboulle, M. (2009). Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. IEEE Trans. Image Process. 18 (11 ):2419–2434. DOI: 10.1109/TIP.2009.2028250.19635705
38 Rockafellar, R. T. (1997). Convex Analysis. Princeton Landmarks in Mathematics. Princeton, NJ: Princeton University Press.
39 Ducotterd, S., Goujon, A., Bohra, P., Perdios, D., Neumayer, S., and Unser, M. (2024). Improving Lipschitz-constrained neural networks by learning activation functions. J. Mach. Learn. Res. 25 :1–30.
40 Bohra, P., Perdios, D., Goujon, A., Emery, S., Unser, M. (2021). Learning Lipschitz-controlled activation functions in neural networks for Plug-and-Play image reconstruction methods. In: NeurIPS 2021 Workshop on Deep Learning and Inverse Problems, Virtual, December 13, 2021.
41 Arbeláez, P., Maire, M., Fowlkes, C., Malik, J. (2011). Contour detection and hierarchical image segmentation. IEEE Trans. Pattern Anal. Mach. Intell. 33 (5 ):898–916. DOI: 10.1109/TPAMI.2010.161.20733228
42 Kingma, D. P., Ba, J. (2015). Adam: A method for stochastic optimization. In: 3rd International Conference of Representation Learning (ICLR 2015), Poster 9, San Diego CA, USA, May 7–9, 2015.
43 Bai, S., Kolter, J. Z., Koltun, V. (2019). Deep equilibrium models. In: Advances in Neural Information Processing Systems, Vol. 32 .
44 Gilton, D., Ongie, G., Willett, R. (2021). Deep equilibrium architectures for inverse problems in imaging. IEEE Trans. Comput. Imaging 7 :1123–113. DOI: 10.1109/TCI.2021.3118944.
45 Schmidt, M., Roux, N., Bach, F. (2011). Convergence rates of inexact proximal-gradient methods for convex optimization. In: Advances in Neural Information Processing Systems, Vol. 24 .
46 Villa, S., Salzo, S., Baldassarre, L., Verri, A. (2013). Accelerated and inexact forward-backward algorithms. SIAM J. Optim. 23 (3 ):1607–1633. DOI: 10.1137/110844805.
47 Dabov, K., Foi, A., Katkovnik, V., Egiazarian, K. (2007). Image denoising by sparse 3-D transform-domain collaborative filtering. IEEE Trans. Image Process. 16 (8 ):2080–2095. DOI: 10.1109/tip.2007.901238.17688213
48 Uecker, M., Virtue, P., Ong, F., Murphy, M. J., Alley, M. T., Vasanawala, S. S., Lustig, M. (2013). Software toolbox and programming library for compressed sensing and parallel imaging. In: ISMRM Workshop on Data Sampling and Image Reconstruction, pp. 41.
49 Uecker, M., Lai, P., Murphy, M. J., Virtue, P., Elad, M., Pauly, J. M., Vasanawala, S. S., Lustig, M. (2014). ESPIRiT-An eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA. Magnet. Reson. Med. 71 (3 ):990–1001. DOI: 10.1002/mrm.24751.
