\newsiamremark

remarkRemark \newsiamremarkhypothesisHypothesis \newsiamthmclaimClaim \newsiamthmproblemProblem \newsiamthmalgorithm2Algorithm \newsiamremarkexampleExample \headersLearning true monotone operatorsY. Belkouchi, J.-C. Pesquet, A. Repetti and H. Talbot

Learning truly monotone operators with applications to nonlinear inverse problems^†^†thanks: Submitted to the editors DATE.

Younes Belkouchi () CentraleSupelec, Inria, Université Paris-Saclay first.last@centralesupelec.fr Jean-Christophe Pesquet²²footnotemark: 2 Audrey Repetti School of Mathematics and Computer Sciences and School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, UK. Maxwell Institute for Mathematical Sciences, Bayes Centre, Edinburgh, UK . a.repetti@hw.ac.uk Hugues Talbot²²footnotemark: 2

Abstract

This article introduces a novel approach to learning monotone neural networks through a newly defined penalization loss. The proposed method is particularly effective in solving classes of variational problems, specifically monotone inclusion problems, commonly encountered in image processing tasks. The Forward-Backward-Forward (FBF) algorithm is employed to address these problems, offering a solution even when the Lipschitz constant of the neural network is unknown. Notably, the FBF algorithm provides convergence guarantees under the condition that the learned operator is monotone.

Building on plug-and-play methodologies, our objective is to apply these newly learned operators to solving non-linear inverse problems. To achieve this, we initially formulate the problem as a variational inclusion problem. Subsequently, we train a monotone neural network to approximate an operator that may not inherently be monotone. Leveraging the FBF algorithm, we then show simulation examples where the non-linear inverse problem is successfully solved.

keywords:

Monotone operator, Optimization, Inverse Problem, Deep learning, Forward-Backward-Forward, Plug and Play (PnP)

{MSCcodes}

47H05, 47H04, 47H10, 49N45

1 Introduction

In many image processing tasks, the objective is to solve a variational problem of the form

(1)

\underset{\begin{subarray}{c}{x\in C}\end{subarray}}{\text{\rm minimize}}\;\;g% (x)+h(x)

where $C$ is a nonempty closed convex subset of $\mathbb{R}^{n}$ , $g\colon\mathbb{R}^{n}\to\left]-\infty,+\infty\right]$ , and $h\colon\mathbb{R}^{n}\to\left]-\infty,+\infty\right]$ are functions which may have different mathematical properties and various interpretations. This problem appears especially in inverse imaging problems when using Bayesian inference methods to define a maximum a posteriori (MAP) estimate from degraded measurements [8]. Basic choices for $g$ and $h$ are

(2)

(\forall x\in\mathbb{R}^{n})\qquad g(x)=\|Wx\|_{p}^{p},\quad h(x)=\frac{1}{2}% \|Hx-y\|^{2},

where $W\in\mathbb{R}^{q\times n}$ and $H\in\mathbb{R}^{m\times n}$ . Here, vector $y\in\mathbb{R}^{m}$ corresponds to measurements corrupted by an additive white Gaussian noise, $H$ represents a linear model for the acquisition process (e.g., convolution, under-sampled Fourier, or Radon transform), $W$ is a well-chosen linear operator mapping the image $x$ to a transformed domain (e.g., wavelet transform), and $\|\cdot\|_{p}$ is an $\ell_{p}$ norm, with $p\in[1,+\infty]$ . When $p=2$ , (2) allows us to recover a least squares formulation with a Tikhonov regularization [5, 8], whereas $p=1$ promotes sparse solutions to the inverse problem. Both $H$ and $W$ may be unknown, or approximately known, and may require learning strategies to achieve high reconstruction quality.

In the last decades, proximal splitting methods [13, 15] have been extensively used to address convex or non-convex variational large-scale problems that encompass constraints, both differentiable and non-differentiable functions, as well as linear operators. Multiple classes of proximal algorithms have been designed to tackle various forms of optimization problems. Among these algorithms, one can cite first-order algorithms such as the proximal point algorithm, Douglas-Rachford algorithm [12], forward-backward (FB) algorithm (also known as proximal gradient algorithm) [4, 13], or the forward-backward-forward (FBF) algorithm (also known as Tseng’s algorithm) [35]. The first two mentioned algorithms allow us to handle non-smooth functions, through their proximity operators, while the last two algorithms mix gradient steps (i.e., explicit or forward steps) with proximal steps (i.e., implicit or backward steps). While both FB and FBF algorithms are capable of minimizing the sum of a differentiable function and a non-differentiable function, the former requires the differentiable function to have a Lipschitz-continuous gradient, which is not needed when using Tseng’s algorithm.

Lately, proximal algorithms have further gained attention for tackling inverse imaging problems, as their efficiency has been enhanced when paired with powerful neural networks (NNs). These hybrid methods consist in replacing some of the steps in the proximal algorithm by a NN that has been trained for a specific task, leading to the so-called “plug-and-play” (PnP) methods. Traditionally, the intuition behind PnPs was to see the proximity operator as a Gaussian denoiser associated with the regularization term to compute the MAP estimate. Such a denoiser can be hand-crafted (e.g., BM3D) [3] or learned (i.e., NNs) [10, 17, 21, 22, 27, 30]. However, recently, a few studies studies began incorporating different types of NNs into proximal algorithms. These NNs may be trained for various tasks (e.g., inpainting [32]) or used to replace other steps in the algorithm than the proximity operator [2, 7, 21, 20]. For example, the NN could be used to replace the gradient step [21]. In this work, the authors learn a denoiser that is not constrained, and use Proximal Gradient Descent (PGD) as the minimization algorithm. The gradient step is modeled by using a denoising NN, which then allowed them to solve various image restoration tasks by plugging the denoiser into PGD, while having a good performance. Although the method involves a modified PGD with a backtracking step search, it does not yield the true minimum of the considered objective function since the modeled prior is non-convex. As such, converging to a critical point is guaranteed (as is the best case in most non-convex settings). On the other hand, the same authors have proposed to learn Bregman proximal operators, which can simplify some computations in order to handle measurement corrupted with Poisson noise [2, 20]. Finally, the variational problem offers a more general setting, namely maximally monotone operator (MMO) theory, to enlarge the parameter space of the NN in a constrained setting, while still offering global convergence guarantees to an optimal solution. Indeed, most proximal algorithms are grounded on MMO theory, in the sense that convergence of the iterates can be investigated for solving monotone inclusions, instead of variational problems [15]. In this context, the authors in [18, 27, 34] proposed to learn the resolvent of a maximally monotone operator, i.e., a firmly-nonexpansive operator.

Motivated by the MMO approach developed in [27], we propose to generalize (1) to monotone inclusion problems. First, it can be obviously noticed that the constrained optimization problem (1) can be rewritten as

(3)

\underset{\begin{subarray}{c}{x\in\mathbb{R}^{n}}\end{subarray}}{\text{\rm minimize% }}\;\;g(x)+h(x)+\iota_{C}(x)

where $\iota_{C}$ denotes the indicator function of $C$ (see definition at the end of this section). Then, assuming that $g$ and $h$ are proper convex and lower semi-continuous functions, under suitable qualification conditions, $\widehat{x}\in\mathbb{R}^{n}$ is a solution to (3) if and only if it is the solution to the variation inclusion

(4)

0\in\partial g(\widehat{x})+\partial h(\widehat{x})+N_{C}(\widehat{x}),

where $N_{C}=\partial\iota_{C}$ is the normal cone to $C$ . Consequently, solving (4) is a special case of the following monotone inclusion problem where we have specified our working assumptions.

Problem 1.1.

Let $C$ be a closed convex set of $\mathbb{R}^{n}$ with a nonempty interior. Let $h\colon\mathbb{R}^{n}\to\mathbb{R}$ be a proper lower-semicontinuous convex function. Assume that $h$ is continuously differentiable on $C\subset\text{\rm int}\,(\text{\rm dom}\,h)$ . Let $A$ be a monotone continuous operator defined on $C$ . We want to

(5)

\text{find }\widehat{x}\in\mathbb{R}^{n}\text{ such that }0\in A(\widehat{x})+% \partial h(\widehat{x})+N_{C}(\widehat{x}).

While (1) is an instance of Problem 1.1 under suitable conditions, the latter provides a more general formulation, and is not restricted to any variational model. For instance, $A$ does not have to be the gradient of any function. We hence move from the Bayesian formulation of inverse problem to a monotone inclusion framework. Note that the main difference between [27] and this work is that the former approach is based on Minty’s theorem which characterizes a maximally monotone operator through its resolvent. In [27], the focus in terms of modeling and learning was thus put on the resolvent operator of $A$ . In contrast, in this work, we model and learn directly $A$ , a “true” monotone operator. We will see that this paradigm shift will also have consequences for the algorithm choice.

Recently, monotone operator theory has found its way in the study of normalizing flows. These generative methods focus on modeling complex probability distributions using simple and known distributions such as Gaussian ones using neural networks [23, 28]. One of the major properties of these neural networks is invertibility, as most models suppose that a tractable diffeomorphism exists between the mappings involving the variables of both densities, and as such, the inverse model can be computed or estimated using different methods. Invertibility was usually imposed by enforcing the non nullity of the determinant of the Jacobian of the mapping, either through re-parametrization or carefully chosen loss functions [24], recent studies have focused on imposing strong monotonicity as a surrogate to invertibility [1, 9]. Ahn et al. imposed monotonicity using the Cayley operator associated to a given operator. They rely on the property stating that the Cayley operator of a monotone mapping is nonexpansive, hence enforce this condition through spectral normalization and newly defined 1-Lipschitz activations. Chaudhari et al. propose two new neural network architectures, Cascaded Network and Modular Network. Both architectures are monotone, and model the gradient field of a convex function. The main difference with our approach is that we aim to learn any monotone operator, which may or may not be a gradient, and thus their proof of monotonicity can be viewed as a special case of the one provided later in this paper. Moreover, [9] present a sufficient condition for the monotonocity imposed by carefully-crafted layers in the proposed model, while our method applies to any architecture, provided that the learning process converges to an acceptable solution. Lastly, we choose to tackle an original nonlinear inverse problem in the context of image restoration, as an application of Problem 1.1, which constitutes a different problem from normalizing flows. Note that, although a lot of works have been dedicated to linear inverse problems, the investigation of nonlinear degradation models is more scarce [11].

In this article, our primary focus lies in scenarios where the operator $A$ is unknown. To address this challenge, our methodological contribution to the field is two-fold. First, we tailor an existing algorithm to effectively tackle the inclusion Problem 1.1. Second, we introduce a comprehensive framework that harnesses neural networks to learn and replace the monotone operator $A$ , leveraging their capabilities of universal approximation. Our contributions encompass both algorithmic adaptation and the establishment of a novel approach for monotone operator learning through neural networks. To demonstrate the potential of our approach, we define and analyze a non-linear inverse problem which is modeled as a monotone inclusion problem, and we show that we are able to solve it using a learned monotone operator.

This article is structured as follows. In Section 2, we introduce Tseng’s algorithm and adapt it to solve Problem 1.1. In Section 3, we propose a PnP framework based on Tseng’s iterations. To this end, we relate monotonicity of neural networks to the study of their Jacobian, and present the training approach for learning monotone neural networks. In Section 4, we formulate a non-linear inverse problem as a variational inclusion problem, and solve it using the presented tools. Lastly, conclusions are drawn in Section 5.

Notation

$(\mathbb{R}^{n},\lVert\cdot\rVert)$ is a Euclidean space, where $\lVert\cdot\rVert$ is the $\ell_{2}$ norm and $\langle\cdot|\cdot\rangle$ is the associated inner product. Let $S\subset\mathbb{R}^{n}$ , the interior of $S$ is denoted by $\text{\rm int}\,S$ . An operator $A:\mathbb{R}^{n}\rightrightarrows 2^{\mathbb{R}^{n}}$ is a set-valued map if it maps every $x\in\mathbb{R}^{n}$ to a subset $A(x)\subset\mathbb{R}^{n}$ . $A$ is single valued if, for every $x\in\mathbb{R}^{n}$ , $\text{\rm card}\,A(x)=1$ , in which case we consider that $A$ can be identified with an application from $\mathbb{R}^{n}$ to $\mathbb{R}^{n}$ . The graph of the set valued operator $A$ is defined as $\operatorname{Gra}A=\{(x,u)\in(\mathbb{R}^{n})^{2}\mid u\in A(x)\}$ . The reflected operator of $A$ is defined as $R_{A}=2A-\operatorname{Id}\,$ .

The Moreau subdifferential of a convex function $f$ is a set valued operator denoted as $\partial f$ . If $f$ is differentiable, then $\partial f$ is single valued, in which case $\nabla f$ refers to its gradient: $(\forall x\in\mathbb{R}^{n})$ $\partial f(x)=\{\nabla f(x)\}$ . In addition, if $f$ is proper and lower-semicontinuous, the proximity operator of $f$ associates to every $x\in\mathbb{R}^{n}$ the unique point $p\in\mathbb{R}^{n}$ such that $x-p\in\partial f(p)$ . The indicator function $\iota_{C}$ of a subset $C$ of $\mathbb{R}^{n}$ is defined as: $(\forall x\in\mathbb{R}^{n})$ $\iota_{C}(x)=0$ if $x\in C$ , and $\iota_{C}(x)=+\infty$ otherwise. The normal cone to a set $C\subset\mathbb{R}^{n}$ at $x\in\mathbb{R}^{n}$ is defined as $N_{C}(x)=\{u\in\mathbb{R}^{n}\mid(\forall y\in\mathbb{R}^{n}){\left\langle{{y-% x}\mid{u}}\right\rangle}\leqslant 0\}$ if $x\in C$ , and $N_{C}(x)={\varnothing}$ otherwise. In particular, if $C$ is a nonempty closed convex set, $\partial\iota_{C}=N_{C}$ and $\text{\rm prox}_{\iota_{C}}=\text{\rm proj}_{C}$ is the projector onto $C$ .

Let $A\colon\mathbb{R}^{n}\rightrightarrows 2^{\mathbb{R}^{n}}$ . Then $A$ is a monotone (resp. strictly monotone) operator if $(\forall(x,u)\in\operatorname{Gra}A)(\forall(y,v)\in\operatorname{Gra}A)\quad{% \left\langle{{u-v}\mid{x-y}}\right\rangle}\geqslant 0$ (resp. ${\left\langle{{u-v}\mid{x-y}}\right\rangle}>0)$ if $x\neq y$ ). Furthermore, $A$ is a $\beta$ -strongly monotone operator, with constant $\beta>0$ , if $(\forall(x,u)\in\operatorname{Gra}A)(\forall(y,v)\in\operatorname{Gra}A)\quad{% \left\langle{{u-v}\mid{x-y}}\right\rangle}\geqslant\beta\|x-y\|^{2}$ . As a limit case, a monotone operator can be considered as $0$ -strongly monotone. We say that $A$ is maximally monotone if there exists no monotone operator $B\neq A$ such that $\operatorname{Gra}A\subset\operatorname{Gra}B$ . When $f\colon\mathbb{R}^{n}\to]-\infty,+\infty]$ is convex, proper, and lower semi-continuous, then $\partial f$ is maximally monotone. A is $\beta$ -cocoercive with constant $\beta>0$ , if $(\forall(x,u)\in\operatorname{Gra}A)(\forall(y,v)\in\operatorname{Gra}A)\quad{% \left\langle{{u-v}\mid{x-y}}\right\rangle}\geqslant\beta\|u-v\|^{2}$ .

Let $T\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a single valued operator, and $x\in\mathbb{R}^{n}$ . If $T$ is Fréchet differentiable, we denote by $\operatorname{J}_{T}(x)\in\mathbb{R}^{n\times n}$ the Jacobian of $T$ at $x$ which is defined as the matrix that verifies $\lim_{\|h\|\to 0,h\neq 0}\frac{\|T(x+h)-T(x)-\operatorname{J}_{T}(x)h\|}{\|h\|% }=0$ . If $T$ is Fréchet differentiable, then $T$ is Gâteaux differentiable and, for every $h\in\mathbb{R}^{n}$ $\operatorname{J}_{T}(x)h=\lim_{t\to 0,t\neq 0}\frac{T(x+th)-T(x)}{t}$ . Further, we define the symmetric part of its Jacobian operator as $\operatorname{J}^{\rm s}_{T}=(\operatorname{J}_{T}+\operatorname{J}_{T}^{\top}% )/2$ .

Let $M\in\mathbb{R}^{n\times n}$ be a symmetric matrix. We denote by $\lambda_{\min}(M)$ the smallest eigenvalue of $M$ , and by $\overline{\lambda}_{\max}(M)$ the maximum absolute value of the eigenvalues of $M$ . Let $(M_{1},M_{2})\in(\mathbb{R}^{n\times n})^{2}$ . The Loewner order $M_{1}\succ M_{2}$ (resp. $M_{1}\succeq M_{2}$ ) is defined as, for every $u\in\mathbb{R}^{n}\setminus\{0\}$ , $u^{\top}M_{1}u>u^{\top}M_{2}u$ (resp. $u^{\top}M_{1}u\geqslant u^{\top}M_{2}u$ ). Then $M$ is positive definite (resp. positive semidefinite) if and only if $M\succ 0$ (resp. $M\succeq 0$ ). An alternative definition is that $M$ is positive definite (resp. positive semidefinite) if and only if $\lambda_{\min}(M)>0$ (resp. $\lambda_{\min}(M)\geqslant 0$ ).

For further background on convex optimization and MMO theory, we refer the reader to [4].

2 Tseng’s algorithms for monotone inclusion problems

The algorithm we propose in this work is grounded on the original work by Tseng in [35]. In this section, we first recall the relevant background, and then give our modified version of Tseng’s algorithm for solving (5).

2.1 Forward-Backward-Forward strategy

The following monotone inclusion problem is considered in [35].

Problem 2.1.

Let $C$ be a nonempty closed convex subset of $\mathbb{R}^{n}$ . Let $f\colon\mathbb{R}^{n}\to]-\infty,+\infty]$ be a proper, lower-semicontinuous, convex function, and let $B$ be a maximally monotone operator, continuous on $\text{\rm dom}\,\partial f\subset C$ . We want to

(6)

\text{find }\widehat{x}\in C\text{ such that }0\in B(\widehat{x})+\partial f(% \widehat{x}).

Tseng proposed to solve Problem 2.1 using the following algorithm. {algorithm2} Let $x_{0}\in\text{\rm dom}\,B$ and $(\gamma_{k})_{k\in\mathbb{N}}$ be a sequence in $[0,+\infty[$ .

(7)

\begin{array}[]{l}\text{For }k=0,1,\ldots,\\ \left\lfloor\begin{array}[]{l}b_{k}=B(x_{k})\\ z_{k}=\text{\rm prox}_{\gamma_{k}f}(x_{k}-\gamma_{k}b_{k})\\ x_{k+1}=\text{\rm proj}_{C}\big{(}z_{k}-\gamma_{k}\big{(}B(z_{k})-b_{k}\big{)}% \big{)}.\end{array}\right.\end{array}

This method, sometimes also called forward-backward-forward (FBF) algorithm, is reminiscent of extragradient methods. The step-size can be determined using Armijo-Goldestein rule, defined below.

Definition 2.2 (Armijo-Goldstein rule).

Let $\sigma\in\left]0,+\infty\right[$ and $(\beta,\theta)\in]0,1[^{2}$ . Let $(x_{k})_{k\in\mathbb{N}}$ be a sequence generated by Section 2.1. At every iteration $k\in\mathbb{N}$ , the stepsize $\gamma_{k}>0$ is chosen such that $\gamma_{k}=\sigma\beta^{i_{k}}$ where

(8)

i_{k}=\inf\left\{i\in\mathbb{N}\,\left|\,\begin{array}[]{l}\gamma=\sigma\beta^% {i}\\ \gamma\|B(z_{k})-B(x_{k})\|\leqslant\theta\|z_{k}-x_{k}\|\\ z_{k}=\text{\rm prox}_{\gamma f}\big{(}x_{k}-\gamma B(x_{k})\big{)}\end{array}% \right.\right\}.

We have the following convergence result [35, Theorem 3.4].

Proposition 2.3.

Consider Problem 2.1. Let $(x_{k})_{k\in\mathbb{N}}$ be a sequence generated by Section 2.1. Assume that

(i)

there exists a solution to the problem;
(ii)

$B+\partial f$ is maximally monotone;

(iii)

for every $x\in C$ ,

(9)

\varrho(x)=\inf_{w\in B(x)+\partial f(x)}\|w\|

is locally bounded;

(iv)

$(\gamma_{k})_{k\in\mathbb{N}}$ satisfies the rule given by Definition 2.2.

Then $(x_{k})_{k\in\mathbb{N}}$ converges to a solution to the problem.

A strength of this algorithm is that it does not require the operator $B$ to be cocoercive as in the standard FB algorithm. In the particular case when $B$ is $\beta$ -Lipschitzian with $\beta\in\left]0,+\infty\right[$ , the stepsize $\gamma_{k}$ can be chosen such that $\inf_{k\in\mathbb{N}}\gamma_{k}>0$ and $\sup_{k\in\mathbb{N}}\gamma_{k}<1/\beta$ , which allows the use of a constant stepsize [35]. Unfortunately, this implies estimating $\beta$ , which may be difficult. For example, when $B$ is a standard neural network employing nonexpansive activation functions, it is indeed Lipschitzian, but computing a tight estimate of its Lischitz constant in an NP-hard problem [14].

2.2 Modified algorithm

In this section, we propose to use the results presented in Section 2.1 to derive an algorithm for solving Problem 1.1. The proposed algorithm is given below. {algorithm2} Let $x_{0}\in C$ and $(\gamma_{k})_{k\in\mathbb{N}}$ be a sequence in $]0,+\infty[$ .

(10)

\begin{array}[]{l}\text{For }k=0,1,\ldots,\\ \left\lfloor\begin{array}[]{l}a_{k}=A(x_{k})+\nabla h(x_{k})\\ z_{k}=\text{\rm proj}_{C}(x_{k}-\gamma_{k}a_{k})\\ x_{k+1}=\text{\rm proj}_{C}\big{(}z_{k}-\gamma_{k}(A(z_{k})+\nabla h(z_{k})-a_% {k})\big{)}.\end{array}\right.\end{array}

The convergence of Section 2.2 can then be deduced from Proposition 2.3, which yields the following result.

Proposition 2.4.

Consider Problem 1.1. Let $(x_{k})_{k\in\mathbb{N}}$ be generated by Section 2.2. Assume that

(i)

there exists a solution to Problem 1.1;
(ii)

for every $k\in\mathbb{N}$ , $\gamma_{k}$ satisfies the Armijo rule given in Definition 2.2, for $B=A+\nabla h$ , and $f=\iota_{C}$ .

Then $(x_{k})_{k\in\mathbb{N}}$ converges to a solution to Problem 1.1.

Proof 2.5.

According to [4, Theorem 20.21], there exists a maximally monotone extension $\widetilde{A}$ of $A$ on $\mathbb{R}^{n}$ , and (5) has the form of (6) with $f=\iota_{C}$ and $B=\widetilde{A}+\partial h$ . Then $f$ is proper, lower semi-continuous, and convex. Further, since $h$ is also proper, lower semi-continuous, and convex, then $\partial h$ is maximally monotone. Since $\widetilde{A}$ is maximally monotone and $\text{\rm dom}\,\widetilde{A}\cap\text{\rm int}\,(\text{\rm dom}\,\partial h)% \supset\text{\rm dom}\,A\cap\text{\rm int}\,(\text{\rm dom}\,\partial h)=C\cap% \text{\rm int}\,(\text{\rm dom}\,\partial h)\supset C\cap\text{\rm int}\,(% \text{\rm dom}\,h)=C\neq{\varnothing}$ , then $B$ is maximally monotone [4, Corollary 25.5(ii)]. In addition, $B=\nabla h+\widetilde{A}=\nabla h+A$ is continuous on $\text{\rm dom}\,\partial f=C$ .

We have $\text{\rm dom}\,B=\text{\rm dom}\,\widetilde{A}\cap\text{\rm dom}\,(\partial h% )\supset\text{\rm dom}\,A\cap\text{\rm int}\,(\text{\rm dom}\,h)=C$ and thus $\text{\rm int}\,(\text{\rm dom}\,B)\cap\text{\rm dom}\,\partial f\supset\text{% \rm int}\,(C)\neq{\varnothing}$ . Consequently, Assumption (ii) in Proposition 2.3 holds [4, Corollary 25.5(ii)].

Let $\varrho$ be defined on $C$ by (9). We have

	$\displaystyle(\forall x\in C)\quad\varrho(x)$	$\displaystyle=\inf_{t\in N_{C}(x)}\\|B(x)+t\\|$
		$\displaystyle\leqslant\\|B(x)\\|+\inf_{t\in N_{C}(x)}\\|t\\|$
		$\displaystyle\leqslant\\|B(x)\\|+\\|x-\text{\rm proj}_{C}x\\|$
(11)			$\displaystyle=\\|B(x)\\|.$

Since we have seen that $B$ is continuous on $C$ , $\varrho$ is locally bounded on $C$ and Assumption (iii) in Proposition 2.3 is satisfied.

The convergence result thus follows from Proposition 2.3 by noticing that $\text{\rm prox}_{\gamma_{k}f}=\text{\rm proj}_{C}$ .

3 Proposed method using a monotone NN

The objective of this work is to develop a PnP version of the proposed Tseng’s algorithm described in Section 2.2. To this aim, we will approximate the operator $A$ in Problem 1.1 by a monotone neural network $F_{\theta}$ , with learned parameters $\theta$ belonging to some set $\Theta$ . Then, the modified form of Section 2.2 is given below. {algorithm2} Let $x_{0}\in C$ and $(\gamma_{k})_{k\in\mathbb{N}}$ be a sequence in $]0,+\infty[$ .

(12)

\begin{array}[]{l}\text{For }k=0,1,\ldots,\\ \left\lfloor\begin{array}[]{l}a_{k}=F_{\theta}(x_{k})+\nabla h(x_{k})\\ z_{k}=\text{\rm proj}_{C}(x_{k}-\gamma_{k}a_{k})\\ x_{k+1}=\text{\rm proj}_{C}\big{(}z_{k}-\gamma_{k}(F_{\theta}(z_{k})+\nabla h(% z_{k})-a_{k})\big{)}.\end{array}\right.\end{array}

In this section, we will describe how to build such a monotone neural network.

3.1 Properties of differentiable monotone operators

As described in Proposition 2.4, to ensure convergence of a sequence $(x_{k})_{k\in\mathbb{N}}$ generated by Section 3, the NN $F_{\theta}$ must be monotone and continuous on $C$ . As most standard neural networks are continuous, especially those that use non-expansive activation function [27, 14], we only need to focus on ensuring monotonicity of $F_{\theta}$ .

Before providing a characterization of (strongly) monotone operators, let us state the following algebraic property.

Lemma 3.1.

Let $M\in\mathbb{R}^{N\times N}$ be a symmetric matrix, let $\rho\in[\overline{\lambda}_{\max}(M),+\infty[$ , and let $M^{\prime}=\rho\operatorname{Id}\,-M$ . Then

(13)

\lambda_{\min}(M)=\rho-\overline{\lambda}_{\max}(M^{\prime}).

Proof 3.2.

Let $(\lambda_{i})_{1\leqslant i\leqslant N}$ be the eigenvalues of $M$ sorted in decreasing order. $M^{\prime}$ is a symmetric matrix with eigenvalues

(14)

(\forall i\in\{1,\ldots,N\})\quad\lambda_{i}^{\prime}=\rho-\lambda_{i}.

Since $\rho\geqslant\overline{\lambda}_{\max}(M)=\max\{\lambda_{1},-\lambda_{N}\}% \geqslant\lambda_{1}$ , we have $0\leqslant\rho-\lambda_{1}\leqslant\cdots\leqslant\rho-\lambda_{N}$ . Thus, $\overline{\lambda}_{\max}(M^{\prime})=\rho-\lambda_{N}$ , which yields (13).

The following conditions will be leveraged later to enforce monotonicity of the network during its training.

Proposition 3.3.

Let $T\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ be Fréchet differentiable, and let, for every $x\in\mathbb{R}^{n}$ , $\operatorname{J}_{T}(x)$ be the Jacobian of $T$ evaluated at $x$ . Let $R_{T}$ be the reflected operator of $T$ , let $\beta\in\left[0,+\infty\right[$ , and let $\rho\in[\overline{\lambda}_{\max}\big{(}\operatorname{J}_{R_{T}}^{\rm s}(x)% \big{)},+\infty[$ . The following properties are equivalent:

(i)

$T$ is $\beta$ -strongly monotone;
(ii)

For every $x\in\mathbb{R}^{n}$ , $\operatorname{J}_{T}^{\rm s}(x)\succeq\beta\operatorname{Id}\,$ ;
(iii)

For every $x\in\mathbb{R}^{n}$ , $J_{R_{T}}^{\rm s}(x)\succeq(2\beta-1)\operatorname{Id}\,$ ;

(iv)

For every $x\in\mathbb{R}^{n}$ ,

(15)

\rho-\overline{\lambda}_{\max}\Big{(}\rho\operatorname{Id}\,-\operatorname{J}_% {R_{T}}^{\rm s}(x)\Big{)}\geqslant 2\beta-1.

In addition, if $\beta>0$ , then $T$ is invertible.

Proof 3.4.

We first show that (i) $\Rightarrow$ (ii). By definition [4, Definition 22.1] of a $\beta$ -strongly monotone operator, we have, for every $(x,h)\in(\mathbb{R}^{n})^{2}$ and $\alpha\in]0,+\infty[$ ,

		$\displaystyle{\left\langle{{h}\mid{T(x+h)-T(x)}}\right\rangle}\geqslant\beta\\|% h\\|^{2}$
(16)		$\displaystyle\Leftrightarrow\quad$	$\displaystyle{\left\langle{{h}\mid{T(x+\alpha h)-T(x)}}\right\rangle}\geqslant% \beta\alpha\\|h\\|^{2}.$

As $T$ is Fréchet differentiable, we deduce that

(17)		$\displaystyle\lim_{\begin{subarray}{c}\alpha\to 0\\ \alpha>0\end{subarray}}{\left\langle{{h}\mid{\frac{T(x+\alpha h)-T(x)}{\alpha}% }}\right\rangle}\geqslant\beta\\|h\\|^{2}$	$\displaystyle\quad\Leftrightarrow\quad{\left\langle{{h}\mid{\operatorname{J}_{% T}(x)h}}\right\rangle}\geqslant\beta\\|h\\|^{2}$
(18)			$\displaystyle\quad\Leftrightarrow\quad{\left\langle{{h}\mid{\operatorname{J}_{% T}^{\rm s}(x)h}}\right\rangle}\geqslant\beta\\|h\\|^{2}.$

Hence, for every $x\in\mathbb{R}^{n}$ , $\operatorname{J}_{T}^{\rm s}(x)\succeq\beta\operatorname{Id}\,$ .

We now show that (ii) $\Rightarrow$ (i). Let $(x,h)\in(\mathbb{R}^{n})^{2}$ , and let $\phi\colon[0,+\infty[\to\mathbb{R}$ be defined as

(19)

(\forall\alpha\in[0,+\infty[)\quad\phi(\alpha)={\left\langle{{h}\mid{T(x+% \alpha h)-T(x)}}\right\rangle}-\beta\alpha\|h\|^{2}.

We notice that $\phi$ is differentiable on $[0,+\infty[$ and its derivative $\phi^{\prime}$ is such that

(20)

\begin{cases}\phi(0)=0\\ (\forall\alpha\in[0,+\infty[)\quad\phi^{{}^{\prime}}(\alpha)={\left\langle{{h}% \mid{\operatorname{J}_{T}(x)h}}\right\rangle}-\beta\|h\|^{2}.\end{cases}

According to (ii) and (17), we have, for every $\alpha\in[0,+\infty[$ , $\phi^{\prime}(\alpha)\geqslant 0$ . Thus $\phi$ is an increasing function and

(21)

(\forall\alpha\in[0,+\infty[)\quad{\left\langle{{h}\mid{T(x+\alpha h)-T(x)}}% \right\rangle}-\beta\alpha\|h\|^{2}=\phi(\alpha)\geqslant\phi(0)=0.

Hence $T$ is $\beta$ -strongly monotone.

By using the definition of positive (semi) definiteness, (ii) is equivalent to

(22)

(\forall x\in\mathbb{R}^{n})\quad\lambda_{\min}\big{(}\operatorname{J}_{T}^{% \rm s}(x)\big{)}\geqslant\beta.

The equivalence with (iv) is a direct consequence of Lemma 3.1.

In addition, we have

	$\displaystyle\quad\Leftrightarrow\quad(\forall u\in\mathbb{R}^{n})\quad u^{% \top}J_{T}^{\rm s}(x)u\geqslant\beta\lVert u\rVert^{2}$
	$\displaystyle\quad\Leftrightarrow\quad(\forall u\in\mathbb{R}^{n})\quad u^{% \top}(2J_{T}^{\rm s}(x)-\operatorname{Id}\,)u\geqslant(2\beta-1)\lVert u\rVert% ^{2}$
	$\displaystyle\quad\Leftrightarrow\quad(\forall u\in\mathbb{R}^{n})\quad u^{% \top}J_{R_{T}}^{\rm s}(x)u\geqslant(2\beta-1)\lVert u\rVert^{2}$
(23)		$\displaystyle\quad\Leftrightarrow\quad\text{\ref{prop:res-reflected:ii}}.$

By using the definition of positive (semi) definiteness, (iii) is equivalent to

(24)

(\forall x\in\mathbb{R}^{n})\quad\lambda_{\min}\big{(}\operatorname{J}_{R_{T}}% ^{\rm s}(x)\big{)}\geqslant 2\beta-1.

The equivalence with (iv) is a direct consequence of Lemma 3.1.

The last statement straightforwarly follows from [4, Corollary 20.28] and [4, Proposition 22.11] (see also [1]). Since $T$ is strongly monotone, it is strictly monotone, hence injective. Moreover, $T$ is single valued, monotone and continuous and is thus maximally monotone. Lastly, a maximally strongly monotone operator is surjective.

Remark 3.5.

(i)

The previous proposition allows us to build more general operators than gradients of convex functions. For example, let $f\colon{\mathcal{H}}\to\mathbb{R}$ be $\beta$ -strongly convex, and twice differentiable. Let $S\colon{\mathcal{H}}\to{\mathcal{H}}$ be Fréchet differentiable and anti-symmetric, i.e., for every $x\in{\mathcal{H}}$ , $S(x)=-S(x)^{\top}$ . Then, as a direct consequence of Proposition 3.3, $T=\nabla f+S$ is $\beta$ -strongly monotone.
(ii)

The Fréchet differentiability assumption required in Proposition 3.3 is fulfilled by neural networks using differentiable activation functions (e.g., sigmoids or softmax). As highlighted in [6], properties satisfied for such activation functions generalize to standard non differentiable ones (e.g., ReLU).
(iii)

Proposition 3.3 (iv) enables to characterize a differentiable (strongly) monotone operator through the maximum magnitude eigenvalue of some matrices, which can be easily computed by power iteration. Further details will be provided in the next section to train a monotone NN.
(iv)

The sufficient strong monotonicity condition for the invertibility of $T$ is restrictive. In [1], a sufficient condition for the strong monotonicity based on the Lipschitzianity of the Cayley operator associated to $T$ is employed. Nevertheless, these conditions are not necessary to guarantee the convergence of Section 3.

3.2 Proposed regularization approach for training monotone NNs

We will now leverage the results given in the previous sections to design a method for training monotone NNs. In a nutshell, we will make the assumption that, for every $\theta\in\Theta$ , $F_{\theta}$ is Fréchet differentiable (see Remark 3.5 (ii)) and resort to the power iteration method to control the eigenvalues of the NN, in accordance with Remark 3.5 (iii).

Subsequently, we will propose an optimization strategy to train monotone NNs. The basic principle for training a NN is to learn its parameters on a specific finite dataset consisting of pairs of input/groundtruth data denoted by $\mathbb{D}_{\rm train}\subset\mathbb{R}^{n}\times\mathbb{R}^{n}$ . Then the objective is to

(25)

\underset{\begin{subarray}{c}{\theta\in\Theta}\end{subarray}}{\text{\rm minimize% }}\;\;\displaystyle\sum_{(x,y)\in\mathbb{D}_{\rm train}}\mathcal{L}(F_{\theta}% (x),y),

where $\mathcal{L}\colon\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ is a loss function. Such a problem can then be solved efficiently using stochastic first-order methods such as Stochastic Gradient Descent, Adagrad, or Adam [31]. In this paper, we are interested in training NNs that are monotone. In this case, (25) becomes a constrained optimization problem:

(26)

\underset{\begin{subarray}{c}{\theta\in\Theta}\end{subarray}}{\text{\rm minimize% }}\;\;\displaystyle\sum_{(x,y)\in\mathbb{D}_{\rm train}}\mathcal{L}(F_{\theta}% (x),y)\quad\text{ such that }F_{\theta}\text{ is monotone}.

Standard gradient-based optimization algorithms are not inherently tailored for solving problems with constraints. While it is possible to pair some of these algorithms with projection techniques (as in [26, 33]), there are cases when identifying an appropriate projection method can be challenging. Moreover, utilizing projection methods may introduce convergence issues, particularly in scenarios involving non-convex constraints or projection operators with no closed form.

In the literature, we can find two deep learning frameworks ensuring constraints to be satisfied by NNs: either using a specific architecture definition [9, 16], or enforcing the constraint iteratively during the training procedure, possibly by adopting a penalized approach rather than a constrained one [1, 27, 33]. The first method consists in modifying the architecture to ensure that the solution always satisfies the desired property (for example the output of a ReLU function is guaranteed to be nonnegative). The second method enables to constrain any NN architecture at the expense of a higher training complexity. In this work, we will adopt the second approach, and propose a training procedure imposing monotonicity to the NN independently of its architecture by solving a penalized problem.

First, according to LABEL:{prop:jacobianmonotony}(iii), we notice that (26) can be rewritten as

(27)

\underset{\begin{subarray}{c}{\theta\in\Theta}\end{subarray}}{\text{\rm minimize% }}\;\;\displaystyle\sum_{(x,y)\in\mathbb{D}_{\rm train}}\mathcal{L}\big{(}F_{% \theta}(x),y\big{)}\quad\text{ such that }(\forall x\in\mathbb{R}^{n})\quad% \lambda_{\text{min}}\big{(}\operatorname{J}_{R_{F_{\theta}}}^{\rm s}(x)\big{)}% \geqslant-1.

Equivalently, we could use the constraint $(\forall x\in\mathbb{R}^{n})$ $\lambda_{\text{min}}\big{(}\operatorname{J}_{F_{\theta}}^{\rm s}(x)\big{)}\geqslant 0$ , but we observed that the former formulation leads numerically to a more stable training. In practice, imposing the constraint for all possible images in $\mathbb{R}^{n}$ is not tractable. Instead, we impose it on a subset of images $\mathbb{D}_{\rm penal}\subset\mathbb{R}^{n}$ .

Further, as suggested earlier, we introduce an external penalization approach for solving (27) [25, Section 13.1]. Hence, we aim to

(28)

\underset{\begin{subarray}{c}{\theta\in\Theta}\end{subarray}}{\text{\rm minimize% }}\;\;\displaystyle\sum_{(x,y)\in\mathbb{D}_{\rm train}}\mathcal{L}(F_{\theta}% (x),y)+\xi\sum_{\widetilde{x}\in\mathbb{D}_{\rm penal}}\mathcal{P}(\theta,% \widetilde{x}),

where $\xi\in\left]0,+\infty\right[$ is the penalization factor and the penalty function $\mathcal{P}$ is given by

(29)

(\forall\theta\in\Theta)(\forall\widetilde{x}\in\mathbb{R}^{n})\quad\mathcal{P% }(\theta,\widetilde{x})=-\min\big{\{}1+\lambda_{\text{min}}(\operatorname{J}_{% R_{F_{\theta}}}^{\rm s}(\widetilde{x})),\epsilon\big{\}}.

In the above definition, $\epsilon\in\left]0,+\infty\right[$ is a severity parameter used to control the enforcement of the constraint. Finally, we define

(30)

\mathbb{D}_{\rm penal}=\left\{{\widetilde{x}\in\mathbb{R}^{n}}~{}\Big{|}~{}{(% \exists(x,y)\in\mathbb{D}_{\rm train})(\exists\nu\in[0,1])\,\widetilde{x}=\nu x% +(1-\nu)y}\right\}.

In practice, we generate points in $\mathbb{D}_{\rm penal}$ by choosing $\nu$ as a realization of a random variable with uniform distribution on $[0,1]$ .

3.3 Penalized training implementation

This section aims to present the practical procedure we employ for solving (28), so as to train a monotone neural network. We will thus provide pseudo-codes usable with deep learning frameworks such as Pytorch. In the following, we first describe how we handle the penalization function $\mathcal{P}$ . We then explain how to leverage this to solve our constrained optimization problem.

Algorithm 1 Computation of

\lambda_{\min}\big{(}\operatorname{J}_{R_{F_{\theta}}}^{\rm s}(\widetilde{x})% \big{)}

1:Input: •

R_{F_{\theta}},\widetilde{x}\in\mathbb{D}_{\rm penal}

\triangleright

Neural network model and data •

N_{\rm iter}

\triangleright

Parameter

2:Disable auto-differentiation

3: Computation of

\rho>\overline{\lambda}_{\max}\big{(}\operatorname{J}_{R_{F_{\theta}}}^{\rm s}% (\widetilde{x})\big{)}

u_{0}\leftarrow

realization of

\mathcal{N}(0,\operatorname{Id}\,)

5: for

k=1,\ldots,N_{\rm iter}

u_{k+1}=\frac{\operatorname{J}_{R_{F_{\theta}}}^{\rm s}(\widetilde{x})u_{k}}{% \|u_{k}\|_{2}^{2}}

7: end for

8: Choose

\widehat{\rho}>\frac{u_{k+1}^{\top}\operatorname{J}_{R_{F_{\theta}}}^{\rm s}(% \widetilde{x})u_{k+1}}{\|u_{k+1}\|_{2}^{2}}

9: Computation of the eigenvector associated with

\chi=\overline{\lambda}_{\max}\big{(}\rho\operatorname{Id}\,-\operatorname{J}_% {R_{F_{\theta}}}^{\rm s}(\widetilde{x})\big{)}

10:

v_{0}\leftarrow

realization of

\mathcal{N}(0,\operatorname{Id}\,)

11: for

k=1,\ldots,N_{\rm iter}

12:

v_{k+1}=\frac{\big{(}\widehat{\rho}\operatorname{Id}\,-\operatorname{J}_{R_{F_% {\theta}}}^{\rm s}(\widetilde{x})\big{)}v_{k}}{\|v_{k}\|_{2}^{2}}

13: end for

14: Enable auto-differentiation

15: Computation of

\chi

16:

\widehat{\chi}=\frac{v_{N_{\rm iter}+1}^{\top}\big{(}\widehat{\rho}% \operatorname{Id}\,-\operatorname{J}_{R_{F_{\theta}}}^{\rm s}(\widetilde{x})% \big{)}v_{N_{\rm iter}+1}}{\|v_{N_{\rm iter}+1}\|_{2}^{2}}

17:return

\widehat{\rho}-\widehat{\chi}\simeq\lambda_{\min}\big{(}\operatorname{J}^{\rm s% }_{R_{F_{\theta}}}(x)\big{)}

Penalization computation

Let $\widetilde{x}\in\mathbb{D}_{\rm penal}$ . We need to evaluate and subdifferentiate $\theta\in\Theta\mapsto\mathcal{P}(\theta,\widetilde{x})$ , and in particular $\lambda_{\min}(\operatorname{J}^{\rm s}_{R_{F_{\theta}}}(\widetilde{x}))$ . As mentioned in Remark 3.5 (iii), (29) can be reexpressed by using maximum absolute eigenvalues instead of a minimum eigenvalue. Similarly to [27], we will use a power iterative method that is designed to provide the largest eigenvalue in absolute value of a given matrix. The power iteration makes the computation tractable since only vector products are necessary. In particular, in our simulations, we will use the Jacobian-vector product (JVP) in Pytorch, using automatic differentiation. It follows from (15) that we actually need to combine two power iterations. The pseudo-code of the resulting method is given in Algorithm 1. In particular, we first compute $\rho\geqslant\overline{\lambda}_{\max}\big{(}\operatorname{J}_{R_{F_{\theta}}}% ^{\rm s}(\widetilde{x})\big{)}$ using $N_{\rm iter}$ power iterations in steps 3-7 of Algorithm 1, followed by a second run of the power iterative method to compute $\chi=\overline{\lambda}_{\max}\big{(}\rho\operatorname{Id}\,-\operatorname{J}_% {R_{F_{\theta}}}^{\rm s}(\widetilde{x})\big{)}$ in steps 9-16 of Algorithm 1. Further, in steps 8, 12 and 16, the dot product between the Jacobian and a vector is computed using JVP as mentioned above. JVP is a function that takes three inputs: an operator $T\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ and two vectors $(\widetilde{x},u)\in(\mathbb{R}^{n})^{2}$ . It returns $\operatorname{J}_{T}(\widetilde{x})u$ . Note that in order to compute $\operatorname{J}^{\rm s}_{T}(\widetilde{x})u$ , we use the double backward trick¹¹1https://j-towns.github.io/2017/06/12/A-new-trick.html and the transpose of the vector-Jacobian product.

Since all the operations used in Algorithm 1 are differentiable operations, a subgradient of $\theta\in\Theta\mapsto\mathcal{P}(\theta,\widetilde{x})$ can be estimated through auto-differentiation to enable its use in standard optimizers (see next paragraph). However, according to [36], the resulting gradient may be noisy and computationally heavy. To overcome these issues, we disable the track of the auto-differentiation of the forward operators (using torch.no_grad() function in Pytorch) to compute the needed eigenvalue and its corresponding eigenvector (steps 3-13). Then, we activate the auto-differentiation tracking to compute the Rayleigh quotient (step 14). Note that this approximation is equivalent to initializing the power method with the correct eigenvector and then perform one step of the algorithm. In Section 4.1.4, we will show that the current approximation is accurate enough to allow us to train networks such that the constraint $\lambda_{\min}(\operatorname{J}^{\rm s}_{R_{F_{\theta}}}(\widetilde{x}))% \geqslant-1$ is properly satisfied.

Remark 3.6.

(i)

The use of the power iteration has the advantage of bypassing the computation of the Jacobian at a given point, which can be huge depending on the input/output dimensions.
(ii)

As a reminder, we cannot decompose our operator into multiple monotone operators as in [29], as the composition of two monotone operators is not necessarily monotone. Thus, being able to enforce monotonicity on a large model in an end to end fashion is required.
(iii)

In order to accelerate the process, we only compute the monotonicity penalization on one random element of each batch, which experimentally yields similar results to using a batch of data, with the benefit of speeding up the training.

Generic training framework

We present the pseudo-code summarizing the overall learning procedure in Algorithm 2. At each epoch $j\in\{1,\ldots,N_{\rm epochs}\}$ , we iteratively select $\mathbb{B}$ from $\mathbb{D}_{\rm train}$ of size $B$ until all the data has been seen. The associated batch loss $\ell_{\mathbb{B}}$ (see step 9 in Algorithm 2) is associated with a Pytorch object containing all the history of the operations (called computational graph)²²2https://pytorch.org/blog/computational-graphs-constructed-in-pytorch/ used to compute the loss value. It is obtained by simply running a forward step with the neural network $F_{\theta}$ , and the associated loss and penalization values. The computational graph is then used to compute the subgradient $g_{\mathbb{B}}$ through auto-differentiation (step (11) in Algorithm 2). Finally, the optimizer step (see step 12 in Algorithm 2) updates the parameters $\theta$ of $F_{\theta}$ with the selected scheme (e.g., gradient descent with momentum, etc) ³³3For example, the Pytorch implementation of Adam can be found in https://pytorch.org/docs/stable/generated/torch.optim.Adam.html..

Algorithm 2 Training a monotone network

F_{\theta}

1:Input: •

F_{\theta}

N_{\rm epochs}

B

\Delta\xi>0

\triangleright

Training parameters •

\mathcal{L}

\mathcal{P}

\mathbb{D}_{\rm train}

\triangleright

Loss, penalization and training set • Optimizer step:

\mathcal{O}

\colon(\theta,g)\mapsto\theta^{+}

\triangleright

e.g., Adam, SGD, etc.

\xi\leftarrow 0

3:for

j=1,\ldots,N_{\rm epochs}

4: for each batch

\mathbb{B}=\{(x_{b},y_{b})\}_{1\leqslant b\leqslant B}\subset\mathbb{D}_{\rm train}

of size

B

5: Computational graph related to the loss and the penalization:

b_{0}\leftarrow

realization of discrete random uniform variable in

\{1,\ldots,B\}

\nu\leftarrow

realization of random uniform variable in

[0,1]

\widetilde{x}_{b_{0}}\leftarrow\nu x_{b_{0}}+(1-\nu)y_{b_{0}}

\ell_{\mathbb{B}}\colon\vartheta\mapsto\frac{1}{B}\sum_{b=1}^{B}\mathcal{L}(F_% {\vartheta}(x_{b}),y_{b})+\xi\,\mathcal{P}(\vartheta,\widetilde{x}_{b_{0}})

\triangleright

Use Algorithm 1

10: Gradient computation and optimizer step:

11:

g_{\mathbb{B}}(\theta)\in\partial\ell_{\mathbb{B}}(\theta)

12:

\theta\leftarrow\mathcal{O}(\theta,g_{\mathbb{B}}(\theta))

13: end for

14:

\xi\leftarrow\xi+\Delta\xi

\triangleright

Increase penalization parameter

15:end for

16:Output:

F_{\theta}

4 Learning non-linear model approximations

In this section, we will study the challenging problem of learning the unknown forward model in a non-linear inverse imaging problem and propose a method to solve this problem.

4.1 Imaging semi-blind non-linear inverse problem

We consider a non-linear inverse problem, where the objective is to find an estimate $\widehat{x}\in C$ , with $C\subset\mathbb{R}^{n}$ , of an original unknown image $\overline{x}\in C$ from degraded measurements $y\in\mathbb{R}^{m}$ obtained as

(31)

y=F(\overline{x})+w,

where $F\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ models a measurement operator (possibly unknown) , and $w\in\mathbb{R}^{n}$ is a realization of an additive i.i.d. white Gaussian random variable with zero-mean and standard deviation $\sigma\geqslant 0$ .

In this section, we assume that the measurement operator $F$ is unknown. In turn, we have access to a collection of pairs of true images and associated measurements $(\overline{x},y)\in\mathbb{D}_{\rm train}$ , that we will use to learn $F$ . Further, since the objective is to estimate $\widehat{x}$ , we will see next that it is judicious to learn a monotone approximation of $F$ . In the following, we will assume that $C$ is closed convex with a nonempty interior.

4.1.1 Measurement operator setting

In the remainder, we assume that the measurement operator $F$ is a sum of $K\in\mathbb{N}$ simple non-linear composite functions, where the inner terms correspond to linear convolution filters. Precisely, we consider

(32)

(\forall x\in\mathbb{R}^{n})\quad F(x)=\frac{1}{K}\sum_{k=1}^{K}S_{\delta}(L_{% k}x),

where $S_{\delta}$ is a saturation function (e.g. Hyperbolic tangent function) with parameter $\delta>0$ and, for every $k\in\{1,\ldots,K\}$ , $L_{k}\in\mathbb{R}^{n\times n}$ (e.g., $L_{k}$ may correspond to a convolution operator).

In our experiments, we consider model (31)-(32) with either $K=1$ or $K=5$ motion blur kernels, of size $9\times 9$ , generated randomly using the motionblur toolbox⁴⁴4https://github.com/LeviBorodenko/motionblur. The $5$ generated kernels are displayed in Figure 1. These kernels are all positive and are normalized so that the sum of the components is equal to $1$ . In addition, we choose $S_{\delta}$ to be a modified $\tanh$ function (more details are given in Appendix A), with $\delta\in\{1,0.6\}$ . Finally, we will consider two noise levels with standard deviation $\sigma\in\{0,10^{-2}\}$ .

Note that there is no guarantee that model (32) is monotone. Solving inverse problems with such an operator can appear to be challenging, especially due to the non-linearity $S$ . A standard practice in such an inverse problem context consists in replacing $F$ by its first-order approximation expressed as

(33)

(\forall x\in\mathbb{R}^{n})\quad{F}^{\rm aff}(x)=\frac{\delta}{K}\sum_{k=1}^{% K}L_{k}x+\dfrac{1-\delta}{2}.

For simplicity, when the parameter $\delta$ is assumed to be unknown, we can use instead the following linear approximation:

(34)

(\forall x\in\mathbb{R}^{n})\quad{F}^{\rm lin}(x)=\frac{\delta}{K}\sum_{k=1}^{% K}L_{k}x.

Such a linearization of the model is standard in inverse problems. Note that it does not introduce any bias when $\delta=1$ .

Since the convolution filters are in practice unknown, ${F}^{\rm lin}$ represents the linear oracle. For this purpose, we will be approximating $F$ using a single linear operator $F_{\theta}^{\rm{lin}}$ , in order to compare our method to a more realistic result. More details about the results of this approximation are provided in Appendix B.

An example is provided in Figure 2, for the case when $K=5$ , $\sigma=0$ , and $S_{\delta}$ is chosen as in LABEL:{s:SaturationFunction}. From left to right, the figure shows the true image $\overline{x}$ , the observed image $y=F(\overline{x})$ , and the approximated linear degradation obtained as $\widetilde{y}={F}^{\rm lin}(\overline{x})$ . It can be observed that the light background behind the penguin is darker when using the true acquisition model $F$ than with the linear version ${F}^{\rm lin}$ , due to the saturation function $S_{\delta}$ .

Refer to caption — Figure 1: Blurring kernels used to model linear operators $(L_{k})_{1\leqslant k\leqslant 5}$ , used in model (32).

4.1.2 Monotone inclusion formulation to solve (31)

The objective of this section is to find an estimate $\widehat{x}$ of $\overline{x}$ from $y$ , solving the inverse problem (31), i.e., inverting $F$ . Note that the approach described below is not restricted to the particular choice of operator $F$ described in Section 4.1.1.

The problem of interest is equivalent to

(35)

\mbox{find}\;\widehat{x}\in C\text{ such that }F(\widehat{x})\approx y.

We propose two approaches to solve this problem.

Direct regularized approach

We propose to define $\widehat{x}$ as a solution to a regularized monotone inclusion problem of the form:

(36)

\text{find }\widehat{x}\in\mathbb{R}^{n}\text{ such that }0\in F_{\theta}(% \widehat{x})-y+\varrho\nabla r(\widehat{x})+N_{C}(\widehat{x}),

where $F_{\theta}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a continuous monotone approximation to $F$ on $C$ , $\varrho\in\left[0,+\infty\right[$ , and $r\colon\mathbb{R}^{n}\to\mathbb{R}$ is a differentiable and convex regularization function.

Interestingly, (36) is an instance of Problem 1.1, where

(37)

(\forall x\in\mathbb{R}^{n})\quad\begin{cases}h(x)=-{\left\langle{{y}\mid{x}}% \right\rangle}\\ A(x)=F_{\theta}(x)+\varrho\nabla r(x).\end{cases}

We have the following result concerning the existence of a solution to (36).

Proposition 4.1.

Assume that one of the following statements holds:

(i)

$\rho>0$ and $r$ is strongly convex,
(ii)

$C$ is bounded.

Then there exists a solution $\widehat{x}$ to (36).

Proof 4.2.

Since $A$ is monotone on $C$ , it admits a maximally monotone extension $\widetilde{A}$ on $\mathbb{R}^{n}$ . Problem (36) is thus equivalent to find a zero of $\widetilde{A}+\nabla h+N_{C}$ . By proceeding similarly to the proof of Proposition 2.4 we can show that this operator is maximally monotone. In addition, if $\rho>0$ and $r$ is strongly convex, it is strongly monotone. The existence of a solution to (36) follows then from [4, Proposition 23.36(ii)] and [4, Corollary 23.37(ii)]

To solve (36) numerically, we can make use of Section 3. The convergence of this algorithm is guaranteed by Proposition 2.4, as soon as condition (ii) on the choice of the step-size is satisfied.

Remark 4.3.

In the particular case when $y$ is replaced by $F_{\theta}(\overline{x})$ in (36), Algorithm (2.2) enables computing a regularized inverse of $F_{\theta}$ delivering an estimate $x^{*}$ to $\overline{x}$ . In particular, if $\rho=0$ and $x^{*}\in\text{\rm int}\,~{}C$ , $x^{*}=F_{\theta}^{-1}(F_{\theta}(\overline{x}))$ .

Least-squares regularized approach

In addition to the direct approach described above, we will investigate an original reformulation of problem (36) as

(38)

\text{Find }\widetilde{x}\in\mathbb{R}^{n}\text{ such that }0\in{{F}^{\rm lin}% }^{\top}F_{\theta}(\widetilde{x})-\widetilde{y}+\varrho\nabla r(\widetilde{x})% +N_{C}(\widetilde{x})

where $\widetilde{y}={{F}^{\rm lin}}^{\top}y$ and $F_{\theta}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a continuous approximation to $F$ on $C$ . This approach is potentially less restrictive than (36) as it does not impose the approximation to the linear model to be monotone, but requires only ${{F}^{\rm lin}}^{\top}F_{\theta}$ to be monotone. Interestingly, when $F_{\theta}={F}^{\rm lin}$ , (38) corresponds to the standard least-squares formulation. Then the existence of a solution to (38) is guaranteed similarly to Proposition 4.1. Section 2.2 can be rewritten for solving (38), similarly to (7).

Remark 4.4.

(i)

Similarly to Remark 4.3, operator $A$ can be inverted using (38) by setting $\widetilde{y}={{F}^{\rm lin}}^{\top}\overline{x}$ and $\varrho=0$ .
(ii)

According to [4, Prop. 20.10], when $K=1$ , ${{F}^{\rm lin}}^{\top}F=L_{1}^{\top}\circ S_{\delta}\circ L_{1}$ is monotone.

Learned approximations to $F$

In this work, we assume that neither the kernels $(L_{k})_{1\leqslant k\leqslant K}$ nor the parameter $\delta$ are known in (32). Then, we leverage both the direct and the least-squares formulations described above to derive our learned approximation strategies for $F$ . The different considered approaches are described below:

(i)

The first classical strategy consists in learning a linear approximation $F_{\theta}^{\rm{lin}}$ of $F$ .
(ii)

We leverage the training approach developed in Section 3.2 to learn a monotone approximation $F_{\theta}^{\rm{mon}}$ of $F$ .
(iii)

To show the necessity of the proposed constrained training approach to obtain a monotone operator, we also consider an approximation $F_{\theta}^{\rm{nom}}$ of $F$ , that is learned without the proposed regularization (i.e., $\xi=0$ in (28)).
(iv)

We finally leverage the least-squares formulation (38) combined with the training approach developed in Section 3.2 to learn a monotone operator ${F_{\theta}^{\rm{lin}}}^{\top}F_{\theta}$ . Since ${F}^{\rm lin}$ is unknown, we use the learned linear operator $F_{\theta}^{\rm{lin}}$ , and learn $F_{\theta}$ under the constraint that ${F_{\theta}^{\rm{lin}}}^{\top}F_{\theta}$ is monotone. In the remainder, we will use the notation $\widetilde{F}^{\rm{mon}}_{\theta}={F_{\theta}^{\rm{lin}}}^{\top}F_{\theta}$ when referring to this approximation strategy.

Remark 4.5.

(i)

In the last restoration approach based on the least-squares strategy (38), ${F}^{\rm lin}$ is replaced by $F_{\theta}^{\rm{lin}}$ and $\widetilde{y}={F_{\theta}^{\rm{lin}}}^{\top}y$ .
(ii)

We emphasize that none of the proposed restoration procedures can be used with $F_{\theta}^{\rm{nom}}$ without losing theoretical guarantees, since neither $F_{\theta}^{\rm{nom}}$ nor $({F_{\theta}^{\rm{lin}}}^{\top}F_{\theta}^{\rm{nom}})$ are monotone.

4.1.3 Model and training procedure

In this section we describe the learning procedures adopted for the considered approximations to $F$ described in the previous section, focusing on non-linear approximations. The learning procedure of the linear operator $F_{\theta}^{\rm{lin}}$ is given in Appendix B.

NN architecture

We use a Residual Unet [37] consisting of $5$ blocks where the output of each down convolution have $32$ , $64$ , $128$ , $256$ , and $512$ feature maps, respectively. LeakyReLU activations were used for the intermediate layers. No activation was used in the last layer, so the learned NN must model the saturation using its weights and biases.

Training dataset

For $\mathbb{D}_{\rm train}$ , we use patches of size $256\times 256$ , randomly extracted from the BSD500 dataset, with pixel-values scaled to the range $[0,1]$ . The couples $(\overline{x},y)$ in $\mathbb{D}_{\rm train}$ are linked through model (31)-(32). To investigate the stability against noise, we consider two cases for the training set: (i) training without noise (i.e., $\sigma_{\rm train}=0$ ), and (ii) training with noise level $\sigma_{\rm train}=0.01$ .
For $\mathbb{D}_{\rm penal}$ , we crop the images of $\mathbb{D}_{\rm train}$ to smaller patches of size $64\times 64$ . We set $68$ images from the dataset as test images (corresponding to BSD68), and use the other images for training.

Training procedure

The NN is trained using Adam optimizer with a learning rate of $2\times 10^{-4}$ for $200$ epochs to solve (28) with an $\ell_{1}$ loss function:

(39)

\mathcal{L}(F_{\theta}(x),y)=\|F_{\theta}(x)-y\|_{1}.

We use a learning rate scheduler that decreases the learning rate by a factor of $0.1$ when the training loss reaches a plateau. Model $F_{\theta}^{\rm{nom}}$ is trained by setting $\xi=0$ in (28) and $\Delta\xi=0$ , and Model $F_{\theta}^{\rm{mon}}$ is obtained by setting $\xi=0.1$ in (28) at start and $\Delta\xi=0.1$ . Finally, $\epsilon$ is chosen equal to $0.01$ in (29).

Filters	Noise	Model	$\text{MAE}(y,F_{\theta}(\overline{x}))$	$\min\lambda_{\text{min}}\big{(}\operatorname{J}_{F_{\theta}}^{\rm s}(\overline% {x})\big{)}$
			( $\times 10^{-2}$ )	( $\times 10^{-2}$ )
$\delta=1$ for $S_{\delta}$ in (32)
$K=1$	$\sigma_{\rm train}=0$	$F_{\theta}^{\rm{mon}}$	1.5658 ( $\pm$ 0.58)	1.67
		$F_{\theta}^{\rm{nom}}$	0.2932 ( $\pm$ 0.11)	-29.99
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.6822 ( $\pm$ 0.25)	$0.69^{*}$
		$F_{\theta}^{\rm{lin}}$	2.1474 ( $\pm$ 1.38)	-24.69
	$\sigma_{\rm train}=0.01$	$F_{\theta}^{\rm{mon}}$	1.1575 ( $\pm$ 0.42)	1.10
		$F_{\theta}^{\rm{nom}}$	0.3020 ( $\pm$ 0.11)	-27.72
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.5351 ( $\pm$ 0.19)	$1.13^{*}$
		$F_{\theta}^{\rm{lin}}$	2.1607 ( $\pm$ 1.37)	-25.87
$K=5$	$\sigma_{\rm train}=0$	$F_{\theta}^{\rm{mon}}$	0.5272 ( $\pm$ 0.20)	1.18
		$F_{\theta}^{\rm{nom}}$	0.2795 ( $\pm$ 0.10)	-22.06
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.6108 ( $\pm$ 0.22)	$1.80^{*}$
		$F_{\theta}^{\rm{lin}}$	2.1473 ( $\pm$ 1.38)	-13.86
	$\sigma_{\rm train}=0.01$	$F_{\theta}^{\rm{mon}}$	0.9414 ( $\pm$ 0.34)	1.80
		$F_{\theta}^{\rm{nom}}$	0.2714 ( $\pm$ 0.09)	-25.48
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.6591 ( $\pm$ 0.25)	$1.64^{*}$
		$F_{\theta}^{\rm{lin}}$	2.1594 ( $\pm$ 1.37)	-16.55
$\delta=0.6$ for $S_{\delta}$ in (32)
$K=1$	$\sigma_{\rm train}=0$	$F_{\theta}^{\rm{mon}}$	1.2816 ( $\pm$ 0.53)	1.91
		$F_{\theta}^{\rm{nom}}$	0.2376 ( $\pm$ 0.09)	-18.73
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.4311 ( $\pm$ 0.14)	$0.37^{*}$
		$F_{\theta}^{\rm{lin}}$	8.2009 ( $\pm$ 2.70)	-42.79
	$\sigma_{\rm train}=0.01$	$F_{\theta}^{\rm{mon}}$	1.1689 ( $\pm$ 0.45)	1.65
		$F_{\theta}^{\rm{nom}}$	0.2275 ( $\pm$ 0.08)	-23.35
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.4679 ( $\pm$ 0.17)	$0.54^{*}$
		$F_{\theta}^{\rm{lin}}$	8.1993 ( $\pm$ 2.69)	-43.18
$K=5$	$\sigma_{\rm train}=0$	$F_{\theta}^{\rm{mon}}$	0.7720 ( $\pm$ 0.30)	1.52
		$F_{\theta}^{\rm{nom}}$	0.1435 ( $\pm$ 0.05)	-17.38
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.4327 ( $\pm$ 0.14)	$0.40^{*}$
		$F_{\theta}^{\rm{lin}}$	8.1920 ( $\pm$ 2.70)	-39.61
	$\sigma_{\rm train}=0.01$	$F_{\theta}^{\rm{mon}}$	0.6867 ( $\pm$ 0.25)	0.66
		$F_{\theta}^{\rm{nom}}$	0.1788 ( $\pm$ 0.06)	-19.04
		$\widetilde{F}^{\rm{mon}}_{\theta}$	0.4809 ( $\pm$ 0.15)	$0.33^{*}$
		$F_{\theta}^{\rm{lin}}$	8.1969 ( $\pm$ 2.70)	-35.39

Table 1: Training results obtained with monotone regularization (

F_{\theta}^{\rm{mon}}

), and without the regularization (

F_{\theta}^{\rm{nom}}

). Results linked to the least squares operator

\widetilde{F}^{\rm{mon}}_{\theta}

and to the learned linear kernel

F_{\theta}^{\rm{lin}}

are also presented. Metrics were computed on cropped images of size

256\times 256

from the BSD68 test set. ^∗For

\widetilde{F}^{\rm{mon}}_{\theta}

, the smallest eigenvalue was computed for the whole operator

{F_{\theta}^{\rm{lin}}}^{\top}F_{\theta}

4.1.4 Training results

In this section, we assess the abilities of the learned NNs in approximating the measurement operator $F$ , as well as satisfying the monotonicity constraint necessary to the reconstruction problem. The results are provided in Table 1. For each model, the average MAE values and the smallest eigenvalue $\min_{\overline{x}}\lambda_{\min}\big{(}\operatorname{J}_{F_{\theta}}^{\rm s}(% \overline{x})\big{)}$ are reported in the table, where $\overline{x}$ are cropped images of size $256\times 256$ from the test set BSD68. Results are provided for the two considered training noise levels $\sigma_{\rm train}\in\{0,0.01\}$ . For $\widetilde{F}^{\rm{mon}}_{\theta}$ , the smallest eigenvalue was computed for the whole operator $\widetilde{F}^{\rm{mon}}_{\theta}={F_{\theta}^{\rm{lin}}}^{\top}F_{\theta}$ . To ensure the monotonicity of the model, the smallest eigenvalue must be non-negative. The MAE values are computed between the non-noisy true output $y=F(\overline{x})$ , and the output obtained through the learned network $F_{\theta}(\overline{x})$ .

On the one hand, we observe that $F_{\theta}^{\rm{nom}}$ slightly performs better in terms of MAE value than $F_{\theta}^{\rm{mon}}$ , in particular when $K=1$ filter is used. This behaviour is expected since the monotonicity constraint reduces the flexibility of the model. We further observe that the noise level added during the training does not appear to impact much the results. On the other hand, the learned model corresponding to $\widetilde{F}^{\rm{mon}}_{\theta}$ seems to better reproduce the true model $F$ , with MAE values closer to the non-monotone model. This is certainly due to the weaker constraint on its structure. Finally, $F_{\theta}^{\rm{lin}}$ is the worst approximation of $F$ , and is not monotone.

Figure 3 and Figure 4 show the output images obtained when using different versions of the measurement operator, for two image examples, for $K=5$ with $\delta=1$ and $\delta=0.6$ , respectively. Output images are provided for the true unknown operator $F$ , as well as for approximated operators ${F}^{\rm lin}$ (true unknown linear approximation), $F_{\theta}^{\rm{lin}}$ , $F_{\theta}^{\rm{nom}}$ , $F_{\theta}^{\rm{mon}}$ , and $\widetilde{F}^{\rm{mon}}_{\theta}$ . Results are shown for models trained without noise (i.e., $\sigma_{\rm train}=0$ ). For each image, we give the MAE between the output from the approximated operator and the output from the true operator. Furthermore, values $(\lambda_{\min},\lambda_{\max})$ are reported. We can observe that $\lambda_{\min}<0$ for all operators apart from the monotone network $F_{\theta}^{\rm{mon}}$ and the relaxed version $\widetilde{F}^{\rm{mon}}_{\theta}$ . The true model $F$ , the linear approximation ${F}^{\rm lin}$ , and the learned linear approximation $F_{\theta}^{\rm{lin}}$ are all non-monotone. This shows that our method enables approximating a non-monotone operator using a monotone network.

	Problem	Operator	$\sigma_{\rm train}=0$		$\sigma_{\rm train}=0.01$
	Problem	Operator	PSNR	SSIM	PSNR	SSIM
$(K,\delta)=(1,1)$	(36)	$F_{\theta}^{\rm{mon}}$	$24.56(\pm 3.96)$	$0.80(\pm 0.11)$	$24.58(\pm 4.26)$	$0.80(\pm 0.11)$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$26.32(\pm 4.14)$	$0.85(\pm 0.04)$	$\boldsymbol{28.31(\pm 3.66)}$	$\boldsymbol{0.89(\pm 0.04)}$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$25.59(\pm 3.14)$	$0.87(\pm 0.07)$	$25.59(\pm 3.11)$	$0.87(\pm 0.07)$
$(K,\delta)=(5,1)$	(36)	$F_{\theta}^{\rm{mon}}$	$27.46(\pm 4.31)$	$0.87(\pm 0.08)$	$26.96(\pm 4.13)$	$0.86(\pm 0.08)$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$28.31(\pm 4.32)$	$0.89(\pm 0.06)$	$\boldsymbol{28.33(\pm 4.33)}$	$\boldsymbol{0.89(\pm 0.06)}$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$25.21(\pm 3.29)$	$0.86(\pm 0.08)$	$25.23(\pm 3.32)$	$0.86(\pm 0.08)$
$(K,\delta)=(1,0.6)$	(36)	$F_{\theta}^{\rm{mon}}$	$25.17(\pm 3.99)$	$0.81(\pm 0.10)$	$25.14(\pm 3.99)$	$0.81(\pm 0.10)$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$25.33(\pm 3.61)$	$0.81(\pm 0.07)$	$\boldsymbol{26.09(\pm 4.02)}$	$\boldsymbol{0.83(\pm 0.07)}$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$18.77(\pm 2.71)$	$0.77(\pm 0.12)$	$18.77(\pm 2.71)$	$0.77(\pm 0.12)$
$(K,\delta)=(5,0.6)$	(36)	$F_{\theta}^{\rm{mon}}$	$26.43(\pm 4.23)$	$\boldsymbol{0.84(\pm 0.09)}$	$\boldsymbol{26.63(\pm 4.32)}$	$\boldsymbol{0.84(\pm 0.09)}$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$24.75(\pm 4.33)$	$0.77(\pm 0.13)$	$24.73(\pm 4.32)$	$0.77(\pm 0.13)$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$18.40(\pm 2.74)$	$0.72(\pm 0.14)$	$18.40(\pm 2.74)$	$0.72(\pm 0.14)$

Table 2: Results for low noise level

\sigma=0.01

: Average PSNR values (and standard deviation), obtained over

10

images of BSD68, for solving the original inverse problem (31), with

K\in\{1,5\}

and

\delta\in\{1,0.6\}

. Results are shown when solving (36) with

F_{\theta}^{\rm{mon}}

, (38) with

\widetilde{F}^{\rm{mon}}_{\theta}

, and (38) with

\widetilde{F}^{\rm{lin}}_{\theta}

. PSNR and SSIM values correspond to the best results obtained when optimizing the regularization parameter

\varrho

	Problem	Operator	$\sigma_{\rm train}=0$		$\sigma_{\rm train}=0.01$
	Problem	Operator	PSNR	SSIM	PSNR	SSIM
$(K,\delta)=(1,1)$	(36)	$F_{\theta}^{\rm{mon}}$	$23.99(\pm 3.54)$	$0.76(\pm 0.13)$	$24.09(\pm 3.56)$	$0.76(\pm 0.13)$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$24.95(\pm 3.29)$	$0.78(\pm 0.08)$	$\boldsymbol{25.44(\pm 3.57)}$	$\boldsymbol{0.80(\pm 0.08)}$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$23.56(\pm 3.04)$	$0.79(\pm 0.11)$	$23.60(\pm 3.07)$	$0.79(\pm 0.11)$
$(K,\delta)=(5,1)$	(36)	$F_{\theta}^{\rm{mon}}$	$24.28(\pm 3.78)$	$0.77(\pm 0.13)$	$24.31(\pm 3.80)$	$0.77(\pm 0.13)$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$25.70(\pm 4.17)$	$0.81(\pm 0.10)$	$\boldsymbol{25.72(\pm 4.21)}$	$\boldsymbol{0.81(\pm 0.11)}$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$23.44(\pm 3.33)$	$0.78(\pm 0.12)$	$23.44(\pm 3.33)$	$0.78(\pm 0.12)$
$(K,\delta)=(1,0.6)$	(36)	$F_{\theta}^{\rm{mon}}$	$21.67(\pm 2.89)$	$0.73(\pm 0.15)$	$21.75(\pm 2.94)$	$0.73(\pm 0.15)$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$24.52(\pm 3.68)$	$0.77(\pm 0.11)$	$\boldsymbol{24.60(\pm 3.72)}$	$\boldsymbol{0.77(\pm 0.11)}$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$18.51(\pm 2.65)$	$0.72(\pm 0.14)$	$18.52(\pm 2.66)$	$0.72(\pm 0.14)$
$(K,\delta)=(5,0.6)$	(36)	$F_{\theta}^{\rm{mon}}$	$21.72(\pm 3.03)$	$0.73(\pm 0.15)$	$21.76(\pm 3.06)$	$0.73(\pm 0.15)$
	(38)	$\widetilde{F}^{\rm{mon}}_{\theta}$	$24.29(\pm 3.99)$	$\boldsymbol{0.76(\pm 0.12)}$	$\boldsymbol{24.30(\pm 4.00)}$	$\boldsymbol{0.76(\pm 0.12)}$
	(38)	$\widetilde{F}^{\rm{lin}}_{\theta}$	$18.36(\pm 2.71)$	$0.70(\pm 0.16)$	$18.36(\pm 2.71)$	$0.70(\pm 0.16)$

Table 3: Results for high noise level

\sigma=0.05

: Average PSNR values (and standard deviation), obtained over

10

images of BSD68, for solving the original inverse problem (31), with

K\in\{1,5\}

and

\delta\in\{1,0.6\}

. Results are shown when solving (36) with

F_{\theta}^{\rm{mon}}

, (38) with

\widetilde{F}^{\rm{mon}}_{\theta}

, and (38) with

\widetilde{F}^{\rm{lin}}_{\theta}

. PSNR and SSIM values correspond to the best results obtained when optimizing the regularization parameter

\varrho

4.1.5 Restoration results

In this section, we consider the original inverse problem (31), with model (32) for four cases: $(K,\delta)\in\{(1,1),(5,1),(1,0.6),(5,0.6)\}$ . We further consider two noise levels on the inverse problem: $\sigma=0.01$ (low noise level) and $\sigma=0.05$ (high noise level). Simulations are run on $10$ images from the BSD68 dataset. We compare the three methods described in Section 4.1.2 that hold convergence guarantees. Precisely, we solve (36) with $F_{\theta}^{\rm{mon}}$ , (38) with $\widetilde{F}^{\rm{mon}}_{\theta}$ , and (38) with $\widetilde{F}^{\rm{lin}}_{\theta}$ . For the regularization $r$ , we choose a smoothed Total Variation term [19] (See Appendix C for details), and we manually choose the regularization parameter $\varrho$ to achieve the best reconstruction quality for each method.

Quantitative results with average PSNR and SSIM values obtained for each method, in each setting, are reported in Tables 2 and 3, for $\sigma=0.01$ and $\sigma=0.05$ , respectively. We can observe overall that the proposed least-squares approach $\widetilde{F}^{\rm{mon}}_{\theta}$ always outperforms its linear counterpart $\widetilde{F}^{\rm{lin}}_{\theta}$ . It also performs better than $F_{\theta}^{\rm{mon}}$ for $K=1$ , and similarly for $K=5$ . Regarding the noise level for training $F_{\theta}^{\rm{mon}}$ and $\widetilde{F}^{\rm{mon}}_{\theta}$ , either choice $\sigma_{\rm train}=0$ or $\sigma_{\rm train}=0.01$ lead to very similar reconstruction results. Hence, the learned monotone approximation of $F$ does not seem affected by the noise level on the training dataset.

Qualitative results are presented in Figures 5 and 6 for the low noise level ( $\sigma=0.01$ ), $K=5$ with $\delta=1$ and $\delta=0.6$ , respectively. Each figure shows results for two images, obtained by selecting the regularization parameter $\varrho$ leading to the best reconstruction quality. Observations for these two images validate the quantitative results reported in Table 2. The linear least squares procedure using $\widetilde{F}^{\rm lin}$ leads to the worst performance due to method failing to correct for the saturation function. For the high saturation case $\delta=0.6$ , we see that $F_{\theta}^{\rm{mon}}$ leads to slightly sharper reconstructions, while both least-squares versions $\widetilde{F}^{\rm{mon}}_{\theta}$ and $\widetilde{F}^{\rm{lin}}_{\theta}$ seem to over-smooth the reconstruction, possibly due to the bias introduced with respect to the true saturation model.

We also present the convergence profiles $\|x_{k+1}-x_{k}\|_{2}/\|y\|_{2}$ with respect to the iterations $k$ of the restored images using the three considered models $F_{\theta}^{\rm{mon}}$ , $\widetilde{F}^{\rm{mon}}_{\theta}$ and $\widetilde{F}^{\rm{lin}}_{\theta}$ , with two training noise levels $\sigma_{\rm train}\in\{0,0.01\}$ . We observe that all methods exhibit a converging behaviour. The oscillations are due to the Goldstein-Armijo rule that introduces a backtracking line search procedure in the algorithms to find the optimal step size. Smoother curves can be obtained by lowering the $\beta$ parameters in (2.2), with the caveat of more processing time due to the additional steps performed within the backtracking line search procedure.

5 Conclusions

In this article, we introduced a novel approach for training monotone neural networks, by designing a penalized training procedure. The resulting monotone networks can then be embedded in the FBF algorithm, within a PnP framework, yet ensuring the convergence of the iterates of the resulting algorithm. This method can be leveraged for addressing a wide range of monotone inclusion problems, including in imaging, as demonstrated in the context of solving non-linear inverse imaging problems.

The proposed PnP-FBF method enables solving generic constrained monotone inclusion problem. Its convergence is ensured even if the involved operators are not cocoercive. Hence the proposed method can be used for a wider class of operators compared to other similar iterative schemes such as the (PnP) forward-backward method. Moreover, combined with the Armijo-Goldstein rule, the proposed FBF-PnP algorithm is guaranteed to converge without needing explicit computation of the Lipschitz constant of the neural network.

The proposed penalized training procedure enables to learn monotone neural networks. To do so, we designed a differentiable penalization function, relying on properties of the network Jacobian, that can be implemented efficiently using auto-differentiation tools. The proposed training approach is very flexible and can be applied to a wide range of network architectures and training paradigms.

We finally illustrated the benefit of the proposed framework in learning monotone operators for semi-blind non-linear deconvolution imaging problems. Our methodology demonstrated an accurate monotone approximation of the true non-monotone degradation function. We show that the monotonocity of the learned network is further instrumental within the proposed PnP-FBF scheme for inverting the network, and solving the image recovery problem.

The versatility of our methodology allows its potential extension to various imaging problems involving a nonlinear model, including super-resolution. The practical utility of monotone operators has already been explored in normalizing flows, offering possibilities for extending our work to monotone normalizing flows [1]. Further investigations into style transfer tasks, leveraging the invertibility property, could yield stable neural networks for image-to-image mapping problems. Finally, the application of non-linear inverse problems in tomography, acknowledging the inherent saturation in tissue absorption of X-rays, offers a promising avenue for future exploration.

Acknowledgments

We would like to thank DATAIA institute from University Paris-Saclay for funding the thesis of YB. The work of JCP was founded by the ANR Research and Teaching Chair in Artificial Intelligence, BRIDGEABLE. The work of AR was partly funded by the Royal Society of Edinburgh; EPSRC grant EP/X028860/1; the CVN, Inria/OPIS, and CentraleSupélec.

References

[1] B. Ahn, C. Kim, Y. Hong, and H. J. Kim, Invertible monotone operators for normalizing flows, Adv. Neural Inf. Process. Syst., 35 (2022), pp. 16836–16848.
[2] A. H. Al-Shabili, X. Xu, I. Selesnick, and U. S. Kamilov, Bregman plug-and-play priors, in Proc. - Int. Conf. Image Process. ICIP, IEEE, 2022, pp. 241–245.
[3] M. S. C. Almeida and M. A. T. Figueiredo, Blind image deblurring with unknown boundaries using the alternating direction method of multipliers, in Proc. - Int. Conf. Image Process. ICIP, Sept. 2013, pp. 586–590.
[4] H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, 2017.
[5] M. Benning and M. Burger, Modern regularization methods for inverse problems, Acta Numer., 27 (2018), pp. 1–111.
[6] J. Bolte and E. Pauwels, Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning, Math. Prog., 188 (2021), pp. 19–51.
[7] C. A. Bouman and G. T. Buzzard, Generative plug and play: Posterior sampling for inverse problems, arXiv preprint arXiv:2306.07233, (2023).
[8] D. Calvetti and E. Somersalo, Inverse problems: From regularization to Bayesian inference, Wiley Interdiscip. Rev.: Comput. Stat., 10 (2018).
[9] S. Chaudhari, S. Pranav, and J. M. Moura, Learning Gradients of Convex Functions with Monotone Gradient Networks, in Proc. - ICASSP IEEE Int. Conf. Acoust. Speech Signal Process., June 2023, pp. 1–5.
[10] R. Cohen, M. Elad, and P. Milanfar, Regularization by denoising via fixed-point projection (RED-PRO), SIAM J. Imaging Sci., 14 (2021), pp. 1374–1406.
[11] F. Colibazzi, D. Lazzaro, S. Morigi, and A. Samoré, Deep-plug-and-play proximal gauss-newton method with applications to nonlinear, ill-posed inverse problems, Inverse Probl. Imaging, (2023), pp. 1226–1248.
[12] P. L. Combettes and J.-C. Pesquet, A Douglas–Rachford Splitting Approach to Nonsmooth Convex Variational Signal Recovery, IEEE J. Sel. Top. Signal Process., 1 (2007), pp. 564–574.
[13] P. L. Combettes and J.-C. Pesquet, Proximal Splitting Methods in Signal Processing, in Fixed-Point Algorithms for Inverse Problems in Science and Engineering, H. H. Bauschke, R. S. Burachik, P. L. Combettes, V. Elser, D. R. Luke, and H. Wolkowicz, eds., Springer Optimization and Its Applications, Springer, New York, NY, 2011, pp. 185–212.
[14] P. L. Combettes and J.-C. Pesquet, Lipschitz Certificates for Layered Network Structures Driven by Averaged Activation Operators, SIAM J. Math. Data Sci., 2 (2020), pp. 529–557.
[15] P. L. Combettes and J.-C. Pesquet, Fixed Point Strategies in Data Science, IEEE Trans. Signal Process., 69 (2021), pp. 3878–3905.
[16] H. Daniels and M. Velikova, Monotone and Partially Monotone Neural Networks, IEEE Trans. Neural Netw. Learn. Syst., 21 (2010), pp. 906–917.
[17] R. Fermanian, M. L. Pendu, and C. Guillemot, Learned gradient of a regularizer for plug-and-play gradient descent, arXiv preprint arXiv:2204.13940, (2022).
[18] C. S. Garcia, M. Larchevêque, S. O’Sullivan, M. Van Waerebeke, R. R. Thomson, A. Repetti, and J.-C. Pesquet, A primal-dual data-driven method for computational optical imaging with a photonic lantern, arXiv preprint arXiv:2306.11679, (2023).
[19] P. Getreuer, Rudin-Osher-Fatemi Total Variation Denoising using Split Bregman, Image Process. Line, 2 (2012), pp. 74–95.
[20] S. Hurault, U. Kamilov, A. Leclaire, and N. Papadakis, Convergent bregman plug-and-play image restoration for poisson inverse problems, Adv. Neural Inf. Process. Syst., 36 (2024).
[21] S. Hurault, A. Leclaire, and N. Papadakis, Gradient step denoiser for convergent plug-and-play, arXiv preprint arXiv:2110.03220, (2021).
[22] S. Hurault, A. Leclaire, and N. Papadakis, Proximal denoiser for convergent plug-and-play optimization with nonconvex regularization, in International Conference on Machine Learning, PMLR, 2022, pp. 9483–9505.
[23] D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling, Improving Variational Inference with Inverse Autoregressive Flow, arXiv:1606.04934 [cs, stat], (2017).
[24] J. Kruse, L. Ardizzone, C. Rother, and U. Köthe, Benchmarking invertible architectures on inverse problems, arXiv preprint arXiv:2101.10763, (2021).
[25] D. G. Luenberger, Y. Ye, et al., Linear and nonlinear programming, vol. 2, Springer, 1984.
[26] A. Neacşu, J.-C. Pesquet, and C. Burileanu, EMG-based automatic gesture recognition using lipschitz-regularized neural networks, ACM Trans. Intell. Syst. Technol., 15 (2024).
[27] J.-C. Pesquet, A. Repetti, M. Terris, and Y. Wiaux, Learning maximally monotone operators for image recovery, SIAM J. Imaging Sci., 14 (2021), pp. 1206–1237.
[28] D. J. Rezende and S. Mohamed, Variational Inference with Normalizing Flows, arXiv:1505.05770 [cs, stat], (2016).
[29] D. Runje and S. M. Shankaranarayana, Constrained monotonic neural networks, in International Conference on Machine Learning, PMLR, 2023, pp. 29338–29353.
[30] E. Ryu, J. Liu, S. Wang, X. Chen, Z. Wang, and W. Yin, Plug-and-play methods provably converge with properly trained denoisers, in International Conference on Machine Learning, PMLR, 2019, pp. 5546–5557.
[31] R. M. Schmidt, F. Schneider, and P. Hennig, Descending through a crowded valley-benchmarking deep learning optimizers, in International Conference on Machine Learning, PMLR, 2021, pp. 9367–9376.
[32] M. Tang and A. Repetti, A data-driven approach for bayesian uncertainty quantification in imaging, arXiv preprint arXiv:2304.11200, (2023).
[33] M. Terris, A. Repetti, J.-C. Pesquet, and Y. Wiaux, Building firmly nonexpansive convolutional neural networks, in Proc. - ICASSP IEEE Int. Conf. Acoust. Speech Signal Process., Proc. - ICASSP IEEE Int. Conf. Acoust. Speech Signal Process., Barcelona, Spain, May 2020, pp. 8658–8662.
[34] M. Terris, A. Repetti, J.-C. Pesquet, and Y. Wiaux, Enhanced Convergent PNP Algorithms For Image Restoration, in Proc. - Int. Conf. Image Process. ICIP, Sept. 2021, pp. 1684–1688.
[35] P. Tseng, A Modified Forward-Backward Splitting Method for Maximal Monotone Mappings, SIAM J. Control Optim., 38 (2000), pp. 431–446.
[36] W. Wang, Z. Dang, Y. Hu, P. Fua, and M. Salzmann, Backpropagation-friendly eigendecomposition, Adv. Neural Inf. Process. Syst., 32 (2019).
[37] Z. Zhang, Q. Liu, and Y. Wang, Road Extraction by Deep Residual U-Net, IEEE Geosci. Remote Sens. Lett., 15 (2018), pp. 749–753.

Appendix A Saturation function

We modify the classical Tanh function to be centered on 0.5 and to take as input images with pixel values in the range $[0,1]$ . We also introduce a scaling factor $\delta$ to accentuate the non-linearity of the tanh (Figure 7). The resulting function is then used for saturation, and is defined as

(40)

S_{\delta}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}\colon\boldsymbol{x}=(x_{i})_{1% \leqslant i\leqslant n}\mapsto(\psi_{\delta}(x_{i}))_{1\leqslant i\leqslant n}

where

(41)

\psi_{\delta}\colon\mathbb{R}\to\mathbb{R}\colon x\mapsto\frac{\text{tanh}(% \delta(2x-1))+1}{2}.

Examples for $\delta\in\{0.6,1\}$ are shown in Figure 7.

Appendix B Linear approximation of a convolution filter

As explained in Section 4.1.2, we compare our model to a learned linear approximation $F_{\theta}^{\rm{lin}}$ of $F$ . To this aim, we define $F_{\theta}^{\rm{lin}}$ as a 2D convolution layer, with convolution kernel $f_{\theta}\in\mathbb{R}^{D}$ . Then, we train $F_{\theta}^{\rm{lin}}$ by solving

(42)

\operatorname*{argmin}_{\theta\in\mathbb{R}^{D}}\displaystyle\sum_{(x,y)\in% \mathbb{D}_{\rm train}}\|\upsilon(f_{\theta})\ast x-y\|_{1},

where

(43)

\displaystyle\upsilon\colon\mathbb{R}^{D}\to\mathbb{R}^{D}\colon f_{\theta}% \mapsto\frac{\text{ReLU}(f_{\theta})}{\sum_{i\in 1}^{D}\max(0,f_{\theta_{i}})}.

This parametrization is introduced as we assume that the true convolution kernels in (32) are normalized and nonnegative.

Problem (42) is then solved using Adam optimizer for 200 epochs, with a learning rate of $0.02$ . The size of the kernels was chosen equal to the size of the true kernels in (32) for simplicity (i.e., $D=9\times 9$ ). The learned kernels are displayed in Figure 8. We can observe that for the normal saturation levels ( $\delta=1$ ), the approximated kernels are almost equal to the true ones, with or without noise. For the high saturation level ( $\delta=0.6$ ) however, the learned kernels are very different from the true ones, certainly to try to compensate for the nonlinear saturation function.

Appendix C Total variation regularization

We provide here the details of the smoothed approximation used in the manuscript for the Total Variation [19] function. We define

(44)

(\forall x\in\mathbb{R}^{n})\quad r(x)=\sum_{i=1}^{n}\sqrt{[\nabla_{h}x]_{i}^{% 2}+[\nabla_{v}x]_{i}^{2}+\varepsilon}

where $\nabla_{h}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ and $\nabla_{v}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ model the linear horizontal and vertical gradient operators, respectively, and $\varepsilon>0$ is a smoothing parameter. Then, the gradient of $r$ is given by

(45)

\displaystyle(\forall x\in\mathbb{R}^{n})\quad\nabla r(x)=\left(\dfrac{\Big{[}% \nabla_{h}^{\top}(\nabla_{h}x)+\nabla_{v}^{\top}(\nabla_{v}x)\Big{]}_{i}}{% \sqrt{[\nabla_{h}x]_{i}^{2}+[\nabla_{v}x]_{i}^{2}+\varepsilon}}\right)_{1% \leqslant i\leqslant n}.

	$\displaystyle(\forall x\in C)\quad\varrho(x)$	$\displaystyle=\inf_{t\in N_{C}(x)}\\|B(x)+t\\|$
		$\displaystyle\leqslant\\|B(x)\\|+\inf_{t\in N_{C}(x)}\\|t\\|$
		$\displaystyle\leqslant\\|B(x)\\|+\\|x-\text{\rm proj}_{C}x\\|$
(11)			$\displaystyle=\\|B(x)\\|.$

(17)		$\displaystyle\lim_{\begin{subarray}{c}\alpha\to 0\\ \alpha>0\end{subarray}}{\left\langle{{h}\mid{\frac{T(x+\alpha h)-T(x)}{\alpha}% }}\right\rangle}\geqslant\beta\\|h\\|^{2}$	$\displaystyle\quad\Leftrightarrow\quad{\left\langle{{h}\mid{\operatorname{J}_{% T}(x)h}}\right\rangle}\geqslant\beta\\|h\\|^{2}$
(18)			$\displaystyle\quad\Leftrightarrow\quad{\left\langle{{h}\mid{\operatorname{J}_{% T}^{\rm s}(x)h}}\right\rangle}\geqslant\beta\\|h\\|^{2}.$


$L_{1}$	$L_{2}$	$L_{3}$	$L_{4}$	$L_{5}$


$y=F(\overline{x})$	$y_{{F}^{\rm lin}}={F}^{\rm lin}(\overline{x})$	$y_{F_{\theta}^{\rm{lin}}}=F_{\theta}^{\rm{lin}}(\overline{x})$
PSNR $=24.07$	MAE( $y,y_{{F}^{\rm lin}}$ ) $=0.035$	MAE( $y,y_{F_{\theta}^{\rm{lin}}}$ ) $=0.035$
$(\lambda_{\min},\lambda_{\max})=(-80.52,81.37)$	$(\lambda_{\min},\lambda_{\max})=(-0.09,1.00)$	$(\lambda_{\min},\lambda_{\max})=(-0.14,0.99)$

$y_{F_{\theta}^{\rm{nom}}}=F_{\theta}^{\rm{nom}}(\overline{x})$	$y_{F_{\theta}^{\rm{mon}}}=F_{\theta}^{\rm{mon}}(\overline{x})$	$y_{\widetilde{F}^{\rm{mon}}_{\theta}}=F_{\theta}(\overline{x})$
MAE( $y,y_{F_{\theta}^{\rm{nom}}}$ ) $=0.002$	MAE( $y,y_{F_{\theta}^{\rm{mon}}}$ ) $=0.004$	MAE( $y,y_{\widetilde{F}^{\rm{mon}}_{\theta}}$ ) $=0.004$
$(\lambda_{\min},\lambda_{\max})=(-0.21,1.01)$	$(\lambda_{\min},\lambda_{\max})=(0.01,1.00)$	$(\lambda_{\min},\lambda_{\max})=(0.00,0.87)$

$y=F(\overline{x})$	$y_{{F}^{\rm lin}}={F}^{\rm lin}(\overline{x})$	$y_{F_{\theta}^{\rm{lin}}}=F_{\theta}^{\rm{lin}}(\overline{x})$
PSNR $=21.17$	MAE( $y,y_{{F}^{\rm lin}}$ ) $=0.037$	MAE( $y,y_{F_{\theta}^{\rm{lin}}}$ ) $=0.037$
$(\lambda_{\min},\lambda_{\max})=(-80.52,81.37)$	$(\lambda_{\min},\lambda_{\max})=(-0.09,1.00)$	$(\lambda_{\min},\lambda_{\max})=(-0.14,0.99)$

$y_{F_{\theta}^{\rm{nom}}}=F_{\theta}^{\rm{nom}}(\overline{x})$	$y_{F_{\theta}^{\rm{mon}}}=F_{\theta}^{\rm{mon}}(\overline{x})$	$y_{\widetilde{F}^{\rm{mon}}_{\theta}}=F_{\theta}(\overline{x})$
MAE( $y,y_{F_{\theta}^{\rm{nom}}}$ ) $=0.003$	MAE( $y,y_{F_{\theta}^{\rm{mon}}}$ ) $=0.005$	MAE( $y,y_{\widetilde{F}^{\rm{mon}}_{\theta}}$ ) $=0.006$
$(\lambda_{\min},\lambda_{\max})=(-0.21,1.01)$	$(\lambda_{\min},\lambda_{\max})=(0.01,1.00)$	$(\lambda_{\min},\lambda_{\max})=(0.00,0.92)$


$y=F(\overline{x})$	$y_{{F}^{\rm lin}}={F}^{\rm lin}(\overline{x})$	$y_{F_{\theta}^{\rm{lin}}}=F_{\theta}^{\rm{lin}}(\overline{x})$
PSNR $=17.33$	MAE( $y,y_{{F}^{\rm lin}}$ ) $=0.110$	MAE( $y,y_{F_{\theta}^{\rm{lin}}}$ ) $=0.108$
$(\lambda_{\min},\lambda_{\max})=(-156.86,157.59)$	$(\lambda_{\min},\lambda_{\max})=(-0.09,1.00)$	$(\lambda_{\min},\lambda_{\max})=(-0.35,0.99)$

$y_{F_{\theta}^{\rm{nom}}}=F_{\theta}^{\rm{nom}}(\overline{x})$	$y_{F_{\theta}^{\rm{mon}}}=F_{\theta}^{\rm{mon}}(\overline{x})$	$y_{\widetilde{F}^{\rm{mon}}_{\theta}}=F_{\theta}(\overline{x})$
MAE( $y,y_{F_{\theta}^{\rm{nom}}}$ ) $=0.004$	MAE( $y,y_{F_{\theta}^{\rm{mon}}}$ ) $=0.006$	MAE( $y,y_{\widetilde{F}^{\rm{mon}}_{\theta}}$ ) $=0.003$
$(\lambda_{\min},\lambda_{\max})=(-0.18,0.63)$	$(\lambda_{\min},\lambda_{\max})=(0.03,0.61)$	$(\lambda_{\min},\lambda_{\max})=(0.00,0.55)$

$y=F(\overline{x})$	$y_{{F}^{\rm lin}}={F}^{\rm lin}(\overline{x})$	$y_{F_{\theta}^{\rm{lin}}}=F_{\theta}^{\rm{lin}}(\overline{x})$
PSNR $=17.33$	MAE( $y,y_{{F}^{\rm lin}}$ ) $=0.110$	MAE( $y,y_{F_{\theta}^{\rm{lin}}}$ ) $=0.111$
$(\lambda_{\min},\lambda_{\max})=(-156.86,157.59)$	$(\lambda_{\min},\lambda_{\max})=(-0.09,1.00)$	$(\lambda_{\min},\lambda_{\max})=(-0.35,1.00)$

$y_{F_{\theta}^{\rm{nom}}}=F_{\theta}^{\rm{nom}}(\overline{x})$	$y_{F_{\theta}^{\rm{mon}}}=F_{\theta}^{\rm{mon}}(\overline{x})$	$y_{\widetilde{F}^{\rm{mon}}_{\theta}}=F_{\theta}(\overline{x})$
MAE( $y,y_{F_{\theta}^{\rm{nom}}}$ ) $=0.009$	MAE( $y,y_{F_{\theta}^{\rm{mon}}}$ ) $=0.009$	MAE( $y,y_{\widetilde{F}^{\rm{mon}}_{\theta}}$ ) $=0.004$
$(\lambda_{\min},\lambda_{\max})=(-0.18,0.63)$	$(\lambda_{\min},\lambda_{\max})=(0.03,0.61)$	$(\lambda_{\min},\lambda_{\max})=(0.00,0.56)$


$\overline{x}$	$y=F(\overline{x})$	$y=F(\overline{x})$	$\widetilde{y}={F}^{\rm lin}(\overline{x})$
	$\delta=1$	$\delta=0.6$


$\overline{x}$	$y$ – $(23.96,0.89)$	$\overline{x}$	$y$ – $(21.11,0.65)$

$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(31.59,0.93)$	$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(30.98,0.93)$	$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(24.95,0.80)$	$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(24.67,0.80)$
$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$

$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(31.84,0.94)$	$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(32.10,0.94)$	$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(25.75,0.83)$	$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(25.72,0.83)$
$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$

$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(24.67,0.93)$	$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(24.69,0.93)$	$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(22.77,0.75)$	$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(22.79,0.75)$
$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$


$\overline{x}$	$y$ – $(17.30,0.85)$	$\overline{x}$	$y$ – $(16.27,0.54)$

$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(30.64,0.92)$	$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(30.59,0.92)$	$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(24.12,0.77)$	$\widehat{x}_{F_{\theta}^{\rm{mon}}}$ – $(24.21,0.78)$
$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$

$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(28.91,0.89)$	$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(28.87,0.89)$	$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(22.55,0.7)$	$\widehat{x}_{\widetilde{F}^{\rm{mon}}_{\theta}}$ – $(22.54,0.7)$
$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$

$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(17.40,0.86)$	$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(17.40,0.86)$	$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(16.42,0.54)$	$\widehat{x}_{\widetilde{F}^{\rm{lin}}_{\theta}}$ – $(16.41,0.54)$
$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$

	True kernel	$\delta=1$		$\delta=0.6$
	True kernel	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$	$\sigma_{\rm train}=0$	$\sigma_{\rm train}=0.01$
$K=1$
$K=5$

Learning truly monotone operators with applications to nonlinear inverse problems††thanks: Submitted to the editors DATE.

Abstract

keywords:

1 Introduction

Problem 1.1.

Notation

2 Tseng’s algorithms for monotone inclusion problems

2.1 Forward-Backward-Forward strategy

Problem 2.1.

Definition 2.2 (Armijo-Goldstein rule).

Proposition 2.3.

2.2 Modified algorithm

Proposition 2.4.

Proof 2.5.

3 Proposed method using a monotone NN

3.1 Properties of differentiable monotone operators

Lemma 3.1.

Proof 3.2.

Proposition 3.3.

Proof 3.4.

Remark 3.5.

3.2 Proposed regularization approach for training monotone NNs

3.3 Penalized training implementation

Penalization computation

Remark 3.6.

Generic training framework

4 Learning non-linear model approximations

4.1 Imaging semi-blind non-linear inverse problem

4.1.1 Measurement operator setting

4.1.2 Monotone inclusion formulation to solve (31)

Direct regularized approach

Proposition 4.1.

Proof 4.2.

Remark 4.3.

Least-squares regularized approach

Remark 4.4.

Learned approximations to F𝐹Fitalic_F

Remark 4.5.

4.1.3 Model and training procedure

NN architecture

Training dataset

Training procedure

4.1.4 Training results

4.1.5 Restoration results

5 Conclusions

Acknowledgments

References

Appendix A Saturation function

Appendix B Linear approximation of a convolution filter

Appendix C Total variation regularization

Learning truly monotone operators with applications to nonlinear inverse problems^†^†thanks: Submitted to the editors DATE.

Learned approximations to $F$