Prof. François Malgouyres

Arthur -Dont Les Gâteaux

Cvn -Alp

Siddhartha Loubna

Evgenios, Guillaume, Hari, Khuê, Stefan Édouard -À La Nouvelle Génération -Matthieu

Je voudrais commencer par remercier mes encadrants Prof. Jean-Christophe Pesquet et Prof. Émilie Chouzenoux pour leurs précieux conseils prodigués durant ces trois années de thèse. Leur rigueur scientifique et l'énergie qu'ils investissent dans leur travail continueront à me servir de modèle et de motivation pour la suite de ma carrière.

Je tiens également à remercier les membres de mon jury, qui ont accepté de consacrer du temps à l'évaluation de mon travail de thèse. Merci au président du jury Prof. Saïd Moussaoui, à

Mots clés : algorithmes proximaux, points intérieurs, algorithme déroulé, imagerie hyperspectrale, imagerie ultrasonore, algorithme de Langevin Résumé : Les problèmes inverses en traitement d'images peuvent être résolus en utilisant des méthodes variationnelles classiques, des approches basées sur l'apprentissage profond, ou encore des stratégies bayésiennes. Bien que différentes, ces approches nécessitent toutes des algorithmes d'optimisation efficaces. L'opérateur proximal est un outil important pour la minimisation de fonctions non lisses. Dans cette thèse, nous illustrons la polyvalence des algorithmes proximaux en les introduisant dans chacune des trois méthodes de résolution susmentionnées. Tout d'abord, nous considérons une formulation variationnelle sous contraintes dont la fonction objectif est composite. Nous développons PIPA, un nouvel algorithme proximal de points intérieurs permettant de résoudre ce problème. Dans le but d'accélérer PIPA, nous y incluons une métrique variable. La convergence de PIPA est prouvée sous certaines conditions et nous montrons que cette méthode est plus rapide que des algorithmes de l'état de l'art au travers de deux exemples numériques en traitement d'images. Dans une deuxième partie, nous étudions iRestNet, une architecture neuronale obtenue en déroulant un algorithme proximal de points intérieurs. iRestNet nécessite l'expression de l'opérateur proximal de la barrière logarithmique et des dérivées premières de cet opérateur. Nous fournissons ces expressions pour trois types de contraintes. Nous montrons ensuite que sous certaines conditions, cette architecture est robuste à une perturbation sur son entrée. Enfin, iRestNet démontre de bonnes performances pratiques en restauration d'images par rapport à une approche variationnelle et à d'autres méthodes d'apprentissage profond. La dernière partie de cette thèse est consacrée à l'étude d'une méthode d'échantillonnage stochastique pour résoudre des problèmes inverses dans un cadre bayésien. Nous proposons une version accélérée de l'algorithme proximal de Langevin non ajusté, baptisée PP-ULA. Cet algorithme est incorporé à un échantillonneur de Gibbs hybride utilisé pour réaliser la déconvolution et la segmentation d'images ultrasonores. PP-ULA utilise le principe de majorationminimisation afin de gérer les distributions non log-concaves. Comme le montrent nos expériences réalisées sur des données ultrasonores simulées et réelles, PP-ULA permet une importante réduction du temps d'exécution tout en produisant des résultats de déconvolution et de segmentation très satisfaisants.

Title : Proximal and Interior Point Optimization Strategies in Image Recovery

Keywords : proximal algorithms, interior points, unfolded algorithm, hyperspectral imaging, ultrasound imaging, Langevin-based schemes Abstract : Inverse problems in image processing can be solved by diverse techniques, such as classical variational methods, recent deep learning approaches, or Bayesian strategies. Although relying on different principles, these methods all require efficient optimization algorithms. The proximity operator appears as a crucial tool in many iterative solvers for nonsmooth optimization problems. In this thesis, we illustrate the versatility of proximal algorithms by incorporating them within each one of the aforementioned resolution methods. First, we consider a variational formulation including a set of constraints and a composite objective function. We present PIPA, a novel proximal interior point algorithm for solving the considered optimization problem. This algorithm includes variable metrics for acceleration purposes. We derive convergence guarantees for PIPA and show in numerical experiments that it compares favorably with state-of-the-art algorithms in two challenging image processing applications. In a second part, we investigate a neural network architecture called iRestNet, obtained by unfolding a proxi-mal interior point algorithm over a fixed number of iterations. iRestNet requires the expression of the logarithmic barrier proximity operator and of its first derivatives, which we provide for three useful types of constraints. Then, we derive conditions under which this optimization-inspired architecture is robust to an input perturbation. We conduct several image deblurring experiments, in which iRestNet performs well with respect to a variational approach and to state-of-theart deep learning methods. The last part of this thesis focuses on a stochastic sampling method for solving inverse problems in a Bayesian setting. We present an accelerated proximal unadjusted Langevin algorithm called PP-ULA. This scheme is incorporated into a hybrid Gibbs sampler used to perform joint deconvolution and segmentation of ultrasound images. PP-ULA employs the majorizeminimize principle to address non log-concave priors. As shown in numerical experiments, PP-ULA leads to a significant time reduction and to very satisfactory deconvolution and segmentation results on both simulated and real ultrasound data.

Université Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

Résumé

Cette thèse de doctorat porte sur le développement de méthodes de résolution pour les problèmes inverses rencontrés en traitement d'images. Ce domaine connaît d'importants changements en raison de l'émergence de méthodes dites d'apprentissage profond, très performantes pour la vision par ordinateur. Ces méthodes récentes concurrencent des approches plus traditionnelles telles que les méthodes variationnelles ou les stratégies bayésiennes. Bien que différentes, les approches mentionnées précédemment ont toutes recours à des algorithmes d'optimisation. Il apparaît donc nécessaire de développer des solveurs itératifs fiables et efficaces, en particulier dans un contexte où la taille des bases de données ainsi que des données ellesmêmes, par exemple dans le cas d'applications 3D, ne cesse de croître. Des outils mathématiques spécifiques sont alors nécessaires. Un exemple d'outil de ce type est l'opérateur proximal, qui joue un rôle important dans la minimisation de fonction non lisses ainsi que dans la résolution de problèmes de grande taille. Cette opérateur a permis le développement d'algorithmes présentant de bonnes performances pratiques et dont la convergence est garantie sous certaines conditions. Dans cette thèse, nous illustrons la polyvalence des algorithmes proximaux au moyen de trois méthodes combinant l'opérateur proximal et des techniques utilisées pour le traitement d'images.

Les méthodes variationnelles sont des techniques classiques qui découlent des modèles utilisés dans la formulation de problèmes inverses; elles se caractérisent par l'étude approfondie de leurs propriétés mathématiques. Ces stratégies consistent à formuler une fonction objectif dont le minimiseur est un bon estimé de l'image ou du signal recherché. Ces méthodes offrent la possibilité d'incorporer au problème d'optimisation des informations préalablement connues sur la solution ou des propriétés que cette dernière doit satisfaire. C'est le cas en particulier des méthodes de points intérieurs, qui sont bien adaptées à la résolution de problèmes d'optimisation sous contraintes présentant une fonction objectif différentiable. Dans le Chapitre 3, nous proposons de recourir à l'opérateur proximal afin d'étendre le champ d'application des méthodes de points intérieurs. Nous considérons une formulation variationnelle dont la fonction objectif est composite et comprend vii viii des contraintes. Nous présentons un nouvel algorithme proximal de points intérieurs pour résoudre ce problème d'optimisation. L'algorithme proposé, appelé PIPA, inclut une barrière logarithmique qui permet aux contraintes d'être toujours strictement satisfaites. Nous utilisons également une métrique variable qui permet d'accélérer la convergence de PIPA. Nous obtenons des garanties de convergence pour l'algorithme proposé, puis nous présentons des résultats numériques démontrant les bonnes performances de PIPA par rapport à d'autre algorithmes dans deux applications en traitement d'image: le démélange hyperspectral et la décomposition en géométrie-texture et reconstruction jointes de données tomographiques.

La création de bases de données de grande taille a permis le développement de méthodes d'apprentissage profond qui atteignent l'état de l'art pour des tâches en vision par ordinateur telles que la segmentation d'images ou la classification. Le développement de ces techniques s'accompagne d'un questionnement concernant leur fiabilité et leur stabilité. Les réseaux de neurones comportent, par nature, de nombreux degrés de liberté, ce qui rend leur analyse difficile et leur interprétation limitée. Malgré des stratégies différentes, les approches variationnelles et celles d'apprentissage profond ont des avantages et des inconvénients complémentaires qui permettent de les combiner. Nous proposons de mettre en oeuvre cette idée dans le Chapitre 4 de cette thèse, où nous étudions un réseau de neurones appelé iRestNet, et dont l'architecture est obtenue en déroulant un algorithme proximal de points intérieurs. iRestNet requiert les expressions de l'opérateur proximal de la barrière logarithmique et de ses dérivées premières. Nous fournissons ces expressions pour trois types de contraintes usuelles. Nous montrons ensuite que sous certaines conditions, cette architecture inspirée de l'optimisation est robuste à une perturbation sur son entrée. Enfin, nous procédons à des expériences numériques en défloutage d'images. iRestNet démontre de bonnes performances par rapport à l'approche variationnelle classique et à des méthodes d'apprentissage profond, y compris d'autre architectures issues d'algorithmes déroulés.

La dernière partie de cette thèse est consacrée à l'étude d'une méthode d'échantillonnage stochastique pour résoudre des problèmes inverses dans un cadre bayésien. De même que les méthodes d'apprentissage profond, les techniques de simulation stochastiques estiment de manière automatique les hyperparamètres d'un problème inverse ainsi que ses variables d'intérêt. Les approches bayésiennes nécessitent la connaissance des distributions de probabilité a priori des différentes variables du modèle considéré. Ainsi, comme pour les méthodes variationnelles, la solution peut être amenée à satisfaire certaines propriétés désirées suivant le choix de ces distributions. Les méthodes de simulation stochastiques ont l'avantage de présenter des garanties de convergence, cependant elles souffrent parfois d'un temps d'exécution très long qui peut rendre impossible leur utilisation pour des applications en temps réel. Dans le Chapitre 5, nous proposons une version accélérée de ix l'algorithme proximal de Langevin non ajusté, baptisée PP-ULA. Cet algorithme est incorporé à un échantillonneur de Gibbs hybride utilisé pour réaliser la déconvolution et la segmentation d'images ultrasonores. PP-ULA utilise le principe de majoration-minimisation afin de gérer les distributions non log-concaves. Comme le montrent nos expériences réalisées sur des données ultrasonores simulées et réelles, PP-ULA permet une importante réduction du temps d'exécution tout en produisant des résultats de déconvolution et de segmentation très satisfaisants.

Dans le Chapitre 6, nous résumons nos principales contributions et nous proposons plusieurs pistes pour étendre les résultats présentés dans cette thèse.

Notation R, R + , R * : sets of real, positive real and non-zero real numbers R m

: set of vectors with m entries R m×n : set of matrices with m rows and n columns N, N * : sets of positive and non-zero positive integers S m : set of symmetric matrices in R m×m S + m : set of symmetric positive-definite matrices in R m×m : Loewner partial order for matrices s : scalars will be denoted by lowercase letters v : vectors will be denoted by lowercase bold letters M : matrices will be denoted by uppercase bold letters M , M -1 : transpose and inverse of M , respectively I m : square identity matrix in R m×m 1 m : vector of R m with all entries equal to 1 0 m×n : matrix in R m×n with all entries equal to 0 0 m : vector of R m with all entries equal to 0 In this work, we develop and investigate resolution methods for inverse problems, and more specifically for applications related to signal and image processing. This field has undergone major changes due to the emergence of very powerful deep learning methods, which compete with more traditional approaches, such as variational strategies or Bayesian methods. At the crossroads of the aforementioned approaches lies optimization. Optimization algorithms occur either in the design or implementation of the methods used to solve inverse problems. Consequently, there is a need for reliable and adequate iterative solvers. In particular, dealing with the evergrowing size of databases and datapoints, for instance in 3D applications, requires powerful mathematical tools. The proximity operator notably plays an important role in large-scale and nondifferentiable optimization. It has led to competitive algorithms with convergence guarantees, available convergence rates and noteworthy practical performance. The goal of this thesis is to illustrate the versatility of proximal algorithms in the context of inverse problems. To that end, we develop three methods combining the proximity operator with different techniques.

v i : ith coefficient of v M i
Variational methods are classical model-based methods with well-known mathematical properties. The variational strategy consists in designing an objective function whose minimizer is an appropriate estimate of the desired signal or image. Hence, prior knowledge about the solution can be incorporated into the resulting optimization problem. This applies to interior point methods, which are well suited for constrained problems in linear programming. In Chapter 3, we show that the scope of the logarithmic barrier method can be extended by combining it with a proximal iterative scheme.

The generalization of data collection and the creation of diverse and large databases have enabled deep learning methods to reach state-of-the-2 Chapter 1. General introduction art results for tasks such as image restoration, compressive sensing, image recovery and classification, to name only a few. Together with the massive development of these approaches, concerns regarding their reliability and stability are growing. Artificial neural networks include, by nature, many degrees of freedom, which makes their study very intricate and their interpretation limited. Although variational and deep learning methods are based on distinct principles, they have complementary benefits and drawbacks, and therefore can greatly benefit from each other. Accordingly, in Chapter 4 we propose an interpretable neural network architecture inspired from a proximal interior point algorithm.

Similarly to machine learning, stochastic simulation methods estimate automatically the variables of interest. The Bayesian framework requires knowledge of prior distributions for the involved variables. Therefore, like in variational strategies, desirable properties for the solution can be promoted based on the choice of these priors. Stochastic simulation methods also offer convergence guarantees. Nonetheless, they may suffer from prohibitively long computational times. In Chapter 5, we propose to accelerate a stochastic sampling method by introducing a novel update strategy based on proximal tools.

We show in numerical experiments that the proposed methods compare favorably with state-of-the-art approaches on several challenging image processing applications. Hence, the outcomes of our work have both a mathematical and an applicative component. § 1.2 Main contributions Chapter 3 includes the following contributions.

(i) We combine the proximity operator with the logarithmic barrier method in a new proximal interior point algorithm called PIPA. One interesting feature of PIPA is its ability to handle variable metrics, which can be used to boost its convergence.

(ii) We also provide convergence results for this algorithm under a mild boundedness condition on the involved variable metrics. Furthermore, we show that the inner loop in PIPA has a linear convergence rate under some suitable assumptions.

(iii) Finally, the performance of the proposed algorithm is evaluated on two applications, namely hyperspectral unmixing, and joint geometrytexture decomposition and reconstruction of computed tomography data. PIPA compares favorably with state-of-the-art algorithms on these two examples.

Our main contributions from Chapter 4 are listed hereafter.

1.3. Collaborations

(i)

We propose a novel neural network architecture called iRestNet, which is obtained by unfolding a proximal interior point algorithm over a finite number of iterations. One key feature of this network is that prior knowledge can be directly incorporated into iRestNet thanks to a logarithmic barrier.

(ii) We derive expressions for the proximity operator of the barrier and for its first derivatives for three useful types of contraints. In particular, this allows our network to be trained using backpropagation and gradient descent.

(iii) We study the robustness of the proposed network with regards to an input perturbation, and derive conditions ensuring its stability in a useful case.

(iv) Numerical experiments performed on several non-blind image deblurring examples show that iRestNet compares favorably with respect to other variational and deep learning methods, including other neural networks based on deep unfolding. Finally, iRestNet benefits from a short execution time per image without any parameter search.

In Chapter 5, we consider a Bayesian framework and contribute to this area as follows.

(i) We propose a new sampling algorithm referred to as PP-ULA. This method is an original accelerated preconditioned version of the proximal unadjusted Langevin algorithm. PP-ULA takes advantage of a splitting strategy. It also includes variable metrics used for acceleration.

(ii) We present an efficient solver based on the majorize-minimize principle to tackle the involved non log-concave priors.

(iii) The proposed PP-ULA is then incorporated within a hybrid Gibbs sampler, yielding a substantial reduction of the computational time needed to perform joint high-quality deconvolution and segmentation of ultrasound images. § 1.3 Collaborations Some parts of this thesis have been done in collaboration with other researchers.

We have worked conjointly with Denis Kouamé and Jean-Yves Tourneret from the IRIT laboratory (UMR CNRS 5505), Toulouse, France, on a Bayesian approach for solving a challenging inverse problem in ultrasound imaging.

More precisely, the topic of this project was to accelerate a hybrid Gibbs sampler using proximal tools, with the aim of performing joint deconvolution and segmentation of ultrasound medical images. I spent two weeks in the TéSA laboratory in Toulouse as part of this collaboration. The results of this work are presented in Chapter 5.

We have also worked with Carla Bertocchi and Marco Prato from the university of Modena, Italy, on an optimization-inspired neural network architecture. Our goal was to apply deep unfolding to a proximal interior point algorithm. During this research collaboration, I spent one month in the university of Modena to work on this project. The corresponding results are presented in Chapter 4. § 1.4 Publications For articles with the * symbol, authors are listed in the alphabetical order, as it is customary in mathematical journals. This manuscript is organized as follows.

Submitted journal articles

In Chapter 2, we provide a general formulation for the study of inverse problems in Section 2.1.1 and describe briefly the three resolution methods mentioned previously, namely variational strategies (Section 2.1.2), Bayesian methods (Section 2.1.3) and deep learning techniques (Section 2.1.4). In particular, we introduce in Section 2.1.4.3 the concept of deep unfolding used in machine learning. Then, we introduce in Section 2.2.1 the main mathematical tools and notation that are used in this thesis. We end this chapter by presenting the optimization algorithms used in the proposed methods or as comparisons in numerical experiments (Section 2.2.2-2.2.5). A special attention is devoted to interior point methods in Section 2.2.5.

Chapter 3 is dedicated to the study of a new proximal interior point algorithm. For clarity, after describing the proposed method in Section 3.3, Chapter 1. General introduction we provide our main convergence results in Section 3.4 before detailing the corresponding proofs in Section 3.5. Finally, the results from two numerical experiments are presented in Sections 3.6 and 3.7. The first application is hyperspectral unmixing, while in the second application two tasks are jointly performed: geometry-texture image decomposition, and reconstruction of computed tomography data.

In Chapter 4, we investigate a novel neural network architecture designed by unfolding a proximal interior point algorithm. We start by describing in Section 4.2 the algorithm which is at the core of our method. Then, in Section 4.3 we derive the expressions of the proximity operator of the barrier and of its first derivatives for three types of constraints. These results are followed by a description of the proposed architecture called iRestNet (Section 4.4), whose robustness to an input perturbation is studied in Section 4.5. Finally, in Section 4.6 we evaluate iRestNet and compare it with state-of-the-art methods in numerical experiments related to non-blind image deblurring.

In Chapter 5, we consider a Bayesian setting and study a new proximal sampling strategy with applications in ultrasound imaging. We first describe the considered inverse problem in Section 5.1, which can be formulated as a joint deconvolution and segmentation of ultrasound images. Next, we detail in Section 5.3 the proposed preconditioned proximal unadjusted Langevin algorithm, referred to as PP-ULA. We then present the investigated hierarchical Bayesian model in Section 5.4, whereupon the proposed hybrid Gibbs sampler including PP-ULA is detailed in Section 5.5. Lastly, numerical experiments on both simulated and real ultrasound data are presented in Section 5.6.

Finally, we draw some conclusions and perspectives in Chapter 6. § 1.6 General notation Throughout this document, R, R + and R * denote the sets of real, positive real and non-zero real numbers, respectively. The set of vectors with m entries is noted R m , while the set of matrices with m rows and n columns is referred to as R m×n . We use the standard notations N and N * for the sets of positive and non-zero positive integers. Scalars will be denoted by lowercase letters, vectors will be denoted by lowercase bold letters, matrices will be denoted by uppercase bold letters and, finally, function will be denoted by lowercase or uppercase letters. For every vector v = (v i ) 1≤i≤m ∈ R m , v i denoted its ith coefficient. For every matrix M = (M i,j ) 1≤i≤m,1≤j≤n ∈ R m×n , M i,j denotes the element in the ith row and jth column of M .

For every matrix M ∈ R m×n , M , M -1 and |||M ||| denote its transpose, inverse and spectral norm, respectively. The square identity matrix in R m×m is referred to as I m , while the vector of R m with all entries equal to 1 is noted 1.6. General notation 1 m . The matrix in R m×n with all entries equal to 0 is noted 0 m×n , and the vector of R m with all entries equal to 0 is denoted as 0 m . For every s ∈ R, the absolute value of s is noted |s|. The Euclidean inner product is denoted as •, • and, for every v ∈ R m , v is the 2 norm of v.

-Chapter 2 -

Background

In many practical cases, the acquired data are an indirect measurement of the desired information. This can be due to physical constraints, for instance in astronomy, or to safety measures, as in medical diagnostics. Furthermore, the data are often degraded or incomplete, either because of the acquisition process or on account of restricted storage capacity. These limitations can induce a significant loss of information. In this context, estimating the sought signal or image constitutes a challenging inverse problem. Optimization plays a critical role in designing efficient algorithms and methods for solving problems of this type.

In this chapter, we provide a mathematical formulation for inverse problems in Section 2.1.1 and introduce three different resolution approaches that are investigated in this thesis: variational strategies (Section 2.1.2), stochastic simulations (Section 2.1.3) and deep learning (Section 2.1.4). Finally, after presenting some key mathematical tools in Section 2.2.1, we present in Sections 2.2.2-2.2.5 useful optimization algorithms that are used throughout this work. A special attention is devoted to interior point methods in Section 2.2.5. § 2.1 Inverse problems Inverse problems arise in a wide range of domains such as astronomy [START_REF] Lucy | Astronomical inverse problems[END_REF][START_REF] Bertero | Image deblurring with Poisson data: from cells to galaxies[END_REF], machine learning [START_REF] Vito | Learning from examples as an inverse problem[END_REF][START_REF] Jenatton | Proximal methods for sparse hierarchical dictionary learning[END_REF], geophysics [START_REF] Menke | Geophysical data analysis: Discrete inverse theory[END_REF], image processing, which includes for instance image restoration [START_REF] Benvenuto | The study of an iterative method for the reconstruction of images corrupted by Poisson and Gaussian noise[END_REF] and image superresolution [START_REF] Yang | Image super-resolution via sparse representation[END_REF]; and medical imaging [START_REF] Bertero | Inverse problems in biomedical imaging: modeling and methods of solution[END_REF], with applications like computed tomography (CT), magnetic resonance imaging, and fluorescence microscopy. 10 Chapter 2. Background

Formulation

In this thesis we consider the following formulation for inverse problems,

y = D(H(x)), (2.1)
where y ∈ R n gathers observations, x ∈ R m is the variable of interest, H : R m → R n is the observation operator, and D : R n → R n is the degradation operator.

The degradation operator represents a loss of information, which can come from the devices used to collect observations. For instance, optical sensors are generally subject to internal fluctuations, which are usually modelled as an additive white Gaussian noise independent from the signal of interest [START_REF] Oliveira | Adaptive total variation image deblurring: a majorization-minimization approach[END_REF]. There exist more realistic models where, for instance, the noise variance is a function of the sought signal [START_REF] Moser | Capacity results of an optical intensity channel with inputdependent Gaussian noise[END_REF]. In some applications such as medical ultrasound [START_REF] Sudha | Speckle noise reduction in ultrasound images by wavelet thresholding based on weighted variance[END_REF] or synthetic aperture radar [START_REF] López-Martínez | Polarimetric SAR speckle noise model[END_REF], measurements are corrupted by speckle noise, which is a multiplicative noise [START_REF] Durand | Multiplicative noise removal using 1 fidelity on frame coefficients[END_REF]. We can also mention impulsive noise [START_REF] Cai | Fast two-phase image deblurring under impulse noise[END_REF], Poisson noise [START_REF] Salmon | Poisson noise reduction with non-local PCA[END_REF], which occurs in photon-counting devices, and noises that are a mixture of the abovementioned noise models [START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF].

Regarding the observation operator, it accounts for the fact that the observation is an indirect measurement of the desired signal. In image deblurring for instance, H is frequently modelled as a linear operator representing the circular convolution of the image with a blur caused by optical defects [START_REF] Bertero | Image deblurring with Poisson data: from cells to galaxies[END_REF] or sensor motion [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF]. This operator is also linear in applications such as hyperspectral unmixing [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF] and compressive sensing [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF], but there exist nonlinear observation models too, e.g. [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF]. Throughout this thesis, we assume that H is known, i.e. we do not address blind inverse problems.

Solving problem (2.1) consists in finding a proper estimate x * of x given y. Assume that the observation operator is a linear operator, modelled by a matrix H ∈ R n×m , and that the noise is Gaussian. Then, an estimate of x can be produced by solving the following least-squares problem,

minimize x∈R m y -Hx 2 , (2.2)
where • is the Euclidean norm. However, (2.2) is generally ill-posed [START_REF] Hadamard | Sur les problèmes aux dérivées partielles et leur signification physique[END_REF], which means that, sometimes, there does not exist a solution, or when it exists, this solution is not necessarily unique, and it can be very sensitive to the presence of noise. Hence, there is a need for more performant approaches for solving inverse problems.

In the next sections, we introduce three different strategies that are investigated in this thesis, namely variational strategies, stochastic simulation methods and neural networks.

Variational strategy

Inverse problems can be efficiently solved by finding a minimizer of a given objective function. We refer to this approach as variational strategy.

Formulation

In order to find an appropriate solution to an inverse problem like (2.1), variational methods incorporate prior information on the sought signal through a penalization term [START_REF] Demoment | Image reconstruction and restoration: Overview of common estimation structures and problems[END_REF]. This leads to the following minimization problem, minimize

x∈C f (H(x), y) + λR(x), (2.3)
where

f : R n × R n →] -∞, +∞] is a data fidelity function, R : R m → ] -∞, +∞] is a regularization function, λ ∈ R + is a regularization param- eter and C is a subset of R m .
The data fidelity function is directly related to the degradation model D. For instance, when the noise is assumed to be Gaussian with variance

σ 2 , then f (•, y) = 1 2σ 2 y -• 2 .
The regularization term allows to enforce some desirable property in the solution like sparsity [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF][START_REF] Zhu | Seismic data denoising through multiscale and sparsity-promoting dictionary learning[END_REF], smoothness [START_REF] Poddar | Dynamic MRI using smoothness regularization on manifolds (SToRM)[END_REF] or grouping [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], while the feasible set C can embed problem-related constraints [START_REF] Musse | Topology preserving deformable image matching using constrained hierarchical parametric models[END_REF][START_REF] Klodt | A convex framework for image segmentation with moment constraints[END_REF]. Some classical choices for R are sparsity-promoting functions [START_REF] Bouman | A unified approach to statistical tomography using coordinate descent optimization[END_REF][START_REF] Bect | A 1 -unified variational framework for image restoration[END_REF][START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF][START_REF] Pustelnik | Waveletbased image deconvolution and reconstruction[END_REF] and the total variation and its various extensions [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF], which promote piecewise constant solutions. The regularization parameter λ is used to weight the prior information with respect to the data fidelity term. When problem (2.3) does not have a closed-form solution, it must be solved by an iterative solver such as the ones presented in Section 2.2.

Limitations

Although useful, variational approaches are sometimes limited by their complexity: solving (2.3) may require advanced algorithms that may be too slow for real-time applications. In addition, λ is a parameter that needs to be set and R is usually parametrized by one or several parameters, whose optimal choice may strongly depend on the data at hand. These parameters are often tuned manually or computed using, for instance, cross validation, the discrepancy principle [START_REF] Scherzer | The use of Morozov's discrepancy principle for Tikhonov regularization for solving nonlinear ill-posed problems[END_REF], or methods based on Stein unbiased risk estimates (SURE) [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF]. However, these methods are often time-consuming and their success is not always guaranteed. Furthermore, despite numerous efforts in designing sophisticated models, the solution to (2.3) could be further away from x than an intermediate iterate produced by a given algorithm used for solving (2.3). Such phenomenon justifies the development of early stopping methods, where the iterative procedure is stopped before convergence [START_REF] Yao | On early stopping in gradient descent learning[END_REF][START_REF] Rosasco | Learning with incremental iterative regularization[END_REF]]. Finding the optimal stopping time depends on the algorithm and requires the use of an oracle such as SURE, which may explain why these techniques are currently restricted to relatively simple objective functions.

In view of these limitations, an alternative to classical variational methods is the Bayesian approach. Accordingly, we present in the next section stochastic simulation methods, which can be used to solve inverse problems in a Bayesian setting.

Bayesian approaches

Bayesian theory leads to useful simulation methods that have demonstrated their efficiency for solving inverse problems such as image restoration [START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics[END_REF], classification [START_REF] Corander | Bayesian model learning based on a parallel MCMC strategy[END_REF], compressive sensing [Xu et al., 2014b] and source separation [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and MCMC sampling[END_REF]. In this section we consider a Bayesian setting, in the sense that y and x from Section 2.1.1 are assumed to be random variables [START_REF] Bernardo | Bayesian theory[END_REF]. After highlighting the links between Bayesian approaches and variational strategies, we are going to briefly introduce the simulation methods used in Chapter 5.

Links with the variational strategy

Although stochastic simulation approaches are quite different from variational methods, Bayesian theory provides an insightful reinterpretation of the data fidelity term and regularization function introduced in Section 2.1.2. Let p(y|x) be the likelihood of the observations and p(x) be the prior density associated to the variable of interest. Assume that p(y|x) ∝ exp(-f (x, y)) and that p(x) ∝ exp(-R(x)). Then, problem (2.3) It is worth noting that Bayesian approaches can be used to obtain other estimators than the MAP, such as the minimum mean squared error (MMSE) estimator [START_REF] Lesieur | MMSE of probabilistic lowrank matrix estimation: Universality with respect to the output channel[END_REF], i.e. the expectation of x with regards to the density p(x|y). In addition, the Bayesian setting can help derive useful uncertainty information about the estimate, for instance its variance.

In the next section, we introduce a class of stochastic methods called Markov chain Monte Carlo methods, which are widely used for solving inverse problems.

Markov chain Monte Carlo methods

Considering an inverse problem as defined in Section 2.1, the goal of Bayesian methods is to estimate x from the posterior probability distribution p(x|y). When this law is not easy to sample from, one can rely on Markov chain Monte Carlo (MCMC) methods [START_REF] Brooks | Handbook of Markov chain Monte Carlo[END_REF][START_REF] Robert | Monte Carlo statistical methods[END_REF][START_REF] Pereyra | A survey of stochastic simulation and optimization methods in signal processing[END_REF], which generate a Markov chain whose stationary distribution is the target law. In Chapter 5 we use two well-known MCMC methods, which are the Metropolis-Hastings random walk algorithm and the Gibbs sampler.

Metropolis-Hastings random walk algorithm

The Metropolis-Hastings (MH) random walk algorithm proposed in [START_REF] Roberts | Weak convergence and optimal scaling of random walk Metropolis algorithms[END_REF]] is a special case of the original MH sampling method from [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]] and [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. For every µ ∈ R m and σ ∈ R * + we denote by N (µ, σ 2 I m ) the Gaussian distribution with mean µ and variance σ 2 I m . In addition, for every (a, b) ∈ R 2 such that a ≤ b, we refer to the uniform distribution in [a, b] as U(a, b). Algorithm 1 draws sample from the target distribution p(x|y) by generating candidates from a proposal law that are accepted with a certain probability. These candidate samples are produced by adding to the current iterate a perturbation which is proportional to a realization of the standard normal distribution. Then, an acceptance test is performed such that, samples that are more probable to be produced by the target law than the current iterate are always accepted, and samples that are less probable states are accepted with a probability p( xk+1 |y)/p(x k |y) < 1.

The scale γ of the perturbation is a key element that drives the speed of convergence of Algorithm 1. In practice, the first iterations of the random walk are used to adjust γ based on some heuristics derived from the acceptance rate [START_REF] Gelman | Efficient Metropolis jumping rules[END_REF]. The iterates produced during this burn-in period are then discarded from the Markov chain. Therefore, the computation of the empirical MMSE estimate x MMSE from T samples (x (t 0 +t) ) 1≤t≤T produced after t 0 burn-in iterations reads

x MMSE = 1 T T t=1
x (t 0 +t) .

(2.6) Chapter 2. Background Algorithm 1: Metropolis-Hastings random walk [START_REF] Roberts | Weak convergence and optimal scaling of random walk Metropolis algorithms[END_REF] Let γ > 0.

Initialization:

x 0 ∈ R m . for k = 0, 1, . . . do xk+1 = x k + γω k with ω k ∼ N (0 m , I m ); t k ∼ U(0, 1); if t k < min 1, p( xk+1 |y) p(x k |y)
then Accept:

x k+1 = xk+1 ; else Reject: x k+1 = x k ; end end

Gibbs sampler

One can take advantage of the specific structure of a variable x ∈ R m that can be decomposed into q ∈ N blocks (x (i) ) 1≤i≤q , such that (∀i ∈ {1, . . . , q}) x (i) ∈ R m i and q i=1 m i = m. The Gibbs sampler [START_REF] Casella | Explaining the Gibbs sampler[END_REF], which is detailed in Algorithm 2, consists in sampling alternatively each block from its conditional distribution given observations y and the latest samples for the other blocks. This strategy is relevant when sampling from the individual conditional densities is simpler than drawing samples from the joint distribution. It can be shown that, after a long enough burnin period, the samples generated by the Gibbs sampler follow the target law p(x|y) [Robert and Casella, 2013, Section 10.2].

Algorithm 2: Gibbs sampler

Initialization: x 0 = (x (i) ) 1≤i≤q ∈ R m . for k = 0, 1, . . . do Generate x (1) k+1 ∼ p x (1) |y, x (2) k , . . . , x (q) k ; Generate x (2) k+1 ∼ p x (2) |y, x (1) k+1 , x (3) k , . . . , x (q) k ; . . . Generate x (q) k+1 ∼ p x (q) |y, x (1)
k+1 , . . . , x (q-1) k+1 ; end When Algorithm 2 is combined with other MCMC methods, it is called a hybrid Gibbs sampler. For instance, when a MH step is used to sample each block, then the convergence of the resulting hybrid Gibbs sampler still holds [Robert and Casella, 2013, Section 10.3].

Although complex inverse problems can be solved using the Bayesian framework, sampling methods can be of prohibitive computational cost. As mentioned previously, Markov chains require burn-in iterations before reaching the target distribution. Then, more iterations are needed to explore this distribution. Furthermore, the fact that the produced samples are not systematically accepted can slow down the process. These different factors can result in a very large computational cost for MCMC methods. Nonetheless, MCMC methods can sometimes be accelerated using specific techniques [START_REF] Orieux | Sampling high-dimensional Gaussian distributions for general linear inverse problems[END_REF][START_REF] Gilavert | Efficient Gaussian sampling for solving large-scale inverse problems using MCMC[END_REF][START_REF] Robert | Accelerating MCMC algorithms[END_REF]. We address this topic in Chapter 5, where we propose an acceleration method related to the proximal sampling scheme from [START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF].

In the next section, we introduce deep learning methods, which are investigated in Chapter 4.

Neural networks for solving inverse problems

As detailed in the review article [START_REF] Mccann | Convolutional neural networks for inverse problems in imaging: A review[END_REF], neural networks provide a popular and efficient alternative to variational methods for solving inverse problems. We clarify in Section 2.1.4.1 the vocabulary which is going to be used regarding neural networks. Then, we discuss in Section 2.1.4.2 some of their characteristics in the context of inverse problems. Lastly, we present in Section 2.1.4.3 the concept of deep-unfolding, which is at the core of Chapter 4.

Scope and vocabulary

In Chapter 4, we will focus on feedforward neural networks [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF]. The latters do not include any cycle, as opposed to recurrent neural networks. A neural network can be thought of as a function h(•, θ), where θ ∈ R p gathers parameters that are usually trained in a supervised manner [Goodfellow et al., 2016, Section 5.1]. This means that a training set of observations (y j ) 1≤j≤n with associated solutions (x j ) 1≤j≤n to problem (2.1) is available. In this context, parameters in θ are trained such that the outputs (h(y j , θ)) 1≤j≤n generated by the network on the training set minimize a given objective function. The latter is referred to as loss function, and is supposed to measure the distance between the solutions to the considered inverse problem and the outputs of the network. This training step may be accelerated by using batch normalization after each layer of the neural network [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF]. A validation set can also help monitoring the training without explicitely acting on the trainable parameters of the network [Goodfellow et al., 2016, Section 5.3]. This thesis does not address issues related to unsupervised training.

Although this training step can be both time and resource consuming, it is generally performed offline. Hence, as opposed to classical variational and MCMC methods, neural networks benefit from very competitive execution times. This has contributed to their considerable development, including for solving inverse problems.

For some inverse problems, it might be easier for the network to learn the residual y -D(H(x)) instead of the solution x itself, because pushing the residual to zero may appear easier than fitting an identity mapping by a stack of layers [START_REF] He | Deep residual learning for image recognition[END_REF]Zhang et al., 2017b]. This strategy is called residual learning and is implemented using skip connections between layers. This strategy may also be used to avoid the problem of vanishing gradient during training [START_REF] Tong | Image super-resolution using dense skip connections[END_REF].

Simple feedforward networks can be expressed as

R K-1 • (W K-1 • +b K-1 ) • • • • • R 0 • (W 0 • +b 0 ), (2.7)
where K ∈ N * is the number of layers, (R k ) 0≤k≤K-1 are nonlinear activation functions, (W k ) 0≤k≤K-1 are weight operators, also called linear layers, and (b k ) 0≤k≤K-1 are bias parameters. Among classical activation functions, one can mention the rectified linear unit (ReLU), x → max(0, x), and Sigmoid, x → 1/(1 + exp(-x)). In Chapter 4, we use two common linear layers that are the fully connected [START_REF] Pal | Multilayer perceptron, fuzzy sets, and classification[END_REF] The aforementioned limitations can be overcome by using deep unfolding, which combines neural networks and classical optimization algorithms. This approach is used in Chapter 4 and presented hereafter.

Deep unfolding

Consider a model-based method relying on an iterative solver A. At each iteration k ∈ N, this algorithm depends on several parameters which can be concatenated in a vector θ k ∈ R p . These hyperparameters are either derived from the model or related to the algorithm itself. Assuming that the convergence of A is proven, several iterates are produced as follows for k = 0, 1, . . .

x k+1 = A(x k , θ k ),
until an appropriate stopping criterion is satisfied. Unfolding or unrolling this method consists in setting the number of iterations to K ∈ N * , and treating the different hyperparameters (θ k ) 0≤k≤K-1 as learnable parameters of a neural network with K layers [START_REF] Hershey | Deep unfolding: Model-based inspiration of novel deep architectures[END_REF]. This leads to the following feed-forward architecture,

A •, L (θ)
K-1 (•) • • • • • A •, L (θ) 0 (•) ,
where every hidden layer L

k is used to infer θ k as follows, L

k (x k ) = θ k . Interestingly, the fact that this technique makes use of a limited number of layers can be viewed as an analogue of early stopping methods [START_REF] Rosasco | Learning with incremental iterative regularization[END_REF]. It is however worth mentioning that, in unfolded algorithms, the number of iterations (i.e., layers) is tuned during the off-line training step and is then fixed for all test images, which differs from early stopping strategies where the iteration number usually differs for each processed image.

Related works apply deep unfolding to probabilistic models, such as Markov random fields [START_REF] Hershey | Deep unfolding: Model-based inspiration of novel deep architectures[END_REF] or topic models [START_REF] Chien | Deep unfolding for topic models[END_REF], and to different algorithms like primal-dual solvers [START_REF] Wang | Proximal deep structured models[END_REF] or the proximal gradient method [START_REF] Mardani | Recurrent generative adversarial networks for proximal learning and automated compressive image recovery[END_REF][START_REF] Diamond | Unrolled optimization with deep priors[END_REF]. Classic optimization algorithms can be unfolded to perform many different tasks in image processing. For instance, FISTA and ISTA can be unfolded to perform sparse coding [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF][START_REF] Kamilov | Learning optimal nonlinearities for iterative thresholding algorithms[END_REF], while the same ISTA and ADMM can be unfolded for image compressive sensing [START_REF] Zhang | ISTA-Net: Interpretable optimization-inspired deep network for image compressive sensing[END_REF][START_REF] Sun | Deep ADMM-Net for compressive sensing MRI[END_REF]. However, in the aforementioned works, some functions and operators are learned, which weakens the link between the resulting network and the original algorithm. Deep unfolding is also used to learn shrinkage functions, which can be viewed as proximity operators of sparsity-promoting functions [Schmidt and

Mathematical analysis tools

We start by introducing some mathematical definitions.

Notation and definitions

The following definitions set up our framework for convex analysis and subdifferential calculus. (i) The domain of f is the set defined by

dom (f ) = {x ∈ R m | f (x) < +∞}.
(ii) The function f is proper if dom (f ) is not empty.

(iii) The function f is convex if for every α ∈]0, 1[ the following holds,

(∀(x, y) ∈ dom (f )) f (αx + (1 -α)y) ≤ αf (x) + (1 -α)f (y). (iv) The function f is coercive if lim x →+∞ f (x) = +∞.
(v) The function f is lower semicontinuous (l.s.c.) if, for every

x 0 ∈ R m , lim inf x→x 0 f (x) ≥ f (x 0 ).
(vi) The conjugate of f is the function f * : R m → [-∞, +∞] defined by

(∀x ∈ R m ) f * (x) = sup y∈R m ( x, y -f (y)) .
We denote by Γ 0 (R m ) the set of functions from R m to ] -∞, +∞] that are proper, l.s.c. and convex. We use the notion of subdifferential as it was introduced by Moreau [START_REF] Moreau | Fonctionnelles convexes[END_REF] and Rockafellar [START_REF] Rockafellar | Convex analysis[END_REF] in the 1960's, and whose definition is recalled below.

Definition 2.2.2 Let f : R m →] -∞, +∞].
The subdifferential of f is the set-valued operator ∂f defined by

(∀x ∈ R m ) ∂f (x) = {u ∈ R m | (∀y ∈ R m ) y -x, u + f (x) ≤ f (y)}.
The subdifferential can be linked to the set of minimizers of a function using Fermat's rule, as stated below.

3] Let f ∈ R m →] -∞, +∞] be proper. Then, Argmin x∈R m f (x) = {x ∈ R m | 0 m ∈ ∂f (x)}.
In the next section, we introduce a fundamental tool in nonconvex optimization.

Kurdyka-Lojasiewicz property

The Kurdyka-Lojasiewicz (KL) inequality was first introduced by Lojasiewicz for real analytic functions [START_REF] Lojasiewicz | Une propriété topologique des sous-ensembles analytiques réels[END_REF], and then extended by Kurdyka [START_REF] Kurdyka | On gradients of functions definable in o-minimal structures[END_REF] to differentiable functions that are definable in a o-minimal structure [van den [START_REF] Van Den Dries | The elementary theory of restricted analytic fields with exponentiation[END_REF]. Following the seminal work of Lojasiewicz and Kurdyka, the KL property has been extensively used for proving the convergence [START_REF] Bolte | The Lojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]Chouzenoux et al., 2014b] and obtaining convergence rates [Attouch and Bolte, 2009;[START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF] of optimization methods for possibly nondifferentiable functions, both in the convex and nonconvex case.

The KL inequality is based on the idea that a large number of functions can be re-parametrized such that they are sharp around their minimizers. This re-parametrization is based on the class of functions specified below. (iii) for every s ∈]0, η[, φ (s) > 0.

We introduce the following notation for the Euclidean distance between a vector v ∈ R m and a set C, dist(v, C) = inf w∈C vw . A general form of the KL inequality has been proposed in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and is recalled hereafter.

f (x) < f (x) < f (x) + η, the following inequality holds φ (f (x) -f (x)) dist(0 m , ∂f (x)) ≥ 1.
If f satisfies the KL property at each point of dom (∂f ), then f is called a KL function.

The KL property is satisfied in most practical optimization applications and, for a wide class of functions, the re-parametrization φ has a specific form which is detailed below and illustrated in Figure 2

(∀s ∈ [0, η[) φ(s) = cs 1-α ,
then f is said to have the KL property at x with exponent α. In other words, there exists c ∈ R * + such that, for every x ∈ Ω such that f (x) < f (x) < f (x) + η, the following inequality is satisfied,

dist(0 m , ∂f (x)) ≥ c(f (x) -f (x)) α .
If f has the same exponent α at any x ∈ dom (∂f ), then f is said to be a KL function with exponent α.

It is worth noting that this particular form of φ encompasses the original inequality introduced by Lojasiewicz in the 1960's. In addition, the KL exponent can be computed explicitly under some conditions, as detailed in [START_REF] Li | Calculus of the exponent of Kurdyka-Lojasiewicz inequality and its applications to linear convergence of first-order methods[END_REF]. This is of particular interest because this exponent can help to derive convergence rates for many optimization methods, as proven for instance in [Attouch and Bolte, 2009, Theorem 2] or [Attouch et al., 2010, Theorem 3.4].

In the next section we introduce the proximity operator, which is a very useful mathematical tool in optimization.

Proximity operators

The proximity operator has been first introduced in [START_REF] Moreau | Proximité et dualité dans un espace Hilbertien[END_REF] as an extension of the projection on a closed convex set. The proximity operator of a convex function f ∈ Γ 0 (R m ) at a point x ∈ R m is then defined as the unique minimizer of y → f (y) + 1 2 yx 2 , namely,

(∀x ∈ R m ) prox f (x) = argmin y∈R m f (y) + 1 2 y -x 2 .
The above definition can be extended to nonconvex functions [Attouch and Bolte, 2009]. Furthermore, it can be made more general by using a weighted norm instead of the Euclidean norm [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF]Combettes and Vũ, 2014;Chouzenoux et al., 2014b]. In this thesis, weighted norms will be noted as follows; for every

v ∈ R m , v M = v, M v 1 2 is the norm of v induced by M ∈ S +
m , with S + m the set of symmetric positivedefinite matrices in R m×m .

These considerations lead to the following general definition, which will be used in the rest of this document. Definition 2.2.7 Let f : R m →] -∞, +∞] be proper and l.s.c., and let M ∈ S + m and x ∈ R m . The proximity operator of f at x with regards to the norm induced by M is defined as the set

prox M f (x) = Argmin y∈R m f (y) + 1 2 y -x 2 M .
If the above set is reduced to a singleton, then, for simplicity, this unique element will also be noted prox M f (x). When the Euclidean norm is used, i.e. M = I m , the proximity operator will be noted prox f (x).

Throughout this thesis, M and its induced norm will be referred to as preconditioner and variable metric, respectively. It is worth noting that, if f ∈ Γ 0 (R m ), then the proximity operator of f is uniquely defined on R m since the objective function to be minimized is strongly convex.

Remark 2.2.8 Let f ∈ Γ 0 (R m ), x ∈ R m and M ∈ S + m .
The proximity operator of f at x in the norm induced by M can be re-written as

prox M f (x) = M -1 2 prox f •M -1 2 M 1 2 x .
This reformulation can be used to compute the proximity operator of a function in a variable metric.

In view of its widespread use in applications related to signal and image processing, the proximity operator has been extensively studied, in particular in the convex case [

) Let f ∈ Γ 0 (R m ), x ∈ R m and γ ∈ R * + . Then, x = prox γf (x) + γprox f * γ x γ . Proposition 2.2.10 (Nonexpansiveness) Let f ∈ Γ 0 (R m ) and M ∈ S + m . Then, prox M f is nonexpansive with respect to • M . Proposition 2.2.11 (Characterization) Let f ∈ Γ 0 (R m ), M ∈ S + m and x ∈ R m . Then, y = prox M f (x) ⇐⇒ x -y ∈ M -1 ∂f (y). Proposition 2.2.12 (Translation) Let f ∈ Γ 0 (R m ), x ∈ R m , y ∈ R m , γ ∈ R * + and µ ∈ R * . Set g = f (µ • -z). Then, prox γg (x) = 1 µ (y + prox γµ 2 f (µx -y)) Proposition 2.2.13 (Composition) Let f ∈ Γ 0 (R), x ∈ R m and y ∈ R m \ {0 m }. Set g = f ( •, y ). Then, prox g (x) = x + prox y 2 f ( x, y ) -x, y y 2 y.
An explicit or closed-form expression of the proximity operator has been derived for many useful functions, e.g. see1 for a collection of known proximity operators with associated codes. (i) If p = 1, we have the 1 -norm, which is widely used in optimization to promote sparsity [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF], and whose proximity operator is known as the soft-thresholding operator defined by

(∀γ ∈ R * + )(∀x ∈ R) prox γ|.| (x) = sign(x) max{|x| -γ}. (ii) If p ∈ { 4 3 , 3
prox γ|•| p (x) =                                            x + 4γ 3×2 1 3 (ξ -x) 1 3 -(ξ + x) 1 3 if p = 4 3 ; with ξ = x 2 + 256 729 γ 3 1 2 x + 9 8 γ 2 sign(x) 1 -1 + 16|x| 9γ 2 1 2 if p = 3 2 ; x 2γ+1 if p = 2; sign(x) (1+12γ|x|) 1 2 -1 6γ if p = 3; ξ+x 8γ 1 3 -ξ-x 8γ 1 3 if p = 4; with ξ = x 2 + 1 27γ 1 2
Example 2.2.15 Our second example is the logarithmic barrier, which is classically used in interior point methods (see Section 2.2.5). (i) The barrier associated with the feasible set [0, +∞[, also called the negative Burg entropy, is the function defined by

(∀x ∈ R) f (x) = -ln(x) if x ∈ R * + , +∞ else.
For every γ ∈ R * + we get the proximity operator of γf from [Bauschke and Combettes, 2017, Example 24.40] as follows,

(∀x ∈ R) prox γf (x) =
x + (x 2 + 4γ)

1 2
2 .

(ii) If we consider the compact feasible set [a, b], with a < b, then the corresponding logarithmic barrier g is defined as follows,

(∀x ∈ R) g(x) = -ln (x -a) -ln (b -x) if x ∈ ]a, b[ , +∞ else.
The proximity operator of g can be found in [START_REF] Chaux | A variational formulation for frame-based inverse problems[END_REF], Example 4.15]: for every γ ∈ R * + and every x ∈ R, prox γg (x) = κ, where κ is the unique solution in ]a, b[ to

κ 3 -(a + b + x) κ 2 + (ab + (a + b) x -2γ) κ -abx + (a + b) γ = 0.
The solution κ to the above cubic equation can be found using Cardano's method. As one can see in Figure 2.3, the smaller γ is, the closer the proximity operator is to the projection on [a, b].

In the next section, we introduce some properties related to operators. These notions are used in the convergence analysis presented in Chapter 3, and in the robustness study discussed in Section 4.5.

Operator properties

The Lipschitz continuity presented hereafter gives an upper bound for the rate of change of a given operator.

(∀(x, y) ∈ R m × R m ) T (x) -T (y) ≤ α x -y .
(ii) The operator T is nonexpansive if it is Lipschitz continuous with constant equal to 1.

(iii) A differentiable function f : R m →] -∞, +∞] is said to be Lipschitz differentiable with constant α ∈ R + if its gradient ∇f is α-Lipschitz continuous.
Next, we introduce the notion of averageness, which can be seen as a generalization of nonexpansiveness. Definition 2.2.17 The following proposition provides an upper bound of the effect of an input perturbation, which depends on the averageness constant. In particular, the smaller the constant is, the more stable the operator is.

(i) If T is averaged, then it is nonexpansive. (ii) The operator T is α-averaged if and only if, for every (x, y) ∈ R m ×R m , T (x) -T (y) 2 ≤ x -y 2 - 1 -α α (I m -T )(x) -(I m -T )(y) 2 .
Now that we have introduced important mathematical tools used in this thesis, we present in the next sections the main algorithms that we use either for designing the proposed methods or as comparisons.

Forward-backward algorithms

First, we consider the case when h is the zero function. Then, (2.8) reads

minimize x∈R m f (x) + g(x).
(2.9)

The well-known forward-backward (FB) algorithm [START_REF] Chen | Convergence rates in forwardbackward splitting[END_REF][START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF] addresses separately the two terms in the composite objective function, which leads to two steps for each iteration, namely a gradient (i.e. forward) step on the differentiable term g, and a proximal (i.e. backward) step on the nonsmooth term f .

Variable metric forward-backward algorithm

The FB algorithm can be accelerated by using a variable metric. This idea was already introduced in the original version of the FB algorithm in [START_REF] Chen | Convergence rates in forwardbackward splitting[END_REF], where the authors obtain the convergence of the algorithm in the convex setting assuming that the preconditioning matrices converge to a given operator [Chen and Rockafellar, 1997, Theorem For every k ∈ N, let

γ k ∈ R * + , λ k ∈ R * + , and M k ∈ S + m . Initialization: Take x 0 ∈ R m . for k = 0, 1, . . . do x k+1 = x k + λ k prox M k γ k f x k -γ k M -1 k ∇g(x k ) -x k ; end 2.2.

Generalized forward-backward algorithm

The FB algorithm can be generalized to the case when f can be written as

f = p i=1 f i . Then, problem (2.8) reads minimize x∈R m p i=1 f i (x) + g(x).
(2.10)

The generalized forward-backward splitting (GFBS) scheme in Algorithm 4 was proposed in [START_REF] Raguet | A generalized forward-backward splitting[END_REF] for solving (2.10) using the proximity operator of each f i . This algorithm will be used as a comparison in Section 3.6 with p = 2 and dom (f

1 ) = R m . The convergence of GFBS in this case is ensured by [Raguet et al., 2013, Theorem 2.1].
Algorithm 4: Generalized forward-backward [START_REF] Raguet | A generalized forward-backward splitting[END_REF] Let ω ∈]0, 1[ p be such that p i=1 ω i = 1. For every k ∈ N, let γ k ∈ R * + and λ k ∈ R * + . Initialization: For every i ∈ {1, . . . , p}, let z

(i) 0 ∈ R m . Set x 0 = p i=1 ω i z (i) 0 . for k = 0, 1, . . . do for i ∈ {1, . . . , p} do z (i) k+1 = z (i) k + λ k prox γ k ω i f i 2x k -z (i) k -γ k ∇g(x k ) -x k ; end x k+1 = p i=1 ω i z (i) k ; end 2.2.

Dual forward-backward algorithm

If f is the zero function and g = 1 2 • -y 2 for some y ∈ R m , then, (2.8) can be re-written as follows,

minimize x∈R m 1 2 x -y 2 + h(Lx), (2.11
Set ρ = |||L||| -2 and let ∈]0, min{1, ρ}[. For every k ∈ N, let γ k ∈ [ , 2ρ -] and λ k ∈ [ , 1]. Initialization: Let u 0 ∈ R n . for k = 0, 1, . . . do x k = y -L u k ; u k+1 = u k + λ k prox γ k h * (u k + γ k Lx k ) -u k ; end It is worth noting that problem (5) is equivalent to computing prox h(L•) (y).
Following this remark, in Section 3.6 and Chapter 5 we will use Algorithm 5 in combination with Remark 2.2.8 for computing the proximity operator of a function in a variable metric.

Primal-dual algorithms

In this section we focus on optimization algorithms that make use of both primal and dual updates.

Primal-dual splitting algorithm

In [Condat, 2013] and [START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF], the authors propose a primal-dual splitting (PDS) algorithm, which relies on a forward-backward scheme to split the differentiable term g and the nonsmooth function f , and which uses duality to decouple h and L. Algorithm 6 is a simplified version of the PDS algorithm where we assume that there is no error in the computation of the proximity operators and of the gradient. Note that, if g is the zero function, then we retrieve the algorithm proposed in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

The convergence result for the primal sequence in Algorithm 6 is given in [Condat, 2013, Theorem 3.1]. We use Algorithm 6 as a comparison method in Section 3.6.

Algorithm 6: Primal-dual splitting [Condat, 2013;[START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF]

(exact version) Let σ ∈ R * + , τ ∈ R * + , and (∀k ∈ N) let λ k ∈]0, 2[. Initialization: Let (x 0 , v 0 ) ∈ R m × R n . for k = 0, 1, . . . do y k = prox τ f x k -τ ∇g(x k ) + L v k ; u k = prox σh * (v k + σL(2y k -x k )); x k+1 = λ k y k + (1 -λ k )x k ; v k+1 = λ k u k + (1 -λ k )v k ; end 2.2.

Accelerated primal-dual splitting algoritm

In [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], the authors proposed an accelerated version of the previous algorithm when g is the zero function and f is strongly convex, i.e. there exists

β f ∈ R * + such that for every x ∈ dom (∂f ), (∀u ∈ ∂f (x))(∀y ∈ R m ) f (x) + u, y -x + β f 2 y -x 2 ≤ f (y).
Then, problem (2.8) reads

minimize x∈R m f (x) + h(Lx).
(2.12)

The accelerated primal-dual splitting (A-PDS) algorithm is given in Algorithm 7. The convergence of this algorithm is proven in [

Initialization: Let τ 0 ∈ R * + and σ 0 = (τ 0 |||L||| 2 ) -1 . Set (x 0 , y 0 ) ∈ R m × R n and x 0 = x 0 . for k = 0, 1, . . . do y k+1 = prox σ k h * (y k + σ k Lx k ); x k+1 = prox τ k f x k -τ k L y k+1 ; θ k = (1 + 2β f τ k ) -1 2 , τ k+1 = θ k τ k , σ k+1 = σ k /θ k ; x k+1 = x k+1 + θ k (x k+1 -x k ); end Since the function f = 1 2 • 2 is
strongly convex, in Section 3.7 we combine Algorithm 7 with Remark 2.2.8 to compute the proximity operator of a function in a variable metric. Chapter 2. Background

Alternating direction method of multipliers

Finally, we consider the following problem,

minimize x∈R m f (J x) + g(Kx) + h(Lx), (2.13)
with where n = p + q + r,

J ∈ R p×m , K ∈ R q×m and L ∈ R r×m . If J = K = I m ,
A =   J K L   , B = -I n ,
and where every y ∈ R n can be decomposed into a triplet (y (p) , y (q) , y

(r) ) ∈ R p × R q × R r , such that y =   y (p) y (q) y (r)   , and f (y) = f y (p) + g y (q) + h y (r) .
Problem (2.14) can be solved with the alternating direction method of multipliers (ADMM) [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF]

∈ R * + . Initialization: Let λ 0 ∈ R n and y 0 ∈ R n . for k = 0, 1, . . . do x k+1 = A A -1 A (λ k -By k ); y k+1 = prox f /µ (Ax k+1 -λ k ); λ k+1 = λ k -Ax k+1 -By k+1 ; end
Although the proximity operator of f /µ is separable with respect to y (p) , y (q) , and y (r) , the computation of the inverse of A A is a serious shortcoming of this algorithm. Nonetheless, reformulation (2.14) is not unique. For instance, another formulation is presented in [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF], where the inverse of A A is easier to compute, but where prox f /µ is not separable. In practice, alternating the minimization on the components of y may perform well [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF], but the convergence is not guaranteed in this case [START_REF] Chen | The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent[END_REF]. We use Algorithm 8 in Sections 3.6 and 3.7 as a comparison method.

When addressing (2.8) directly is too difficult, an interesting strategy consists in replacing the original problem by a sequence of intermediate problems that are easier to solve. This idea is at the core of the approaches presented in the next sections, namely the majorize-minimize principle and interior point methods.

Majorize-Minimize Principle

The majorize-minimize principle (MM) was first introduced in [START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF]. This principle is built on the notion of majorant, which is defined hereafter.

Definition 2.2.19 (Majorant) Let f : R m →] -∞, +∞]. The function h : R m × R m →] -∞, +∞] is said to be a majorant of f at y ∈ R m if h(y, y) = f (y) and (∀x ∈ dom (f )) h(x, y) ≥ f (x).
There exist different strategies for designing a majorant. For instance, if the function f is differentiable, then a quadratic majorant of f can be created as follows [START_REF] Böhning | Monotonicity of quadratic-approximation algorithms[END_REF]]

(∀x ∈ R m )(∀y ∈ R m ) h(x, y) = f (y) + ∇f (y) (x -y) + 1 2 (x -y) A(y)(x -y),
where A(y) ∈ S + m is such that the conditions in Definition 2.2.19 are satisfied. In addition, if

f is Lipschitz differentiable with constant L ∈ R + , then A(•) can be set equal to αI m with α ≥ L. More generally, A can be chosen such that A(y) -∇ 2 f (y) ∈ S +
m [START_REF] Hunter | A tutorial on MM algorithms[END_REF]. As shown in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF], this strategy can be used for designing efficient preconditioning matrices in the VMFB algorithm. We will follow this approach in Chapters 3 and 4 for choosing the involved variable metrics.

Finding a minimizer of the majorant can be easier than computing a solution to the original optimization problem. For instance, if the majorant is quadratic, then its minimizer has an explicit expression. Therefore, the MM principle consists in the following approach, instead of minimizing f directly, a sequence of iterates is produced where each iterate is a minimizer of a majorant of f taken at the previous iterate. This leads to the following Chapter 2. Background scheme, where for every iteration k ∈ N the update rule reads

x k+1 = argmin x∈dom (f ) h(x, x k ),
with h a majorant of f at x k . The convergence of this algorithm was established under some assumptions [START_REF] Jacobson | An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms[END_REF]. Furthermore, this approach has demonstrated good performance in various applications, such as image recovery [START_REF] Sotthivirat | Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms[END_REF][START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF], denoising [START_REF] Selesnick | Total variation denoising (an MM algorithm[END_REF] or quantile regression [START_REF] Hunter | Quantile regression via an MM algorithm[END_REF]. It is worth noting that this algorithm can be applied in a nonconvex setting [Chouzenoux et al., 2013a]. In Chapter 5, we derive an algorithm based on the MM principle to solve a nonconvex minimization problem.

In the next section, we present a specific class of methods used to address constrained optimization problems.

Interior point methods

Interior point methods (IPMs) encompass a wide range of methods which share the ability to handle constraints [Boyd and Vandenberghe, 2004, Chapter 11]. In this section we introduce the fundamental ideas and interior point methods that are used in Chapters 3 and 4.

IPMs became very popular from 1980 onwards due to the polynomialtime projection algorithm proposed in [START_REF] Karmarkar | A new polynomial-time algorithm for linear programming[END_REF], and which was later reinterpreted as an IPM in [START_REF] Gill | On projected Newton barrier methods for linear programming and an equivalence to KarmarkarâẠỎs projective method[END_REF]]. These approaches then demonstrated excellent performance in linear programming [START_REF] Wright | Interior methods for constrained optimization[END_REF]. They were extended to nonlinear problems [START_REF] Forsgren | Interior methods for nonlinear optimization[END_REF] thanks to the notion of self-concordance [START_REF] Nesterov | Interior-point polynomial algorithms in convex programming[END_REF]], which will be explained hereafter. IPMs have been applied to a wide variety of applications, including large-scale portfolio optimization [START_REF] Gondzio | Direct solution of linear systems of size 10 9 arising in optimization with interior point methods[END_REF], semidefinite programming [START_REF] Alizadeh | Interior point methods in semidefinite programming with applications to combinatorial optimization[END_REF] and image reconstruction [START_REF] Kim | An interiorpoint method for large-scale 1 -regularized least squares[END_REF]. Although IPMs are well-established methods, there are still new developments to be made regarding these approaches [START_REF] Gondzio | Interior point methods 25 years later[END_REF][START_REF] Armand | Uniform boundedness of the inverse of a Jacobian matrix arising in regularized interior-point methods[END_REF][START_REF] Ge | Interior-point methods strike back: Solving the Wasserstein barycenter problem[END_REF]].

Principle

Let us consider the following variational formulation,

minimize x∈C f (x), (2.15)
where the feasible set C is defined using p inequality constraints,

C = {x ∈ R m | (∀i ∈ {1, . . . , p}) c i (x) ≤ 0}.
(2.16)

Figure 2.4: Illustration of the merit function and influence of the barrier parameter when the barrier is chosen as the logarithmic one. In this example, the feasible set is taken as

C = {x ∈ R | 0.2 ≤ x ≤ 0.9}.
Barrier methods are standard IPMs in which the original problem (2.15) is replaced by a sequence of intermediate unconstrained problems (P µ k ) k∈N parametrized by a sequence of barrier coefficients (µ k ) k∈N . For every µ ∈ R * + , P µ is defined as

P µ : minimize x∈R m f (x) + µB(x), (2.17)
where B is a barrier function. The objective function in (2.17) is called the merit function and, for every

µ ∈ R * + , it is noted Ψ µ = f + µB.
The barrier can be seen as an approximation of the indicator function of the feasible set; its purpose is to ensure that the constraints involved in C are always strictly satisfied. As illustrated in Figure 2.4, it should be chosen such that the continuity properties of the constraints are preserved and such that it goes to infinity at the border of the feasible set [START_REF] Forsgren | Interior methods for nonlinear optimization[END_REF]. There exist various barrier functions in the literature [START_REF] Carroll | The created response surface technique for optimizing nonlinear, restrained systems[END_REF][START_REF] Kowalik | Nonlinear programming procedures and design optimization[END_REF], the most popular one being the logarithmic barrier, defined as follows

B : R m → ] -∞, +∞] x →      - p i=1 ln(-c i (x)) if (∀i ∈ {1, . . . , p}) c i (x) < 0, +∞ otherwise.
The logarithmic barrier has interesting features, in particular it is selfconcordant for every affine or quadratic constraints [Boyd and Vandenberghe, 2004, Chapter 9]. This means that B is three times continuously differentiable on its domain and that

(∀x ∈ dom (B)) |B (x)| ≤ 2B (x) 3/2 .
In

) ∈ R * + × R, the τ -lower level set of Ψ µ = f + µB, defined by lev ≤τ (Ψ µ ) = {x ∈ dom (B) | Ψ µ (x) ≤ τ }, is compact.
From now on, B will refer to the logarithmic barrier. The path-following method presented hereafter is a specific barrier method where each intermediate problem P µ k is solved before decreasing the barrier coefficient µ k .

Path-following method

If there exists a solution to every intermediate problem (2.17), then one can define a trajectory x * where, for every µ ∈ R * + , x * (µ) is a solution to P µ . The standard path-following interior point method detailed in Algorithm 9 follows this trajectory called central path. A general convergence result regarding this algorithm is given in Theorem 2.2.22.

Algorithm 9: Path-following barrier method

Let ρ > 1 and ∈ R * + . Initialization: Take x 0 ∈ dom (B) and µ 0 ∈ R * + . for k = 0, 1, . . . do Find x * (µ k ) a solution to P µ k using an iterative solver initialized with x k ;

x k+1 = x * (µ k ); Stop if pµ k ≤ ; µ k+1 = µ k /ρ; end Theorem 2.2.22 [Wright, 1992, Theorem 5 (iv)] 2 Under Assumption 2.2.20, if (x k ) k∈N
is a sequence generated by Algorithm 9, then there exists at least one convergent subsequence, and every of its cluster point is a solution to problem (2.15). We illustrate in Figure 2.5 the central path trajectory given by Algorithm 9 applied to the following linear programming problem,

minimize x∈R m c x s.t. Ax -b ∈] -∞, 0] p , (2.18)
where p = 5, m = 2,

c = 1 2 , A =       -1 -1 -1 1 0.7 1 3 -1 0.5 -1       , and b =       2 2 1 3 1      
. This problem will be used in Section 3.7 and serves as an illustration in the remaining of this section.

The Lagrangian theory provides an interesting reinterpretation of the path-following method.

Link with the Lagrangian theory

We consider problem (2.18) and make the following assumption. Assumption 2.2.23 Consider problem (2.18) and assume that there exists x ∈ R m such that Axb ∈] -∞, 0[ p . In addition, assume that for every

τ ∈ R, the sublevel set {x ∈ R m | Ax -b ∈] -∞, 0] p and c x ≤ τ } is bounded.
Under Assumption 2.2.23, Slater's constraint qualification is satisfied, and finding a solution x * to problem (2.18) is equivalent to finding a saddle point of the Lagrangian L, which is defined as follows,

(∀x ∈ R m )(∀λ ∈ [0, +∞[ p ) L(x, λ) = c x + λ (Ax -b).
Saddle points of L are characterized by the Karush-Kuhn-Tucker (KKT) conditions; for every solution x * to problem (2.18), there exists λ * ∈ R p such that the following system of equations is satisfied,

     c + A λ * = 0 p , (2.19) (∀i ∈ {1, . . . , p}) λ * i ≥ 0 and (Ax * -b) i ≤ 0, (2.20) (∀i ∈ {1, . . . , p}) λ * i (Ax * -b) i = 0. (2.21)
For every µ ∈ R * + , P µ can be rewritten as follows, minimize

(x,z)∈R m ×R p c x + µ B(z) s.t. z = Ax -b,
where

B(z) = -p i=1 ln(-z i ) if z ∈] -∞, 0[ p , +∞ otherwise.
Finding the solution x * (µ) to P µ is equivalent to identifying a saddle point of the associated Lagrangian L µ defined by

(∀x ∈ R m )(∀z ∈ R p )(∀λ ∈ R p ) L µ (x, z, λ) = c x + µ B(z) + λ (Ax -b -z).
There exist z * ∈ R p and λ * ∈ R p such that the following KKT conditions are satisfied

         c + A λ * = 0 p , (2.22) (∀i ∈ {1, . . . , p}) (Ax * (µ) -b) i < 0, (2.23) (∀i ∈ {1, . . . , p}) λ * i (Ax * (µ) -b) i = -µ, (2.24)
z * = Ax * (µ) -b. ( 2
i ∈ {1, . . . , p}, if (Ax * -b) i = 0 then λ * i > 0.
Remark 2.2.24 (Strict complementarity) In a more general context than that of linear programming, strict complementarity is not always necessary, see for instance [START_REF] Jittorntrum | Solution point differentiability without strict complementarity in nonlinear programming[END_REF], [Bonnans et al., 2006, Chapter 18] and the concept of linear monotone complementarity used in [Bonnans et al., 2006, Chapter 20] for quadratic programming. However, this assumption is often made for deriving convergence results regarding IPMs [START_REF] Wright | Interior methods for constrained optimization[END_REF].

Note that this assumption will also be used in Chapter 3.

Algorithm 9 assumes that one can easily compute a solution to each intermediate problem. We present hereafter the Newton barrier method, which uses Newton's method at each iteration to find a solution to P µ k . The Newton barrier method is used in particular in Section 3.7.

Newton barrier method

Algorithm 10: Newton barrier method

Let ρ > 1, ρ n > 1, ∈ R * + , n ∈ R * + and α ∈]0, 1 2 [. Initialization: Take x 0 ∈ dom (B) and µ 0 ∈ R * + . Set µ ← µ 0 . for k = 0, 1, . . . do ∆x k = -∇ 2 Ψ µ (x k ) -1 ∇Ψ µ (x k ) (Newton's direction); d k = -∇Ψ µ (x k ) ∆x k (Newton's decrement); if d k /2 ≤ n then Stop if pµ ≤ ; µ ← µ/ρ; else γ ← 1; while Ψ µ (x k + γ∆x k ) > Ψ µ (x k ) -αγd k do γ ← γ/ρ n (Backtracking); end x k+1 = x k + γ∆x k ; end end
The Newton barrier method is presented in Algorithm 10. In problem (2.18), the cost is linear and the constraints are affine, hence the merit function is self-concordant. This important property leads to the following result regarding the convergence rate of Algorithm 10.

(∀k ≥ k ) c x k -c x * ≤ ,
where k is defined by

k = ln(pµ 0 / ) ln(ρ) p(ρ -1 -ln(ρ)) η 1 + η 2 ,
with

η 1 = α(1 -2α) 2 ρ n (20 -8α) and η 2 = log 2 (log 2 (1/ n )).
It is worth noting that the Newton direction can also be derived from the modified KKT conditions by replacing λ * in (2.22) by its expression given in (2.24), and then applying Newton's method on the resulting equation. Following this remark, it is also possible to apply Newton's method directly on the modified KKT conditions without first eliminating λ * . This approach leads to primal-dual barrier methods with search directions on both the primal and dual variables. In particular, one can mention the infeasible primal-dual path-following algorithm presented in [START_REF] Gondzio | Interior point methods 25 years later[END_REF] which has polynomial complexity.

We introduce one last concept related to IPMs, which is the analytic center. Consider problem (2.15) and let J P be the set of indices i ∈ {1, . . . , p} such that there exists a solution to (2.15) for which the ith constraint is not active. In addition, let C and K be defined as follows,

Analytic center

C = {x ∈ R m | (∀i ∈ J P ) c i (x) < 0},
K = {x ∈ R m | (∀i ∈ {1, . . . , p} \ J P ) c i (x) = 0 and f (x) = f * },
where f * is the minimum value of f in C. Then, if it exists, the analytic center of the set C ∩ K is a solution to problem (2.15). When it exists, this specific solution is referred to as the analytic center of the solution set to (2.15). This insightful characterization will be used in Chapter 3.

In this section, we have presented classical IPMs which have proven their effectiveness, in particular for linear and quadratic programming. However, when dealing with more general optimization problems, these methods have some notable drawbacks. For instance, they require the inversion of a linear system at each iteration, and they are usually restricted to twice-differentiable objective functions. In addition, results about their pointwise convergence might be limited, as illustrated in Theorem 2.2.22. We address these issues in Chapter 3, where we propose to combine the barrier method with proximal tools. § 2.3 Summary In this chapter, we have introduced a mathematical formulation for inverse problems, as well as three resolution methods investigated in the next chapters. The notation and mathematical definitions used throughout this document have been presented in Section 2.2.1. As mentioned previously, the goal of this thesis is to illustrate the versatility of the proximity operator in the context of inverse problems, by combining it with each one of the three aforementioned approaches.

In Chapter 3, we propose a variational method based on a new algorithm combining the proximity operator with the barrier method, which has been presented in Section 2.2.5. Chapter 4 is dedicated to the study of an optimization-inspired neural network. The proposed architecture is obtained by applying the concept of deep unfolding presented in Section 2.1. 4 , 2011]. In this chapter, we propose to combine a proximal algorithm with the logarithmic barrier method, leading to a proximal interior point algorithm referred to as PIPA. One advantage of PIPA is that it can handle nonsmooth functions while ensuring that the constraints are always satisfied. In addition, PIPA includes a variable metric that can be used to boost its speed of convergence. This chapter is organized as follows. After describing the challenges involved and motivating our approach in Section 3.1, we introduce the variational formulation and our assumptions in Section 3.2. The proposed method is then presented in Section 3.3 and summarized in Algorithms 11 and 12.

Our main theoretical results are provided in Section 3.4, while the proofs for these results are given in Section 3.5. Numerical experiments are presented in Sections 3.6 and 3.7. § 3.1 Challenges and motivation Combining a proximal algorithm with the classical barrier method requires to address two main challenges that are discussed in this section: (i) handling the barrier term while maintaining acceptable computational complexity and speed, (ii) ensuring the convergence of the algorithm without the gradient-Lipschitz property. The first difficulty will be addressed thanks to a splitting strategy and a variable metric, while the second point will be tackled by including an appropriate line search in the proposed algorithm.

In most applications, the objective function is composite, in that it can be split into a nonsmooth term and a differentiable term. In particular, intermediate problems (2.17) arising in the barrier method may have this composite nature. As detailed in Section 2.2.2, some proximal resolution methods take advantage of this decomposition, like the VMFB algorithm, summarized in Algorithm 3, where a variable metric can be used for acceleration or simplification. The KL inequality presented in Section 2.2.1.2 proves to be a central argument for proving the convergence of this algorithm [Chouzenoux et Proximal and interior point methods can be combined to produce efficient solvers, as illustrated in [START_REF] Kaplan | Proximal methods in view of interior-point strategies[END_REF]], where the authors proposed an algorithm that minimizes a convex differentiable function over convex inequality constraints. However, this framework does not make use of any splitting strategy, hence it assumes that the proximity operator of the merit function is easy to compute. More recently, a preconditioning strategy based on the logarithmic barrier was proposed in [START_REF] Valkonen | Interior-proximal primal-dual methods[END_REF] to modify the dual update in a proximal primal-dual algorithm.

A major challenge, when dealing with IPMs is that the logarithmic barrier does not satisfy the gradient-Lipschitz property. Therefore, specific line search strategies, like the majorize-minimize line search [START_REF] Chouzenoux | Majorize-minimize linesearch for inversion methods involving barrier function optimization[END_REF], have to be designed in order to preserve the convergence properties of the methods used to solve the intermediate problems. This question has also been addressed in the context of VMFB algorithms, when the gradient of the smooth term in the objective function is not globally Lipschitz-continuous. Following the work of [START_REF] Bibliography | A coordinate gradient descent method for nonsmooth separable minimization[END_REF], Armijo-type line searches were proposed in [Bello Cruz and Nghia, 2016; [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF], where the convergence of the algorithm is obtained in a convex setting under suitable assumptions on the variable metrics. However, the line search strategy in [Bello Cruz and Nghia, 2016] requires multiple gradient computations, 3.2. Optimization problem and assumptions while [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF] requires the domain of the nondifferentiable function to be closed. Other line searches have also been studied in [START_REF] Bibliography | The variable metric forward-backward splitting algorithm under mild differentiability assumptions[END_REF], where convergence guarantees and convergence rate in terms of function values are provided for the convex case. It is worth noting that, similarly to [Bonettini and Prato, 2015b], [START_REF] Bibliography | The variable metric forward-backward splitting algorithm under mild differentiability assumptions[END_REF] considers a relaxed version of the monotonicity condition on the variable metrics from [Combettes and Vũ, 2014], where the metrics converge to a multiple of the identity operator, with a multiplicative factor which is allowed to vary along iterations. This assumption still remains restrictive compared to the conditions required in [Chouzenoux et al., 2014b], where the proof of convergence is carried out using the KL property.

In this chapter, we propose to combine the VMFB algorithm with the logarithmic barrier method, leading to our proposed algorithm PIPA. We will make use of one of the line searches investigated in [START_REF] Bibliography | The variable metric forward-backward splitting algorithm under mild differentiability assumptions[END_REF] to determine the stepsize value for each iteration of the proposed algorithm. Although our assumptions on the function domains are different from those in [START_REF] Bibliography | The variable metric forward-backward splitting algorithm under mild differentiability assumptions[END_REF], we prove that the line search remains valid in our context. Furthermore, we carry out the convergence analysis of the proposed algorithm under a mild boundedness condition on the involved variable metrics. Under some additional assumptions, we derive a linear convergence rate for the inner loop involved in PIPA. § 3.

Optimization problem and assumptions

In this chapter we consider the following constrained minimization problem,

P 0 : minimize x∈C f (x) + g(x) (3.1)
where the feasible set C is defined as follows,

C = {x ∈ R m | (∀i ∈ {1, . . . , p}) c i (x) ≤ 0}.
We also define the set D as

D = {x ∈ R m | (∀i ∈ {1, . . . , p}) c i (x) < 0},
and assume that Assumption 3. [START_REF] Musse | Topology preserving deformable image matching using constrained hierarchical parametric models[END_REF] to ensure that the estimated image deformation is injective and preserves the topology. Constraints can also serve to enforce some a priori knowledge about the solution, and act as regularization terms, as in the image segmentation approach in [START_REF] Klodt | A convex framework for image segmentation with moment constraints[END_REF], where bound constraints are imposed on the segmented areas and their barycenters.

Following the framework of IPMs, we propose to reformulate the constrained problem P 0 as a sequence (P µ j ) j∈N of modified subproblems parametrized for every j ∈ N by a barrier coefficient µ j > 0 and defined as

P µ j : minimize x∈R m f (x) + g(x) + µ j B(x) (3.2)
where B is the logarithmic barrier associated with the constraints:

B : R m → ] -∞, +∞] x →      - p i=1 ln(-c i (x)) if x ∈ D +∞ otherwise.
For simplicity, for every x ∈ R m we introduce the shorter notation c(x) = (c i (x)) 1≤i≤p ∈ R p , and the following functions, (∀µ > 0) ϕ µ = g + µB, and

Ψ µ = f + ϕ µ ,
where Ψ µ is designated as the merit function. § 3.3 PIPA algorithm

The proposed method, PIPA, is made of two interlocked loops. These are detailed in Algorithms 11 and 12, where denotes the Loewner partial order for matrices. Given j ∈ N, Algorithm 11 produces an approximate solution to P µ j via VMFB iterations consisting in a gradient step on the smooth term ϕ µ j , and a proximal step on the nondifferentiable term f . The proximity operators are computed within the metric induced by symmetric definite positive preconditioning matrices, whose eigenvalues are bounded from below Algorithm 11: A µ (x 0 , δ, θ, γ, , ν, ν)

Inputs are such that (δ, θ) ∈]0, 1[ 2 , (γ, µ) ∈ (R * + ) 2 , ∈ R + and 0 < ν ≤ ν; Initialization: x 0 ∈ D; for k = 0, 1, . . . do Choose M k ∈ S + m such that νI m M k νI m ; for l = 0, 1, . . . do xk,l = prox M k γθ l f x k -γθ l M -1 k ∇ϕ µ j (x k ) ; Exit loop if (3.3) is satisfied; end x k+1 = xk,l ; γ k = γθ l ; v k+1 = 1 γ k M k (x k -x k+1 ) -∇ϕ µ (x k ) + ∇ϕ µ (x k+1 ); Stop if v k+1 < ; end Return (x k+1 , v k+1 ); Algorithm 12: Proximal Interior Point Algorithm (PIPA) Let (δ, θ) ∈]0, 1[ 2 , (γ, µ 0 ) ∈ (R * + ) 2 ,
x 0 ∈ D, and let (µ j ) j∈N and ( j ) j∈N satisfy Assumption 3.3.1;

for j = 0, 1, . . . do Let 0 < ν j ≤ ν j ; (x j+1 , v j+1 ) = A µ j (x j , δ, θ, γ, j , ν j , ν j ); λ j+1 = - µ j c i (x j+1 ) 1≤i≤p ; end Return x j+1 ;
and from above (using ν and ν positive constants). For well-chosen matrices, this variable metric strategy can significantly improve the convergence speed. Preconditioning matrices can be determined for instance through a majorize-minimize procedure [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF].

It must be emphasized that, since the barrier is logarithmic, the gradient of ϕ µ j is not Lipschitz-continuous on R m . Thus, the VMFB algorithm must be associated with a line search to find an appropriate value for the stepsize so as to guarantee convergence of the scheme. In Algorithm 11, we use the backtracking line search method investigated in [Salzo,

ϕ µ ( xk,l ) -ϕ µ (x k ) -xk,l -x k , ∇ϕ µ (x k ) ≤ δ γθ l xk,l -x k 2 M k , (3.3)
for some l ∈ N. Such line search ensures both sufficient decrease of the criterion and feasibility of the next iterate. It is interesting to note that if the whole cost function in P 0 is smooth (f = 0), then it reduces to the standard Armijo line search along the steepest direction. Applications related to p -norms, with 1 < p < 2, and Bregman distances where considered in [START_REF] Bibliography | The variable metric forward-backward splitting algorithm under mild differentiability assumptions[END_REF], but, to the best of our knowledge, it is the first time that this line search is applied in the context of an interior point approach.

The resolution of the inner subproblem through Algorithm 11 is stopped once a certain accuracy is reached. In practice, we propose to stop the iterations once the norm of one element of the sequence (v k+1 ) k∈N , where (∀k ∈ N) v k+1 ∈ ∂Ψ µ (x k+1 ), is sufficiently small. As stated in Assumption 3.3.1 below, the sequence ( j ) j∈N and the barrier parameter (µ j ) j∈N , have to be chosen properly to secure the convergence of the sequence of iterates produced by Algorithm 12 to a solution to the initial problem P 0 . Hence, the decrease of the barrier parameter in Algorithm 12 is chosen such that Assumption 3.3.1 is satisfied ; examples of such decrease strategies can be found in Sections 3.6 and 3.7. Assumption 3.3.1 (Hyperparameters) For every j ∈ N, µ j > 0 and j > 0. In addition, lim j→+∞ µ j = 0 and lim j→+∞ j /µ j = 0. Remark 3.3.2 The sequence (λ j+1 ) j∈N , produced by Algorithm 12, can be linked to the Lagrangian parameters associated with the constrained problem P 0 . This is used in proof of convergence in Section 3.5.3.

Related works

It can be noted that there exist links between the proposed PIPA algorithm and a different class of methods called diagonal or penalization methods [Attouch et al., 2011b;[START_REF] Garrigos | Iterative regularization via dual diagonal descent[END_REF], for which a general study was recently provided for the continuous setting in [START_REF] Attouch | Asymptotic behavior of nonautonomous monotone and subgradient evolution equations[END_REF]. In [Attouch et al., 2011a], the authors proposed a similar approach in the discrete setting based on the FB algorithm. However, they assume that the gradient of the penalization function, which is the equivalent of our barrier, is Lipschitz continuous, and they let the penalization parameter tend to infinity, whereas our barrier parameter decreases to zero. In addition, their method does not solve a problem of the form (3.1), instead, it performs a hierarchical minimization [Attouch et al., 2011a, Eq. ( 6)]. This key difference is also highlighted in [START_REF] Alvarez | Asymptotic selection of viscosity equilibria of semilinear evolution equations by the introduction of a slowly vanishing term[END_REF]] and [START_REF] Cabot | Proximal point algorithm controlled by a slowly vanishing term: applications to hierarchical minimization[END_REF], in the continuous and discrete settings, respectively. In [START_REF] Cabot | Proximal point algorithm controlled by a slowly vanishing term: applications to hierarchical minimization[END_REF] the author studied the proximal point algorithm, without any forward-backward step, and considered a penalization parameter that, similarly to our barrier parameter, vanishes to zero. The optimization problem that is solved by this algorithm then depends on the rate of reduction chosen for this coefficient. In addition, it is assumed that the penalization function is bounded from below, which is not necessarily satisfied by the logarithmic barrier.

Related works also include Bregman distance approaches and entropylike proximal algorithms [START_REF] Iusem | Entropy-like proximal methods in convex programming[END_REF][START_REF] Brito | Interior proximal algorithm for quasiconvex programming problems and variational inequalities with linear constraints[END_REF][START_REF] Quiroz | An inexact proximal method for quasiconvex minimization[END_REF], where the Euclidean norm in the definition of the proximity operator is replaced by a divergence measure. The latter can be chosen such that feasibility is ensured at each iteration. However, the computation of the modified proximity operator in such methods is usually not straightforward.

In the next sections, we provide a convergence study of PIPA. § 3.

Main convergence results

In this section we state our main theoretical results.

Well-definedness of Algorithm 11

First, Theorem 3.4.1 guarantees that the stopping criterion in Algorithm 11 is well-defined.

Theorem 3.4.1 Under Assumption 3.2.1, for every (δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν] and x 0 ∈ D, for every ∈ R * + there exists k ∈ N such that v k+1 < .
Proof. See Section 3.5.1.4.

Convergence rate of Algorithm 11

An important feature of Algorithm 11 is the decay rate of the sequence ( v k+1 ) k∈N . This rate can be made explicit for the particular instance of linear inequality constraints. More precisely, we focus on the case when the constrained problem takes the form:

minimize x∈R m g(Hx) + f (x) subject to Ax -b ∈] -∞, 0] p , (3.4)
where g : (ii) g is l.s.c. with an open domain, it is strongly convex on any compact subset of dom ( g) and it is twice continuously differentiable on dom ( g).

R n →]-∞, +∞], H ∈ R n×m , A ∈ R p×m , b ∈ R p ,
(iii) f + g • H is proper and it is continuous on the domain of its subdifferential.

Remark 3.4.3 Assumption 3.4.2 holds for instance if g = 1 2 • -y 2 with y ∈ R n and if f = κ W • 1 where κ ≥ 0 and W ∈ R q×m is a linear transform (e.g., a wavelet analysis operator [START_REF] Pustelnik | Waveletbased image deconvolution and reconstruction[END_REF]Chaux et al., 2010]). This corresponds to an 1 -regularized least-squares problem, at the core of many applications such as denoising [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], image restoration [START_REF] Fu | Efficient minimization methods of mixed 2-1 and 1-1 norms for image restoration[END_REF], machine learning, or biological data analysis [Hastie et al., 2009, Chapter 18.4].

The following result can be deduced from existing results concerning the use of the KL inequality in optimization. Proof. See Section 3.5.2.2.

Convergence of Algorithm 12

We finally present results regarding the convergence of the proposed method PIPA, i.e. Algorithm 12.

Theorem 3.4.5 Suppose that Assumptions 3.2.1 and 3.3.1 hold for every (δ, θ) ∈]0, 1[ 2 , (γ, µ 0 ) ∈ (R * + ) 2 and x 0 ∈ D. Then, any sequence (x j+1 , λ j+1 ) j∈N generated by Algorithm 12 is bounded. In addition, every of its cluster point (x * , λ * ) is a primal-dual solution to P 0 , i.e. (x * , λ * ) is a saddle point for the Lagrangian defined in (3.35).

Proof. See Section 3.5.3.1.

A stronger convergence result can be obtained under additional assumptions. In particular, the following condition will turn out to play an important role. Assumption 3.4.6 Either the constraints are affine, i.e. c : x → Axb where A ∈ R p×m and b ∈ R p , and A has full column rank (i.e. A is injective), or there exists i ∈ {1, . . . , p} such that c i is strictly convex.

Let S P be the set of solutions to P 0 (primal solutions), and let S D be the set of solutions to the Lagrange dual problem associated with P 0 (dual solutions), whose definition is recalled below [Boyd and Vandenberghe, 2004, Section 5.2].

maximize

λ∈[0,+∞[ p inf x∈R m f (x) + g(x) + λ c(x)
In addition, let

J P = {i ∈ {1, . . . , p} | (∃x ∈ S P ) c i (x) < 0}, (3.5)
and (ii) Suppose that there exist at least two distinct elements in S P (resp. S D ), that Assumption 3.4.6 holds, and that P 0 has the strict complementarity property, i.e. J P ∪ J D = {1, . . . , p}. Then the sequence (x j+1 ) j∈N (resp. (λ j+1 ) j∈N ) produced by Algorithm 12 converges to a primal (resp. dual) solution to P 0 , which is the analytic center of S P (resp. S D ).

J D = {i ∈ {1, . . . , p} | (∃λ ∈ S D ) λ i > 0}. ( 3
Proof. See Section 3.5.3.3.

The next sections providing the proofs for the aforementioned theorems are organized as follows. First, we show in Section 3.5.1.2 that, under the considered assumptions, the chosen line search is well-defined. Then, we derive Lemma 3.5.6, Corollary 3.5.8 and Lemma 3.5.10 in Section 3.5.1.3, which lead to the proof of Theorem 3.4.1 in Section 3.5.1.4. Section 3.5.2 is dedicated to the convergence analysis and convergence rate of Algorithm 11. In Section 3.5.2.1 we start by deriving Lemma 3.5.12 which, together with Lemmas 3.5.6 and 3.5.10, ensures that the sufficient decrease, relative error and continuity conditions required in [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forwardbackward splitting, and regularized Gauss-Seidel methods[END_REF], Theorem 2.9] are satisfied. This leads to Proposition 3.5.13. We then derive Proposition 3.5.14 which directly leads to the proof of Theorem 3.4.4 in Section 3.5.2.2. Finally, we study the convergence of Algorithm 12 in Section 3.5.3 based on a Lagrangian approach. The proof for Theorem 3.4.5 is given in Section 3.5.3.1 and the proof of Theorem 3.4.7 is provided in 3.5.3.3. § 3.5 Proofs Proof. Let µ ∈ R * + . By assumption, there exists

x 0 ∈ D such that x 0 ∈ dom (Ψ µ ). The set lev ≤Ψµ(x 0 ) (Ψ µ ) = {x ∈ D | Ψ µ (x) ≤ Ψ µ (x 0 )} is nonempty since it includes x 0 , it is convex since D is convex and Ψ µ ∈ Γ 0 (R m
), and it is compact in view of Theorem 2.2.21. Solving P µ amounts to minimizing Ψ µ over lev ≤Ψµ(x 0 ) (Ψ µ ). Hence, the solution set to P µ is nonempty, convex, closed and bounded, as a subset of lev ≤Ψµ(x 0 ) (Ψ µ ) which is compact.

Proofs 51

Next, we derive a preliminary result about functions ϕ µ and Ψ µ with µ > 0.

Lemma 3.5.2 Under Assumption 3.2.1(iii), for every µ > 0, ∇ϕ µ is Lipschitz-continuous on every compact subset of D.

Proof. Let K be a compact subset of D. By assumption, for every i ∈ {1, . . . , p}, c i is differentiable on K, so it is continuous on K and, according to the extreme value theorem, it is bounded on K and it attains its bounds. Thus, there exist

(c, c) ∈] -∞, 0[ 2 such that (∀i ∈ {1, . . . , p})(∀x ∈ K) c ≤ c i (x) < c. Hence, for every (x, y) ∈ K 2 , ∇B(x) -∇B(y) ≤ p i=1 c i (y)∇c i (x) -c i (x)∇c i (y) |c i (x)c i (y)| ≤ p i=1 |c i (y)| c 2 ∇c i (x) -∇c i (y) + |c i (y) -c i (x)| c 2 ∇c i (y) . (3.7)
In addition, by assumption, for every i ∈ {1, . . . , p}, ∇c i is L i -Lipschitz continuous on K for some L i > 0; in particular, it is bounded by some constant K i > 0. Hence, for every i ∈ {1, . . . , p}, c i is K i -Lipschitz continuous on K and we deduce from (3.7) that

∇B(x) -∇B(y) ≤ p i=1 cL i + K 2 i c 2 x -y .
Therefore, for every µ > 0, ∇ϕ µ = ∇g + µ∇B is Lipschitz continuous on K.

Line search

We show in this section that the chosen line search is well-defined given our assumptions. Let (γ, µ)

∈ (R * + ) 2 , θ ∈]0, 1[ and M ∈ S + m .
Let h be such that

(∀x ∈ D)(∀l ∈ N) h(x, l) = prox M γθ l f x -γθ l M -1 ∇ϕ µ (x) .
Note that, from Proposition 2.2.11, for every x ∈ D and l ∈ N,

M (x -h(x, l)) -γθ l ∇ϕ µ (x) ∈ γθ l ∂f (h(x, l)). (3.8)
First, we check that, in the backtracking procedure, if the stepsize tends to zero then the expression for the next iterate in Algorithm 11 converges to the current iterate. Proof. Let l ∈ N and x ∈ D. From Proposition 2.2.10 we have

h(x, l) -prox M γθ l f (x) M ≤ γθ l M -1 ∇ϕ µ (x) M .
(3.10)

Taking the limit in (3.10) we deduce that,

lim l→+∞ h(x, l) -prox M γθ l f (x) M = 0 (3.11)
In addition, M -1 ∂f is a maximally monotone operator with respect to • M . From [Bauschke and Combettes, 2017, Proposition 16.27] and since

∅ = intdom (f ) ⊂ dom (∂f ), we have D ⊂ C ⊂ dom (f ) ⊂ dom (f ) = intdom (f ) ⊂ dom (∂f ).
Thus, [Bauschke and Combettes, 2017, Theorem 23.48] leads to

lim l→+∞ prox M γθ l f (x) -x M = 0. (3.12)
Finally, from (3.11) and (3.12) and the triangular inequality it follows that lim l→+∞ h(x, l)x M = 0, hence the result since M ∈ S + m .

We now show that, [Salzo, 2017, Lemma 3.6(ii)] holds for the line search.

Lemma 3.5.4 Under Assumption 3.2.1, for every θ ∈]0, 1[, (γ, µ)

∈ (R * + ) 2 and M ∈ S + m , if x ∈ D \ Argmin(f + ϕ µ ), then lim l→+∞ γθ l (ϕ µ (h(x, l)) -ϕ µ (x) -h(x, l) -x, ∇ϕ µ (x) ) h(x, l) -x 2 M = 0.
Proof. Take l ∈ N and x ∈ D \ Argmin(f + ϕ µ ). If x = h(x, l) then, from (3.8), it follows that 0 ∈ ∂f (x) + ∇ϕ µ (x) and x is a minimizer of f + ϕ µ , which leads to a contradiction. Hence, h(x, l)x M = 0. Since M -1 ∂f is a maximally monotone operator with respect to • M , [Huang and Dong, 2014, Lemma 1] leads to

x -h(x, l) M γθ l ≤ x -h(x, l + 1) M γθ l+1
3.5. Proofs 53 and (γθ l / h(x, l)x M ) l∈N is a decreasing sequence. Hence, there exists l 0 ∈ N and a ∈ R * + such that

(∀l > l 0 ) γθ l x -h(x, l) M ≤ a. (3.13)
In addition, from Lemma 3.5.3 we deduce that (3.9) holds. According to Assumption 3.2.1(iii), D is an open set, so there exist l 1 ≥ l 0 and a convex subset K of D such that x ∈ K and for every l ≥ l 1 , h(x, l) ∈ K. From Lemma 3.5.2 it follows that ∇ϕ µ is uniformly continuous on any compact subset of K. Thus, [Salzo, 2017, Corollary 3.4 (ii)] and the norm equivalence

λ min (M ) 1/2 • ≤ • M ≤ λ max (M ) 1/2
• , where λ min (M ) and λ max (M ) are the minimal and maximal eigen values of M , lead to Proof. Let us prove the result by induction. First, note that x 0 ∈ D. Assume that the property is satisfied at iteration k -1 if k ≥ 1. Under Assumption 3.2.1, if x k ∈ D is not a minimizer of f + ϕ µ then we can apply Lemma 3.5.4, which implies that the line search will stop: for any δ > 0, (3.3) will be satisfied for a finite l. If x k ∈ Argmin(f + ϕ µ ), then xk,0 = x k and the line search is satisfied for l = 0. Moreover, if x k+1 ∈ D then ϕ µ (x k+1 ) = +∞ and the inequality (3.3) could not be satisfied. Hence, x k+1 ∈ D.

lim l→+∞ ϕ µ (h(x, l)) -ϕ µ (x) -h(x, l) -x, ∇ϕ µ (x) h(x, l) -x M = 0. ( 3
Ψ µ (x k+1 ) ≤ Ψ µ (x k ) - ν(1 -δ) γ x k+1 -x k 2 . (3.15)
Proof. Let k ∈ N. The stopping criterion (3.3) for the backtracking procedure on γ k leads to

Ψ µ (x k+1 ) ≤ ϕ µ (x k ) + x k+1 -x k , ∇ϕ µ (x k ) + δ γ k x k+1 -x k 2 M k + f (x k+1 ). (3.16)
In addition, we have

M k (x k -x k+1 ) -γ k ∇ϕ µ (x k ) ∈ γ k ∂f (x k+1 ),
and it follows from the definition of the subdifferential that

γ k f (x k ) ≥ M k (x k -x k+1 ) -γ k ∇ϕ µ (x k ), x k -x k+1 + γ k f (x k+1 ). (3.17)
Re-writing (3.17) in a more convenient form yields

x k+1 -x k , ∇ϕ µ (x k ) + f (x k+1 ) ≤ f (x k ) - 1 γ k x k -x k+1 2 M k . (3.18)
Plugging (3.18) into (3.16) and using νI n M k completes the proof.

Remark 3.5.7 It is worth noting that, without the assumption of existence of bounds (ν, ν) on matrices (M k ) k∈N , the proof of Lemma 3.5.6 still allows us to conclude that the sequence (Ψ µ (x k )) k∈N is decreasing. Thus, in view of Theorem 2.2.21, there exists a compact K ⊂ D such that (∀k ∈ N) x k ∈ K. Therefore, by the continuity of function c on K we deduce that there exist

c and c in ] -∞, 0[ such that (∀k ∈ N)(∀i ∈ {1, . . . , p}) c ≤ c i (x k ) ≤ c < 0.
This remark will be useful in Sections 3.6 and 3.7 to prove that the chosen variable metrics satisfy the boundedness condition.

Before deriving a lowerbound for (γ k ) k∈N in Lemma 3.5.9, we show that the distance between two iterates produced by Algorithm 11 tends to zero and that the iterates are bounded. Corollary 3.5.8 Under Assumption 3.2.1, for every (δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν], and x 0 ∈ D, the sequence (x k ) k∈N produced by Algorithm 11 with = 0 satisfies the following properties:

(i) lim k→+∞ x k+1 -x k = 0; (ii) there exists a compact K ⊂ D such that (∀k ∈ N) x k ∈ K. Proof. (i) Summing (3.15) for k = 0 to N -1 ≥ 0 gives N -1 k=0 x k+1 -x k 2 ≤ γ ν(1 -δ) (Ψ µ (x 0 ) -Ψ µ (x N )) (3.19) ≤ γ ν(1 -δ) Ψ µ (x 0 ) -Ψ µ , (3.20)
where

Ψ µ = min x∈R m Ψ µ (x).
The existence of Ψ µ is ensured by Corollary 3.5.1. Letting N tend to infinity gives ∞ k=0 x k+1x k 2 < +∞, which leads directly to the result. (ii) From Theorem 2.2.21, the set {x ∈ D | Ψ µ (x) ≤ Ψ µ (x 0 )} is compact, and from Lemma 3.5.6, for every k ∈ N, x k belongs to this set, which completes the proof.

Before deriving relative error and continuity conditions, we show that, for every barrier problem, the stepsize computed with the line search is bounded below from zero. Lemma 3.5.9 Under Assumption 3.2.1, for every (δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν], and x 0 ∈ D, there exists γ µ > 0 such that the sequence (γ k ) k∈N generated by Algorithm 11 with = 0 is bounded below by γ µ .

Proof. Let I be the set of iterations in Algorithm 11 for which the stepsize value produced by the backtracking is strictly smaller than γ, i.e. I = {k ∈ N | γ k < γ}. For every k ∈ I there exists an integer l k > 0 such that γ k = γθ l k . By applying [Huang and Dong, 2014, Lemma 1], we have

(∀k ∈ I) x k -xk,l k -1 M k γθ l k -1 ≤ x k -xk,l k M k γθ l k ,
which leads to

(∀k ∈ I) x k -xk,l k -1 ≤ 1 θ ν ν 1 2 x k -x k+1 . (3.21)
From Corollary 3.5.8(ii), there exists a compact subset K of D such that, for every k ∈ N, On the other hand, since ∪ z∈K B(z, η/4) is a cover of K, it follows from the compactness of this latter set that there exists (z j ) 1≤j≤J in K such that On the other hand, according to Corollary 3.5.8(i), there exists k 0 ∈ N such that

x k ∈ K. Let ϑ : K → [0, +∞[ : x → dist(x, R m \ D).
K ⊂ J j=1 B(z j , η/4). ( 3
(∀k ≥ k 0 ) x k -x k+1 ≤ θ ν ν 1 2 η 4 . ( 3
(∀k ∈ I 0 ) xk,l k -1 -z j k < η 2 ,
which shows that xk,l k -1 ∈ S. Since ϕ µ is convex, the following inequality holds for every k ∈ I 0 :

xk,l k -1 -x k , ∇ϕ µ ( xk,l k -1 ) -∇ϕ µ (x k ) ≥ ϕ µ ( xk,l k -1 ) -ϕ µ (x k ) -xk,l k -1 -x k , ∇ϕ µ (x k ) . (3.26)
In addition, l k is the smallest integer such that (3.3) is satisfied. Hence, (3.3) is not satisfied for xk,l k -1 and, for every k ∈ I 0 , the following holds,

ϕ µ ( xk,l k -1 ) -ϕ µ (x k ) -xk,l k -1 -x k , ∇ϕ µ (x k ) > νθδ γ k xk,l k -1 -x k 2 .
(3.27) Necessarily, x k = xk,l k -1 . From (3.26) and (3.27), it follows that (∀k

∈ I 0 ), xk,l k -1 -x k ∇ϕ µ ( xk,l k -1 ) -∇ϕ µ (x k ) > νθδ γ k xk,l k -1 -x k 2 .
Moreover, according to Lemma 3.5.2, ∇ϕ µ is Lipschitz continuous on S.

Hence, there exists L S ∈ R * + such that

(∀k ∈ I 0 ) γ k > νθδ L S .

Proofs 57

In addition, (∀k ∈ I) γ k = γ, and the set I \ I 0 has a finite number of elements. Hence, the proof is complete by setting γ µ = min γ, νθδ L S , (γ k ) k∈I\I 0 .

We are now ready to identify a sequence of subgradients of Ψ µ converging to zero.

Lemma 3.5.10 (Relative error condition) Under Assumption 3.2.1, for every

(δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν]
, and x 0 ∈ D, the sequence (v k+1 ) k∈N produced by Algorithm 11 when = 0 is such that (∀k ∈ N) v k+1 ∈ ∂Ψ µ (x k+1 ) and

(∃K µ > 0)(∀k ∈ N) v k+1 ≤ K µ x k+1 -x k . Proof. Let k ∈ N. By definition of x k+1 and v k+1 , v k+1 -∇ϕ µ (x k+1 ) = 1 γ k M k (x k -x k+1 ) -∇ϕ µ (x k ) ∈ ∂f (x k+1 ). (3.28)
By definition of Ψ µ , it follows that v k+1 ∈ ∂Ψ µ (x k+1 ). In addition, the triangle inequality and Lemma 3.5.9 lead to

v k+1 ≤ ν γ µ x k -x k+1 + ∇ϕ µ (x k+1 ) -∇ϕ µ (x k ) .
From Corollary 3.5.8(ii) and Lemma 3.5.2, we know that (x k ) k∈N belongs to a compact subset K of D, on which ∇ϕ µ is Lipschitz-continuous for some constant L K > 0. Setting K µ = ν/γ µ + L K completes the proof.

Remark 3.5.11 It can also be deduced from (3.28) that, for every k ∈ N, In view of Lemma 3.5.10 and Corollary 3.5.8(i), for every ∈ R * + , there exists k ∈ N such that v k+1 < .

v k+1 -µ∇B(x k+1 ) ∈ ∂(f + g)(x k+1 ).

Convergence analysis of Algorithm 11

Preliminary results

We first derive the following continuity condition. (3.29)

Proof. From Corollary 3.5.8(ii) we know that (x k ) k∈N belongs to a compact subset K of D. Hence, there exists a subsequence (x kq ) q∈N converging to an element x ∈ D. Since Ψ µ is l.s.c., we have

Ψ µ (x) ≤ lim inf q→+∞ Ψ µ (x kq ). (3.30)
Without loss of generality one can assume that k 0 > 0. From Lemma 3.5.10, for every q ∈ N, v kq belongs to ∂Ψ µ (x kq ) and v kq ≤ K µ x kqx kq-1 .

For every q ∈ N, we have

Ψ µ (x kq ) ≤ -v kq , x -x kq + Ψ µ (x) ≤ K µ x kq -x kq-1 x -x kq + Ψ µ (x). (3.31)
From Corollary 3.5.8(i), x kqx kq-1 → 0 as q → +∞. Hence, taking the limit in (3.31) yields the following inequality The next result guarantees that, in the absence of stopping rule, Algorithm 11 converges to a solution to the barrier problem.

lim sup q→+∞ Ψ µ (x kq ) ≤ Ψ µ (x). ( 3
Proposition 3.5.13 Let (δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν]
, and x 0 ∈ D. Suppose that Ψ µ is a KL function and that Assumption 3.2.1 holds. Then the sequence (x k ) k∈N produced by Algorithm 11 with = 0 converges to a solution to P µ and has a finite length, i.e.

+∞ k=0

x k+1x k < +∞. The linear convergence properties follow from the fact that exp(-cK µ ) < 1 and exp(-cK µ /2) < 1.

Proof of Theorem 3.4.4

Proposition 3.5.13 ensures the convergence of (x k ) k∈N to some x * . According to Lemma 3.5.10, there exists K µ > 0 such that

(∀k ∈ N) v k+1 ≤ K µ x k+1 -x k ≤ K µ ( x k+1 -x * + x * -x k ).
We then deduce from Proposition 3.5.14 that there exists c ∈ R * + such that v k+1 = O exp(-cK µ /2) k .

Convergence analysis of Algorithm 12

We are now ready to establish the convergence of Algorithm 12 to a solution to problem P 0 described in (3.1). Under Assumption 3.2.1, finding a solution to P 0 is equivalent to finding a saddle point of the associated Lagrangian L 0 , which is defined, for every x ∈ R m and λ ∈ [0, +∞[ p , as

L 0 (x, λ) = f (x) + g(x) + λ c(x).
(3.35)

For every x ∈ C, let ∇c(x) be the matrix in R m×p such that (∀i ∈ {1, . . . , p}) its ith column is equal to ∇c i (x). For every v ∈ R p let Diag(v) denotes the diagonal matrix whose elements are given by v. A point (x * , λ * ) is a saddle point of L 0 if and only if the following optimality conditions are satisfied.

     0 ∈ ∂f (x * ) + ∇g(x * ) + ∇c(x * )λ * (3.36) Diag(λ * )c(x * ) = 0 p (3.37) λ * ∈ [0, +∞[ p and c(x * ) ∈] -∞, 0] p (3.38)
We prove in the following section that the sequences (x j ) j∈N and (λ j+1 ) j∈N , produced by Algorithm 12, converge to a saddle point of L 0 . We remind that, for every j ∈ N, the dual variable λ j+1 is defined as follows,

λ j+1 = - µ j c i (x j+1 ) 1≤i≤p . ( 3

.39)

As it is detailed in Section 2.2.5.2, (3.39) can be seen as a perturbation of condition (3.37) and, in the case of affine inequality constraints, (3.39) can be directly derived from the Lagrangian formulation of the barrier problem. We first show that the primal and dual sequences produced by Algorithm 12 are bounded. Lemma 3.5.15 Under Assumptions 3.2.1 and 3.3.1, for every (δ, θ) ∈]0, 1[ 2 , (γ, µ 0 ) ∈ (R * + ) 2 , and x 0 ∈ D, the sequences (x j+1 ) j∈N and (λ j+1 ) j∈N produced by Algorithm 12 are bounded.

Proof. Let j ∈ N. According to Lemma 3.5.10, v j+1 belongs to ∂Ψ µ j (x j+1 ). Using the definition of the subdifferential, the stopping criterion v j+1 ≤ j , and the Cauchy-Schwarz inequality leads to

Ψ µ j (x j+1 ) ≤ v j+1 , x j+1 -x 0 + Ψ µ j (x 0 ) ≤ j x j+1 -x 0 + Ψ µ j (x 0 ). (3.40)
Eq. (3.40) is equivalent to

(f + g)(x j+1 ) ≤ µ j (B(x 0 ) -B(x j+1 )) + j x j+1 -x 0 + (f + g)(x 0 ). (3.41)

Proofs 61

Moreover, since B is convex,

(f + g)(x j+1 ) ≤ -µ j ∇B(x 0 ), x j+1 -x 0 + j x j+1 -x 0 + (f + g)(x 0 ) ≤ ( j + µ j ∇B(x 0 ) ) x j+1 -x 0 + (f + g)(x 0 ). (3.42)
Assume that (x j+1 ) j∈N is unbounded. Then there exists a subsequence (d q ) q∈N = (x jq+1 ) q∈N of (x j+1 ) j∈N such that the sequence (t q ) q∈N = ( x jq+1 ) q∈N has only strictly positive elements and satisfies

lim q→+∞ t q = +∞, lim q→+∞ d q t q = d ∈ R m and d = 1. (3.43)
The last two equalities are derived from the compactness of the unit ball. Let x * be a solution to P 0 and let t ∈ R * + . For every q 0 ∈ N, let τ q 0 = min{t q | q ≥ q 0 }. Since (∀i ∈ {1, . . . , p}) c i (x * ) ≤ 0, (∀q ∈ N) c i (d q ) < 0, and c i is convex, we have

(∀q ≥ q 0 ) 1 - τ q 0 t q x * + τ q 0 t q d q ∈ D.
By taking the limit in the above inclusion as q → ∞, we obtain x * +τ q 0 d ∈ C. In addition, for every q ≥ q 0 , (f + g) 1 -τ q 0 t q x * + τ q 0 t q d q ≤ 1 -τ q 0 t q (f + g)(x * ) + τ q 0 t q (f + g)(d q ).

(3.44) We deduce from (3.42) that

1 t q (f + g)(d q ) ≤ ( jq + µ jq ∇B(x 0 ) ) d q t q - x 0 t q + 1 t q (f + g)(x 0 ). (3.45)
As q → +∞, we have jq → 0, µ jq → 0, t q → +∞ and d q /t qx 0 /t q → 1. Hence, taking the limit in (3.45) leads to lim sup q→∞ (f + g)(d q )/t q ≤ 0.

Using now the lower-semicontinuity of f + g and letting q tend to +∞ in (3.44) lead to

(∀q 0 ∈ N) (f + g) x * + τ q 0 d ≤ lim inf q→+∞ (f + g) 1 - τ q 0 t q x * + τ q 0 t q d q ≤ (f + g)(x * ) + lim sup q→+∞ τ q 0 t q (f + g)(d q ) ≤ (f + g)(x * ).
Therefore, for every q 0 ∈ N, x * + τ q 0 d is a solution to P 0 . Since τ q 0 → +∞ as q 0 → +∞, the set of solution to P 0 is unbounded. This is however in contradiction with Assumption 3.2.1(i), thus showing that (x j+1 ) j∈N is bounded.

A fast proximal IPM for constrained variational formulations

Similarly, we prove that (λ j ) j∈N is bounded. Let j ∈ N and i ∈ {1, . . . , p}. Since c i is convex, c i (x j+1 ) < 0 and µ j > 0, the following inequality holds

µ j c i (x 0 ) c i (x j+1 ) ≤ µ j + µ j c i (x j+1 ) ∇c i (x j+1 ), x 0 -x j+1 . (3.46)
Summing (3.46) for all i ∈ {1, . . . , p} leads to

-c(x 0 ), λ j+1 ≤ µ j p -µ j ∇B(x j+1 ), x 0 -x j+1 . (3.47)
In addition, from Remark 3.5.11 and the definition of the subdifferential of f + g, we deduce that

(f + g)(x j+1 ) ≤ -v j+1 -µ j ∇B(x j+1 ), x 0 -x j+1 + (f + g)(x 0 ). (3.48)
Combining (3.47) and (3.48) yields

-c(x 0 ), λ j+1 ≤ µ j p + (f + g)(x 0 ) -(f + g)(x j+1 ) + j x 0 -x j+1 .
Moreover, every component of λ j+1 and of -c(x 0 ) is strictly positive, hence

0 < -c(x 0 ), λ j+1 ≤ µ j p + (f + g)(x 0 ) -(f + g)(x * ) + j x 0 -x j+1 ,
where x * is a solution to P 0 . Since (µ j , j ) → (0, 0) as j → +∞, and since (x j+1 ) j∈N has been shown to be bounded, we conclude thatc(x 0 ), λ j+1 is bounded and so is (λ j+1 ) j∈N .

Proof of Theorem 3.4.5

According to Lemma 3.5.15, the sequences (x j+1 ) j∈N and (λ j+1 ) j∈N are bounded. Hence, there exists a subsequence (x jq+1 , λ jq+1 ) q∈N converging to some point (x * , λ * ). By construction, for every q ∈ N, c(x jq+1 ) ∈] -∞, 0[ p and λ jq+1 ∈]0, +∞[ p . Since (∀i ∈ {1, . . . , p}) c i is l.s.c., taking the limit as q → +∞ yields (3.38). Moreover, by definition, for every q ∈ N we have Diag(λ jq+1 )c(x jq+1 ) = -µ jq 1 p . Since (∀i ∈ {1, . . . , p}) c i is continuous on C, letting q → +∞ in the previous equality leads to (3.37). Let q ∈ N. From Remark 3.5.11, it follows that x jq+1 = prox f +g (x jq+1 + v jq+1µ jq ∇B(x jq+1 )). In addition, prox f +g is nonexpansive. Hence,

x jq+1 -prox f +g (x * -∇c(x * )λ * ) ≤ x jq+1 + v jq+1 -µ jq ∇B(x jq+1 ) -x * + ∇c(x * )λ * . (3.49)

Proofs 63

By applying the triangle inequality, we deduce from (3.49) that

x * -prox f +g (x * -∇c(x * )λ * ) ≤ 2 x * -x jq+1 + j + p i=1 (λ jq+1 ) i ∇c i (x jq+1 ) -λ * i ∇c i (x * ) . (3.50)
The sequence (x j+1 ) j∈N is a bounded sequence in D. Therefore, all its cluster points belong to a compact subset of C. In view of Assumption 3.2.1(iii), ∇c i is continuous at x * . Thus, taking the limit in (3.50) as q → +∞ leads to

x * = prox f +g (x * -∇c(x * )λ * ) , (3.51)
which is equivalent to (3.36). Finally, (x * , λ * ) is a saddle-point for the Lagrangian (3.35), which completes the proof.

Analytic center

When there are several primal or dual solutions to the constrained problem (3.1), assumptions are needed to prove the uniqueness of the cluster point exhibited in Theorem 3.4.5. Under these assumptions, the analytic center introduced in Section 2.2.5.4 provides an insightful characterization for the limit point. We show in this section the existence of a primal and a dual analytic centers.

Let us recall that under Assumption 3.2.1, strong duality holds and the set of saddle points for the Lagrangian L 0 is equal to S P × S D where S P and S D are the so-called sets of primal and dual solutions to P 0 [Boyd and Vandenberghe, 2004, Section 5.4.2]. We derive the following result for these sets.

Proposition 3.5.16 Under Assumptions 3.2.1(i)-(ii), S P and S D are nonempty bounded convex sets.

Proof. The results for S P directly follows from Assumptions 3.2.1(i)-(ii). The convexity of S D follows from standard results [Boyd and Vandenberghe, 2004, Section 5.2]. The fact that S D is not empty also follows from Assumption 3.2.1(i)-(ii). For every (x * , λ * ) ∈ S P × S D , the inequality (f + g)(x * ) ≤ (f + g)(x 0 ) + λ * c(x 0 ) holds for every x 0 ∈ D. If S D is unbounded, then the right-hand side can tend to -∞ which is in contradiction with (f + g)(x * ) being a finite number. So S D is bounded.

We also show that under an additional assumption, the sets J P , defined in (3.5), and J D , defined in (3.6), are nonempty. The number of elements in J P and J D are noted #J P and #J D , respectively. Proposition 3.5.17 If Assumption 3.4.6 holds and if the set S P (resp. S D ) contains at least two distinct elements, then the set J P (resp. J D ) is nonempty.

Proof. Assume that there are at least two distinct elements x 1 and x 2 in S P . Then, for every i ∈ {1, . . . , p}, c i (x 1 ) ≤ 0 and c i (x 2 ) ≤ 0. Assume that Assumption 3.4.6 holds. If the constraints are affine, i.e. c : x → Axb, with A an injective matrix, then there exists i 0 ∈ {1, . . . , p} such that (Ax 1 ) i 0 = (Ax 2 ) i 0 and c i 0 (x 1 ) = c i 0 (x 2 ). The same conclusion obviously holds if i 0 ∈ {1, . . . , p} is such that c i 0 is strictly convex. Hence, we have either c i 0 (x 1 ) < 0 or c i 0 (x 2 ) < 0, that is i 0 ∈ J P . Assume that there are at least two distinct elements λ 1 and λ 2 in S D . Both λ 1 and λ 2 belong to [0, +∞[ p . Since λ 1 and λ 2 are distinct there exists i 0 ∈ {1, . . . , p} such that (λ 1 ) (i 0 ) = (λ 2 ) (i 0 ) . Hence, we have either

(λ 1 ) (i 0 ) > 0 or (λ 2 ) (i 0 ) > 0, that is i 0 ∈ J D .
If J P is nonempty, we define the following quantities: for every x ∈ R m , c J P (x) = (c i (x)) i∈J P and B J P (x) = -i∈J P ln(-c i (x)) if c J P (x) ∈] -∞, 0[ #J P , +∞ otherwise. We also consider the following problem.

P P : minimize x∈R n B J P (x)
subject to x ∈ S P and c J P (x) ∈] -∞, 0[ #J P .

(3.52)

Similarly, if J D is nonempty, we define the following quantities: for every

λ ∈ R p , λ J D = (λ i ) i∈J D and b J D (λ) = -i∈J D ln(λ i ) if λ J D ∈]0, +∞[ #J D ,
+∞ otherwise. We also consider the following problem.

P D : minimize λ∈R p b J D (λ) subject to λ ∈ S D and λ J D ∈]0, +∞[ #J D (3.53)
Lemma 3.5.18 Under Assumptions 3.2.1 and 3.4.6, if S P (resp. S D ) does not reduce to a singleton, then there exists a unique solution to P P (resp. P D ) called the analytic center of S P (resp. S D ).

Proof. Assume that S P does not reduce to a singleton. According to Proposition 3.5.16, S P is nonempty and it thus contains at least two distinct elements. It then follows from Proposition 3.5.17 that the set J P is nonempty.

In addition, by invoking again Proposition 3.5.16, S P is bounded and convex. Since the constraint functions (c i ) 1≤i≤p are convex, we deduce that the 3.5. Proofs 65 feasible set of P P is bounded and convex. Because of the convexity of the set S P and the functions (c i ) 1≤i≤p , it can be checked that this feasible set is nonempty. In addition, under Assumptions 3.2.1 and 3.4.6, B J P is l.s.c. and strictly convex, and it is finite-valued on {x ∈ R m | c J P (x) ∈] -∞, 0[ #J P }.

Hence, there exists a unique solution to P P . Assume that S D does not reduce to a singleton. It then follows from Propositions 3.5.16 and 3.5.17 that J D = ∅. By using a similar reasoning as for P P we deduce that there exists a unique solution to P D .

The complementary slackness property in (3.37) ensures that J P ∩ J D = ∅. We say that P 0 has the strict complementarity property if J P ∪ J D = {1, . . . , p}. It follows from [Bonnans and Shapiro, 2013, Theorem 3.133] that strict complementarity always holds in linear programming. We refer to Remark 2.2.24 for a discussion regarding this assumption. (ii) Let (x * , λ * ) be a primal-dual solution to P 0 . Let (x, λ) be a cluster point of (x j+1 , λ j+1 ) j∈N and let (x jq+1 , λ jq+1 ) q∈N be a subsequence converging to this point. Pick q ∈ N. In view of Remark 3.5.11 and (3.36) we have that

v jq+1 -∇c(x jq+1 )λ jq+1 ∈ ∂(f + g)(x jq+1 ),
and

-∇c(x * )λ * ∈ ∂(f + g)(x * ).
Since f + g is convex, we deduce from the monotonicity of its subdifferential that

0 ≤ x jq+1 -x * , v jq+1 -x jq+1 -x * , ∇c(x jq+1 )λ jq+1 -∇c(x * )λ * .
(3.54) In addition, v jq+1 ≤ jq and (∀i ∈ {1, . . . , p}) c i is convex. Hence, we deduce from (3.54) that

0 ≤ x jq+1 -x * jq + c(x jq+1 ) -c(x * ) λ * + c(x * ) -c(x jq+1 ) λ jq+1 .
(3.55) From (3.37) and (3.39), c(x * ) λ * = 0 and, for every i ∈ {1, . . . , p}, c i (x jq+1 )(λ jq+1 ) i = -µ jq . Hence,

0 ≤ x jq+1 -x * jq µ jq + p - p i=1 λ * i (λ jq+1 ) i + c i (x * ) c i (x jq+1 ) . ( 3

.56)

Note that (∀i ∈ J P ) c i (x * ) = 0 and (∀i ∈ J D ) λ * i = 0. If J P (resp. J D ) is nonempty, we can then choose x * (resp. λ * ) such that c J P (x * ) ∈]-∞, 0[ #J P 66 A fast proximal IPM for constrained variational formulations (resp. λ * J D ∈]0, +∞[ #J D ). Consequently, by using Assumption 3.3.1, as q → +∞ ,(3.56) becomes

i∈J D λ * i λ i + i∈J P c i (x * ) c i (x) ≤ p, (3.57)
where we necessarily have

(∀i ∈ J D ) λ i > 0 and (∀i ∈ J P ) c i (x) < 0. (3.58)
Because of the strict complementarity, there are exactly p positive terms in the left-hand side of (3.57). Therefore, we can apply the arithmeticgeometric mean inequality which leads to

  i∈J D λ * i λ i     i∈J P c i (x * ) c i (x)   ≤ 1, (3.59)
with the convention that, if J P (resp. J D ) is empty, the corresponding product is equal to 1. From Theorem 3.4.5, we deduce that (x, λ) is a primal-dual solution to P 0 . Hence, (x * , λ) and (x, λ * ) are also primal-dual solutions to P 0 . Therefore, because (x, λ) satisfies (3.58), (3.59) also holds when either λ * = λ or x * = x. Consequently,

i∈J P (-c i (x * )) ≤ i∈J P (-c i (x)) and i∈J D λ * i ≤ i∈J D λ i . (3.60)
If there exist at least two distinct elements in S P (resp. S D ) then, from Propositions 3.5.16 and 3.5.17, J P (resp. J D ) is nonempty. It follows from (3.60) that x (resp. λ) is a solution to P P (resp. P D ). In turn, Lemma 3.5.18 guarantees that P P (resp. P D ) has a unique solution. Thus, there exists a unique cluster point for the primal (resp. dual) sequence and (x j+1 ) j∈N (resp. (λ j+1 ) j∈N ) converges to the analytic center of S P (resp. S D ).

Now that we have established some theoretical guarantees regarding the proposed algorithm, we show that PIPA performs well with respect to stateof-the-art methods on two applications in image processing, namely hyperspectral unmixing and joint geometry-texture decomposition and reconstruction of computed tomography data. Our numerical experiments demonstrate in addition the benefits of using a variable metric to accelerate the convergence. § 3.6 Application: Hyperspectral unmixing Hyperspectral imaging devices are remote sensing systems that acquire the emitting light spectrum of a distant scene, here modeled as a 2D image [Bioucas-Dias et al., 2012]. Let s and r be respectively the number of acquired spectral bands and pixels in the image, and let Y ∈ R s×r denote the measured hyperspectral cube. Assume that we have access to a library S ∈ R s×q , where each column of S contains the spectral signatures of one material (or endmember) among q that are expected to be present in the scene. The proportion or abundance of every material in every pixel is described through the abundance matrix X ∈ R q×r . The following linear model is frequently used to relate the data, the endmembers and the abundances, when there is no microscopic interaction between the materials:

Y = SX + Ω,
with Ω ∈ R s×r a realization of an additive white Gaussian noise. The estimation of X from Y and S is an inverse problem called unmixing [START_REF] Chan | A two-stage method for spectral-spatial classification of hyperspectral images[END_REF]].

Problem formulation

minimize X∈R q×r 1 2 Y -SX 2 2 + κ q i=1 (W X i ) d 1 subject to (∀j ∈ {1, . . . , r}) q i=1 X i,j ≤ 1 (∀i ∈ {1, . . . , q})(∀j ∈ {1, . . . , r}) X i,j ≥ 0, (3.61)
where • 2 denotes the Frobenius norm, (∀i ∈ {1, . . . , q}) X i ∈ R r is the ith line of the abundance matrix X, W ∈ R r×r is a wavelet decomposition operator, (•) d 1 is the 1 -norm of the detail wavelet coefficients, and κ ∈ R + is a regularization parameter. It is worth noting that the linear constraints account for the atmospheric absorption [START_REF] Keshava | Spectral unmixing[END_REF] since, for every pixel, the sum of all fractional abundances may be less than one. Moreover, the wavelet-based penalization [START_REF] Pustelnik | Waveletbased image deconvolution and reconstruction[END_REF] allows us to enforce useful spatial regularity on the sought abundance maps. From this point forward, the vectorizations of X ∈ R q×r and Y ∈ R s×r , in lexicographic order, are denoted by x ∈ R m with m = qr and y ∈ R sr , respectively, and ⊗ denotes the Kronecker product. Problem (3.61) can thus be re-written as in (3.1), with p = m + r and (∀x ∈ R m ) g(x) = 1/2 y -(I r ⊗ S)x 2 , f (x) = κ q i=1 (W P i x) d 1 , (∀i ∈ {1, . . . , q}) P i ∈ R r×m is a decimation matrix such that P i x = X i , and c(x) = Axb with

A = I r ⊗ 1 q -I m and b = 1 r 0 m .
The resulting minimization problem satisfies Assumptions 3.2.1, 3.4.2 and 3.4.6. Hence, Theorems 3.4.1, 3.4.4 and 3.4.5 regarding convergence and convergence rate of Algorithm 11 hold. In addition, in the considered example, the rank of S is equal to q, so there exists a unique solution to (3.61) and Theorem 3.4.7(i) holds.

Realistic data simulation and test configuration

In order to simulate Y , we make use of the Urban 1 dataset, which provides the spectral signatures and abundance maps for q = 6 materials in s = 162 spectral bands. We consider images of size r = 256 × 256. The product of the spectral library and attenuated abundance map is corrupted with an additive white Gaussian noise with a standard deviation of 0.06. Let the signal-to-noise ratio be defined as

SNR = 20 log 10 ( x / x -x ),
where x is the ground-truth for x. In addition, for each material i ∈ {1, . . . , q}, the signal-to-noise ratio of its associated abundance map

X i ∈ R r is SNR i = 20 log 10 X i / X i -X i .
The regularization weight κ is tuned by a grid search so as to reach the largest SNR, in that case κ = 0.01. Regarding the operator W , we selected an orthogonal Daubechies 4 wavelet decomposition performed over 2 resolution levels. As for the variable metric, we consider two cases: the proposed method without variable metric (taken as the identity matrix), which is referred to as PIPA; and PIPA-VM, which denotes the case when, following the strategy in [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF], for every j ∈ N, the variable metric is chosen as the Hessian of ϕ µ j . The proximity operator of the regularization term in the variable metric is computed numerically using Algorithm 5.

Let us now discuss the boundedness condition required for the variable metrics in PIPA-VM. For every x ∈ D, we have

ϕ µ (x) = 1 2 y -(I r ⊗ S)x 2 -µ p i=1 ln(b i -A i x),
1 www.escience.cn/people/feiyunZHU/Dataset_GT.html

where A i ∈ R m (resp. b i ∈ R) is the ith row (resp. component) of A (resp. of b). For every µ ∈ R * + and x ∈ D, the Hessian of ϕ µ at x is equal to

∇ 2 ϕ µ (x) = I r ⊗ (S S) + µ p i=1 A i A i (A i x -b i ) 2 .
Finally, in view of Remark 3.5.7, we deduce that there exist c and c in ]-∞, 0[ such that, for every k ∈ N and every i ∈ {1, . . . , p}, c

≤ A i x k -b i ≤ c < 0.
Since the rank of A is equal to m, for every µ ∈ R * + , the aforementioned variable metrics are bounded from below and above by strictly positive constants, as required in Algorithm 11.

According to Assumption 3.3.1, the precision with which the subproblem is solved must decrease faster to zero than the barrier parameter. Hence, in order to satisfy Assumption 3.3.1, we choose the barrier parameter and precision sequences as follows,

(∀j ∈ N) j = µ j ζ j , and µ j+1 = µ j ρ j , (3.62)
where ρ j ≥ ρ > 1, ζ > 1 and > 0. Regarding PIPA, we take = 10 3 , µ 0 = 1, ρ j = 1.5 for every j ∈ N and ζ = 1 + 10 -5 . For PIPA-VM we choose = 10 5 , µ 0 = 0.01, ζ = 1 + 10 -5 and, to avoid numerical instabilities when µ j is very small, we pick ρ j = 1.5 for every j ∈ N such that µ j ≥ 10 -6 , and decrease it gradually: if 4 × 10 -9 ≤ µ j < 10 -6 then ρ j = 1.1, if 10 -12 ≤ µ j < 4 × 10 -9 then ρ j = 1.01, and finally, if µ j < 10 -12 , then ρ j = 1.001.

We compare PIPA and PIPA-VM with three state-of-the-art convex optimization algorithms: ADMM [Setzer et We also implement preconditioned versions of ADMM and GFBS, which are referred to as ADMM-VM and GFBS-VM, respectively. ADMM-VM is based on [

Results

The solution to (3.61) with κ = 0, i.e. without regularization, can be obtained with the primal-dual interior point method from [Chouzenoux et al., 2014a]. It is referred to as IPLS and yields SNR = 11.02 dB after about 12 seconds, whereas solving the same problem with κ = 0.01 leads to a better reconstruction with SNR = 12.45 dB, illustrating the benefits of regularizing in this example. Figure 3.1 shows that the SNR increases faster with PIPA-VM than with the other algorithms. Moreover, it can be clearly seen in this figure that PIPA-VM exhibits a much faster convergence than PIPA, which demonstrates the advantage of using a variable metric in this example. The SNR obtained for each material after running the different methods for 11 seconds can be found in Table 3.1. For 5 out of 6 endmembers, the SNR of the abundance maps obtained with PIPA-VM after 11 seconds is better than for all other methods. In addition, for all materials PIPA-VM gives better results after 11 seconds than the non-regularized solution IPLS.

Visual results for Asphalt and Dirt materials are displayed in Figure 3.2, where we only show the results for IPLS, ADMM and PIPA-VM, since they outperformed the other methods. One can notice that, after running all algorithms for 11 sec, the abundance maps produced by PIPA-VM for these two materials are visually more satisfactory than the ones obtained with ADMM, while the non-regularized solution IPLS is significantly noisy.

In order to evaluate the algorithms based on their pointwise convergence, we let them run for a very large number of iterations and compute the relative distance between the current iterate and the solution x ∞ . As one can see on Figure 3.3, the sequence generated by PIPA-VM converges faster to the solution than the iterates produced by the other algorithms. It is finally worth noticing that, although an inexact computation of the proximity operator is performed in PIPA-VM, the method appears to be robust to the error generated by this approximation. In the following, we will denote by x t ∈ R r and x g ∈ R r , the texture and geometry components of a natural image x t+g ∈ R r , so that x t+g = x t +x g . The geometry x g represents a piecewise smooth version of the image, and can be extracted by using the total variation semi-norm [START_REF] Osher | Image decomposition and restoration using total variation minimization and the H -1 norm[END_REF]. The texture highlights local components with higher spatial frequencies. Depending on the considered application, different texture models can be found in the literature, based on wavelet decompositions [START_REF] Briceño-Arias | Proximal algorithms for multicomponent image recovery problems[END_REF] or on the Hölder exponent [START_REF] Pustelnik | Local regularity for texture segmentation: Combining wavelet leaders and proximal minimization[END_REF], to name only a few. Here, we will focus on material images in which the texture is located near the boundaries of different objects, as it can happen for instance in material image analysis, where the samples are subject to erosion and microporosity. Therefore, we will rely on the Laplacian detector for texture extraction, as the latter is known to be useful for edge and blob detection [START_REF] Haralick | Statistical and structural approaches to texture[END_REF]. Another issue is the acquisition procedure used for this type of images.

X-ray computed tomography

X-Ray Computed Tomography (CT) is a fast non-destructive scanning technique [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF], which is frequently used to acquire images from material samples. The acquisition process in CT consists in measuring the absorption of an object along a sampled grid of size L for N θ angular positions. It can be modeled through the discrete Radon projection operator H ∈ R n×r , with r the number of pixels and n = L × N θ the number of measurements. Although matrix H is sparse, it is also high-dimensional and ill-conditioned. Hence, reconstructing the image from the measured data y ∈ R n (also called sinogram) is a challenging inverse problem [START_REF] Chouzenoux | A majorizeminimize memory gradient algorithm applied to X-ray tomography[END_REF], which involves ringing artifacts, contrast issues due to beam hardening, and noise caused by sensor motion [START_REF] Ketcham | Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences[END_REF]. These artifacts can be partially removed by improving scanning techniques or by increasing the number of measurements, but this comes at a price and lots of applications can benefit from a faster subsampled CT acquisition.

Reconstructing the image from the acquired data is a heavy processing step which can introduce a bias in subsequent image processing tasks, such as the classification of material components [START_REF] Gouillart | Belief propagation reconstruction for discrete tomography[END_REF]. Therefore, we propose to perform jointly two tasks: the reconstruction and the geometry-texture decomposition. We show that this decomposition can be performed in a reasonable time with PIPA.

Problem formulation

We consider the following variational formulation, minimize

(x t ,x g )∈R r ×R r 1 2 F x t 2 + κTV(x g ) subject to x t + x g ∈ [x min , x max ] r x t ∈ [-α, α] r H(x t + x g ) -y ∞ ≤ χ (3.63)
where x min = 0 and x max = 1 are the minimal and maximal pixel intensity values, α > 0 is a range value parameter for the texture, κ > 0 is a regularization parameter, χ > 0 is an upperbound on the measurement uncertainty, TV denotes the isotropic total variation semi-norm with (zero) Dirichlet boundary conditions. Moreover, F = I r -∆ ∈ R r×r where ∆ ∈ R r×r is the Laplacian associated with the following 2D kernel padded with circulant assumption:

  0 1 0 1 -4 1 0 1 0   .
The first term in the objective function enforces edge detection in the texture, while the geometry is made piecewise smooth thanks to the total variation regularization. The first set of constraints represents bounds on the pixel values of the natural image. The texture, which is supposed to capture small variations in the image, is modeled as a zero-centered component in the second set of constraints, where we take α = x max /3. The last constraint is the data-fit term, which can be decomposed into 2 × n linear inequalities. Hence, the constraints can be reformulated as

Ax -b ∈] -∞, 0] p where x = [(x t ) , (x g ) ] , M ∈ R p×m , b ∈ R p , p = 2(n + 2r), m = 2r, M =         I r I r -I r -I r I r 0 r×r -I r 0 r×r H H -H -H         and b =         x max 1 r -x min 1 r α1 r α1 r y + χ1 n -y + χ1 n         .
(3.64) As in the previous example, (3.63) is an instance of Problem (3.1) where (∀x ∈ R m ) f (x) = κTV(P g x) and g(x) = 1/2 F P t x2 , with P t = (I r 0 r×r ) and P g = (0 r×r I r ). It can be noted that Assumptions 3.2.1 and 3.4.6 are satisfied. Thus, Theorems 3.4.1 and 3.4.5 hold. Since Assumption 3.4.2 does not hold, Theorem 3.4.4 does not apply here. The solution to (3.63) is not necessarily unique. Although the strict complementarity required to apply Theorem 3.4.7 is difficult to check, the convergence of PIPA to a single cluster point was observed in practice.

Initialization

In order to find an initial point that satisfies strictly the constraints, we set x t to zero. Following the method in [Boyd and Vandenberghe, 2004, Chap. 11.4], we solve the minimization problem below to initialize x g , minimize

(s,x g )∈R×R r s subject to s ≥ 0, x g ∈ [x min , x max ] r Hx g -y ∞ ≤ χ + s, (3.65)
where s ≥ 0 is the maximal infeasibility. The ground-truth natural image x t+g satisfies Hx t+gy ∞ < χ so that the solution to (3.65) is reached for s = 0. Problem (3.65) is a linear programming problem just like problem (2.18). In addition, it is straightforward to check that Assumption 2.2.23 is satisfied. Hence, we can use the Newton barrier method summarized in Algorithm 10 for solving problem (3.65). We use the code from [Boyd and Vandenberghe, 2004, Chap. 11.4] which is available online 2 This algorithm generates iterates that belong to ]x min , x max [ r so that, in our numerical experiments, we are able to find a strictly feasible initial point in a reasonable time.

Variable metric and hyperparameters

Let µ j ∈ R * + and k ∈ N. For every x ∈ D, the Hessian of ϕ µ j at x ∈ R m is equal to (∀i ∈ {1, . . . , r})

∇ 2 ϕ µ j (x) = F F + µ j Diag(d 3 (x)) + G(x) G(x) G(x) G(x) , ( 3
(d 1 (x)) i = x t+g i -x min -2 + x max -x t+g i -2 (∀j ∈ {1, . . . , n}) (d 2 (x)) j = Hx t+g -y j + χ -2 + y -Hx t+g j + χ -2 (∀i ∈ {1, . . . , r}) (d 3 (x)) i = x t i + α -2 + α -x t i -2
Given the huge size and ill-conditioning of H, the inversion of ∇ 2 ϕ µ j (x) is hardly feasible. Hence, instead of using the full Hessian of ϕ µ j for the variable metric as in Section 3.6, we propose to use an upper bound of it, i.e. M k ∈ S + m such that M k -∇ 2 ϕ µ j (x k ) also belongs to S + m . We propose to majorize µ j Diag(d 3 (x)) by β(x)I r where

β(x) = max 1≤i≤r µ j (d 3 (x)) i .
For H Diag(d 2 (x))H, we propose to follow the strategy in [Chouzenoux et al., 2014c] and upper-bound it by the diagonal matrix Diag(P d 2 (x)) ∈ R r×r , where P ∈ R n×r is such that for every i ∈ {1, . . . , r} and j ∈ {1, . . . , n},

P j,i = H j,i r s=1 H j,s .
This leads to the following variable metric in Algorithm 11,

(∀k ∈ N) M k = F F + β(x k )I r + D(x k ) D(x k ) D(x k ) D(x k ) (3.68) where D(x k ) = µ j Diag(d 1 (x k )) + Diag(P d 2 (x k )) .
Since D(x k ) is diagonal, the operator M k is straightforward to invert using the Schur formula.

In addition, similarly to Section 3.6, we deduce from Remark 3.5.7 that matrix (3.68) satisfies the boundedness condition required in Algorithm 11. In order to compute the proximity operator of f in such variable metric, we use Algorithm 7. Regarding the hyperparameters in the proposed method, in order to satisfy Assumption 3.3.1, as in the previous example, we take sequences of the form (3.62) with µ 0 = 10 -3 , ζ = 1 + 10 -5 , = 8.3 × 10 3 and ρ j = 1.1 for every j ∈ N.

Test settings

We perform the joint reconstruction and decomposition of two high-quality scans, referred to as Glass and Agaricus, which are displayed in Figure 3.4. These images are of size r = 128×128. The discrete Radon operator H models parallel projections along N θ = 180 angular positions on a detector grid of size L = 128, so that n = 180 × 128. To account for measurement uncertainty, the sinograms are degraded with a uniform noise with an amplitude χ equal to 2% of the maximal entry of y. We set manually the regularization parameter κ so that it leads to a visually satisfatory decomposition: it is set to 0.25 for Glass and to 0.5 for Agaricus. The proposed algorithm PIPA-VM is compared to ADMM summarized in Algorithm 8, which was the most competitive method in Section 3.6. Remark that, in order to make the implementation of ADMM feasible, we follow the same strategy as in [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF], and alternate the minimization on the splitting variables. In our example, we need seven splitting variables. Since ADMM does not require a feasible starting point, we run it with two different initializations: ADMM1 refers to ADMM initialized like PIPA-VM, and ADMM2 refers to ADMM initialized with x t taken as the zero vector and x g set to 1/2(x min +x max )1 r . All computational times are given for experiments run on Matlab 2018b on an Intel Xeon CPU E5-1650 at frequency 3.20 GHz.

Results

To compare the convergence speed of the different methods, we plot for each of them the relative distance between the current iterate x and the final solution x ∞ , obtained after running the algorithms for 12 hours. As one can see in Figures 3.5 and 3.6, PIPA-VM converges faster to its limit point than ADMM for both initializations. Remark that the time necessary to solve (3.65) and to find a feasible point is taken into account in the graphs. The results clearly show the advantage of using a feasible starting point over a simple initial guess.

To assess the visual quality of the geometry-texture decomposition we consider the solution obtained after reaching the stopping criterion xx ∞ / x ∞ ≤ 10 -2 . This accuracy is reached first for PIPA-VM after 14 min for Glass and 18 min for Agaricus. The corresponding visual decomposition and reconstruction after these durations are given in Figure 3.7. As one can see in this figure, PIPA-VM identifies correctly the geometry as an almost piecewise-constant image, free from locally-fast varying components like the gills in the Agaricus mushroom. Moreover, the texture obtained for Glass image captures well the elements on the borders between the two species in presence, and the Agaricus gills can be found in the texture (Figure 3.7 bottom left).

Finally, we evaluate the reconstruction quality based on the signal-tonoise ratio: SNR = 20 log 10 ( x t+g / x t+gx t+g ), (3.69) where x t+g denotes the ground-truth image. The SNR values obtained with the three methods for the reconstructions x t+g after the same durations are summarized in Table 3.2. This table shows that the SNR obtained with PIPA-VM is better than with the other methods after the same durations. It can be further observed in Figures 3.8 and 3.9, that even if the SNR converges to the same value for the three methods, PIPA-VM follows a path which would lead to a better SNR if stopped before convergence. § 3.8 Summary

In this chapter, we have shown that it is possible to combine efficiently two powerful optimization frameworks: proximal splitting methods and interior point algorithms. One interesting feature of the resulting iterative method is the use of a variable metric, which can boost the convergence, as illustrated in our hyperspectral unmixing application. The convergence of the proposed method, as well as a convergence rate for the inner loop, have been obtained under suitable assumptions. Regarding numerical performance, as shown on two large-scale image processing applications, our method compares favorably in terms of speed of convergence with state-of-the-art algorithms.

Although we are able to produce fast algorithms for solving optimization problems, they are often sensitive to hyperparameters that need to be set carefully. Moreover, variational strategies for addressing inverse problems may suffer from downsides, some of which are mentioned in Section 2.1.2.2. In that regard, neural networks appear as a powerful alternative, which is studied in the next chapter.

-Chapter 4 -

An optimization-inspired neural network architecture for image deblurring

As we mentioned in Section 2.1.4.2, neural networks demonstrate state-ofthe-art performance in several applications linked to inverse problems. However, these methods suffer from important downsides including a lack of interpretability, few mathematical results, and no real possibility for imposing constraints on their parameters. On the other hand, variational formulations, whose limitations are discussed in Section 2.1.2.2, allow to incorporate some a priori knowledge or desirable properties for the solution, and the convergence results regarding classical optimization algoritms are well established. Hence, DNNs and variational-based methods can be considered as complementary. One straightforward way to combine the benefits of both approaches is to unfold an iterative method and untie the parameters of both the model and the algorithm across the layers of the network [START_REF] Hershey | Deep unfolding: Model-based inspiration of novel deep architectures[END_REF].

In this chapter, we propose a novel neural network architecture called iRestNet, which is obtained by unfolding a proximal interior point algorithm over a finite number of iterations. More information about deep unfolding methods can be found in Section 2.1.4.3. One key feature of this algorithm is that it produces only feasible iterates thanks to a logarithmic barrier. This barrier enables prior knowledge to be directly incorporated into iRestNet and, as opposed to a projection onto the feasible set, it allows differentiation and gradient backpropagation throughout the network. Hence, gradient descent can be used for training. The stepsize, barrier parameter, and regularization weight are untied across the network and learned for each layer. Thus, once the network has been trained, its application on test images requires only a short execution time per image without any parameter search, as opposed to traditional variational methods.

This chapter is organized as follows. First, we present our contributions with regards to related works in Section 4.1, then, in Section 4.2, we describe 82 An optimization-inspired neural network for image deblurring the algorithm which is at the core of our method. Since the proposed algorithm requires the computation of the proximity operator of the barrier, we provide the latter for three useful cases in Section 4.3. The proposed neural network architecture and its associated backpropagation method are presented in Section 4.4. Then, in Section 4.5, we conduct a stability analysis of the proposed network when the data fidelity term and the regularization function are quadratic. Section 4.6 is dedicated to numerical experiments and comparison to state-of-the-art methods for image deblurring; finally, some conclusions are drawn in Section 4.7. § 4.1 Link with related works Several recent works consider replacing handcrafted algorithms by learned iterative methods [START_REF] Andrychowicz | Learning to learn by gradient descent by gradient descent[END_REF][START_REF] Li | Learning to optimize[END_REF]. In these approaches, the goal is to find the minimizer of a given objective function, whereas, in the proposed method, the architecture is inspired by an optimization strategy applied to the minimization of an objective function, but a better indicator of perceptual quality is optimized during the training step.

Only a few works so far have considered combining IPMs with deep learning. Every layer of the network from [Amos and Kolter, 2017] solves a small quadratic problem using an IPM, while in [Trafalis et al., 1997], hard constraints are enforced on weights by using the logarithmic barrier function during training. More recently, an interior point strategy was used to design a recurrent network, whose purpose is to solve a specific convex constrained problem [START_REF] Krasopoulos | An Interior Point Recurrent Neural Network for Convex Optimization Problems[END_REF]]. In our case, however, we have two distinct objective functions. The first one leads to a constrained problem from which the proposed architecture is inspired, while the second one is used during training as a loss function. It is worth noting that the output of the trained network is not necessarily a minimizer of the first objective. Moreover, the second objective could not be a substitute to the first one since it requires the knowledge of the ground-truth, which is available for training time but not in testing conditions. In addition, iRestNet appears to have more flexibility since the regularization weight can vary among layers.

To the best of our knowledge, this chapter presents the first architecture corresponding to a deep unfolded version of an interior point algorithm with untied stepsize and regularization parameter. As opposed to other unfolding methods like [START_REF] Mardani | Recurrent generative adversarial networks for proximal learning and automated compressive image recovery[END_REF][START_REF] Diamond | Unrolled optimization with deep priors[END_REF], the proximity operator and the regularization term are kept explicit, which establishes a direct relation between the original algorithm and the network. Other contributions of this work include the expression of the required proximity operator, and of its corresponding gradient, for three standard variational formulations, along with numerical experiments demonstrating the benefit of using the proposed approach over other machine learning and variational methods for image deblurring. § 4.2 Proposed algorithm

Variational formulation and assumptions

Let us focus on inverse problems related to model (2.1), where the observation model is linear. Under this assumption, (2.1) becomes

y = D(Hx), (4.1)
where y ∈ R n is the observed data, x ∈ R m is the sought signal or image, H ∈ R n×m is the observation operator, and D is the noise perturbation operator. The linear operator H is assumed to be known from a physical model or prior identification step [START_REF] Lagendijk | The Handbook of Image and Video Processing[END_REF]; Xu and Jia, 2010]. As detailed in Section 2.1.2, the sought image x can be classically approximated by the minimizer of a penalized cost function expressed as the sum of a data-fitting term, which measures the fidelity of the solution to the observation model (4.1), and a regularization term, which is introduced so as to avoid meaningless solutions and improve stability relative to noise. This leads to problem (2.3), which is rewritten below in the case of a linear observation model.

minimize x∈C f (Hx, y) + λR(x) (4.2)
In the remaining of this chapter, we will assume that, for every y ∈ R n , f (•, y) ∈ Γ 0 (R n ) and R ∈ Γ 0 (R m ) are twice-differentiable functions. Note that such assumption is necessary to define the derivative steps involved in the backpropagation procedure for the training of our network. The feasible set C is defined by p inequality constraints, which enforce the fulfillment of some properties that are expected to be satisfied a priori by the image:

C = {x ∈ R m | (∀i ∈ {1, . . . , p}) c i (x) ≤ 0}, (4.3)
where, for every i ∈ {1, . . . , p}, c i ∈ Γ 0 (R m ). The strict interior of the feasible domain is equal to

D = {x ∈ R m | (∀i ∈ {1, . . . , p}) c i (x) < 0},
and it is assumed to be nonempty. The interest of such constraints in the resolution of inverse problems has been discussed in the previous chapter in Remark 3.2.2. Finally, we will assume that either f (H•, y) + λR is coercive, or C is bounded. Then the existence of solutions for (4.2) is guaranteed. It is worthy to emphasize that a large class of penalized formulations encountered in the literature of image restoration fulfills the above requirements, see e.g. [START_REF] Durand | Stability of minimizers of regularized least squares objective functions I: study of the local behaviour[END_REF]] and references therein. For simplicity, for all

(x, y, λ) ∈ R m × R n × R + , we define h(x, y, λ) = f (Hx, y) + λR(x),
and

∇ 1 h(x, y, λ) = H ∇ 1 f (Hx, y) + λ∇R(x),
where ∇ 1 f is the partial gradient of f with respect to its first variable.

In general, problem (4.2) does not have a closed-form solution on account of the inequality constraints, even for simple regularizations, hence an iterative solver must be used. Several resolution approaches are available, either based on projected gradient strategies [START_REF] Iusem | On the convergence properties of the projected gradient method for convex optimization[END_REF]Bonettini and Prato, 2015a], ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], primal-dual schemes [START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF], or interior point techniques [Bonettini and Serafini, 2009b]. As it is shown in Chapter 3, combining the interior point framework with a proximal forward-backward strategy leads to very competitive solvers for inverse problems. As detailed in Section 2.2.5, the idea behind IPMs is to replace the initial constrained optimization problem by a sequence of unconstrained subproblems of the form:

min x∈R m f (Hx, y) + λR(x) + µB(x) (4.4)
where µ ∈ R * + is the barrier parameter and B : R m →] -∞, +∞] is the logarithmic barrier function whose definition is recalled below.

(∀x ∈ R m ) B(x) =      - p i=1 ln(-c i (x)) if x ∈ D, +∞ otherwise. (4.5)
We have made the assumption that either f (H•, y) + λR is coercive, or C is bounded, hence, the set of solutions to (4.2) is bounded. Since D is not empty, we can apply Corollary 3.5.1, which ensures the existence of solutions to (4.4).

Proposed iterative schemes

Thanks to the proximity operator, the IPM from [START_REF] Kaplan | Proximal methods in view of interior-point strategies[END_REF]] does not require any matrix inversion, as opposed to the classical Newton barrier method. When the proximity operator is computed in an exact manner, the proposed IPM can be rewritten as Algorithm 13, whose convergence has been proven under some assumptions [Kaplan and Tichatschke, 1998, Theorem 4.1].

Algorithm 13: Exact version of the proximal IPM in [START_REF] Kaplan | Proximal methods in view of interior-point strategies[END_REF]] applied to problem (4.2).

Let γ ∈ R * + and (γ k ) k∈N be a sequence such that

(∀k ∈ N) γ ≤ γ k ; Initialization: Let x 0 ∈ D; for k = 0, 1, . . . do x k+1 = prox γ k (h(•,y,λ)+µ k B) (x k ) end
Algorithm 13 requires evaluating the proximity operator of the sum of the barrier and the regularized cost function, which can be an issue since, in most of the cases, this operator does not have a closed-form expression. This is the reason why we propose to modify it by introducing a forward step, which leads to Algorithm 14.

Algorithm 14: Proposed forward-backward proximal IPM.

Let γ ∈ R * + and (γ k ) k∈N be a sequence such that (∀k ∈ N) γ ≤ γ k ; Initialization: Let x 0 ∈ D; for k = 0, 1, . . . do x k+1 = prox γ k µ k B (x k -γ k ∇ 1 h (x k , y, λ)) end
To the best of our knowledge, there is no available convergence study for Algorithm 14 among the literature of interior-point methods. There exist links between the above algorithm and the diagonal or penalization method introduced in [START_REF] Czarnecki | Splitting forward-backward penalty scheme for constrained variational problems[END_REF]. Indeed, taking A ≡ 0 and Ψ 1 ≡ 0 in [START_REF] Czarnecki | Splitting forward-backward penalty scheme for constrained variational problems[END_REF] leads to Algorithm 14, whose convergence is proven. However, there are some key differences between both approaches, namely i) in [START_REF] Czarnecki | Splitting forward-backward penalty scheme for constrained variational problems[END_REF], the barrier parameter tends to infinity while it goes to zero in our case, and ii) the algorithm in [START_REF] Czarnecki | Splitting forward-backward penalty scheme for constrained variational problems[END_REF] solves a hierarchical minimization problem instead of the constrained optimization problem (4.2). It is worth noting that Algorithm 14 only requires computing the proximity operator of the logarithmic barrier. We will provide its expression in Section 4.3 for three different types of constraints.

Limitations

In IPMs, the barrier parameter and stepsize sequences, (µ k ) k∈N and (γ k ) k∈N , are usually set by following some heuristic rules, which ensure the convergence of the method to a minimizer of the considered objective function. However, handcrafted variational formulations do not necessarily capture perceptual image quality well. These heuristics can thus lead to a loss in terms of efficiency and versatility of the resulting restoration schemes. Moreover, as already mentioned, an accurate setting of the regularization weights is particularly critical in order to obtain a satisfactory image quality when using such penalized restoration approaches. Existing approaches for selecting λ, which are based on statistical considerations, are usually associated with a substantial increase of the computational cost.

To overcome these limitations, we propose to unfold Algorithm 14 over a given number of iterations and to learn the stepsize, the barrier and the regularization parameters for every iteration in a supervised fashion. Our machine learning method will make use of gradient backpropagation for its training step. The latter requires the derivatives of the proximity operator in Algorithm 14 with respect to its input and to the aforementioned parameters which are to be learned. Therefore, we first conduct an analysis of the proximity operator of the barrier and of its derivatives, for three examples of interest in Section 4.3. § 4.3 Proximity operator of the barrier Let B be defined as in (4.5) and for all µ ∈ R * + , γ ∈ R * + and x ∈ R m , let ϕ be defined as follows:

ϕ(x, µ, γ) = prox γµB (x).
We provide in this section expressions of ϕ and of its derivatives with respect to its input variable x and the involved barrier and stepsize parameters (µ, γ), for three common types of constraints. The latter will be necessary for training the proposed neural network using a gradient backpropagation scheme.

Affine constraints

Let us first consider the following half-space constraint:

C = {x ∈ R m | a x ≤ b}, (4.6)
with a ∈ R m \ {0 m } and b ∈ R. Proposition 4.3.1 Let γ ∈ R * + , µ ∈ R * +
, and let B be the function associated to (4.6), defined as

(∀x ∈ R m ) B(x) = -ln(b -a x) if a x < b, +∞ otherwise.
Then, for every x ∈ R m , the proximity operator of γµB at x is given by

ϕ(x, µ, γ) = x + b -a x -(b -a x) 2 + 4γµ a 2 2 a 2 a. (4.7)
In addition, the Jacobian matrix of ϕ with respect to x and the gradients of ϕ with respect to µ and γ are given by

J (x) ϕ (x, µ, γ) = I m - 1 2 a 2 1 + a x -b (b -a x) 2 + 4γµ a 2
aa , (4.8)

∇ (µ) ϕ (x, µ, γ) = -γ (b -a x) 2 + 4γµ a 2 a, (4.9)
and

∇ (γ) ϕ (x, µ, γ) = -µ (b -a x) 2 + 4γµ a 2 a. (4.10)
Proof. The expression for the proximity operator (4.7) directly follows from Example 2.2.15(i), and Propositions 2.2.12 and 2.2.13. Taking the derivative of (4.7) with respect to x, µ and γ leads to (4.8)-(4.10).

Hyperslab constraints

We now consider the following hyperslab set:

C = {x ∈ R m | b ≤ a x ≤ b}, (4.11)
where

a ∈ R m \ {0 m }, b ∈ R and b ∈ R with b < b. Proposition 4.3.2 Let γ ∈ R * + , µ ∈ R * +
, and let B be the barrier function associated to (4.11), defined as

(∀x ∈ R m ) B(x) = -ln(b -a x) -ln(a x -b) if b < a x < b, +∞ otherwise.
Then, for every x ∈ R m , the proximity operator of γµB at x is given by

ϕ(x, µ, γ) = x + κ(x, µ, γ) -a x a 2 a, (4.12)
where κ(x, µ, γ) is the unique solution in ]b, b[, of the following cubic equation:

0 = z 3 -(b + b + a x)z 2 + (bb + a x(b + b) -2γµ a 2 )z -bba x + γµ(b + b) a 2 . (4.13)
In addition, the Jacobian matrix of ϕ with respect to x and the gradients of ϕ with respect to µ and γ are given by

J (x) ϕ (x, µ, γ) = I m + 1 a 2 (b -κ(x, µ, γ))(b -κ(x, µ, γ)) η(x, µ, γ) -1 aa , (4.14)
∇ (µ) ϕ (x, µ, γ) = -γ(b + b -2κ(x, µ, γ)) η(x, µ, γ) a, (4.15)
and

∇ (γ) ϕ (x, µ, γ) = -µ(b + b -2κ(x, µ, γ)) η(x, µ, γ) a, (4.16)
where

η(x, µ, γ) = (b -κ(x, µ, γ))(b -κ(x, µ, γ)) -(b + b -2κ(x, µ, γ))(κ(x, µ, γ) -a x) -2γµ a 2 . Proof. Let x ∈ R m , γ ∈ R * + , and µ ∈ R * + .
The expression for the proximity operator (4.12) follows from Example 2.2.15(ii) and Proposition 2.2.13. Let F be defined as follows:

F (x, µ, γ, z) = (b -z)(b -z)(z -a x) + γµ(b + b -2z) a 2 , (4.17)
for z ∈]b, b[. Expanding (4.17) gives the following:

F (x, µ, γ, z) = z 3 -(a x + b + b)z 2 + (bb + a x(b + b) -2γµ a 2 )z -bba x + γµ(b + b) a 2 .
Hence, by definition of κ(x, µ, γ), we have F (x, µ, γ, κ(x, µ, γ)) = 0. In addition, the derivative of F with respect to its last variable is equal to

∇F (z) (x, µ, γ, z) = (b -z)(b -z) -(b + b -2z)(z -a x) -2γµ a 2 .
By construction, (b-κ(x, µ, γ))(b-κ(x, µ, γ)) < 0. Moreover, -2γµ a 2 < 0 and, since F (x, µ, γ, κ(x, µ, γ)) = 0, it follows that (b + b -2κ(x, µ, γ)) and κ(x, µ, γ)a x share the same sign. Hence,

η(x, µ, γ) = ∇F (z) (x, µ, γ, κ(x, µ, γ)) = 0.
From the implicit function theorem [START_REF] Dontchev | Implicit functions and solution mappings[END_REF], Theorem 1B.1], we deduce that the gradient of κ with respect to x and the partial derivatives of κ with respect to µ and γ exist and are equal to

∇κ (x) (x, µ, γ) = (b -κ(x, µ, γ))(b -κ(x, µ, γ)) η(x, µ, γ) a, (4.18)
∇κ (µ) (x, µ, γ) = -γ a 2 (b + b -2κ(x, µ, γ)) η(x, µ, γ) , (4.19)
and

∇κ (γ) (x, µ, γ) = -µ a 2 (b + b -2κ(x, µ, γ)) η(x, µ, γ) . (4.20)
Differentiating (4.12) with respect to x, µ and γ and using It can be noted that Example 2.2.15(ii) is a special case of Proposition 4.3.2. The three roots of (4.13) can easily be computed using the Cardano formula.

Bounded 2 -norm

We now consider the case when the feasible set in (4.2) is a Euclidean ball,

C = {x ∈ R m | x -c 2 ≤ α}, (4.21)
with α ∈ R * + and c ∈ R m .

(∀x ∈ R m ) B(x) = -ln(α -x -c 2 ) if x -c 2 < α, +∞ otherwise.
Then, for every x ∈ R m , the proximity operator of γµB at x is given by

ϕ(x, µ, γ) = c + α -κ(x, µ, γ) 2 α -κ(x, µ, γ) 2 + 2γµ (x -c), (4.22)
where κ(x, µ, γ) is the unique solution in [0, √ α[ of the cubic equation:

0 = z 3 -x -c z 2 -(α + 2γµ)z + α x -c . (4.23)
In addition, the Jacobian matrix of ϕ with respect to x and the gradients of ϕ with respect to µ and γ are given by

J (x) ϕ (x, µ, γ) = α -ϕ(x, µ, γ) -c 2 α -ϕ(x, µ, γ) -c 2 + 2γµ M (x, µ, γ), (4.24)
∇ (µ) ϕ (x, µ, γ) = -2γ α -ϕ(x, µ, γ) -c 2 + 2γµ M (x, µ, γ)(ϕ(x, µ, γ) -c), (4.25) 90
An optimization-inspired neural network for image deblurring and c) .

∇ (γ) ϕ (x, µ, γ) = -2µ α -ϕ(x, µ, γ) -c 2 + 2γµ M (x, µ, γ)(ϕ(x, µ, γ) -c), (4.26) where
M (x, µ, γ) = I m - 2(x -ϕ(x, µ, γ))(ϕ(x, µ, γ) -c) α -3 ϕ(x, µ, γ) -c 2 + 2γµ + 2(ϕ(x, µ, γ) -c) (x
(4.27)

Proof. Let x ∈ R m , γ ∈ R * + , µ ∈ R * + .
Let us first consider the case when c = 0. We denote with ϕ 0 the following proximity operator:

ϕ 0 (x, µ, γ) = argmin u∈D 1 2 x -u 2 -γµ ln(α -u 2 ).
Hence, ϕ 0 (x, µ, γ) 2 < α and ϕ 0 (x, µ, γ) is a solution to the following equation:

0 = ϕ 0 (x, µ, γ) -x + 2γµ α -ϕ 0 (x, µ, γ) 2 ϕ 0 (x, µ, γ). (4.28)
Since αϕ 0 (x, µ, γ) 2 + 2γµ > 0, (4.28) becomes

ϕ 0 (x, µ, γ) = α -ϕ 0 (x, µ, γ) 2 α -ϕ 0 (x, µ, γ) 2 + 2γµ
x.

(4.29)

By taking the norm in both sides of (4.29), we deduce that ϕ 0 (x, µ, γ) = κ(x, µ, γ) is a solution to the cubic equation (4.23). Since the proximity operator at a given x is uniquely defined, there exists only one real solution to (4.23) which belongs to [0, √ α[. Plugging the latter into (4.29) leads to (4.22). The analysis when c = 0 is deduced from the case c = 0 by using Proposition 2.2.12: the proximity operator of γµB at x is given by

ϕ(x, µ, γ) = c + ϕ 0 (x -c, µ, γ). (4.30)
Now we can study the derivatives of ϕ 0 . For every v ∈ R m , let F be defined as

F (x, µ, γ, v) = (α -v 2 )(v -x) + 2γµv.
The Jacobian of F with respect to its last variable is equal to

J (v) F (x, µ, γ, v) = (α -v 2 + 2γµ)I m + 2(x -v)v .
Since αϕ 0 (x, µ, γ) 2 > 0, according to the Sherman-Morrison Lemma [START_REF] Bartlett | An inverse matrix adjustment arising in discriminant analysis[END_REF],

J (v) F (x, µ, γ, ϕ 0 (x, µ, γ)) is invertible if and only if α -ϕ 0 (x, µ, γ) 2 + 2γµ + 2ϕ 0 (x, µ, γ) (x -ϕ 0 (x, µ, γ)) = 0.

Proximity operator of the barrier 91

Furthermore, it follows from (4.28) that

F (x, µ, γ, ϕ 0 (x, µ, γ)) = 0 m . (4.31)
Applying ϕ 0 (x, µ, γ) on (4.31) leads to ϕ 0 (x, µ, γ) (xϕ 0 (x, µ, γ)) ≥ 0.

In addition, αϕ 0 (x, µ, γ) 2 + 2γµ > 0. Hence, J

F (x, µ, γ, ϕ 0 (x, µ, γ)) is invertible and its inverse is given by the Sherman-Morrison formula:

J (v) F (x, µ, γ, ϕ 0 (x, µ, γ)) -1 = 1 α -ϕ 0 (x, µ, γ) 2 + 2γµ × I m - 2(x -ϕ 0 (x, µ, γ))ϕ 0 (x, µ, γ) α -3 ϕ 0 (x, µ, γ) 2 + 2γµ + 2ϕ 0 (x, µ, γ) x .
From the implicit function theorem [START_REF] Dontchev | Implicit functions and solution mappings[END_REF], Theorem 1B.1] we deduce that the Jacobian of ϕ 0 with respect to x and the gradients of ϕ 0 with respect to µ and γ exist and are equal to

J (x) ϕ 0 (x, µ, γ) = -J (v) F (x, µ, γ, ϕ 0 (x, µ, γ)) -1 J (x) F (x, µ, γ, ϕ 0 (x, µ, γ)), ∇ (µ) ϕ 0 (x, µ, γ) = -J (v) F (x, µ, γ, ϕ 0 (x, µ, γ)) -1 ∇ (µ) F (x, µ, γ, ϕ 0 (x, µ, γ)), and
∇ (γ) ϕ 0 (x, µ, γ) = -J (v) F (x, µ, γ, ϕ 0 (x, µ, γ)) -1 ∇ (γ)
F (x, µ, γ, ϕ 0 (x, µ, γ)).

When c = 0, the derivatives of ϕ are deduced from those of ϕ 0 using (4.30):

J (x) ϕ (x, µ, γ) = -J (v) F (x -c, µ, γ, ϕ(x, µ, γ) -c) -1 × J (x) F (x -c, µ, γ, ϕ(x, µ, γ) -c), ∇ (µ) ϕ (x, µ, γ) = -J (v) F (x -c, µ, γ, ϕ(x, µ, γ) -c) -1 × ∇ (µ) F (x -c, µ, γ, ϕ(x, µ, γ) -c), and
∇ (γ) ϕ (x, µ, γ) = -J (v) F (x -c, µ, γ, ϕ(x, µ, γ) -c) -1 × ∇ (γ) F (x -c, µ, γ, ϕ(x, µ, γ) -c),
which lead to (4.24)-(4.26).

Similarly to the previous case, the three solutions to (4.23) can be obtained thanks to the Cardano formula. The form of the resulting proximity As shown in this section, the proximity operator of the barrier is easily computable and differentiable for several classic types of constraints. Next, we detail the proposed approach in Section 4.4. § 4.4 iRestNet architecture

Overview

Our proposal is to adopt a supervised learning strategy in order to determine, from a training set of images, an optimal setting for the parameters of Algorithm 14, which should lead to an optimal image restoration quality. To this aim, Algorithm 14 is unfolded over K iterations and the regularization parameter λ is untied across the network, so as to provide more flexibility to the approach [START_REF] Hershey | Deep unfolding: Model-based inspiration of novel deep architectures[END_REF]. The update rule at a given iteration k ∈ {0, . . . , K -1} reads

x k+1 = A (x k , µ k , γ k , λ k ) with A (x k , µ k , γ k , λ k ) = prox γ k µ k B (x k -γ k ∇ 1 h (x k , y, λ k )) . (4.32)
For every k ∈ {0, . . . , K -1}, we build the kth layer L k as the association of three hidden structures, L

k , L

k and L (λ) k , followed by the update A. aim at inferring the barrier parameter µ k , the stepsize γ k and the regularization weight λ k , respectively. Since a finite number K of layers (i.e., updates) is used, the convergence of the resulting scheme is not an issue. Note that we also allow in our framework the use of a post-processing step after going through the K layers, that will be denoted as L pp . The resulting architecture is depicted in Figure 4.2.

Hidden structures

Let us now provide more details about the hidden structures. For every k ∈ {0, . . . , K -1}, the outputs

(µ k , γ k , λ k ) of the structures L (µ) k , L (γ) k , and L (λ)
k must be positive. To enforce such constraint, we use the Softplus function [START_REF] Dugas | Incorporating second-order functional knowledge for better option pricing[END_REF], defined below, which can be viewed as a smooth approximation of the ReLU activation function:

(∀z ∈ R) Softplus(z) = ln(1 + exp(z)).
Unlike the ReLU, the gradient of Softplus is never strictly equal to zero, which, given our architecture, helps propagate the gradient through the network. The stepsize is estimated as follows,

(∀k ∈ {0, . . . , K -1}) γ k = L (γ) k = Softplus (a k ) , (4.33)
where (a k ) 0≤k≤K-1 are scalar parameters of the network learned during training. The barrier parameter is obtained using two convolutional and average pooling layers followed by a fully connected layer. The detailed architecture of L

k is depicted in Figure 4.3. Traditional methods for estimating the regularization parameter generally depend on the signal-to-noise ratio and on the image statistics [START_REF] Vogel | Computational methods for inverse problems[END_REF]. For most applications the noise level is unknown and can be estimated, for instance, by applying a median filter over the wavelet diagonal coefficients of the image [START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Ramadhan | Image denoising by median filter in wavelet domain[END_REF]. This strategy is used in the numerical experiments presented in Section 4.6. The advantage is to yield a network which can handle datasets for which the signal-to-noise

k .

ratio is unknown and can vary within a reasonable range. The expression of L

k is then problem-dependent since its expression depends on the regularization function R. A specific example is given in Section 4.6 for the total variation regularization function.

Regarding the post-processing step L pp , its detailed architecture also depends on the task to be performed. An example is provided in Section 4.6 for the case of deblurring: the purpose of L pp is then to remove remaining artifacts using convolutional layers, residual learning, batch normalization, and dilation (the aforementioned notions are defined in Section 2.1.4.1).

Differential calculus

To train the neural network presented in Figure 4.2 using gradient descent, one needs to compute the gradient of x K with respect to the different parameters of the network. The chain rule can be applied since most of the steps in the network correspond to operators having straightforward derivatives. However, particular care should be taken when differentiating A. Since f and R are assumed to be twice differentiable, the only area of concern is related to prox γµB . If prox γµB is simple enough, automatic differentiation [START_REF] Paszke | Automatic differentiation in Py-Torch[END_REF]] can be used. Otherwise, as shown in Section 4.3, for common examples of barrier functions, the differential of this term is well-defined. The corresponding expressions for the derivatives are provided in Propositions 4.3.1-4.3.3. § 4

.5 Network stability

There is a growing interest for interpretable and stable deep learning architectures [Zhang et al., 2018;[START_REF] Malgouyres | Multilinear compressive sensing and an application to convolutional linear networks[END_REF]. One critical issue concerning neural networks is to guarantee that their performance re-mains acceptable when the input is perturbed. For example, the authors of [START_REF] Szegedy | Intriguing properties of neural networks[END_REF] show that the class prediction made by AlexNet can be arbitrarily changed by using small nonrandom perturbations on the test image. For some applications involving high risk and legal responsibility, for instance in medical image processing, the lack of theoretical guarantees is a significant curb on the utilization of deep learning approaches.

A recent work [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF] provides a theoretical framework which enables to evaluate the robustness of a network. In this section, we will focus on a subclass of problem (4.2) where both f (•, y) and R are quadratic functions. After highlighting the similarities between the proposed architecture and generic feedforward networks in that case, we will give explicit conditions under which the robustness of the proposed architecture is ensured.

Relation to generic deep neural networks

Although the proposed architecture may seem specific to Algorithm 14, it is actually very similar to generic feedforward neural networks, which are presented in Section 2.1.4.1. For the sake of simplicity, we will consider the variational problem,

minimize x∈C 1 2 Hx -y 2 + λ 2 Dx 2 , (4.34)
where y ∈ R m , H ∈ R m×m , D ∈ R m×m , and C is defined as in (4.3). Moreover, we assume that no post-processing layer L pp is used. Following the notation of Section 4.4, (∀k ∈ {0, . . . , K -1}) (µ k , γ k , λ k ) are given positive real numbers, K being the number of layers of the network. Then, for every k ∈ {0, . . . , K -1}, layer L k corresponds to the following update,

x k+1 = prox γ k µ k B x k -γ k H (Hx k -y) + λ k D Dx k = prox γ k µ k B I m -γ k H H + λ k D D x k + γ k H y ,
where B is defined as in (4.5). For every k ∈ {0, . . . , K -1}, we set

W k = I m -γ k H H + λ k D D , b k = γ k H y, R k = prox γ k µ k B . (4.35) Then, the K-layer network L K-1 • • • • • L 0 is equivalent to (2.7): R K-1 • (W K-1 • +b K-1 ) • • • • • R 0 • (W 0 • +b 0 ),
where (W k ) 0≤k≤K-1 and (b k ) 0≤k≤K-1 are interpreted as weight operators and bias parameters, respectively. The operators (R k ) 0≤k≤K-1 defined in (4.35) can be viewed as specific activation functions since, as shown in [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF], every standard activation function can be derived from a proximity operator. In addition, using [Bauschke and Combettes, 2017, Proposition 24.8(iii)], for every k ∈ {0, . . . , K -1}, R k can be re-written as the sum of a proximal activation operator [Combettes and Pesquet, 2018, Definition 2.20] and a bias.

Preliminary results

Before stating our main stability theorem, we recall the result from [Combettes and Pesquet, 2018, Lemma 3.3] in Proposition 4.5.1 below. We then derive Proposition 4.5.2, which will appear useful when addressing the robustness of the global network. Here, S m denotes the set of symmetric matrices in R m×m .

Proposition 4.5.1 [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF] Let K ≥ 1 be an integer and set θ -1 = 1. For every k ∈ {0, . . . , K -1}, let W k ∈ R m×m and let θ k be defined by

θ k = |||W k • • • • • W 0 ||| + k-1 =0 0≤j 0 <•••<j ≤k-1 |||W k • • • • • W j +1 |||× W j • • • • • W j -1 +1 • • • |||W j 0 • • • • • W 0 |||.
Then, for every k ∈ {0, . . . , K -1},

θ k = k =0 θ -1 |||W k • . . . • W l |||.
W -2 K (1 -α)I m -|||W ||| + 2θ ≤ 2 K α (4.36)
is satisfied if and only if one of the following conditions holds:

(i) β + + β -≤ 0 and θ ≤ 2 K-1 (2α -1); (ii) 0 ≤ β + + β -≤ 2 K+1 (1 -α) and 2θ ≤ β + + β -+ 2 K (2α -1); (iii) 2 K+1 (1 -α) ≤ β + + β -and θ ≤ 2 K-1 . Proof. Let α ∈ [1/2, 1]. Since W ∈ S m , we have, |||W ||| = max{β + , -β -}, and
W -2 K (1 -α)I m = max β + -2 K (1 -α), -β -+ 2 K (1 -α) .
(4.37) Three different cases arise that we review below. (ii

(i) If β + + β -≤ 0 then |||W ||| = -β -and β + -2 K (1 -α) ≤ -β -+ 2 K (1 -α).
) If 0 ≤ β + + β -≤ 2 K+1 (1 -α) then |||W ||| = β + and (4.38) is satisfied. Hence, W -2 K (1 -α)I m = -β -+ 2 K (1 -α).
Replacing |||W ||| and W -2 K (1α)I m by their value in (4.36) leads to Proposition 4.5.2(ii).

(iii) If 2 K+1 (1 -α) ≤ β + + β -then |||W ||| = β + and β + -2 K (1 -α) ≥ -β -+ 2 K (1 -α). (4.39)
From (4.37) and (4.39), we deduce that

W -2 K (1 -α)I m = β + - 2 K (1 -α).
Replacing |||W ||| and W -2 K (1α)I m by their value in (4.36) leads to Proposition 4.5.2(iii), which completes the proof.

Robustness of iRestNet to an input perturbation

Nonexpansiveness can be used to measure the stability of an operator. For instance, let T : R m → R m be a nonexpansive operator and let x and δ x be in R m , where δ x can be seen as an input perturbation. From the definition of a nonexpansive operator (Section 2.2.1.4), we deduce that the distance between the perturbed and the original outputs is bounded by the norm of the perturbation; namely, ∆T (x) ≤ δ x , with ∆T (x) = T (x+δ x )-T (x).

An even stronger indicator of robustness is given by the notion of averageness. If we assume that T is an averaged operator, then we deduce from Proposition 2.2.18 that, as δ x tends to 0, ∆T (x) tends to δ x . Hence, the notion of averageness, in addition to providing a tighter bound on the norm of the ouput perturbation, also provides an insight regarding its direction.

Let us consider problem (4.34), where we assume additionally that H H and D D are diagonalizable in a same basis denoted P. The latter is satisfied for instance if H and D are the results of cyclic convolutive operators. Theorem 4.5.3 below gives sufficient conditions under which the proposed network applied to problem (4.34) is averaged.

β -= min 1≤p≤m K-1 k=0 β (p) k , β + = max 1≤p≤m K-1 k=0 β (p) k , (4.40)
and

θ k = k l=0 θ l-1 max 1≤q l ≤m β (q l ) k . . . β (q l ) l . (4.41)
In addition, if one of the following conditions is satisfied

(i) β + + β -≤ 0 and θ K-1 ≤ 2 K-1 (2α -1); (ii) 0 ≤ β + + β -≤ 2 K+1 (1 -α) and 2θ K-1 ≤ β + + β -+ 2 K (2α -1); (iii) 2 K+1 (1 -α) ≤ β + + β -and θ K-1 ≤ 2 K-1 , then the operator R K-1 •(W K-1 •+b K-1 )•• • ••R 0 •(W 0 •+b 0 ) is α-averaged.
Proof In this section, we present numerical experiments on a set of problems of image restoration, demonstrating that in many cases the proposed approach yields a better reconstruction quality than standard variational and machine learning methods.

Image deblurring

We consider the non-blind color image deblurring problem, whose degradation model reads

y = Hx + ω,
where n is the number of pixels, y = (y (j) ) 1≤j≤3 ∈ R 3n is the blurred RGB image, x = (x (j) ) 1≤j≤3 ∈ R 3n is the ground-truth, H ∈ R 3n×3n is a linear operator that models the circular convolution of a known blur kernel with each channel of the color image, and ω ∈ R 3n is a realization of an additive white Gaussian noise with standard deviation σ. An estimate of x can be derived from the following penalized formulation, which includes a smoothed total variation regularization, also called 2 -1 norm or edge-preserving regularization [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF],

minimize x∈C 1 2 Hx -y 2 + λ 3n i=1 (D v x) 2 i + (D h x) 2 i δ 2 + 1, (4.42)
where the feasible set C is the hypercube [x min , x max ] 3n , x min and x max are a lower and an upper bound on the pixel intensity, respectively, D v ∈ R 3n×3n and D h ∈ R 3n×3n are the vertical and horizontal gradient operators, respectively, δ ∈ R * + is a smoothing parameter and λ ∈ R + is the regularization parameter. Here, x min = 0, x max = 1 and we set δ = 0.01 in all experiments, which appears as an appropriate order of magnitude. To find this value for δ, we solved problem (4.42) for a small set of images of the database and used the simplex method to find the best values for δ and λ in terms of image quality. It is worth noting that the value for δ has not been fine-tuned, but that the proposed architecture could also be easily modified to include the inference of δ. The update A, defined in (4.32), is derived from (4.42), and is unfolded over K iterations, as it is described in Section 4.4. The bound constraints in problem (4.42) fall under the framework studied in Section 4.3.2, which provides the expression for the proximity operator of the barrier and its gradient.

Network characteristics

The tuning of the number of unfolded iterations K must achieve a compromise between training time, memory requirement, and performance. In order to determine a suitable setting for K, we trained networks with different numbers of layers and increased the number of layers until the performance of the network did not improve significantly. Using this procedure, the depth of iRestNet is taken equal to K = 40.

(∀k ∈ {0, . . . , K -1}) λ k = L (λ) k (x k ) = Softplus (b k ) σ(y) η(x k ) + Softplus (c k ) , (4.43)
where (b k , c k ) is a pair of scalars learned by the network, η(x k ) is the standard deviation of

[(D v x k ) (D h x k ) ],
which are the concatenated spatial gradients of x k , and σ(y) is an approximation of the noise level in the blurred image. The noise level is estimated as in [Mallat, 1999, Section 11.3.1]

σ(y) = median(|W H y|)/0.6745,
where |W H y| is the vector gathering the absolute value of the diagonal coefficients of the first level Haar wavelet decomposition of y. It is worth noticing that the proposed architecture does not require any prior knowledge about the noise level, in particular the noise standard deviation does not have to be the same for all input images. The architecture of the post-processing layer L pp is inspired from [START_REF] Zhang | Learning deep CNN denoiser prior for image restoration[END_REF], it is made of 9 convolutional layers with filters of size 3 × 3. The dilation factor changes from one layer to another, so as to widen the receptive field without creating memory issues. There is little correlation between the artifacts that remain in the image after going through the 40 blocks of iRestNet and the ground-truth image. Hence, it is easier for the network to learn the residual mapping instead of the image itself. Therefore, we add a skip connection between the input of L pp and its output. Finally, a ReLU activation function is used after each convolution, the final activation function is chosen as the Sigmoid function, and residual learning is combined with batch normalization, a technique which is widely used in deep learning to accelerate and stabilize the training process [START_REF] Zhang | Learning deep CNN denoiser prior for image restoration[END_REF]. The final architecture of L pp can be found in Figure 4.4.

Dataset and experimental settings

The training set is made of 1200 RGB images: 200 images stem from the Berkeley segmentation (BSD500) training set, while the remaining 1000 images are taken from the COCO training set. We use the BSD500 validation set, which is made of 100 images, to monitor the training and check if there is overfitting. The performance of the proposed method is evaluated on two different test sets: the BSD500 test set, which is made of 200 RGB images, and the Flickr30 test set used in [START_REF] Xu | Deep convolutional neural network for image deconvolution[END_REF], which is made of 30 RGB images. The test images have been center-cropped using a window of size 256 × 256. Blurry images are produced using the following 25 × 25 blur kernels and noise levels:

-A Gaussian kernel, which models atmospheric turbulence, with a standard deviation of 1.6 pixels, and a Gaussian noise standard deviation of σ = 0.008. This configuration is denoted as GaussianA. To evaluate the robustness of the proposed method with respect to the noise level, the same kernel is used with a Gaussian noise whose standard deviation is uniformly distributed between 0.01 and 0.05. The latter is denoted as GaussianB.

-The Gaussian kernel with a standard deviation of 3 pixels, and a Gaussian noise standard deviation of σ = 0.04, denoted as GaussianC.

-The eighth and third motion test kernels from [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF], which are real-world camera shake kernels, with a Gaussian noise standard deviation of σ = 0.01. These settings are denoted as MotionA and MotionB, respectively.

-The square uniform kernel of size 7×7, with a Gaussian noise standard deviation of σ = 0.01. This configuration is referred to as Square.

Training

For each degradation model, one iRestNet network is trained. We use a greedy approach for training the first 30 layers. For L 0 , a minibatch of 10 images is selected at every iteration, randomly cropped using a window of size 256 × 256, blurred with the given kernel, and degraded with Gaussian noise; the training of L 0 stops after a fixed number of epochs. Then, for each image of the training set, a random crop of size 256 × 256 is selected, blurred, corrupted with noise and passed through L 0 , the output is saved and used as an input to train L 1 . When the training of L 1 is complete, its output is used to train the next layer, etc... This training strategy is chosen with regards to its low memory requirement: the number of layers is not limited by the hardware. The rest of the network, L pp • L 39 • . . . • L 30 , is trained as one block and the learning rate is multiplied by 0.9 every 50 epochs. To accelerate the training, for every k ∈ {1, . . . , K -1}, the weights of L k are initialized with those of L k-1 . Detailed information about learning rates and number of epochs can be found in Tables 4.1 and 4.2 below.

The validation set is used to monitor this last step of the training. In particular, the configuration of network parameters that gives the best performance on the validation set during the training is the one saved and used for the tests. Note that for the first 30 layers, after each layer the quality of the restored training images should improve. This property comes from the training strategy, it is not encoded in the network: if memory was not an issue, then iRestNet could be trained in an end-to-end fashion.

We use the Adam optimizer [Kingma and Ba, 2014] to minimize the training loss, which is taken as the negative of the structural similarity measure (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]] defined below

SSIM(x, x) = (2µ x µ x + c 1 )(2σ x σ x + c 2 )(2cov xx + c 3 ) (µ 2 x + µ 2 x + c 1 )(σ 2 x + σ 2 x + c 2 )(σ x σ x + c 3 )
, where x is the ground truth, x is the restored image, (µ x , σ x ) and (µ x , σ x ) are mean and standard deviation of x and x, respectively, cov xx is the cross-covariance of x and x, and c 1 , c 2 and c 3 are constants. As explained in [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], the SSIM is a good measure of perceived visual quality, since it is based on how the human eye extracts structural information from an image. Hence, it is more discriminative with regards to artifacts than the mean square error for instance. The gradient of the SSIM loss with respect to the trainable parameters of the network is computed using a code available online1 and based on [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], the chain rule, automatic differentiation [START_REF] Paszke | Automatic differentiation in Py-Torch[END_REF], and the expression given in Section 4.3.2 for the derivatives of the barrier proximity operator. Codes are implemented in Pytorch. Some hidden layers in the postprocessing part make use of ReLU, which is not differentiable everywhere. Since this nondifferentiability happens only at specific points for which the left and right derivatives are well-defined, Pytorch can handle it as explained in [START_REF] Goodfellow | Deep Learning[END_REF]. All trainings are conducted using a GeForce GTX 1080 GPU or a Tesla V100 GPU. The training, which can be performed off-line, takes approximately 3 to 4 days for each blur kernel, while the time taken per test image is only about 1.4 sec on a GeForce GTX 1080 GPU.

Evaluation metrics and competitors

The restoration is evaluated in terms of the SSIM metric. The reconstruction given by the proposed approach is compared with a solution to problem (4.42) obtained using the projected gradient algorithm [START_REF] Iusem | On the convergence properties of the projected gradient method for convex optimization[END_REF]. For every blurred image, the pair (λ, δ) which leads to the best SSIM is selected using the simplex method. The solution given by this variational approach is referred to as VAR. The latter is an unrealistic scenario since it assumes that there is a perfect estimator of the error for selecting the hyperparameters, but it gives an upper bound on the image quality that one can expect by solving (4.42). We also use the following deep learning image restoration methods for comparison.

(i) EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] corresponds to a Bayesian approach where the authors divide images into patches and assign a prior to each patch. This results in a Gaussian mixture model whose parameters are learnt. The deblurred image is then a maximum a posteriori estimate of the image.

(ii) MLP refers to the Multi-Layer Perceptron network proposed in [START_REF] Schuler | A machine learning approach for non-blind image deconvolution[END_REF], which is fed with a pre-deconvolved image produced by a Wiener deconvolution filter. The network is trained for each blurring kernel.

Finally, we include comparisons with three unfolded-based methods.

(iii) In IRCNN [START_REF] Zhang | Learning deep CNN denoiser prior for image restoration[END_REF], an empirical algorithm derived from an augmented Lagrangian formulation is unfolded over 30 iterations 2 , and a CNN is used as a denoiser to update the splitting variable.

(iv) In FCNN [Zhang et al., 2017a], the authors unfold the same algorithm as in the previous reference, and use a network to learn an effective regularization function.

(v) The method from [START_REF] Meinhardt | Learning proximal operators: Using denoising networks for regularizing inverse imaging problems[END_REF] is referred to as PDHG.

The authors perform a maximum of 30 iterations of a primal dual hybrid gradient algorithm, and the proximity operator of the second regularization function is replaced by a neural network.

For FCNN, we use the code that is available online, in which the authors provide a model that has only been trained for motion blurs. Hence, for a fair comparison, we only provide the results of FCNN on MotionA and MotionB, and we specify that this method is not applicable (n/a) to the other configurations. Similarly, for MLP and PDHG, the authors do not provide models that were trained specifically for MotionB and Square, so we do not test these methods on these two configurations.

Since MLP, EPLL and IRCNN require the knowledge of the noise level, for the GaussianB degradation model, we make use of the estimation of the noise standard deviation given by the method in [Mallat, 1999, Section 11.3.1]. In addition, since some comparison methods, like EPLL for instance, do not estimate well the borders of the images, the SSIM index is computed excluding a 6-pixel-wide frame for all images and all tested methods.

Results and discussion

The average SSIM obtained with the different methods for the various blur kernels and noise levels on the BSD500 test set can be found in Table 4.3. The mean SSIM achieved with iRestNet on this test set is greater than those obtained with the other methods for all degradation models except MotionA. For this kernel, the average SSIM achieved with iRestNet is the second highest value after IRCNN, which appears as the most competitive method. IRCNN involves two steps: first, a Wiener filter is applied to the blurred image, then, a neural network is used to predict the residual and denoise the image. These two steps are repeated 30 times, for 30 different manually tuned regularization parameters. In contrast, iRestNet does not require any tuning from the user regarding the regularization parameters during training. For completeness, the SSIM of all images of the BSD500 test set are plotted in Figure 4.5 for the 6 different degradation models. As one can see, iRestNet performs well in terms of SSIM on most of the images.

Since no image was taken from Flickr for training iRestNet, the results on the Flickr30 test set show how well the performance of the trained networks are transferable on test sets with statistics that are different from those of Examples of visual results obtained with the different methods can be found in Figures 4.6 and 4.7 for two images from the BSD500 test set and the blur kernels GaussianB and Square, respectively. We also provide the results obtained for one image from the Flickr30 test set that has been degraded with MotionB. As one can see from inspecting these pictures, details from the snake's and caterpillar's skin patterns are better retrieved with iRestNet, which provides more visually-satisfactory results than competitors. Similarly, on Figure 4.8, competitors tend to smooth too much the details on the leaves as it can be seen in the top left-hand corner. Regarding Figure 4.6, which belongs to the test set with a level-varying noise, it is worth noting that, on the result obtained with the proposed method, the green background is free from artifacts, which is not the case for the other methods, in particular for PDHG and IRCNN. This suggests that those two competitors are not robust to a small change in the noise level. interior point algorithm. It can be noted that the proposed approach can be extended to a set of regularization functions, or to penalizations which are parametrized by several variables. Useful constraints on the sought solution can be enforced thanks to a logarithmic barrier, so providing more control over the output of the network. We have shown for three standard types of constraints that the involved proximity operator can easily be computed, and that its derivatives are well-defined and computable. In the case of a quadratic cost function, the theoretical result of Section 4.5 regarding the robustness of the network with respect to an input perturbation, ensures the reliability of the proposed method, which is crucial for many applications.

It would be interesting to extend the scope of this study to a wider class of 4.7. Summary 109 problems, and to illustrate this stability result by numerical experiments on different applications like classification. As demonstrated by our experiments in image restoration, iRestNet performs favorably compared to state-of-theart variational and machine learning methods. An advantage of the proposed approach is that, in contrast with its evaluated competitors, it does not require any knowledge about the noise level and it does not involve any hand-selection of the regularization parameters. One limitation of iRestNet is that the network needs to be trained for a given blur kernel. A direction for future works is to extend the method to situations in which the observation model is not fully known, so as to address blind or semi-blind deconvolution problems.

The unfolded approach proposed in this chapter overcomes a limitation of classical variational strategies, which is the setting of their hyperparameters. Another alternative to address this issue is to use the Bayesian framework, which leads to useful stochastic simulation methods. In addition, similarly to iRestNet whose training is based on a better indicator of the error than the original objective function, Bayesian approaches allow to derive insightful estimators other than the MAP, like the MMSE for instance (please refer to Section 2.1.3 for more details). In the next chapter, we propose to incorporate a proximal algorithm in a Bayesian approach.

-Chapter 5 -A preconditioned proximal sampling algorithm with applications in ultrasound imaging Ultrasound (US) imaging is a non-invasive procedure which is widely used in medical applications, such as obstetrics or cardiology. Although the quasi absense of side-effects caused by this procedure, and its low cost and aquisition time, make it the norm for numerous diagnostic exams [START_REF] Szabo | Diagnostic ultrasound imaging: inside out[END_REF], ultrasonic measurements suffer from some downsides, including a low signalto-noise ratio, attenutation, speckle, shadows and a limited resolution [START_REF] Noble | Ultrasound image segmentation: a survey[END_REF][START_REF] Ploquin | Resolution enhancement in medical ultrasound imaging[END_REF]. The echo signals recorded during an ultrasound exam are processed to produce radio-frequency (RF) images. While image processing tasks are commonly performed directly on the resulting RF image [START_REF] Nair | Coronary plaque classification with intravascular ultrasound radiofrequency data analysis[END_REF], or on its B-mode [START_REF] Noble | Ultrasound image segmentation: a survey[END_REF], which is a log-compression of its envelope, useful information can also be drawn from the statistics of the scatterer field, also called tissue reflectivity function (TRF) [START_REF] Jensen | A model for the propagation and scattering of ultrasound in tissue[END_REF][START_REF] Ng | Modeling ultrasound imaging as a linear, shift-variant system[END_REF]. Since the TRF models well tissue inhomogeneity, it can be used to perform efficiently segmentation [START_REF] Pereyra | Segmentation of skin lesions in 2D and 3D ultrasound images using a spatially coherent generalized Rayleigh mixture model[END_REF], tissue characterization [START_REF] Bernard | Statistics of the radio-frequency signal based on K distribution with application to echocardiography[END_REF], or classification [START_REF] Alessandrini | A restoration framework for ultrasonic tissue characterization[END_REF]. Estimating the TRF given the RF image is a challenging inverse problem, which we are going to tackle jointly with the segmentation task in this chapter.

Models that are studied in US imaging are derived from the physical laws governing the involved aquisition process. Such models may be significantly complex depending of the task to be performed. The Bayesian framework and stochastic simulation methods, like the ones described in Section 2.1.3, have proven to be very useful for solving inverse problems arising in US imaging [START_REF] Comer | The EM/MPM algorithm for segmentation of textured images: analysis and further experimental results[END_REF][START_REF] Husby | Estimating blood vessel areas in ultrasound images using a deformable template model[END_REF][START_REF] Dobigeon | Regularized Bayesian compressed sensing in ultrasound imaging[END_REF]]. In the context of Bayesian theory, we propose in this chapter a preconditioned proximal sampling algorithm for the TRF, which is incorporated in a hybrid Gibbs sampler with the aim of performing joint deconvolution and

111 every x ∈ R m (∀i ∈ {1, . . . , m}) (b(x)) i = 1 2 m j=1 (A(x)) i,j ∂ log π(x) ∂x j +det(A(x)) 1 2 m k=1 ∂ ∂x j (A(x)) i,k det(A(x)) -1 2 .
In the above definition, (∀x ∈ R m ) A(x) = V (x)V (x) ∈ R m×m is a symmetric positive definite matrix, and det(A(x)) denotes its determinant. In addition, π is the density of the stationary distribution of the diffusion. Euler's discretization scheme can be used on (5.2) to generate a Langevin Markov chain (x (t) ) t∈N as follows

(∀t ∈ N) x (t+1) = x (t) + 2γb(x (t) ) + 2γ(A(x (t) )) 1 2 ω (t) .
Hereabove, (∀t ∈ N) ω (t) is generated using a zero-mean Gaussian distribution N (0 m , I m ), and γ ∈ R * + is the discretization stepsize that controls the length of the jumps, while the scale matrix A(•) drives their direction. Instead of letting A(•) be the identity operator as in the standard Metropolis adjusted Langevin algorithm [START_REF] Roberts | Exponential convergence of Langevin distributions and their discrete approximations[END_REF], we follow [START_REF] Stuart | Conditional path sampling of SDEs and the Langevin MCMC method[END_REF][START_REF] Marnissi | Majorize-minimize adapted Metropolis-Hastings algorithm[END_REF] and use a preconditioning matrix A ∈ S + m to accelerate the Langevin scheme, which leads to

(∀t ∈ N) x (t+1) = x (t) + γA∇ log π(x) + 2γA 1 2 ω (t+1) .
(5.3)

Approximation of the target diffusion

Let f : R m →] -∞, +∞] and let g ∈ Γ 0 (R m ) be differentiable on R m with a Lipschitz-continuous gradient. Assume that the target distribution π satisfies the following relation,

(∀x ∈ R m ) π(x) ∝ exp(-(f + g)(x)).
Let γ ∈ R * + and M ∈ S + m . Following [START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF], we replace π by its Moreau approximation π M γ defined by

(∀x ∈ R m ) π M γ (x) = sup u∈R m π(u) exp - u -x 2 M 2γ .
Note that we dropped the normalization constant and that, for convergence acceleration purposes, we have introduced the preconditioning matrix M , which was not included in [START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF]. When M is not specified, the identity matrix is used, i.e. M = I m . Hence, the approximated version of (5.3) reads

(∀t ∈ N) x (t+1) = x (t) + γA∇ log π M γ (x) + 2γA 1 2 ω (t+1) . (5.4)
We can then deduce the following result when f + g is convex.

∈ R * + , M ∈ S + m and x ∈ R m , if (f + g) ∈ Γ 0 (R m ), then we have ∇ log π M γ (x) = M prox M γ(f +g) (x) -x γ .
Proof. By definition of π M γ , we have for all

x ∈ R m log π M γ (x) = -(f + g) prox M γ(f +g) (x) - 1 2γ prox M γ(f +g) (x) -x 2 M .
Hence, applying [Combettes and Wajs, 2005, Lemma 2.5] in the metric induced by M directly leads to the result.

From Proposition 5.3.1, (5.4) becomes

(∀t ∈ N) x (t+1) = x (t) + γAM prox M γ(f +g) (x (t) ) -x (t) γ + 2γA 1 2 ω (t+1) .
Finally, we set A = M -1 , which leads to

(∀t ∈ N) x (t+1) = prox M γ(f +g) (x (t) ) + 2γM -1 2 ω (t+1) . (5.5)
Next, we introduce a splitting strategy in (5.5).

Forward-backward approximation

By assumption, g is differentiable on R m and its gradient ∇g is Lipschitzcontinuous. It is worth noting that the computation of the proximity operator of the sum of two functions is generally intractable [START_REF] Pustelnik | Proximity operator of a sum of functions; application to depth map estimation[END_REF]. Hence, as suggested in [START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF], we use a first-order Taylor expansion to approximate the proximity operator of f + g and introduce a forward step in the proposed PP-ULA iteration. Let o denotes Landau's notation. 2Let x ∈ R m and u ∈ R m . Using

g(u) = g(x) + (u -x) ∇g(x) + o( u -x ),
we have

(f + g)(u) + 1 2γ u -x 2 M = f (u) + g(x) + 1 2γ u -x 2 M + (u -x) ∇g(x) + o( u -x ),
which can be re-written as

(f + g)(u) + 1 2γ u -x 2 M = f (u) + g(x) + 1 2γ u -x + γM -1 ∇g(x) 2 M - γ 2 M -1 2 ∇g(x) 2 + o( u -x ).
Hence, the proximity operator of f + g can be expressed as follows,

prox M γ(f +g) (x) = argmin u∈R m (f + g)(u) + 1 2γ u -x 2 M = argmin u∈R m f (u) + 1 2γ u -x + γM -1 ∇g(x) 2 M +o( u -x ) .
In addition, we have

prox M γf (x -γM -1 ∇g(x)) = argmin u∈R m f (u) + 1 2γ u -x + γM -1 ∇g(x) 2 M .
Therefore, when γ is small, prox M γf (x -γM -1 ∇g(x)) is a good approximation of prox M γ(f +g) (x). Plugging this preconditioned forward-backward scheme [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]] in (5.5) leads to the proposed sampling method PP-ULA,

(∀t ∈ N) x (t+1) = prox M γf (x (t) -γM -1 ∇g(x (t) )) + 2γM -1 2 ω (t+1) . (5.6)
Now that we have presented our proposed sampling strategy PP-ULA, we describe in the next section the Bayesian model to which it is going to be applied. § 5.4 Bayesian model

In this section, we describe the hierarchical Bayesian model attached to the variables introduced in Section 5.1. This model follows the same strategy as in [START_REF] Zhao | Joint segmentation and deconvolution of ultrasound images using a hierarchical Bayesian model based on generalized Gaussian priors[END_REF] and is illustrated in Figure 5.1.

Priors

The likelihood corresponding to the RF image y, which follows from the Gaussian nature of the noise, is derived from (5.1) and expressed below. The TRF is a mixture of GGDs which, under the assumption that the pixel values are independent given the segmentation labels z, leads to the following prior

p(y|x, σ 2 ) = 1 (2πσ 2 ) m/2 exp - y -Hx 2 2σ 2
p(x|α, β, z) = m i=1 1 2β 1/αz i z i Γ(1 + 1/α z i ) exp - |x i | αz i β z i , (5.7)
where Γ is the Gamma function. For every set C, let ı C denote the indicator function of this set. Uninformative Jeffreys priors are assigned to the noise variance σ 2 and scale parameters β, while the shape parameters α are assumed to be uniformly distributed between 0 and 3. These assumptions lead to the following priors for these variables,

p(σ 2 ) = 1 σ 2 ı R + (σ 2 ), p(α) = K k=1 1 3 ı [0,3] (α k ), and p(β) = K k=1 1 β k ı R + (β k ).
The segmentation labels z are modeled by a Potts Markov random field with prior

p(z) = 1 C(θ) exp   m i=1 j∈V(i) θδ(z i -z j )   , (5.8)
with δ the Kronecker function, C(θ) ∈ R * + a normalizing constant, θ ∈ R * + a granularity coefficient, and V(i) the set of four closest neighbours of the ith pixel. The different variables are sampled according to their conditional distributions, which are provided in the next section.

Conditional distributions

The conditional distribution of the noise variance is derived from the Bayes theorem and can be expressed as follows,

p(σ 2 |y, x) ∝ IG m 2 , y -Hx 2 2 , (5.9)
where IG denotes the inverse gamma distribution. Assuming that the different regions have independent shape and scale parameters, for every region k ∈ {1, . . . , K}, we obtain

p(α k |x, β, z) ∝ i∈I k ı [0,3] (α k ) 2β 1/α k k Γ (1 + 1/α k ) exp - |x i | α k β k , (5.10) and p
(β k |x, α, z) ∝ IG   m k α k , i∈I k |x i | α k   , (5.11)
with I k = {i ∈ {1, . . . , m} | z i = k} and m k the number of elements in I k .

For every pixel i ∈ {1, . . . , m} and every region k ∈ {1, . . . , K}, the Bayes rule applied to the segmentation labels leads to the following probability distribution,

p(z i = k|x, α, β, z V(i) ) ∝ exp j∈V(i) θδ(z j -k) -|x i | α k β k 2β 1/α k k Γ(1 + 1/α k ) , (5.12)
where z V(i) denotes the vector made of the labels in the neighborhood of the ith pixel. Finally, the conditional distribution of the TRF is equal to

π(x) = p(x|y, σ 2 , α, β, z) ∝ exp (-f (x) -g(x)) , (5.13)
where, for every x ∈ R m ,

f (x) = m i=1 |x i | β z i αz i and g(x) = y -Hx 2 2σ 2 .
(5.14) § 5.5 Hybrid Gibbs sampler

We can now describe the proposed hybrid Gibbs sampler used to perform joint deconvolution and segmentation of ultrasound images. A special attention is devoted to the sampling scheme of the TRF, which is our main contribution.

Sampling from the conditional distributions

Samples for the noise variance σ 2 and scale parameters β can easily be drawn from their respective conditional distribution (5.9) and (5.11). Regarding the shape parameters α, we use the Metropolis-Hastings random walk Algorithm 1, presented in Section 2.1.3.2, to generate samples according to distribution (5.10). For the segmentation labels z, we first suitably 5.5. Hybrid Gibbs sampler 119 normalize the distributions given by (5.12). The resulting normalized conditional probabilities are given below

(∀i ∈ {1, . . . , m})(∀k ∈ {1, . . . , K}) pi,k = p(z i = k|x, α, β, z V(i) ) K l=1 p(z i = l|x, α, β, z V(i) )
.

(5.15) Hence, for every pixel i ∈ {1, . . . , m}, label z i is drawn from {1, . . . , K} using the probabilities (p i,k ) 1≤k≤K defined above. Finally, we propose to apply our PP-ULA scheme (5.6) to the conditional distribution (5.13) in order to generate samples for the TRF x. This leads to the following scheme,

(∀t ∈ N) x (t+1) = prox M γf ( xt ) + 2γM -1 2 ω (t+1) , (5.16)
where f is defined as in (5.14), M ∈ S + m is a preconditioning matrix used for acceleration, ω (t+1) ∼ N (0 m , I m ), and

x(t) = x (t) - γ σ 2 M -1 H (Hx (t) -y).
(5.17)

It can be noted that, in Proposition 5.3.1, f is assumed to be convex, which is not necessarily satisfied in our case. However, for simplicity, we take the discrete scheme (5.16) even in the nonconvex case. Since the proposed sampling strategy is unadjusted, (5.16) is not followed by an acceptance test.

The bias with respect to π increases with γ, as the speed of convergence of the algorithm. A compromise must be found when setting γ.

As described in the next section, when there exist shape parameters that are below 1, we use the MM principle [START_REF] Schifano | Majorizationminimization algorithms for nonsmoothly penalized objective functions[END_REF] introduced in Section 2.2.4 to replace the nonconvex minimization problem involved in the computation of prox M γf with a sequence of convex surrogate problems.

MM principle for computing the proximity operator

The function f defined in (5.14) can be written as the sum of a convex term and a nonconvex term. We propose to approximate the nonconvex part using the MM principle.

Let σ ∈ R * + , α ∈ [0, 3] K , β ∈ (R * + ) K and z ∈ {1, . . . , K} m . We define the sets I <1 = {i ∈ {1, . . . , m} | α z i < 1} and I ≥1 = {1, . . . , m} \ I <1 . Let J ⊂ I <1 . We define h J at every (u, v) ∈ R m × (R * + ) m by h J (u, v) = i∈I ≥1 |u i | αz i β z i + j∈J (1 -α z j )v αz j j + α z j v αz j -1 j |u j | β z j .
From concavity, we deduce that, for every v ∈ (R * + ) m and u ∈ R m such that J ⊂ {i ∈ I <1 | |u i | > 0}, the following majoration property holds The dual forward-backward Algorithm 5 is then combined with Remark 2.2.8 in order to compute the proximity operator of the convex function h J (•, v) in any metric. This leads to Algorithm 15 below, which generates a sequence (s (p) ) p∈N converging to prox M h J (•,v) (x) for any x ∈ R m and M ∈ S + m .

h J (u, v) ≥ i∈I ≥1 ∪J |u i | αz i β z i = h J (u, (|u i |) 1≤i≤m ). Since h J (•, v) is convex
Algorithm 15: DFB algorithm to compute prox

M h J (•,v) (x) Set ρ = 1/ M -1 , ∈]0, min{1, ρ}[, η ∈ [ , 2ρ -]; Initialization: dual variable w (1) ∈ R m ; for p = 1, ... do s (p) = x -M -1 w (p) ; w (p+1) = w (p) + ηs (p) -ηprox η -1 h J (•,v) (η -1 w (p) + s (p) ) end
Finally, Algorithm 16 generates a sequence (u (q) ) q∈N estimating prox M γf ( x(t) ).

Algorithm 16: MM principle to compute prox M γf . Initialize u (1) ∈ R m ; for q = 1, ... do

J (q) = {i ∈ I <1 | |u (q) i | > 0}; v (q) = (|u (q) i |) 1≤i≤m ; u (q+1) = prox M γh J (q) (•,v (q)
) ( x(t) ) (using Algorithm 15) end

Proposed sampler

The resulting hybrid Gibbs sampler is summarized in Algorithm 17. It is evaluated on six different experiments in the next section. obtained with a Wiener filter, where the noise level has been estimated as in [START_REF] Mallat | A wavelet tour of signal processing[END_REF], and with the solution to the Lasso problem, where the regularization weight is set i ) manually when the ground-truth is not available, or ii ) using a golden-section search to maximize the peak signalto-noise ratio (PSNR) defined as PSNR = 10 log 10 m max 1≤i≤m (x i , x * i ) 2 xx * 2 , (5.18) with x ∈ R m the true TRF and x * ∈ R m the estimated one. We also compare our results with the segmentation given by Otsu's method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF] applied to the Wiener-deconvolved image, and with the SLaT method [START_REF] Cai | A three-stage approach for segmenting degraded color images: Smoothing, lifting and thresholding (SLaT)[END_REF] applied to the Lasso-deconvolved image. The proposed method, PP-ULA, is used with γ = 0.09 and M an approximation of the Hessian of the differentiable term in (5.13) [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF], M = (H H + λI m )/σ 2 , with λ = 0.1 so that the inverse of M is well-defined.

In order to evaluate the deconvolution results, we compute the PSNR and the structural similarity measure (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] between the restored TRF and its ground-truth. The contrast in the restored TRF is assessed according to the contrast-to-noise ratio (CNR) [START_REF] Krishnan | Improved estimation of phase aberration profiles[END_REF] between two windows from different regions of the B-mode TRF images. The CNR is defined as follows, CNR = |µ 1µ 2 |/(ν 1 + ν 2 ) 1/2 , where (µ 1 , µ 2 ) and (ν 1 , ν 2 ) are, respectively, the means and variances of the two windows. The segmentation is evaluated according to the percentage of correctly predicted labels, which is referred to as overall accuracy (OA). The minimum mean square error (MMSE) estimators of all parameters in HMC and PP-ULA are computed after the burn-in regime. Moreover, to evaluate the mixing property of the Markov chain after convergence, we compute the mean square jump (MSJ) per second, which is the ratio of the MSJ to the time per iteration. The MSJ is obtained using T samples of the TRF (x (t 0 +1) , . . . , x (t 0 +T ) ) generated after the burn-in period, i.e.

MSJ = 1 T -1 T -1 t=1 x (t 0 +t) -x (t 0 +t+1) 2 1/2
. First, we present the results obtained on the proposed simulated data.

Results on simulated data

The convergence speed of Algorithm 17 is empirically observed for Simu1 and Simu2, as illustrated in Figure 5.3, where we also display the results of the non-preconditioned P-ULA, for which M = I m and γ = 1.99σ Comparing P-ULA and PP-ULA on these simulated data allows us to study the effect of adding a preconditioner in the proposed sampling scheme. As reported in Table 5.2, P-ULA needs more iterations and more time to converge than PP-ULA: the proposed method is 12.2 and 4.8 times faster than P-ULA on Simu1 and Simu2, respectively. In addition, from Table 5.3 and Figure 5.4, we deduce that P-ULA is more biased than PP-ULA, which samples correctly the target distributions. Finally, as one can see in Figure 5.3 and Table 5.4, P-ULA leads to lower PSNR, SSIM and OA values than PP-ULA. These results clearly emphasize the benefits of preconditioning in this example.

From Table 5.2, PP-ULA is 5.7 and 6.6 times faster than HMC on Simu1 and Simu2 and has better mixing properties, as shown by the MSJ per second. Visual results from Figures 5.5 and 5.6, and CNR values in Table 5. 4 show that the contrast obtained with PP-ULA is better than with competitors on Simu2, and is second best after P-ULA on Simu1. However, it should be noted that the PSNR and SSIM obtained on Simu1 with P-ULA are much lower than with the other methods. In addition, the PSNR and

Iterations

Results on a tissue-mimicking phantom and on real data

The convergence of Algorithm 17 is also empirically observed for the experiments on the tissue-mimicking phantom and on real data, i.e. Kidney, Thyroid, Bladder and KidneyReal. As mentioned in Table 5.5, the proposed method leads to a significant acceleration since it is between 3.7 and 6.3 times faster than HMC on these experiments. Visual results from Figures 5.9 and 5.10, and CNR values in Table 5.6 show that the contrast obtained with PP-ULA is better than with competitors on all these test images. In addition, the PSNR and SSIM values from Although the ground-truth of the segmentation is not available for these experiments, one can see from the visual segmentation results shown in In this chapter, we have investigated a new method based on a preconditioned proximal unadjusted Langevin algorithm for the joint restoration and segmentation of ultrasound images, which showed faster convergence than an existing Hamiltonian Monte Carlo algorithm. Hence, the proposed method has the potential to speed-up the approach proposed in [START_REF] Pereyra | Segmentation of skin lesions in 2D and 3D ultrasound images using a spatially coherent generalized Rayleigh mixture model[END_REF] for the segmentation of ultrasound images. Another direction for future work is to extend this framework to a spatially variant, possibly unknown, PSF. -Chapter 6 -Conclusion § 6.1 Summary

During this thesis, our main axes of research aimed at studying the use of the proximity operator in different contexts, and at building bridges between distinct methods through optimization. Accordingly, we have investigated three new methods for solving challenging inverse problems.

First, we have focused on the traditional variational strategy, which consists in finding a solution to an inverse problem by minimizing an objective function. Interior point methods are a specific class of variational methods that are well suited for constrained optimization problems. Despite their impressive performance in linear and quadratic programming, they progressively lost popularity to give rise to more flexible techniques, such as proximal algorithms.

In Chapter 3, we have combined the logarithmic barrier method with the proximal variable metric forward-backward algorithm, leading to the proposed PIPA algorithm. In the case of composite objective functions, PIPA benefits from a splitting strategy. In addition, it includes variable metrics, which can be used to boost its convergence as proven in our numerical experiments. Last but not least, the logarithmic barrier in PIPA ensures that every iterate satisfies strictly the constraints, which might be critical if, for instance, these constraints model a safety margin or if they have a physical meaning.

We have then carried out a mathematical analysis of the proposed algorithm. The convergence of PIPA was obtained under mild boundedness conditions on the variable metrics, and we have shown that, under some assumptions, the inner loop in PIPA converges linearly. Finally, we have evaluated the performance of PIPA on two applications in image processing. The first example concerns hyperspectral unmixing and the second one 134 Chapter 6. Conclusion addresses joint geometry-texture decomposition and reconstruction of CT data. In both experiments, we have formulated a constrained optimization problem, with a nondifferentiable term in the objective function, for solving the corresponding inverse problem. PIPA compares favorably with widelyused state-of-the-art algorithms on these two applications.

Second, we have considered the framework offered by recent deep learning methods for solving inverse problems. Neural networks have reached outstanding performances in applications related to computer vision. They may however suffer from a lack of interpretability and mathematical guarantees.

In Chapter 4, we have applied deep-unfolding to a proximal interior point algorithm. This algorithm has been unfolded over a fixed number of iterations identified with the network layers. As opposed to other unfolded methods, we have maintained the link between the original algorithm and the resulting architecture by keeping the proximal update explicit, so that only hyperparameters are inferred by the network.

In order to train the proposed network, referred to as iRestNet, using gradient descent and backpropagation, we have derived explicit expressions for the proximity operator of the barrier and for its first derivatives for three common types of constraints. It is worth noting that the logarithmic barrier allows to incorporate prior knowledge or desired properties about the solution within iRestNet. Moreover, when the optimization problem is quadratic, we have derived sufficient conditions under which iRestNet is robust with regards to an input perturbation. Such stability property is very important to promote the use of deep neural networks in fields involving high responsibility, such as medical imaging.

Lastly, we have demonstrated the good performance of iRestNet on image deblurring for several types of blurring kernels and noise levels. It is worth noting that iRestNet compares favorably in terms of image restoration with the classical variational approach and state-of-the-art neural network approaches for image deblurring, including architectures derived from other types of unfolded algorithms. Third, we have combined proximal tools with a Bayesian approach. The Bayesian framework provides useful stochastic simulation methods, that can be used as alternatives to variational methods which may require setting hyperparameters. In addition, Bayesian methods can help derive useful estimators, other than the MAP, such as the MMSE. However, Bayesian techniques can suffer from high computational costs.

In Chapter 5, we have introduced a new accelerated proximal unadjusted Langevin algorithm called PP-ULA. The preconditioning matrices and splitting scheme used in PP-ULA constitute two original features of this algorithm. Since PP-ULA is an unadjusted sampling method, it is not followed by an acceptance test. Then, we have considered a hybrid Gibbs sampler which has been derived for a specific application, the deconvolution and segmentation of ultrasound images. We have incorporated PP-ULA into this Gibbs sampler so as to accelerate the method. Finally, we have conducted tests on six different images, half of which are synthetic images, and the remaining are real in vivo ultrasound data. Using the proposed PP-ULA has led to a significant reduction of the computational time compared to an Hamiltonian Monte Carlo algorithm known for its good performance. In addition, the proposed method also produced better deconvolution and segmentation results than several variational approaches used as comparisons.

In the next section, we propose several extensions of the aforementioned methods that could be investigated for future works. § 6.2 Perspectives

The proximal interior point algorithm studied in Chapter 3 could be further improved as suggested hereafter.

Study the convergence of PIPA for nonconvex functions One interesting lead to extend the scope of application of PIPA would be to conduct its convergence analysis in a nonconvex setting. This could be done for instance by better relying on the KL property, as in the following works [Attouch and Bolte, 2009;[START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Lojasiewicz inequality[END_REF]Chouzenoux et al., 2014b].

Consider inexact proximity operators in PIPA In our numerical experiments from Sections 3.6 and 3.7, we use an iterative scheme to compute the proximity operator in the variable metric. Despite this inexact computation, we still observe the practical convergence of PIPA in both applications. Therefore, another improvement in the analysis of this algorithm would be to prove mathematically its convergence assuming that the proximity operator is computed in an inexact manner. Inexact proximity operators have been addressed for example in [START_REF] Salzo | Inexact and accelerated proximal point algorithms[END_REF], where the authors consider several types of errors. It is worth noting that, depending on the error type, the sequence of errors does not necessarily need to be summable for convergence to hold.

Derive accelerated versions of PIPA Based on the seminal work [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF], accelerated versions of proximal algorithms have been recently proposed, for example in [START_REF] Calatroni | Backtracking strategies for accelerated descent methods with smooth composite objectives[END_REF].

It would be interesting to integrate these acceleration techniques into PIPA and to perform numerical experiments to test if the resulting algorithm performs better than the simple version of PIPA.

Conduct experiments with nonlinear constraints It should be noted that the convergence results presented in Section 3.4 do not require the constraints to be linear. However, in the numerical experiments presented in Sections 3.6 and 3.7, constraints are taken linear. Hence, it would be insightful to use PIPA for an application that includes nonlinear constraints, like quadratic constraints for instance. This new setting could influence the choice of the variable metrics since it would modify the Hessian of the logarithmic barrier.

Building on the neural network architecture presented in Chapter 4, we propose the following research directions related to deep learning methods.

Extend the stability result of iRestNet The result obtained in Section 4.5 regarding the robustness of iRestNet holds under the assumption that the objective function is quadratic. Providing mathematical guarantees regarding the behavior of a neural network under an input perturbation is a prerequisite in many real-world applications. Therefore, we could consider studying the robustness of iRestNet in a more general setting so as to widen the scope of this method. In addition, this result could be improved by considering the hidden and post-processing layers in the study.

Improve the architecture of the hidden layers The proposed neural network makes use of very simple hidden layers to predict the hyperparameters involved in the unfolded algorithm, namely the stepsize, the barrier parameter, and the regularization coefficient. Using more sophisticated and deeper hidden layers could improve the restoration results. For instance, we could make the stepsize depend on the input image, as it is done for the barrier parameter. Moreover, the regularization parameter could depend on additional features in the image such as the mean, the variance, or second order spatial derivatives. We could also treat the parameter δ used to smooth the total variation as a learnable parameter.

Perform experiments with more complex models The numerical examples presented in Chapter 4 could be supplemented by additional experiments with different noise models. For instance, we could consider a Poisson noise or a Poisson-Gaussian mixture noise, which would lead to a different data fidelity function f as detailed in [START_REF] Chouzenoux | A convex approach for image restoration with exact Poisson-Gaussian likelihood[END_REF]. Furthermore, as opposed to classical variational methods, iRestNet is not limited by the number of hyperparameters. Hence, considering the local total variation used in Section 4.6, it would be interesting to introduce one regularization weight per pixel, or per group of pixels, and to measure the influence of this modification on the quality of the deblurred image. More sophisticated regularization functions could also be used, such as nonlocal regularization functions [START_REF] Peyré | Non-local regularization of inverse problems[END_REF] Extend iRestNet to blind deconvolution One limitation of the proposed architecture is that we need to train one network per blurring kernel. It would be interesting to include the inference of the observation operator H in iRestNet for the purposes of bringing more flexibility to the method and reducing the training time.

Finally, we suggest the following improvements and future leads related to our study of a Bayesian approach in Chapter 5.

Include the estimation of the granularity The hierarchical Bayesian model introduced in Section 5.4 includes a granularity coefficient in the segmentation Potts model. The segmentation result highy depends on the value of this parameter. In our numerical experiment we set it in a semi-automatic way that requires launching several Markov chains.

The proposed method could be improved by considering the granualrity as a random variable and estimating it jointly with the other variables in the Gibbs sampler. This could however slow down the simulation.

Extend the Bayesian approach to nonblind deconvolution For real in vivo ultrasound data, the PSF of the probe is not known. In our experiments, we have estimated the PSF from the RF image of a thin wire. It would be interesting to include the estimation of the PSF in the proposed model, so as to improve the performance of the method on real data.

Study the convergence of PP-ULA Studying the convergence of the full hybrid Gibbs sampler is very challenging since PP-ULA is an unadjusted algorithm. Nonetheless, we could start by studying the convergence of PP-ULA alone by building on the results derived for the classical proximal unadjusted Langevin algorithm in [START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF].

Investigate parallelization Although PP-ULA leads to a significant reduction of the time needed to perform joint deconvolution and segmentation of US images, the computational cost remains too long for online applications. We have implemented our method using Matlab. The proposed approach would greatly benefit from an implementation in a language with better performance. In addition, we have not made use of any parallelization or GPU implementation, which could be promising research directions for improving the efficiency of PP-ULA.

  ,j : element in the ith row and jth column of M |s| : absolute value of s •, • : Euclidean inner product v : 2 norm of v v M : norm induced by M ∈ S + m and equal to v, M v 1 2 |||M ||| : spectral norm of M Diag(v) : diagonal matrix whose elements are given by v ∇f (v) : gradient vector of the function f at v ∇ 2 f (v) : Hessian matrix of f at v ı C : indicator function of the set C dist(v, C) : inf w∈C vw , Euclidean distance between v and the set C N (µ, σ 2 ) : Gaussian distribution with mean µ ∈ R and variance σ 2 ∈ R * + U(a, b) : uniform distribution in [a, b] xi Acronyms ADMM : alternating direction method of multipliers A-PDS : accelerated primal-dual splitting CNN : convolutional neural network CNR : contrast-to-noise ratio CT : computed tomography DFB : dual forward-backward DNN : deep neural network FB : forward-backward GFBS : generalized forward-backward splitting HMC : Hamiltonian Monte Carlo IPM : interior point method KL : Kurdyka-Lojasiewicz l.s.c. : lower semicontinuous MAP : maximum a posteriori MCMC : Markov Chain Monte Carlo MH : Metropolis-Hasting MM : majorize-minimize MMSE : minimum mean squared error OA : overall accuracy PDS : primal-dual splitting PSF : point-spread function PSNR : peak signal-to-noise ratio P-ULA : proximal unadjusted Langevin algorithm ReLU : rectified linear unit RF : radio-frequency SSIM : structural similarity measure TRF : tissue reflectivity function US . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Collaborations . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 General notation . . . . . . . . . . . . . . . . . . . . . . . . . -Chapter 1 -General introduction § 1.1 Context

Definition 2 . 2 . 1

 221 Let f : R m →] -∞, +∞].

Theorem 2 . 2 . 3 (

 223 Fermat's rule) [Bauschke and Combettes, 2017, Theorem 16.

Definition 2 . 2 . 4

 224 For every η ∈ R * + , Φ η is the set of continuous concave functions φ : [0, η[→ [0, +∞[ such that (i) φ(0) = 0, (ii) φ is C 1 on ]0, η[ and continuous at 0,

Figure 2 . 1 :

 21 Figure 2.1: Illustration of the KL re-parametrization of the nonconvex function f : x ∈ R → -|x 2 -1|, which is a KL function with exponent 1 2 [Li and Pong, 2018, Corollary 5.2], with φ : s ∈ [0, 1[ → s 1 2 at the point x = 0.

Definition 2 . 2 . 5 [

 225 Bolte et al., 2014, Definition 3] Let f : R m →] -∞, +∞] be proper and l.s.c. The function is said to have the Kurdyka-Lojasiewicz property at x ∈ dom (∂f ) = {x ∈ R m | ∂f (x) = ∅} if there exist η ∈ R * + , a neighborhood Ω of x, and a function φ ∈ Φ η such that, for all x such that

2 , 2 , 3 ,

 223 4}, then we can use the results from [Chaux et al., 2007, Examples 4.3 and 4.4] which are recalled below for every γ ∈ R * + and every x ∈ R.

Figure 2 . 2 :

 22 Figure 2.2: Function | • | p (left) and its proximity operator (right) for different values of p.

Figure 2 . 3 :

 23 Figure 2.3: Logarithmic barrier function (left) and its proximity operator (right) for different values of multiplicative factor γ.

Definition 2 . 2 .

 22 16 Let T : R m → R n .(i) The operator T is Lipschitz continuous with constant α ∈ R + , or α-Lipschitz continuous, if

  [Bauschke and Combettes, 2017, Definition 4.33] Let T : R m → R m and α ∈]0, 1[. The operator T is averaged with constant α, or α-averaged, if there exists a nonexpansive operator R : R m → R m such that T = (1α)I m + αR.

Proposition 2 . 2 .

 22 18 [Bauschke and Combettes, 2017, Remark 4.34, Proposition 4.35] Let T : R m → R m and α ∈]0, 1].

Figure 2 . 5 :

 25 Figure 2.5: Illustration of Algorithm 9 applied to problem (2.18) with c = (1, 2) , and where the feasible set is represented as a blue polytope. Hyperparameters are chosen such that µ 0 = 10, ρ = 4 and = 0.02. Brown dotted lines show the iso-B contours.

Theorem 2 .

 2 2.25 [Boyd and Vandenberghe, 2004, Section 11.5.3] Consider problem (2.18) and assume that Assumption 2.2.23 is satisfied. Let (x k ) k∈N be generated by Algorithm 10. Then,

Definition 2 . 2 . 26 (

 2226 Analytic center) Let C by defined as in(2.16) and let K be defined by a set of equalities,K = {x ∈ R m | (∀i ∈ {1, . . . , q}) h i (x) = 0},where (∀i ∈ {1, . . . , q}) h i : R m →] -∞, +∞]. Let B be the logarithmic barrier associated to C. If there exists a unique solution to the following problem, minimize x∈K B(x), then this point is called the analytic center of C ∩ K. This notion has led to a specific class of IPMs called analytic center methods [Sonnevend, 1986] [Jarre, 1992, Section 2.1]. The analytic center of a set of inequalities can be thought of as the feasible point which maximizes the product of its margins with regards to the boundaries of this set [Boyd and Vandenberghe, 2004, Section 8.5.3]. It is worth noting that the analytic center depends on how the feasible set is specified. Adding a redundant inequality in the definition of C will change the position of the analytic center.

Figure 2 .

 2 6 illustrates this observation.

2. 3

 3 Figure 2.6: Left: analytic center x a of the feasible set defined by the affine blue constraints. Right: influence of a redundant constraints in the definition of the feasible set on the position of the analytic center. Colors represent the logarithmic barrier value.

  and the involved functionals satisfy the following assumption. Assumption 3.4.2 (i) f is a polyhedral function, i.e. its epigraph a finite intersection of closed halfspaces.

Theorem 3 . 4 . 4

 344 Let µ > 0 and consider the barrier problem P µ associated to problem(3.4). Under Assumptions 3.2.1 and 3.4.2, for every (δ, θ) ∈ ]0, 1[ 2 , (γ, ν) ∈ (R * + ) 2 , ν ∈]0, ν], and x 0 ∈ D, the sequence (v k+1 ) k∈N generated by Algorithm 11 converges linearly to 0 when = 0.

. 6 ) 2 .

 62 Our main convergence result, summarized in Theorem 3.4.7, provides a useful characterization for the limit point of Algorithm 12 using the notion of analytic center introduced in Section 2.2.5.4. The definition and the proof of existence of the analytic center are given in Section 3.5.3.Theorem 3.4.7 Under Assumptions 3.2.1 and 3.3.1, the following statements hold. (i) If there exists only one element in S P (resp. S D ), then the sequence (x j+1 ) j∈N (resp. (λ j+1 ) j∈N ) produced by Algorithm 12 converges to this unique primal (resp. dual) solution to P 0 .

3. 5 .

 5 1 Well-definedness of Algorithm 113.5.1.1 Preliminary resultsFirst, we show that there exists a solution to every intermediate problem (3.2). Corollary 3.5.1 Under Assumption 3.2.1, for every µ ∈ R * + , the solution set to P µ is a nonempty convex compact subset of D.

Lemma 3 . 5 . 3

 353 Under Assumption 3.2.1, for every θ ∈]0, 1[, (γ, µ) ∈ (R * + ) 2 and M ∈ S + m , if x ∈ D, then lim l→+∞ h(x, l) = x.(3.9)

. 14 )

 14 Combining (3.13) and (3.14) completes the proof. Finally, we derive Corollary 3.5.5 below which states that the line search performed at each iteration of Algorithm 11 is properly defined. Corollary 3.5.5 Let (δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν], and x 0 ∈ D. Suppose that Assumption 3.2.1 holds and that Algorithm 11 is run at iteration k ∈ N. Then, x k+1 ∈ D and condition (3.3) is met for some l ∈ N.

3. 5

 5 .1.3 Key elements for the proof of Theorem 3.4.1 Let us first show the following sufficient decrease property regarding the values of the merit function.Lemma 3.5.6 (Sufficient decrease) For every (δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν], and x 0 ∈ D, if Assumption 3.2.1 holds, then the sequence (x k ) k∈N produced by Algorithm 11 with = 0 satisfies, for every k ∈ N, the inequality

  ϑ is a continuous function defined on a compact set and, since D is open, it is positive valued. It follows from the extreme value theorem, that there exists η ∈]0, +∞[ such that η = min x∈K ϑ(x). For every z ∈ K, let B(z, η/2) be the open ball with center z and radius η/2. For every y ∈ B(z, η/2), yz ≤ η/2 < η ≤ dist(z, R m \ D), which implies that y ∈ R m \ D, that is y ∈ D. This shows that (∀z ∈ K) B(z, η/2) ⊂ D.(3.22)

  .25) Set I 0 = {k ∈ I | k ≥ k 0 }. By applying the triangle inequality, we deduce from (3.21), (3.24), and (3.25) that

3. 5 . 1 . 4

 514 Proof of Theorem 3.4.1

Lemma 3 . 5 . 12 (

 3512 Continuity condition) Under Assumptions 3.2.1, let (δ, θ) ∈]0, 1[ 2 , (γ, µ, ν) ∈ (R * + ) 3 , ν ∈]0, ν]and x 0 ∈ D. If the sequence (x k ) k∈N is produced by Algorithm 11 with = 0, then there exists a subsequence (x kq ) q∈N and there exists x ∈ D such that lim q→+∞ x kq = x and lim q→+∞ Ψ µ (x kq ) = Ψ µ (x).

. 32 )

 32 Altogether (3.30) and (3.32) lead to lim q→+∞ Ψ µ (x kq ) = Ψ µ (x).

  Under Assumption 3.2.1, Lemmas 3.5.6, 3.5.10 and 3.5.12 hold. If, in addition, Ψ µ is a KL function, then we can apply[Attouch et al., 2013, Theorem 2.9]. Thus, (x k ) k∈N converges to a critical point of Ψ µ and has finite length. By convexity, every critical point of Ψ µ is a global minimizer of Ψ µ , and a solution to P µ .We now show that, for a useful special case, Algorithm 11 converges linearly in terms of iterate and objective function value. Proposition 3.5.14 Let µ ∈ R * + and consider the barrier problem P µ associated to problem (3.4). Under Assumptions 3.2.1 and 3.4.2, for every (δ, θ) ∈]0, 1[ 2 , (γ, ν) ∈ (R * + ) 2 , ν ∈]0, ν], and x 0 ∈ D, the sequence (x k ) k∈N generated by Algorithm 11 with = 0 converges linearly to a solution x * to P µ , and (Ψ µ (x k )) k∈N converges linearly to Ψ(x * ). Proof. Let µ ∈ R * + . Under Assumptions 3.2.1 and 3.4.2 we can apply [Li and Pong, 2018, Corollary 5.1] which states that Ψ µ is a KL function with exponent 1/2. The convergence of (x k ) k∈N to a solution x * ∈ D to P µ is guaranteed by Proposition 3.5.13. From Lemmas 3.5.6, 3.5.10 and 3.5.12, we can apply [Frankel et al., 2015, Theorem 4(ii)]. Hence, there exist c ∈ R * + and k 0 ∈ N such that for every k ≥ k 0 , Ψ µ (x k ) -Ψ µ (x * ) = O exp(-cK µ ) k (3.33) and x *x k = O exp(-cK µ /2) k . (3.34)

3. 5 . 3 . 3

 533 Proof of Theorem 3.4.7 (i) The result follows from Theorem 3.4.5.

Following[

  Iordache et al., 2012; Moussaoui et al., 2012; Chouzenoux et al., 2014a], we propose to formulate the following constrained minimization problem to perform the unmixing task,

  Shefi and Teboulle, 2014, Algorithm 2], where the metrics are taken constant as in [Shefi and Teboulle, 2014, Example 3.4]. Regarding GFBS-VM, we implement [Raguet and Landrieu, 2015, Algorithm 1] with a modified metric based on the Hessian of the data-fitting term in (3.61). All computational times are given for experiments run on Matlab 2018b on an Intel Xeon CPU E5-1650 at frequency 3.20 GHz.

Figure 3 . 1 :

 31 Figure 3.1: SNR as a function of time.

Figure 3 . 2 :

 32 Figure 3.2: Abundance map of Asphalt road (first row) and Dirt (second row): (a) ground-truth, (b) IPLS solution, visual results after running (c) ADMM and (d) PIPA-VM for 11 seconds.

Figure 3 . 3 :

 33 Figure 3.3: Relative distance from current iterate to limit point as a function of time.

  .66) where G(x) = µ j Diag(d 1 (x)) + H Diag(d 2 (x))H . (3.67) Hereabove, Diag(d 1 (x)), Diag(d 2 (x)) and Diag(d 3 (x)) are the diagonal matrices whose elements are given by d 1 (x) ∈ R r , d 2 (x) ∈ R n , and d 3 (x) ∈ R r , respectively. These vectors are defined as follows.

Figure 3 . 4 :

 34 Figure 3.4: Natural images: (left) phase-separated barium borosilicate glass sample, imaged at the ESRF synchrotron (courtesy of David Bouttes), (right) mushroom Agaricus bisporus (courtesy of DigiMorph.org, The University of Texas High-Resolution X-ray CT Facility (UTCT), and NSF grant IIS-0208675).

Figure 3 . 5 :

 35 Figure 3.5: Relative distance from the iterates to the limit point as a function of time for Glass.

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Relative distance from the iterates to the limit point as a function of time for Agaricus.

Figure 3 . 8 :

 38 Figure 3.8: SNR (dB) for x t+g as a function of time for Glass.

Figure 3 . 9 :

 39 Figure 3.9: SNR (dB) for x t+g as a function of time for Agaricus.

  (4.18)-(4.20) yields (4.14)-(4.16).

Proposition 4 . 3 . 3

 433 Let γ ∈ R * + and let µ ∈ R * + . Let B be the barrier function associated to (4.21), defined as

Figure 4 . 1 :

 41 Figure 4.1: Proximity operator of the logarithmic barrier: prox γµB (x) 1 for a constraint on the 2 -norm as in Section 4.3.3 with α = 0.7.

Figure 4 . 2 :

 42 Figure 4.2: iRestNet global architecture.

Figure 4 . 3 :

 43 Figure 4.3: Architecture of L

Proposition 4 . 5 . 2

 452 Let K ∈ N * , θ ∈ R * + , and α ∈ [1/2, 1]. Let W ∈ S m and let β -and β + denote the smallest and largest eigenvalues of W , respectively. Then, the condition

(4. 38 )From ( 4 .

 384 37) and(4.38), we deduce thatW -2 K (1α)I m = -β -+ 2 K (1α).Replacing |||W ||| and W -2 K (1α)I m by their value in (4.36) leads to Proposition 4.5.2(i).

Theorem 4 . 5 . 3 D

 453 Let α ∈ [1/2, 1], (W k , b k , R k ) 0≤k≤K-1 be defined by (4.35), and (θ k ) -1≤k≤K-1 be defined as in Proposition 4.5.1. Let β -and β + be the smallest and largest eigenvalues ofW = W K-1 • • • • • W 0 , respectively.For every p ∈ {1, . . . , m} and every k ∈ {0, . . . , K -1}, let β denote the pth eigenvalue of H H and D D in P, respectively. Then, β -, β + , and (∀k ∈ {0, . . . , K -1}) θ k can be computed as follows:

Figure 4 . 4 :

 44 Figure 4.4: Architecture of L pp . BN: batch normalization.

Figure 4 . 5 :

 45 Figure 4.5: Sorted improvement of iRestNet with regards to other methods on the BSD500 test set using the SSIM metric: a negative value indicates a better performance of iRestNet. (a): GaussianA, (b): GaussianB, (c): Gaus-sianC, (d): MotionA, (e): MotionB, (f): Square.

Figure 4 .Figure 4 . 8 :

 448 Figure 4.9 shows the stepsize, barrier parameter and regularization weight sequences obtained by passing the image from Figure 4.6 through the 40 layers of iRestNet.

Figure 4 . 9 :

 49 Figure 4.9: Left to right: estimated stepsize (γ k ) 0≤k≤K-1 , barrier parameter (µ k ) 0≤k≤K-1 and regularization weight (λ k ) 0≤k≤K-1 for the image from Figure 4.6 passed through the network layers.

5. 3 .Proposition 5 . 3 . 1

 3531 Preconditioned P-ULA sampling algorithm 115 For every γ

2 Figure 5 . 1 :

 251 Figure 5.1: Hierarchical Bayesian model. Parameters in boxes are fixed in advance.

Simu1 Simu2 σ 2 α 1 β 1 α 2 β 2 σ 2 α 1 β 1 α 2 β 2 α 3 β 3 Table 5 . 3 :

 353 True 0.013 1.5 1.0 0.60 1.0 33 1.5 100 1.0 50 0.50 4.0 P-ULA 0.041 2.0 0.5 0.59 1.0 122 2.0 330 2.0 3186 0.48 3.4 HMC 0.013 1.8 1.2 0.61 1.0 34 1.4 66 1.1 111 0.54 5.2 PP-ULA 0.013 1.4 0.9 0.62 1.1 35 2.3 2676 1.2 122 0.55 5.8 MMSE Estimates of the noise variance and GGD parameters.

Figure 5 . 3 :

 53 Figure 5.3: PSNR along time for Simu1. Dotted lines indicate the PSNR of the MMSE estimator of the TRF after the burn-in regime.

Figure 5 . 4 :

 54 Figure 5.4: Simu1, GGD distributions (5.7) of regions 1 (left) and 2 (right).

Figure 5 . 5 :Figure 5 Figure 5 . 7 :

 55557 Figure 5.5: B-mode visualization of the TRF for Simu1: (a) ground-truth, (b) Wiener, (c) Lasso, (d) P-ULA, (e) HMC, (f) PP-ULA.

Figure 5 . 8 :

 58 Figure 5.8: Segmentation results for Simu2: (a) ground-truth, (b) Otsu, (c) SLaT, (d) P-ULA, (e) HMC, (f) PP-ULA. Main differences are circled in green.

Figure 5 .Figure 5 . 9 :

 559 Figure 5.11, that the segmentation based on the Potts model (PP-ULA and HMC) gives more homogeneous areas than Otsu, and recovers more details than SLaT.

Figure 5 . 10 :Figure 5 . 11 :

 510511 Figure 5.10: B-mode visualization of the TRF results for real data: (a) Wiener, (b) Lasso, (c) HMC, (d) PP-ULA. From top to bottom: Thyroid, Bladder and KidneyReal.

  .1.

	Definition 2.2.6 [Li and Pong, 2018, Definition 2.3] Let f : R m →] -∞, +∞] be a proper l.s.c. function satisfying the KL property at x ∈ dom (∂f ) within a neighborhood Ω for η ∈ R * + and + . If there exist c ∈ R * α ∈ [0, 1[ such that the corresponding re-parametrization φ ∈ Φ η can be chosen as

  Hereafter, we list two examples that are relevant to this thesis.

	Example 2.2.14 In Chapter 5, we consider functions of the type • p p for p ∈ R * p is separable with regards + . Since the proximity operator of • p to its entries, we consider here functions of the form | • | p . When p ≥ 1, | • | p is convex and we can obtain an closed-form expression of the proximity operator for some specific values of p. The resulting operators are illustrated
	in Figure 2.2.

  5.2].The variable metric forward-backward (VMFB) algorithm has recently been extended in[Combettes and Vũ, 2014], where the authors consider summable errors in the computation of the gradient and of the proximity operator. The convergence of this algorithm is obtained in[Combettes and Vũ, 2014, Theorem 4.1] under a monotonicity assumption on the variable metrics. Some recent works, relying on the KL inequality, extend this convergence result to nonconvex problems under milder boundedness conditions on the variable metrics[Chouzenoux et al., 2014b;[START_REF] Frankel | Splitting methods with variable metric for Kurdyka-Lojasiewicz functions and general convergence rates[END_REF].Algorithm 3 corresponds to a simplified version of the VMFB algorithm without the summable errors. By taking (∀k ∈ N) M k = I m in Algorithm 3, we retrieve the classical FB algorithm. Algorithm 3: Variable metric forward-backward [Combettes and Vũ, 2014] (exact version)

)

  Chapter 2. BackgroundIn the dual problem associated to(2.11), h * and L are decoupled. Hence, using duality can lead to iterative solvers that only require the proximity operator of h * , which can easily be derived from the one of h using Moreau's decomposition from Theorem 2.2.9. We present in Algorithm 5 a simplified version of the dual forward-backward (DFB) algorithm proposed in[Combettes et al., 2011] applied to problem(2.11), where we assume that there is no error in the computation of the proximity operator of h * . The convergence result for Algorithm 5 is given in[Combettes et al., 2011, Theorem 2.2].

	Algorithm 5: Dual forward-backward [Combettes et al., 2011] (exact
	version) applied to problem (2.11)

  then we recover our starting problem (2.8). Problem (2.13) can be reformulated as follows,

	minimize (x,y)∈R m ×R n	f (y)	(2.14)
	s.t.	Ax + By = 0 n

  addition, if Assumption 2.2.20 below holds and if the barrier is chosen as the logarithmic one, then the sublevel sets of the merit function are compact. Theorem 2.2.21 [Wright, 1992, Theorem 4] 2 Under Assumption 2.2.20, for every (µ, τ

This result is detailed in Theorem 2.2.21 hereafter. Assumption 2.2.20 Consider problem

(2.15)

, and assume that f ∈ Γ 0 (R m ) and (∀i ∈ {1, . . . , p}) c i ∈ Γ 0 (R m ). Let B be the logarithmic barrier associated with the set C defined in

(2.16)

. Assume in addition that dom (B) is nonempty, C ⊂ dom (f ), and that the set of solutions to problem (2.15) is nonempty and bounded.

  .25) In the above system of equations, the complementary slackness from (2.21) is replaced by condition (2.24). Hence, the path-following method consists in finding solutions to modified KKT conditions such that the iterates produced by Algorithm 9 almost satisfy (2.19)-(2.21) when µ is small. It is worth noting that strict complementarity always holds in linear programming [Bonnans and Shapiro, 2013, Theorem 3.133], i.e. there exists (x * , λ * ) a solution to the KKT conditions (2.19)-(2.21) such that, for every

  .3, on a proximal interior point iterative scheme. Lastly, in Chapter 5 we consider a Bayesian setting and propose to accelerate a hybrid Gibbs sampler with a new sampling algorithm based on the proximity operator.

	-Chapter 3 -
	A fast proximal interior point algorithm for
	constrained variational formulations

As pointed out in Section 2.1.2, many problems in image processing, such as segmentation

[START_REF] Chan | Active contours without edges[END_REF]

, classification

[START_REF] Briceño-Arias | A random block-coordinate Douglas-Rachford splitting method with low computational complexity for binary logistic regression[END_REF]

, or restoration

[START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF][START_REF] Nikolova | A variational approach to remove outliers and impulse noise[END_REF]

, can be formulated as the minimization of a convex objective function under convex constraints. Such problem can be successfully addressed by interior point methods, which have been presented in Section 2.2.5. From a numerical perspective, IPMs have demonstrated very good performance on several challenging applications, such as image reconstruction and multispectral image unmixing

[START_REF] Johnson | Interior-point methodology for 3-D PET reconstruction[END_REF] Chouzenoux et al., 2014a

]. However, it is worth noting that lots of interior point approaches rely on second-order methods and, therefore, assume that the objective function is at least twice-differentiable

[START_REF] Armand | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF] Bonettini and Serafini, 2009a]

. This characteristic restricts their scope of application since the quality of the solution to an inverse problem and its robustness to noise, can often be improved by including a nondifferentiable regularization term in the objective function. Although IPMs can handle the 1 norm

[START_REF] Fu | Efficient minimization methods of mixed 2-1 and 1-1 norms for image restoration[END_REF][START_REF] Kim | An interiorpoint method for large-scale 1 -regularized least squares[END_REF][START_REF] Bibliography | Performance of first-and second-order methods for 1 -regularized least squares problems[END_REF]

, for more general non-smooth penalizations, optimization approaches relying on the proximity operator seem more appropriate

[Combettes and Pesquet

  al., 2014b; Frankel et al., 2015], and for deriving convergence rates [Attouch and Bolte, 2009; Attouch et al., 2013].

  2.1 below is satisfied. (iii) Functions g and (c i ) 1≤i≤p are differentiable on C; ∇g and (∇c i ) 1≤i≤p are Lipschitz-continuous on any compact subset of C. Since the functions (c i ) 1≤i≤p are l.s.c., the closure of D is equal to C.

	Remark 3.2.2 (Role of the constraints) In image processing, the con-
	straints can be derived from the underlying geometry of the problem [Hariza-
	nov et al., 2013]. For instance, inequality constraints are used in a problem
	of deformable image matching in
	Assumption 3.2.1
	(i) The set of solutions to P 0 is nonempty and bounded.

(ii) Functions f , g and (c i ) 1≤i≤p belong to Γ 0 (R m ), and f + g is bounded from below. The set D is assumed to be nonempty, open, and C ⊂ dom (f ) ⊂ dom (g).

  2017, LS1], itself a generalization of [Bello Cruz and Nghia, 2016]. At iteration k ∈ N, this backtracking procedure stops if

  .23) Let S = ∪ J j=1 B(z j , η/2). It follows from (3.22) and (3.23) that this set is a compact subset of D including K. From (3.23), for every k ∈ I, there exists j k ∈ {1, . . . , J} such that x kz j k <

	η 4	.	(3.24)

  al., 2010; Iordache et al., 2012], which is presented in Algorithm 8, PDS [Komodakis and Pesquet, 2015; Combettes et al., 2014], summarized in Algorithm 6, and GFBS [Raguet et al., 2013], which is presented in Algorithm 4.

Table 3 . 1 :

 31 VM 11.31 12.25 13.04 15.27 7.12 14.52 Signal-to-noise ratio (SNR i ) 1≤i≤6 (in dB) for each material after 11 seconds.

	3.7. Joint geometry-texture decomposition and reconstruction	71
		Asphalt Grass	Tree	Roof Metal Dirt
	IPLS	10.12	11.21 11.86 14.91	4.90	13.68
	PDS	2.23	3.65	4.83	8.56	7.12	10.30
	GFBS	2.31	3.82	4.88	9.38	0.40	10.24
	GFBS-VM	2.50	3.58	4.86	7.46	-0.42 10.71
	ADMM	7.40	11.37 12.45 15.08	7.25	12.34
	ADMM-VM	2.78	4.90	4.01	8.02	0.70	6.79
	PIPA	1.71	2.17	2.90	2.20	-0.94	8.51
	PIPA-

Table 3 . 2 :

 32 SNR (dB) of the reconstruction x t+g obtained after running the algorithms for the same duration.

	3.8. Summary

  . If H H and D D are diagonalizable in the same basis then W ∈

	S m , which, combined with Proposition 4.5.1, leads to (4.40) and (4.41).
	If one of the conditions (i)-(iii) is satisfied, then we deduce from Proposi-
	tion 4.5.2 that W satisfies [Combettes and Pesquet, 2018, Proposition 3.6(iii)]
	and [Combettes and Pesquet, 2018, Condition 3.1]. In addition, for every
	k ∈ {0, . . . , K -1}, R k (• + b k ) is firmly nonexpansive [Bauschke and Com-bettes, 2017, Proposition 12.28]. Finally, [Combettes and Pesquet, 2018,
	Theorem 3.8] completes the proof.
	The conditions provided by Theorem 4.5.3 can be easily checked using (4.40).
	Theorem 4.5.3 provides a framework under which iRestNet is robust to a
	perturbation of its input: the upper bound of the output perturbation can
	then be derived from Proposition 2.2.18.
	§ 4.6 Numerical experiments

Table 4 .

 4 1: Training information for Gaussian blur kernels: initial learning rates and number of epochs. For every couple, the first and second numbers correspond to the training of (L k ) 0≤k≤29 and L pp •L 39 •. . .•L 30 , respectively.

		GaussianA	GaussianB	GaussianC
	Rates (0.01,0.001) (0.01,0.001) (0.001,0.001)
	Epochs	(40,393)	(40,340)	(40,300)
		MotionA	MotionB	Square
	Rates (0.01,0.002) (0.01,0.001) (0.01,0.005)
	Epochs	(40,1200)	(40,1250)	(40,740)

Table 4 . 2 :

 42 Training information for motion and square blur kernels: initial learning rates and number of epochs. For every couple, the first and second numbers correspond to the training of (L k ) 0≤k≤29 and L pp • L 39 • . . . • L 30 , respectively.

Table 4 .

 4 3: SSIM results on the BSD500 test set.

		GaussianA GaussianB GaussianC MotionA MotionB Square
	Blurred	0.676	0.526	0.326	0.383	0.549	0.544
	VAR	0.804	0.723	0.587	0.819	0.829	0.756
	EPLL	0.800	0.708	0.565	0.816	0.839	0.755
	MLP	0.821	0.734	0.608	0.854	n/a	n/a
	PDHG	0.796	0.716	0.563	0.801	n/a	n/a
	IRCNN	0.841	0.768	0.619	0.902	0.907	0.834
	FCNN	n/a	n/a	n/a	0.794	0.847	n/a
	iRestNet	0.853	0.787	0.641	0.898	0.910 0.840
		GaussianA GaussianB GaussianC MotionA MotionB Square
	Blurred	0.723	0.545	0.355	0.376	0.590	0.579
	VAR	0.857	0.776	0.639	0.856	0.869	0.818
	EPLL	0.860	0.770	0.616	0.857	0.887	0.827
	MLP	0.874	0.798	0.668	0.891	n/a	n/a
	PDHG	0.853	0.781	0.623	0.855	n/a	n/a
	IRCNN	0.885	0.819	0.676	0.927	0.930 0.886
	FCNN	n/a	n/a	n/a	0.801	0.890	n/a
	iRestNet	0.892	0.833	0.696	0.919	0.930 0.886

Table 4 . 4 :

 44 SSIM results on the Flickr30 test set.the training set. Table4.4 contains the average SSIM obtained with the different methods on the Flickr30 test set. Similarly to the BSD500 test set, iRestNet compares favorably with the other approaches on the Flickr30 test set.

  and separable, its proximity operator in the Euclidean metric is straightforward to compute. More precisely, for every i ∈ I ≥1 , η ∈ R * + and s ∈ R, prox η -1 |•| αz i (s) has either a closed form, given in Example 2.2.14, or can be found using a bisection search in [0, |s|].

Table 5 . 2 :

 52 2 /|||H||| 2 . Number of iterations, computational time and MSJ per sec for experiments Simu1 and Simu2.

			Iterations	Time	Mixing
			Burn-in Total	Duration	PP-ULA speed gain (per sec) MSJ
		P-ULA	70000 140000 2 h 27 min	12.2	665
	Simu1	HMC	4000	8000	1 h 08 min	5.7	173
		PP-ULA	4000	8000	12 min	1	970
		P-ULA	70000 140000 3 h 06 min	4.8	590
	Simu2	HMC	10000 20000 4 h 14 min	6.6	22
		PP-ULA	10000 20000	39 min	1	793

Table 5 . 4 :

 54 ULA 40.3 0.62 1.51 99.7 38.6 0.71 1.64 98.7 PSNR, SSIM, CNR and segmentation OA for simulated data. SSIM values from Table 5.4 obtained with PP-ULA are equivalent or higher than all competitors for these two experiments. Visual segmentation results are shown in Figures 5.7 and 5.8, and OA values can be found in Table 5.4. For these simulated images, more pixels are correctly labeled with PP-ULA than with competitors.

	Simu1	Simu2
	PSNR SSIM CNR OA PSNR SSIM CNR OA
	Wiener -Otsu 37.1 0.57 1.26 99.5	35.4 0.63 0.97 96.0
	Lasso -SLaT 39.2 0.60 1.15 99.6	37.8 0.70 0.99 98.3
	P-ULA 38.9 0.45 1.82 98.7	37.1 0.57 1.59 94.9
	HMC 40.0 0.62 1.47 99.7	36.4 0.64 1.59 98.5
	PP-

Table 5 . 5 :

 55 Number of iterations, computational time and MSJ per sec for experiments on the tissue-mimicking phantom and on real data.

			Time	Mixing
		Burn-in Total	Duration	PP-ULA speed gain (per sec) MSJ
	Kidney	HMC 7000 14000 4 h 23 min PP-ULA 7000 14000 42 min	6.3 1	167 657
	Thyroid	HMC 3000 6000 2 h 09 min PP-ULA 3000 6000 35 min	3.7 1	175 950
	Bladder	HMC 5000 10000 2 h 45 min PP-ULA 5000 10000 32 min	5.2 1	13 1396
	KidneyReal	HMC 5000 10000 1 h 49 min PP-ULA 5000 10000 19 min	5.8 1	11 1361

Table 5

 5 .6 obtained with PP-ULA on the Kidney experiment are equivalent or higher than all competitors.

		Kidney	Thyroid Bladder KidneyReal
		PSNR SSIM CNR	CNR	CNR	CNR
	Wiener	27.6 0.58 0.66	0.56	1.66	1.61
	Lasso	28.5 0.59 0.67	0.99	1.76	1.76
	HMC 29.5 0.62 1.10	1.52	2.23	1.88
	PP-ULA	29.3 0.62 1.14	1.56	2.48	1.93

Table 5 .

 5 6: PSNR, SSIM and CNR results.

  , which have demonstrated better performances than local strategies for several examples in image restoration [Kindermann et al., 2005; Li et al., 2017].

http://proximity-operator.net/

It is straightforward to extend the original theorem, which was established assuming continuous functions, to l.s.c. functions.

https://web.stanford.edu/~boyd/cvxbook/cvxbook_examples/chap11/

https://github.com/Po-Hsun-Su/pytorch-ssim

In[START_REF] Zhang | Learning deep CNN denoiser prior for image restoration[END_REF], this algorithm is improperly called half-quadratic splitting, but it does not correspond to usual half-quadratic optimization methods described for instance in[START_REF] Allain | On global and local convergence of half-quadratic algorithms[END_REF]. Actually, the algorithm unfolded in[START_REF] Zhang | Learning deep CNN denoiser prior for image restoration[END_REF] can be interpreted as a preconditioned forward-backward algorithm.

Following Landau's notation, we will write that F (u) = o( u-x ), where F : R m → R and x ∈ R m , if F (u)/ ux → 0 as u → x.

https://github.com/mccorbineau/PP-ULA

Remerciements

An optimization-inspired neural network for image deblurring

From a variational formulation of an inverse problem, we have derived in this chapter a novel neural network architecture by unfolding a proximal A preconditioned proximal sampling algorithm for US imaging segmentation of US images. This chapter is organized as follows. After formulating the studied inverse problem in Section 5.1 and motivating our approach in Section 5.2, we detail the proposed strategy used to sample the TRF in Section 5.3. Section 5.4 then describes the investigated hierarchical Bayesian model, while the proposed hybrid Gibbs sampler is detailed in Section 5.5. Finally, numerical experiments on both simulated and real US data are presented in Section 5.6 and some conclusions are drawn in Section 5.7. § 5.1 Problem statement Let x ∈ R m and y ∈ R m be the vectorized TRF and radio-frequency image, respectively. The following simplified model is used [START_REF] Jensen | Deconvolution of in-vivo ultrasound B-mode images[END_REF][START_REF] Ng | Modeling ultrasound imaging as a linear, shift-variant system[END_REF] y = Hx + ω, (

where H ∈ R m×m is a linear operator that models the convolution with the point spread function (PSF) of the probe, and ω ∼ N (0 m , σ 2 I m ). In this chapter, we assume that the PSF is known, while the variance of the Gaussian noise σ 2 ∈ R * + is an unknown parameter to be estimated. The TRF is comprised of K different tissues, which are identified by a hidden label field z = (z i ) 1≤i≤m ∈ {1, . . . , K} m . For every k ∈ {1, . . . , K}, the kth region is modeled by a generalized Gaussian distribution (GGD) [START_REF] Alessandrini | A restoration framework for ultrasonic tissue characterization[END_REF][START_REF] Zhao | Joint segmentation and deconvolution of ultrasound images using a hierarchical Bayesian model based on generalized Gaussian priors[END_REF], which is parametrized by a shape parameter α k ∈ [0, 3], related to the scatterer concentration, and a scale parameter β k ∈ R * + , linked to the signal energy. Given y and H, the aim is to estimate a deblurred image x [START_REF] Jensen | Deconvolution of ultrasound images[END_REF][START_REF] Michailovich | Blind deconvolution of medical ultrasound images: A parametric inverse filtering approach[END_REF], as well as the noise variance σ 2 , the GGD parameters α = (α k ) 1≤k≤K and β = (β k ) 1≤k≤K , and the label field z. § 5.2 Motivation Due to the interdependence of the unknowns introduced in the previous section, it is beneficial to perform the deconvolution and segmentation tasks in a joint manner [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF][START_REF] Pirayre | HOGMep: Variational Bayes and higher-order graphical models applied to joint image recovery and segmentation[END_REF]. This is achieved in [START_REF] Zhao | Joint segmentation and deconvolution of ultrasound images using a hierarchical Bayesian model based on generalized Gaussian priors[END_REF] by considering a hierarchical Bayesian model, which is used within an MCMC method 1 to sample x, σ 2 , α, β, and z according to the full conditional distribution. Despite promising results in image restoration and segmentation, the method in [START_REF] Zhao | Joint segmentation and deconvolution of ultrasound images using a hierarchical Bayesian model based on generalized Gaussian priors[END_REF] is of significant computational complexity, in particular due to the adjusted 1 Please refer to Section 2.1.3.2 for more information about MCMC methods.

Preconditioned P-ULA sampling algorithm 113

Hamiltonian Monte Carlo (HMC) method [START_REF] Neal | MCMC using Hamiltonian dynamics[END_REF][START_REF] Robert | Accelerating MCMC algorithms[END_REF] used to sample the TRF.

Recently, efficient and reliable stochastic sampling strategies have been devised [START_REF] Durmus | Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau[END_REF][START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF][START_REF] Schreck | A shrinkage-thresholding Metropolis adjusted Langevin algorithm for Bayesian variable selection[END_REF] using the proximity operator presented in Section 2.2.1.3, which is known as a useful tool for large-scale nonsmooth optimization [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. In this work, we investigate an MCMC algorithm to perform the joint deconvolution and segmentation of US images, where the TRF is sampled with a scheme inspired from the proximal unadjusted Langevin algorithm (P-ULA) [START_REF] Pereyra | Proximal Markov chain Monte Carlo algorithms[END_REF]. P-ULA generates samples according to an approximation of the target distribution without acceptance test, while being geometrically ergodic, whereas classical unadjusted Langevin algorithms may have convergence issues.

In the next section, we propose an original accelerated preconditioned version of P-ULA (PP-ULA), which relies on the use of a variable metric forward-backward strategy [START_REF] Stuart | Conditional path sampling of SDEs and the Langevin MCMC method[END_REF]Chouzenoux et al., 2014b]. We also introduce in Section 5.5.2 an efficient solver based on the majorizeminimize (MM) principle to tackle the involved non log-concave priors. The proposed PP-ULA is then incorporated within a new hybrid Gibbs sampler summarized in Algorithm 17, yielding a substantial reduction of the computational time needed to perform joint high-quality deconvolution and segmentation of US images. § 5.3 Preconditioned P-ULA sampling algorithm

In this section, after reminding results about the Langevin diffusion and its discretization using Euler's scheme, we provide details about the derivation of the proposed sampling algorithm, called PP-ULA. This algorithm will be used to sample the TRF in Section 5.5.

Discrete Langevin diffusion

An m-dimensional Langevin diffusion is a continuous time Markov process x : R + → R m , which is the solution to the following stochastic differential equation [START_REF] Roberts | Langevin diffusions and Metropolis-Hastings algorithms[END_REF],

where B : R + → R m is a Brownian motion, V : R m → R m×m is the volatility matrix-valued function, and the drift term b : R m → R m is such that, for Algorithm 17: Hybrid Gibbs sampler 1 Sample the noise variance σ 2 according to (5.9); 2 Sample the shape parameter α using Algorithm 1 applied to (5.10);

3 Sample the scale parameter β according to (5.11);

4 Sample the hidden label field z using (5.15) ; 5 Sample the TRF x using PP-ULA (5.16)-(5.17). § 5.6 Numerical experiments

Experimental settings

Simu1 and Simu2 refer to simulated images with two and three regions, respectively. Kidney denotes the tissue-mimicking phantom produced from 10 6 scatterers uniformly distributed over a digital image of human kidney tissue provided with the Field II ultrasound simulator [START_REF] Jensen | Simulation of advanced ultrasound systems using Field II[END_REF]. The amplitude of each scatterer is produced using a zero-mean Gaussian distribution, whose variance is linked to the amplitude of the corresponding point on the digital image. The PSF for the aforementioned simulations is obtained with Field II and corresponds to a 3.5 MHz linear probe. The number of regions K is set to 2 for Simu1, and it is set to 3 for Simu2 and Kidney. We also perform tests on three real ultrasound images. Thyroid denotes a real RF image of thyroidal flux obtained in vivo with a 7.8 MHz probe. The unknown PSF is identified using the RF image of a wire cross-section which was acquired with the same probe. Since the diameter of the wire is of the order of a few µm, its cross-section can almost be viewed as a point. Thus, its RF image provides a good approximation of the PSF. Finally, Bladder and KidneyReal refer to the RF images of a mouse bladder and mouse kidney, respectively. Both images were obtained in vivo with a 20 MHz probe. The PSF for these two real images is estimated using the same method as for Thyroid. The number of regions K is set to 2 for KidneyReal, and it is set to 3 for Thyroid and Bladder.

Test settings are summarized Table 5.1, while the simulated and real RF images used in these numerical experiments are shown in Figure 5.2.

The TRF is initialized using a pre-deconvolved image obtained with a Wiener filter, while the segmentation is initialized by applying a 7×7 median filter, followed by the Otsu method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF], to the B-mode of the initial TRF. Shape and scale parameters are randomly selected in [0.5, 1.5], and [1, 200], respectively. The granularity parameter θ for the Potts model (5.8) is adjusted to ensure that the percentage of isolated points in the segmentation, obtained with a 3×3 median filter, is close to 0.05, 0.1, 0.8, 0.08, 0.08 and 0.08 for Simu1, Simu2, Kidney, Thyroid, Bladder and KidneyReal, respectively.

Comparisons and evaluation metrics

All computational times are given for simulations run on Matlab 2018b on an Intel Xeon CPU E5-1650 3.20 GHz. The code for the proposed method is available online 3 . In addition to comparing Algorithm 17 with HMC [START_REF] Zhao | Joint segmentation and deconvolution of ultrasound images using a hierarchical Bayesian model based on generalized Gaussian priors[END_REF], the quality of the deconvolution is compared with the one

List of Algorithms

List of Figures