# Optimization Methods for Signal Processing

Emilie Chouzenoux, Jean-Christophe Pesquet

► To cite this version:

Emilie Chouzenoux, Jean-Christophe Pesquet. Optimization Methods for Signal Processing. C.Jutten, L.T. Duarte, S. Moussaoui. Source Separation in Physical-Chemical Sensing, Chapter 2, IEEE Press; JOHN WILEY & SONS, 2023, 9781119137276. hal-04250055

HAL Id: hal-04250055

https://inria.hal.science/hal-04250055v1

Submitted on 12 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Contents

Notations iii

## 1 Optimization Methods for Signal Processing 1

E. Chouzenoux and J.-C. Pesquet

1.1 Introduction to optimization problems 2

1.1.1 Problem formulation 2

1.1.2 Theoretical background 3

1.1.3 Examples in the context of source separation 5

1.1.4 Chapter outline 7

1.2 Majorization-Minimization approaches 8

1.2.1 Majorization-Minimization principle 8

1.2.2 Majorization techniques 9

1.2.3 MM quadratic methods 14

1.2.4 Variable metric forward-backward algorithm 20

1.2.5 Block-coordinate MM algorithms 23

1.3 Primal-dual methods 28

1.3.1 Lagrange duality 28

1.3.2 Alternating direction method of multipliers 29

1.3.3 Primal-dual proximal algorithms 31

1.3.4 Primal-dual interior point algorithm 32

1.4 Illustration in the context of DOSY NMR signal restoration 38

1.4.1 Quadratic penalization 39

1.4.2 Sparsity prior in the signal domain 40

1.4.3 Sparsity prior in a transformed domain 41

1.4.4 Sparsity prior and range constraints 42

1.4.5 Concluding remarks 43

1.5 Conclusion 44
Notations

<table><tr><td>x</td><td>column vector of components x<sub>p</sub>, 1 ≤ p ≤ P</td></tr><tr><td>s, x, y</td><td>sources, observations, separator outputs</td></tr><tr><td>R</td><td>number of sources</td></tr><tr><td>P</td><td>number of sensors</td></tr><tr><td>T</td><td>number of observed samples</td></tr><tr><td>*</td><td>convolution</td></tr><tr><td>A</td><td>matrix with components A<sub>ij</sub></td></tr><tr><td>A, B</td><td>mixing and separation matrices</td></tr><tr><td>G, W, Q</td><td>global, whitening, and separating unitary matrices</td></tr><tr><td>g̃</td><td>Fourier transform of g</td></tr><tr><td>ŝ</td><td>estimate of quantity s</td></tr><tr><td>p<sub>x</sub></td><td>probability density of x</td></tr><tr><td>ψ</td><td>joint score function</td></tr><tr><td>φ<sub>i</sub></td><td>marginal score function of source s<sub>i</sub></td></tr><tr><td>E x, E{x}</td><td>mathematical expectation of x</td></tr><tr><td>I{y} or I(p<sub>y</sub>)</td><td>mutual information of y</td></tr><tr><td>K{x; y} or K(p<sub>x</sub>; p<sub>y</sub>)</td><td>Kullback divergence between p<sub>x</sub> and p<sub>y</sub></td></tr><tr><td>Hx or H(p<sub>x</sub>)</td><td>Shannon entropy x</td></tr><tr><td>L</td><td>likelihood</td></tr><tr><td>A, B</td><td>mixing, and separating (non linear) operators</td></tr><tr><td>cum{x<sub>1</sub>, ..., x<sub>P</sub>}</td><td>joint cumulant of variables {x<sub>1</sub>, ..., x<sub>P</sub>}</td></tr><tr><td>cum<sub>R</sub>{y}</td><td>marginal cumulant of order R of variable y</td></tr></table>
<table>
  <tr>
    <td>Q<sup>T</sup></td>
    <td>transposition</td>
  </tr>
  <tr>
    <td>Q<sup>H</sup></td>
    <td>conjugate transposition</td>
  </tr>
  <tr>
    <td>Q<sup>*</sup></td>
    <td>complex conjugation</td>
  </tr>
  <tr>
    <td>Q<sup>†</sup></td>
    <td>pseudo-inverse</td>
  </tr>
  <tr>
    <td>Υ</td>
    <td>contrast function</td>
  </tr>
  <tr>
    <td>ℝ</td>
    <td>real field</td>
  </tr>
  <tr>
    <td>ℝ</td>
    <td>complex field</td>
  </tr>
  <tr>
    <td>Λ̂</td>
    <td>estimator of mixing matrix</td>
  </tr>
  <tr>
    <td>rank{A}</td>
    <td>rank of matrix A</td>
  </tr>
  <tr>
    <td>krank{A}</td>
    <td>Kruskal's k-rank of matrix A</td>
  </tr>
  <tr>
    <td>diag{A}</td>
    <td>vector whose components are the diagonal of matrix A</td>
  </tr>
  <tr>
    <td>Diag{α}</td>
    <td>diagonal matrix whose entries are those of vector α</td>
  </tr>
  <tr>
    <td>trace{A}</td>
    <td>trace of matrix A</td>
  </tr>
  <tr>
    <td>det A</td>
    <td>determinant of matrix A</td>
  </tr>
  <tr>
    <td>š(ν)</td>
    <td>Fourier transform of process s(t)</td>
  </tr>
  <tr>
    <td>☒</td>
    <td>Kronecker product between matrices</td>
  </tr>
  <tr>
    <td>⊙</td>
    <td>Khatri-Rao (column-wise Kronecker) product between matrices</td>
  </tr>
  <tr>
    <td>⊠</td>
    <td>Hadamard (entry-wise) product between arrays</td>
  </tr>
  <tr>
    <td>⊗</td>
    <td>tensor product</td>
  </tr>
  <tr>
    <td>•<sub>j</sub></td>
    <td>contraction over index j</td>
  </tr>
  <tr>
    <td>∫, ∬</td>
    <td>sets</td>
  </tr>
  <tr>
    <td>T</td>
    <td>tensor or tensor array (representing multilinear operators)</td>
  </tr>
  <tr>
    <td>T<sub>ijk</sub></td>
    <td>entries of a tensor array</td>
  </tr>
  <tr>
    <td>vec(·)</td>
    <td>vectorization of an array</td>
  </tr>
  <tr>
    <td>Unvec<sub>n</sub>(·)</td>
    <td>unvectorization into a matrix with n rows</td>
  </tr>
  <tr>
    <td>T<sup>(d)</sup></td>
    <td>dth mode unfolding of a tensor array into a matrix (also called matricization or flattening)</td>
  </tr>
  <tr>
    <td>aff{a<sub>1</sub>, ..., a<sub>p</sub>}</td>
    <td>affine hull of the set {a<sub>1</sub>, ..., a<sub>p</sub>}</td>
  </tr>
  <tr>
    <td>conv{a<sub>1</sub>, ..., a<sub>p</sub>}</td>
    <td>convex hull of the set {a<sub>1</sub>, ..., a<sub>p</sub>}</td>
  </tr>
  <tr>
    <td>vol{a<sub>1</sub>, ..., a<sub>p</sub>}</td>
    <td>volume of conv{a<sub>1</sub>, ..., a<sub>p</sub>}</td>
  </tr>
</table>
dom *f* domain of function *f*

prox<sup>A</sup><sub>f</sub>(*x*) proximity operator of function *f* within the metric induced by **A** computed at *x*

ι<sub>C</sub> indicator function of set *C*

P<sub>C</sub> projector on set *C*

epi *f* epigraph of function *f*

argmin<sub>C</sub> *f* minimum argument of *f* over set *C*

argmax<sub>C</sub> *f* maximum argument of *f* over set *C*

sup<sub>C</sub> *f* supremum of *f* over set *C*

inf<sub>C</sub> *f* infimum of *f* over set *C*

∇*f* gradient of *f*

∇<sup>2</sup>*f* Hessian of *f*

∂*f*(*u*) subdifferential set of *f* at *u*

*f*<sup>*</sup> conjugate of function *f*

γ *f* Moreau envelope of *f* of parameter γ

□ infimum-convolution

||·||<sub>F</sub> Frobenius norm

||||·||| spectral norm
vi
# 1 Optimization Methods for Signal Processing

E. Chouzenoux and J.-C. Pesquet

## Abstract

In a wide range of problems arising in source separation, large scale optimization problems need to be solved. The objective of this chapter is to introduce the theoretical background which makes it possible to develop efficient algorithms to successfully address these problems. This chapter will be mainly focused on nonlinear optimization tools for dealing with convex and nonconvex problems. Proximal tools, parallel splitting techniques, Majorization-Minimization strategies and alternating minimization approaches will be presented. Illustrations of these methods on various problems will be provided.

## Keywords

Optimization, Convex analysis, Proximal algorithms, Splitting techniques, Primal-dual algorithms, Majorization-Minimization, Half-Quadratic strategies, Block coordinate minimization, Interior point algorithms.

## Affiliations

E. Chouzenoux, and J.-C. Pesquet are with OPIS of Inria Saclay (Palaiseau, France) and CVN at CentraleSupélec, Université Paris-Saclay (Saclay, France).
1.1

Introduction to optimization problems

1.1.1

Problem formulation

In this chapter, we are interested in general problems involving the minimization of a cost function $f$ with respect to a vector $\mathbf{u} \in \mathbb{R}^N$ with components $(u_n)_{1 \le n \le N}$ representing variables to estimate. In the context of source separation, the unknown vector can often be block-decomposed into $J \ge 1$ blocks having various roles and properties. More specifically, we will decompose $\mathbf{u}$ as $(\mathbf{u}^{(j)})_{1 \le j \le J} \in \mathbb{R}^{N_1} \times \dots \times \mathbb{R}^{N_J}$, where, for every $j \in \{1, \dots, J\}$, $\mathbf{u}^{(j)} = (u_n)_{n \in \mathcal{J}_j} \in \mathbb{R}^{N_j}$ and $\mathcal{J}_j$ is the index set of the components corresponding to the $j$-th block. With this notation, the optimization problem reads

$$
\begin{equation}
\begin{array}{ll}
\underset{\mathbf{u}^{(1)} \in \mathbb{R}^{N_1}, \dots, \mathbf{u}^{(J)} \in \mathbb{R}^{N_J}}{\text{minimize}} & f(\mathbf{u}). \hfill (1.1) \\
& \mathbf{u}^{(1)} \in \mathbb{R}^{N_1}, \dots, \mathbf{u}^{(J)} \in \mathbb{R}^{N_J}
\end{array}
\end{equation}
$$

Function $f$ is usually expressed as a composite objective function consisting of a sum of terms: a data fidelity term describing the existing relations existing between the observed data and the target variables, and penalization (or regularization) terms modeling the available prior information on the sought variables. Various examples of possible choices for these functions in the context of source separation will be given in Section 1.1.3. In particular, one may want to introduce some constraints on the domain of each block, namely that, for every $j \in \{1, \dots, J\}$, $\mathbf{u}^{(j)}$ belongs to some set $C_j \subset \mathbb{R}^{N_j}$. Problem (1.1) then becomes a constrained optimization problem

$$
\begin{equation}
\begin{array}{ll}
\underset{\mathbf{u}^{(1)} \in C_1, \dots, \mathbf{u}^{(J)} \in C_J}{\text{minimize}} & f(\mathbf{u}). \hfill (1.2) \\
& \mathbf{u}^{(1)} \in C_1, \dots, \mathbf{u}^{(J)} \in C_J
\end{array}
\end{equation}
$$

It is important to note that Problem (1.2) actually is an instance of Problem (1.1). Indeed, it suffices to replace the cost function $f$ in (1.1) by the sum $f + \iota_C$, where $\iota_C$ is the indicator function of set $C$:

$$
(\forall \mathbf{u} \in \mathbb{R}^N) \quad \iota_C(\mathbf{u}) = \begin{cases} 0 & \text{if } \mathbf{u} \in C \\ +\infty & \text{elsewhere,} \end{cases} \tag{1.3}
$$

and $C$ gathers all the constraints on the block variables, i.e. $C$ is the Cartesian product $C_1 \times \dots \times C_J$ of the $J$ constraint spaces. In a nutshell, allowing a cost function value to be equal to $+\infty$ over some set is a convenient way for constraining its minimizer to be out of this set.

With the exception of very specific cases (typically, fully quadratic problems), solving Problem (1.1) requires the design of an algorithm that builds a sequence of iterates $(\mathbf{u}_k)_{k \in \mathbb{N}}$ converging to a minimizer $\hat{\mathbf{u}}$ of $f$. There is in the literature a bunch of optimization methods, which differ by their convergence properties, their
<figure><img src="image_2.png" /><figcaption>Figure 1.1 Indicator function of the interval $[\delta_1, \delta_2]$</figcaption></figure>

computational cost, their memory requirements, their convergence speed, and their sensitivity to numerical errors. The goal of this chapter is to help the reader to bring answers to the following questions:

* How to exploit the mathematical properties of each term involved in $f$? How to handle constraints efficiently? How to deal with non differentiable terms in $f$? Which convergence result can be expected if $f$ is non convex?

* How to reduce the memory requirements of an optimization algorithm? How to avoid large-size matrix inversion?

* What are the benefits of block alternating strategies? What are their convergence guarantees?

* How to accelerate the convergence speed of a first-order (gradient-like) optimization method?

### 1.1.2

#### Theoretical background

##### 1.1.2.1 Convex functions

Throughout this chapter, all the functions will be assumed to take their values in $\mathbb{R}$, possibly extended to $\mathbb{R} \cup \{+\infty\}$. The domain dom $f$ of a function $f$ from $\mathbb{R}^N$ to $\mathbb{R} \cup \{+\infty\}$ is the subset of $\mathbb{R}^N$ over which this function is finite-valued. When this domain is nonempty, it is said that the function is proper. In optimization, one is able to handle discontinuous functions, but most of the time, the involved functions must be assumed to be lower-semicontinuous (lsc). A function $f$ from $\mathbb{R}^N$ to $\mathbb{R} \cup \{+\infty\}$ is lower-semicontinuous if its epigraph

$$ \mathrm{epi}(f) = \{(u, \zeta) \in \mathbb{R}^N \times \mathbb{R} \mid f(u) \le \zeta\}, $$

which basically represents the area above the graph of the function, is a closed set. It is then easy to check that continuous functions are lower-semicontinuous. In the rest of this chapter, we will assume, without recalling it explicitly, that all the considered functions are proper and lower-semicontinuous.

Another desirable property of a cost function is its convexity. Recall that a set $C$ is convex if any line segment linking two arbitrary points of $C$ is included in $C$. A
<figure><img src="image_3.png" /><figcaption>Figure 1.2 Epigraph of the absolute value function.</figcaption></figure>

function $f$ from $\mathbb{R}^N$ to $\mathbb{R} \cup \{+\infty\}$ is convex if its epigraph is convex, which is also equivalent to the following inequality:

$$ (\forall (u, v) \in (\mathbb{R}^N)^2) (\forall \lambda \in ]0, 1[) \quad f(\lambda u + (1-\lambda)v) \le \lambda f(u) + (1-\lambda)f(v). $$

In addition, the function is strictly convex if the above inequality is strict whenever $(u, v) \in (\text{dom } f)^2$ and $u \ne v$. An important property of a convex function is that any of its local minimizers i.e., minimizer on a local neighborhood in its domain, is a (global) minimizer over $\mathbb{R}^N$. If the function is strictly convex, such a minimizer, whenever it exists, is unique. These properties are important since many optimization algorithms are only guaranteed to converge to a local minimizer, which may lead to suboptimal performance in practice, when the associated cost function is non convex.

#### 1.1.2.2 Differentiability and subdifferentiability

The subdifferential of a function $f$ from $\mathbb{R}^N$ to $\mathbb{R} \cup \{+\infty\}$, at $u \in \mathbb{R}^N$, is defined as

$$ \partial f(u) = \{t \in \mathbb{R}^N \mid (\forall v \in \mathbb{R}^N) \ f(v) \ge f(u) + t^T (v - u)\}. \quad (1.4) $$

Any vector $t \in \partial f(u)$ is called a subgradient of $f$ at $u \in \mathbb{R}^N$. The above definition actually corresponds to the Moreau subdifferential, but other definitions are possible [1], which may be more relevant in the non convex case. It is worth noting that the affine term $f(u) + t^T(v - u)$, $t \in \partial f(u)$ actually correspond to a supporting hyperplane crossing the graph of $f$ at $u$, thus generalizing the notion of tangent plane to non necessarily differentiable functions. In particular, a function may have several gradients at a given point $u \in \mathbb{R}^N$ or none (in which case, $\partial f(u)$ is empty). But, if $f$ is convex and continuous at $u \in \mathbb{R}^N$, then $\partial f(u)$ is nonempty. Furthermore, if $f$ is differentiable at $u \in \mathbb{R}^N$ and it is convex, its only subgradient at $u \in \mathbb{R}^N$ is its gradient $\nabla f(u)$.

A key property of the Moreau subdifferential is that it allows the characterization of solutions to an optimization problem. Indeed, according to Fermat's rule, $\hat{u}$ is a (global) minimizer of $f$ if and only if $0$ is a subgradient of $f$ at $\hat{u}$. However, if $f$
<figure><img src="image_4.png" /><figcaption>**Figure 1.3** Function *f* and a supporting line of it at *u* = 0 (left). Graph of the subdifferential (right).</figcaption></figure>

is differentiable at $\hat{u}$ and non convex, the condition $\nabla f(\hat{u}) = 0$ only constitutes a necessary condition for $\hat{u}$ to be a minimizer of $f$.

1.1.3

**Examples in the context of source separation**

1.1.3.1 **Non-negative matrix factorization**

The problem in non-negative matrix factorization is to decompose a given , possibly noisy, data matrix $\mathbf{X} \in \mathbb{R}^{P \times T}$ into a product $\mathbf{AS}$ where $\mathbf{A} \in \mathbb{R}^{P \times R}$ and $\mathbf{S} \in \mathbb{R}^{R \times T}$ are constrained to have non-negative entries. A common strategy to solve this algebraic problem is to define $\mathbf{A}$ and $\mathbf{S}$ as the solutions of a constrained optimization problem of the form:

$$
\underset{\mathbf{A} \ge 0, \mathbf{S} \ge 0}{\text{minimize}} \Phi(\mathbf{A}\mathbf{S}, \mathbf{X}) + \Theta_{\mathrm{a}}(\mathbf{A}) + \Theta_{\mathrm{s}}(\mathbf{S}), \qquad (1.5)
$$

where $\Phi$ measures the discrepancy between the data $\mathbf{X}$ and the sought product $\mathbf{AS}$
(see Table 1.1 for some examples), and $\Theta_a$ (resp. $\Theta_s$) introduces some a priori
information on the unknown matrices.

<table><caption>**Table 1.1** Examples of discrepancy measures in NMF (with $\beta \in \mathbb{R} \setminus \{0, 1\}$).</caption>
  <thead>
    <tr>
      <th>Name</th>
      <th>Expression</th>
      <th>Reference</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Euclidean</td>
      <td><sup>1</sup>&frasl;<sub>2</sub> ||<b>X</b> &minus; <b>AS</b>||<sub>F</sub><sup>2</sup></td>
      <td>[2]</td>
    </tr>
    <tr>
      <td>Kullback-Leibler</td>
      <td>&sum;<sub>p=1</sub><sup>P</sup> &sum;<sub>t=1</sub><sup>T</sup> (<b>AS</b>)<sub>p,t</sub> &minus; X<sub>p,t</sub> log ((<b>AS</b>)<sub>p,t</sub>)</td>
      <td>[2]</td>
    </tr>
    <tr>
      <td>Itakura-Saito</td>
      <td>&sum;<sub>p=1</sub><sup>P</sup> &sum;<sub>t=1</sub><sup>T</sup> X<sub>p,t</sub> / (<b>AS</b>)<sub>p,t</sub> + log ((<b>AS</b>)<sub>p,t</sub>)</td>
      <td>[3]</td>
    </tr>
    <tr>
      <td>&beta;-divergence</td>
      <td>&sum;<sub>p=1</sub><sup>P</sup> &sum;<sub>t=1</sub><sup>T</sup> &beta;<sup>&minus;1</sup>(<b>AS</b>)<sub>p,t</sub><sup>&beta;</sup> &minus; (&beta; &minus; 1)<sup>&minus;1</sup> X<sub>p,t</sub>(<b>AS</b>)<sub>p,t</sub><sup>&beta;&minus;1</sup></td>
      <td>[4]</td>
    </tr>
  </tbody>
</table>

1.1.3.2 **Independent component analysis**

The problem of ICA consists in retrieving unknown statistically independent com-
ponents of a vector $s(t)$ which are mixed by an unknown linear operator $\mathbf{A}$, leading
to an observed vector $\boldsymbol{x}(t)$, at time $t = 1, \dots, T$. If $\mathbf{A}$ is an $R \times R$ invertible matrix, and $\mathbf{B} = \mathbf{A}^{-1}$ denotes the separating matrix, the likelihood of $\mathbf{B}$ is

$$ \mathcal{L}(\mathbf{B}) = \prod_{t=1}^{T} p_{s(t)}(\mathbf{B}\boldsymbol{x}(t)) |\det \mathbf{B}|, \qquad (1.6) $$

where $p_{s(t)}$ is the probability density function of the source vector at time $t$. For simplicity, let us now assume that the vectors $(s(t))_{1\le t \le T}$ are independent and identically distributed. The maximum likelihood estimator of $\mathbf{B}$ is then obtained by minimizing the neg-log-likelihood of $\mathbf{B}$, that is

$$ f(\mathbf{B}) = \sum_{t=1}^{T} \sum_{r=1}^{R} \varphi_r(\mathbf{B}, \boldsymbol{x}(t)) - T \log |\det \mathbf{B}|, \qquad (1.7) $$

where, for every $r = 1, \dots, R$, $\varphi_r(\mathbf{B}, \cdot) = -\log p_{s_r}([\mathbf{B}]_r)$ and $p_{s_r}$ is the probability density function of the $r$-th component of the source vector. Note that, even if these probability density functions are log-concave, $f$ is a non-convex function due to the last logarithmic term. In addition, several global minimizers may exist, especially when some of the source component probability distributions are equal (permutation ambiguity) or even (sign ambiguity). When the source distributions are unknown, a criterion similar to (1.7) can be minimized where some approximations to the probability distributions of the sources are employed [5].

A commonly used technique in ICA is to perform a prewhitening of the observations which aims at computing a matrix $\mathbf{W}$ such that, for every $t = 1, \dots, T$, $\mathbf{y}(t) = \mathbf{W}\mathbf{x}(t)$ has an identity correlation matrix $\mathbb{E}\{\mathbf{y}(t)\mathbf{y}(t)^{\top}\}$. We have then $\mathbf{s}(t) = \mathbf{Q}\mathbf{y}(t)$ where $\mathbf{Q}$ is a matrix which remains to be estimated. If the sources are assumed to be zero-mean with unit variance, then this matrix is orthogonal and the neg-log-likelihood reduces to

$$ \mathbf{Q} \mapsto \sum_{t=1}^{T} \sum_{r=1}^{R} \varphi_r(\mathbf{Q}, \mathbf{y}(t)) \qquad (1.8) $$

This simplification comes at the expense of the orthogonality constraint which needs to be imposed on the sought matrix. The resulting cost function to be minimized on the space of $R \times R$ matrices reads

$$ f(\mathbf{Q}) = \sum_{t=1}^{T} \sum_{r=1}^{R} \varphi_r(\mathbf{Q}, \mathbf{y}(t)) + \iota_{\mathcal{O}}(\mathbf{Q}), \qquad (1.9) $$

where $\mathcal{O}$ is the non-convex closed set of orthogonal matrices.

#### 1.1.3.3 Tensor decomposition

Source separation problems can often be recast either directly, or after some preprocessing (e.g. by computing higher-order cumulants of the observations), under a tensorial form [6]. One important problem in such tensor analysis is the Canonical
Polyadic Decomposition (CPD) which consists in decomposing a tensor $\mathcal{T}$ of order $D$ in a sum of $R$ one-rank tensors, i.e.

$$ \mathcal{T} = \sum_{r=1}^{R} \mathbf{a}^{(r,1)} \otimes \dots \otimes \mathbf{a}^{(r,d)}, \qquad (1.10) $$

where $R$ is the tensor rank and, for every $r = 1, \dots, R$ and $d = 1, \dots, D$, $\mathbf{a}^{(r,d)}$ is a $P_d$-dimensional vector. More generally, one may be interested in approximating a tensor $\mathcal{T}$ by a sum corresponding to the right-hand side summation in (1.10) where $R$ is then less than or equal to the tensor rank. Then, it may be useful to introduce constraints on the sought vectors $(\mathbf{a}^{(r,d)})_{1 \le r \le R, 1 \le d \le D}$, which are related to properties of the underlying physical model. For example, the components of these vectors can sometimes be assumed to belong to $\mathbb{R}_+$. Due to possible symmetry properties, some of these vectors can also be constrained to be equal. Finally, to limit indeterminacies, some of these vectors can be normalized by setting their norms to 1. All these constraints, mean that $\mathbf{a} = (\mathbf{a}^{(r,d)})_{1 \le r \le R, 1 \le d \le D}$ has to belong to a given constraint set $C$. The resulting constrained CPD problem can therefore be reformulated as the minimization of

$$ f(\mathbf{a}) = \Phi(\mathcal{T} - \sum_{r=1}^{R} \mathbf{a}^{(r,1)} \otimes \dots \otimes \mathbf{a}^{(r,d)}) + \iota_C(\mathbf{a}), \qquad (1.11) $$

where $\Phi$ is a given distance measure. The latter function is often chosen equal to the squared Frobenius norm. Because of the multilinearity of model (1.10), even if $\Phi$ is a convex function and $C$ is a convex set, the resulting optimization problem is non convex.

### 1.1.4

#### Chapter outline

In Section 1.2, we will introduce the Majorization-Minimization (MM) principle which will be our main guideline for all the approaches presented in this chapter. General recipes for building a good surrogate majorant to a given cost function will be discussed. Then, we will focus on quadratic MM methods, which will be shown to cover a wide range of existing algorithms for minimizing both smooth functions (e.g. MM subspace algorithms) and non smooth ones (e.g. the Variable Metric Forward-Backward algorithm). In the last part of this section, we will present block-coordinate approaches which now play a prominent role in the solution of large-scale optimization problems. Section 1.3 will be dedicated to primal-dual algorithms which have gained much popularity in the last years. These methods will be introduced from the viewpoint of Lagrange duality. Both proximal approaches and interior point methods will be described. We conclude this chapter by presenting in Section 1.4 an example of signal restoration, arising in spectroscopy, and various optimization approaches to resolve it.
<figure><img src="image_5.png" /><figcaption>**Figure 1.4** Graph of function *f* and a majorant function *h*(*v*) of *f* at *v*</figcaption></figure>

1.2

Majorization-Minimization approaches

Majorization-Minimization (MM) [7], also known as optimization transfer [8], aux-
iliary function method [9], quadratic approximation [10], surrogate minimization
[11], successive upper-bound minimization [12], is at the core of many optimiza-
tion algorithms for large scale signal/image processing. In particular, it appears that
most of the optimization methods that are useful to solve source separation problems
(gradient descent schemes, NMF algorithms, proximal methods, ...) can be rewrit-
ten under this generic framework, which motivates us for making a detailed focus
on MM algorithms, along with their convergence guarantees, and their applicability
on specific examples inspired from the context of source separation.

1.2.1

Majorization-Minimization principle

The main idea behind MM methods lies in the concept of majorant function, defined below.

DEFINITION 1.1 Let $f: \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$. Let $v \in \mathbb{R}^N$.
$h(\cdot|v): \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$ is a majorant function of $f$ at $v$ if

$$
(\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) \leq h(\mathbf{u}|\mathbf{v}), \quad \text{and} \quad f(\mathbf{v}) = h(\mathbf{v}|\mathbf{v}).
$$

The MM algorithm consists of solving the minimization problem

$$
\[
\underset{\mathbf{u} \in \mathbb{R}^N}{\text{minimize}} \ f(\mathbf{u}), \tag{1.12}
\]
$$

by alternating between two steps:

1) Majorize the criterion at current iterate with a majorant function,

2) Minimize the majorant function to define the next iterate.
A simple instance of MM algorithm for minimizing function $f$ over $\mathbb{R}^N$ reads:

$$ \text{For } k = 0, 1, \dots \left[ \mathbf{u}_{k+1} \in \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{Argmin}} h(\mathbf{u}|\mathbf{u}_k) \right. \qquad (1.13) $$

where $h(\cdot|\mathbf{u}_k)$ is a majorant function for $f$ at $\mathbf{u}_k$. By construction, we have the following inequalities:

$$ (\forall k \in \mathbb{N}) \quad f(\mathbf{u}_{k+1}) \le \underset{\text{Maj.}}{h(\mathbf{u}_{k+1}|\mathbf{u}_k)} \le \underset{\text{Min.}}{h(\mathbf{u}_k|\mathbf{u}_k)} = f(\mathbf{u}_k). $$

We thus see that the sequence $(f(\mathbf{u}_k))_{k \in \mathbb{N}}$ is decreasing. The construction of an MM algorithm requires to define a strategy for building majorant functions, as well as a strategy for minimizing (or at least, decreasing) them, with the aim to guarantee the convergence of the sequence $(\mathbf{u}_k)_{k \in \mathbb{N}}$ to a solution to (1.12).

### 1.2.2 Majorization techniques

In this section, we recall some important properties from [8, 11] that may be useful to construct a suitable majorant function for a given criterion $f$.

**PROPERTY 1.1** Let $f_1: \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$, $f_2: \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$ and $\mathbf{v} \in \mathbb{R}^N$. Assume that $h_1(\cdot|\mathbf{v}): \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$ is a majorant function of $f_1$ at $\mathbf{v}$, and $h_2(\cdot|\mathbf{v}): \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$ is a majorant function of $f_2$ at $\mathbf{v}$. Then,

(i) $h_1(\cdot|\mathbf{v}) + h_2(\cdot|\mathbf{v})$ is a majorant function of $f_1 + f_2$ at $\mathbf{v}$.

(ii) If, for all $\mathbf{u} \in \mathbb{R}^N$, $f_1(\mathbf{u}) \ge 0$ and $f_2(\mathbf{u}) \ge 0$, then $h_1(\cdot|\mathbf{v})h_2(\cdot|\mathbf{v})$ is a majorant function of $f_1 f_2$ at $\mathbf{v}$.

(iii) If $\phi: \mathbb{R} \to \mathbb{R} \cup \{+\infty\}$ is an increasing function, then $\phi(h_1(\cdot|\mathbf{v}))$ is a majorant function of $\phi(f_1)$ at $\mathbf{v}$.

The definition of the subdifferential in (1.4) yields the following useful property.

**PROPERTY 1.2** Let $f: \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$ and $\mathbf{v} \in \mathbb{R}^N$. If $f$ is concave (i.e., $-f$ is convex) on $\mathbb{R}^N$, then a majorant function for $f$ at $\mathbf{v}$ is

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = f(\mathbf{v}) + \mathbf{t}^\mathrm{T}(\mathbf{u} - \mathbf{v}), $$

where $(-\mathbf{t}) \in \partial(-f)(\mathbf{v})$, where $\partial f$ denotes the subdifferential of $f$.

Note that the above Property 1.2 is at the core of DC (difference of convex) programming algorithms for minimizing the difference of two convex functions [13]. At each iteration of such methods, the concave part of the criterion is majorized at the current point by a linear function while the convex part is majorized by a simpler convex (typically quadratic) function. The resulting majorant function is then minimized, according to an MM scheme.

Quadratic majorant functions can be derived, by making use of second-order properties of $f$.
**PROPERTY 1.3 (i)** Assume that $f$ is $\beta$-Lipschitz differentiable on $\mathbb{R}^N$, i.e. there exists $\beta > 0$ such that

$$ (\forall (\mathbf{u}, \mathbf{v}) \in \mathbb{R}^N \times \mathbb{R}^N) \quad \|\nabla f(\mathbf{v}) - \nabla f(\mathbf{u})\| \le \beta \|\mathbf{u} - \mathbf{v}\|. $$

Then a majorant function for $f$ at $\mathbf{v}$ is

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = f(\mathbf{v}) + (\nabla f(\mathbf{v}))^\top (\mathbf{u} - \mathbf{v}) + \frac{\mu}{2} \|\mathbf{u} - \mathbf{v}\|^2, $$

where $\mu \in [\beta, +\infty[$.

(ii) If $f$ is twice differentiable on $\mathbb{R}^N$ with Hessian $\nabla^2 f$, then a majorant function for $f$ at $\mathbf{v}$ is

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = f(\mathbf{v}) + (\nabla f(\mathbf{v}))^\top (\mathbf{u} - \mathbf{v}) + \frac{1}{2} \|\mathbf{u} - \mathbf{v}\|_{\mathbf{A}}^2, $$

where $\mathbf{A} \in \mathbb{R}^{N \times N}$ is a positive semidefinite matrix such that, for every $\mathbf{u} \in \mathbb{R}^N$, $\mathbf{A} - \nabla^2 f(\mathbf{u})$ is positive semidefinite and $\|\mathbf{u}\|_{\mathbf{A}} = \sqrt{\mathbf{u}^\top \mathbf{A} \mathbf{u}}$ denotes the weighted Euclidean norm of $\mathbf{u}$ within the metric induced by $\mathbf{A}$.

**EXAMPLE 1.4** Consider Problem (1.12), where $f$ is $\beta$-Lipschitz differentiable on $\mathbb{R}^N$. Then, Property 1.3(i) holds, and, for this particular choice of majorizing function (by setting $\mu = \gamma^{-1}$), the MM algorithm (1.13) reads:

$$ \begin{array}{l} \mathbf{u}_0 \in \mathbb{R}^N \\ \gamma \in [0, \beta^{-1}] \\ \text{For } k = 0, 1, \dots \\ \quad \lfloor \mathbf{u}_{k+1} = \mathbf{u}_k - \gamma \nabla f(\mathbf{u}_k) \rfloor \end{array} \tag{1.14} $$

The above algorithm identifies with the steepest descent algorithm (or gradient algorithm) with fixed stepsize. Its convergence properties can be found in [14], for the over-relaxed case when $\gamma \in [0, 2\beta^{-1}]$. The advantage of this method is its simplicity, but this may come at the expense of a slow convergence rate. The convergence speed of the steepest descent algorithm can be much improved by introducing a variable stepsize associated to a linesearch strategy [14, Chap.3]. More sophisticated strategies, based on preconditioning, or subspace acceleration will be discussed in the next section.

**EXAMPLE 1.5** Consider Problem (1.12), where $f$ is a convex lsc function on $\mathbb{R}^N$. For all $\mathbf{v} \in \mathbb{R}^N$, a strongly convex majorant function of $f$ can be derived by setting

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = f(\mathbf{u}) + \frac{1}{2} \|\mathbf{u} - \mathbf{v}\|_{\mathbf{A}}^2, $$

with $\mathbf{A} \in \mathbb{R}^{N \times N}$ a positive definite matrix. The minimizer of the later majorant function is unique, and defines the so-called proximity operator of $f$ at $\mathbf{v}$ within the metric induced by $\mathbf{A}$:

$$ \operatorname{prox}_{f}^{\mathbf{A}}(\mathbf{v}) = \operatorname{argmin}_{\mathbf{u} \in \mathbb{R}^{N}} f(\mathbf{u}) + \frac{1}{2} \|\mathbf{u} - \mathbf{v}\|_{\mathbf{A}}^{2} \tag{1.15} $$
with the simplified notation $\text{prox}_f^I \equiv \text{prox}_f$. The proximity operator benefits from many interesting properties [15]$^{1)}$. In particular, if $\hat{u}$ is a minimizer of $f$, we have the fixed-point relation:

$$ \hat{u} = \text{prox}_f^A(\hat{u}), $$

which is at the core of the convergence analysis of the proximal point algorithm given below:

$$ \begin{array}{l} \mathbf{u}_0 \in \mathbb{R}^N \\ \text{For } k = 0, 1, \dots \\ \quad \lfloor \mathbf{u}_{k+1} = \text{prox}_f^A(\mathbf{u}_k) \rfloor. \end{array} \qquad (1.16) $$

Let $\mathbf{A} = \gamma^{-1}\mathbf{I}$, with $\gamma \in ]0, +\infty[$. According to the definition of the proximity operator, the main iteration (1.16) can be rewritten as:

$$ \mathbf{u}_{k+1} = \mathbf{u}_k - \gamma \mathbf{t}_k, \quad \text{with} \quad \mathbf{t}_k \in \partial f(\mathbf{u}_{k+1}), \qquad (1.17) $$

which (1.17) can be viewed as an implicit subgradient scheme. The main advantage of this implicit scheme is that it allows the stepsize $\gamma$ to be constant along the iterations, whereas the step-size is required to decrease to zero in the case of the standard (explicit) subgradient algorithm [16].

EXAMPLE 1.6 Let us consider the resolution of Problem (1.12) where

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \frac{1}{2} \|\mathbf{H}\mathbf{u} - \mathbf{z}\|_2^2 + \lambda \sum_{n=1}^N \log(\eta + |u_n|) \qquad (1.18) $$

with $\mathbf{H} \in \mathbb{R}^{M \times N}$, $\mathbf{z} \in ]0, +\infty[$, $\mathbf{M}$, $(\lambda, \eta) \in ]0, +\infty[$. Function $u \mapsto \log(\eta + u)$ is concave on $[0, +\infty[$, so according to Property 1.2, for every $\mathbf{v} \in \mathbb{R}^N$, a tangent majorant function of (1.18) is

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = \frac{1}{2} \|\mathbf{H}\mathbf{u} - \mathbf{z}\|_2^2 + \lambda \sum_{n=1}^N \frac{|u_n|}{\eta + |v_n|} + C(\mathbf{v}) \qquad (1.19) $$

where $C(\mathbf{v})$ is a constant such that $h(\mathbf{v}|\mathbf{v}) = f(\mathbf{v})$. The application of the MM algorithm to the majorizing approximation (1.19) leads to the so-called iterative reweighted $\ell_1$ (IRL1) algorithm for sparse signal reconstruction [17, 18, 19, 20]. This name comes from the fact that function (1.19) reads as a least squares term, penalized by a weighted $\ell_1$ norm, whose weights depend on the entries of $\mathbf{v}$.

The following property is at the core of iterative reweighted least-squares algorithms [21, 22, 23, 24, 25], and of half quadratic methods [26, 27, 28, 29, 30, 31, 32] for the minimization of penalized quadratic functions.

PROPERTY 1.7 Let $f: \mathbb{R} \to \mathbb{R}$ be a continuous function such that, for every $u \in \mathbb{R}$, $f(u) = \psi(|u|)$, where:

<footnote n="1">1) See also http://proximity-operator.net</footnote>
(i) $\psi$ is differentiable on $]0, +\infty[$,

(ii) $\psi(\sqrt{\cdot})$ is concave on $]0, +\infty[$,

(iii) $\psi$ is increasing on $]0, +\infty[$.

Then, for all $v \in \mathbb{R}^*$,

$$ (\forall u \in \mathbb{R}) \quad f(u) \le f(v) + \dot{f}(v)(u-v) + \frac{1}{2}\omega(|v|)(u-v)^2 $$

where, for every $\xi \in ]0, +\infty[$, $\omega(\xi) = \dot{\psi}(\xi)/\xi$ and $\dot{f}$ denotes the derivative of the scalar function $f$.

**EXAMPLE 1.8** Let us focus on the search for a minimizer on $\mathbb{R}^N$ of the non smooth function $f: \mathbb{R}^N \to \mathbb{R}$ defined as:

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \sum_{s=1}^{S} \| \mathbf{L}_s \mathbf{u} - \mathbf{c}_s \| \qquad (1.20) $$

with $(\forall s \in \{1, \dots, S\}) \mathbf{L}_s \in \mathbb{R}^{P_s \times N}, \mathbf{c}_s \in \mathbb{R}^{P_s}$. Applying Property 1.7 with $f$ the absolute value function, yields, for every $v \in \mathbb{R}^*$,

$$ (\forall u \in \mathbb{R}) \quad |u| \le |v| + \mathrm{sign}(v)(u-v) + \frac{1}{2|v|}(u-v)^2 = \frac{u^2}{2|v|} + \frac{|v|}{2}. $$

A direct consequence is that, for all $v$ belonging to the set

$$ C = \{v \in \mathbb{R}^N \mid (\forall s \in \{1, \dots, S\}) \mathbf{L}_s v - \mathbf{c}_s \neq 0 \}, $$

the quadratic function

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = \frac{1}{2} \sum_{s=1}^{S} \frac{\|\mathbf{L}_s \mathbf{u} - \mathbf{c}_s\|^2}{\|\mathbf{L}_s \mathbf{v} - \mathbf{c}_s\|} + \|\mathbf{L}_s \mathbf{v} - \mathbf{c}_s\|, $$

is a majorant function for $f$ at $v$. This leads to an MM algorithm for minimizing $f$ usually known as the iterative reweigthed least-squares (IRLS) algorithm. A difficulty is that IRLS needs to be stopped as soon as an iterate does not belong to $C$ (since the majorant function would not be defined anymore). This constitutes a difficulty to prove the convergence of the method. In particular, one must ensure that the algorithm does not stopped before reaching the convergence. An analysis of IRLS convergence properties is available in specific cases, namely if for every $s \in \{1, \dots, S\}$, $\mathbf{L}_s = \mathbf{I}$, IRLS identifying with the Weiszfeld algorithm [23, 33], or when for every $s \in \{1, \dots, S\}$, $P_s = 1$, $\mathbf{L}_s = \mathbf{e}_s^\top$, with $\mathbf{e}_s$ the s-th canonical basis vector of $\mathbb{R}^N$, and $f$ contains an additional quadratic term [21]. Note that, even if the iterates $(\mathbf{u}_k)_{k \in \mathbb{N}}$ stay in $C$, numerical issues may arise if, for some $k$, some component of $(\|\mathbf{L}_s \mathbf{u}_k - \mathbf{c}_s\|)_{1 \le s \le S}$ becomes close to zero. A common strategy adopted in the literature to avoid this problem is to replace the expression of the
majorant function, by a smooth approximation of it, defined on the whole space $\mathbb{R}^N$.
For instance, [34, 35, 12] define, for every $\mathbf{v} \in \mathbb{R}^N$, the smooth majorant function

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad \tilde{h}(\mathbf{u}|\mathbf{v}) = \frac{1}{2} \sum_{s=1}^{S} \frac{\|\mathbf{L}_s \mathbf{u} - \mathbf{c}_s\|^2}{\sqrt{\|\mathbf{L}_s \mathbf{v} - \mathbf{c}_s\|^2 + \eta^2}} + C(\mathbf{v}), $$

where $C(\mathbf{v})$ is constant with respect to $\mathbf{u}$ such that $f(\mathbf{v}) = \tilde{h}(\mathbf{v}|\mathbf{v})$, and $\eta > 0$ is a smoothing parameter. It can be easily shown that this modified IRLS algorithm identifies with an MM algorithm for minimizing a smooth approximation to (1.20), given by

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad \tilde{f}(\mathbf{u}) = \sum_{s=1}^{S} \sqrt{\|\mathbf{L}_s \mathbf{u} - \mathbf{c}_s\|^2 + \eta^2 - \eta}. $$

In practice, parameter $\eta$ is chosen close to zero, so that the minimizers of $f$ and $\tilde{f}$ almost coincide, while numerical issues are avoided.

We finish this section by stating a property resulting from Jensen's inequality which is recalled below:

**LEMMA 1.9** Let $\psi: \mathbb{R} \to \mathbb{R} \cup \{+\infty\}$ be a convex function and let $\omega = (\omega_n)_{1 \le n \le N} \in [0, +\infty)^N$ be such that $\sum_{n=1}^N \omega_n = 1$. Then,

$$ (\forall (u_1, \dots, u_N) \in \mathbb{R}^N) \quad \psi\left(\sum_{n=1}^N \omega_n u_n\right) \le \sum_{n=1}^N \omega_n \psi(u_n). $$

**PROPERTY 1.10** Let $\psi: \mathbb{R} \to \mathbb{R} \cup \{+\infty\}$ be a convex function.

(i) $(\forall (\mathbf{u}, \mathbf{v}, \mathbf{c}) \in (]0, +\infty)^3) \quad \psi(\mathbf{c}^\top \mathbf{u}) \le \sum_{n=1}^N \frac{c_n v_n}{c^\top \mathbf{v}} \psi(\frac{\mathbf{c}^\top \mathbf{v}}{v_n} u_n)$.

(ii) Let $\omega \in ]0, +\infty]^N$ such that $\sum_{n=1}^N \omega_n = 1$. We have

$$ (\forall (\mathbf{u}, \mathbf{v}, \mathbf{c}) \in (] -\infty, +\infty]^3) \quad \psi(\mathbf{c}^\top \mathbf{u}) \le \sum_{n=1}^N \omega_n \psi\left(\frac{c_n}{\omega_n}(u_n - v_n) + \mathbf{c}^\top \mathbf{v}\right). $$

This property is useful in the derivation of several optimization algorithms in the context of signal/image restoration.

**EXAMPLE 1.11** Let us consider the resolution of Problem (1.12) where

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \sum_{m=1}^{M} \psi((\mathbf{H}\mathbf{u})_m; z_m), \quad (1.21) $$
$\mathbf{H} \in \mathbb{R}^{M \times N}$, $\mathbf{z} \in [0, +\infty]^{[M]}$, and $\psi$ is the Kullback-Leibler divergence defined as

$$ (\forall u \in \mathbb{R})(\forall z \in [0, +\infty]) \quad \psi(u; z) = \begin{cases} u - z \log(u), & \text{if } u > 0 \text{ and } z > 0 \\ 0, & \text{if } u = 0 \text{ and } z = 0 \\ +\infty, & \text{elsewhere.} \end{cases} \tag{1.22} $$

Function $f$ is convex but its minimizers have no closed form expression. A strategy is then to build a majorant function of $f$, which is simple to minimize. The Richardson-Lucy algorithm [36] consists in applying the MM method (1.13) where the majorant function is given by

$$ (\forall \mathbf{u} \in ]0, +\infty[) \quad h(\mathbf{u}|\mathbf{v}) = \sum_{n=1}^{N} a_n u_n - b_n(\mathbf{v}) \log(u_n) + c_n(\mathbf{v}), $$

with, for every $n \in \{1, \dots, N\}$,

$$ a_n = [\mathbf{H}^\mathrm{T}\mathbf{1}]_n \\ b_n(\mathbf{v}) = [\mathbf{H}^\mathrm{T}(\mathbf{z} \otimes \mathbf{H}\mathbf{v})]_n v_n \tag{1.23} $$

$$ c_n(\mathbf{v}) = [\mathbf{H}^\mathrm{T}(\mathbf{z} \odot \log(\mathbf{H}\mathbf{v}) \otimes \mathbf{H}\mathbf{v})]_n v_n - [\mathbf{H}^\mathrm{T}(\mathbf{z} \otimes \mathbf{H}\mathbf{v})]_n v_n \log(v_n) $$

(consequence of Property 1.10(i)). This majorant function is easy to optimize since it is separable with respect to $\mathbf{u}$. Its minimization leads to the following multiplicative update rule:

$$ \begin{array}{l} \mathbf{u}_0 \in ]0, +\infty[ \\ \text{For } k = 0, 1, \dots \\ \quad \left[ \begin{array}{l} \text{For } n = 1, 2, \dots, N \\ \quad \left[ u_{n,k+1} = u_{n,k} \frac{[\mathbf{H}^\mathrm{T}(\mathbf{z} \otimes \mathbf{H}\mathbf{u}_k)]_n}{a_n} \right. \\ \quad \left. \mathbf{u}_{n,k} \in ]0, +\infty[ \right. \\ \end{array} \right. \end{array} \tag{1.24} $$

The convergence of the iterations (1.24) has been studied for instance in [37]. Note that extensions of the Richardson-Lucy algorithm to the case of a penalized KL criterion are available (see, for instance, [38]).

### 1.2.3

#### MM quadratic methods

The numerical efficiency of the MM algorithm relies on the use of simple majorant functions providing good approximations to $f$ and whose minimizer is simple to compute. An intuitive choice is thus to choose quadratic majorant functions, which leads to the MM quadratic algorithm [10]. Half-quadratic methods [29, 30, 32, 31] are a particular case of this method.
<figure><img src="image_6.png" /><figcaption>**Figure 1.5** MM quadratic algorithm. Function *f* and a quadratic majorizing function *h*(*v*|*u*<sub>*k*</sub>). The gray interval depicts associated values for *u*<sub>*k*+1</sub> for *θ*<sub>*k*</sub> ∈ [0, 2[. The choice *θ*<sub>*k*</sub> = 1 corresponds to choose *u*<sub>*k*+1</sub> as the minimizer of *h*(*v*|*u*<sub>*k*</sub>)</figcaption></figure>

1.2.3.1 MM quadratic algorithm

Assume that, for every $\mathbf{v} \in \mathbb{R}^N$, there exists a positive definite symmetric matrix
$\mathbf{A}(\mathbf{v}) \in \mathbb{R}^{N \times N}$ such that the quadratic function

$$
(\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = f(\mathbf{v}) + (\mathbf{u}-\mathbf{v})^\top \nabla f(\mathbf{v}) + \frac{1}{2} \|\mathbf{u}-\mathbf{v}\|_{\mathbf{A}(\mathbf{v})}^2 \quad (1.25)
$$

is a majorant function of *f* at *v*. The MM quadratic algorithm then reads:

$$
\begin{equation}
\begin{gathered}
\mathbf{u}_0 \in \mathbb{R}^N \\
\text{For } k = 0, 1, \dots \\
\quad \left[ \begin{array}{@{}l@{}}
\theta_k \in [0, 2[ \\
\mathbf{u}_{k+1} = \mathbf{u}_k - \theta_k \mathbf{A}(\mathbf{u}_k)^{-1} \nabla f(\mathbf{u}_k).
\end{array} \right.
\end{gathered}
\tag{1.26}
\end{equation}
$$

The parameters $(\theta_k)_{k \in \mathbb{N}}$ act as stepsizes. For $\theta_k \equiv 1$, we recover the basic MM algorithm. An illustration of the role of the stepsize is illustrated in the scalar case in Figure 1.5. The simple Algorithm (1.26), which can be viewed as a preconditioned gradient algorithm, for the specific class of preconditioners $(\mathbf{A}(\mathbf{u}_k)^{-1})_{k \in \mathbb{N}}$ [39, Sec.1.3], enjoys the following convergence properties:

THEOREM 1.12 Let $f: \mathbb{R}^N \to \mathbb{R}$ be a differentiable function such that $\lim_{\|\mathbf{u}\| \to +\infty} f(\mathbf{u}) = +\infty$. Assume that there exists $(\underline{v}, \bar{v}) \in [0, +\infty)^2$ such that $(\forall k \in \mathbb{N}) \underline{v}I \preceq \mathbf{A}(\mathbf{u}_k) \preceq \bar{v}I, ^{2)}$ and $(\theta, \bar{\theta}) \in [0, +\infty)^2$ such that, $(\forall k \in \mathbb{N})$ $\underline{\theta} \le \theta_k \le 2 - \bar{\theta}$. Then, the following statements hold:

1) As $k \to +\infty$, $\nabla f(\mathbf{u}_k) \to 0$ and $f(\mathbf{u}_k) \searrow f(\tilde{\mathbf{u}})^{3)}$, where $\tilde{\mathbf{u}}$ is a critical point of $f$.

2) If $f$ is convex, any sequential cluster point of $(\mathbf{u}_k)_{k \in \mathbb{N}}$ is a minimizer of $f$.

3) If $f$ is strictly convex, then $\mathbf{u}_k \to \hat{\mathbf{u}}$ where $\hat{\mathbf{u}}$ is the unique minimizer of $f$.

2) We use $\preceq$ to define the Loewner partial ordering on $\mathbb{R}^{N \times N}$: For every $(\mathbf{A}, \mathbf{B}) \in \mathbb{R}^{N \times N} \times \mathbb{R}^{N \times N}$, $\mathbf{A} \preceq \mathbf{B}$ if and only if, for every $\mathbf{u} \in \mathbb{R}^N$, $\mathbf{u}^\top \mathbf{A} \mathbf{u} \le \mathbf{u}^\top \mathbf{B} \mathbf{u}$.

3) We use the shorter notation $\searrow$ to indicate monotonically decreasing convergence.
EXAMPLE 1.13 Algorithm (1.26) has been applied to the problem of chromatogram baseline correction and noise reduction in [40]. Chromatogram measurements $\mathbf{y} \in \mathbb{R}^N$ are modeled as $\mathbf{y} = \mathbf{u} + \mathbf{v} + \mathbf{w}$ with $\mathbf{u} \in \mathbb{R}^N$ the unknown signal of interest, $\mathbf{v} \in \mathbb{R}^N$ an unknown smooth background also called baseline and $\mathbf{w} \in \mathbb{R}^N$ a noise term. The goal is to retrieve $\mathbf{u}$ under some prior assumptions on $\mathbf{u}$ and $\mathbf{v}$. In [40], the authors proposed the following penalized cost function

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \frac{1}{2} \|\mathbf{H}(\mathbf{y} - \mathbf{u})\|^2 + \sum_{i=0}^{M} \lambda_i \phi(D_i \mathbf{u}), \quad (1.27) $$

with $\mathbf{H}$ a suitable high-pass filter, and, for every $i = 0, \dots, M$, $D_i$ the i-th order difference operator, $\lambda_i > 0$ and $\phi$ a smooth even potential function. The minimization of $f$ is performed by Algorithm (1.26), leading to so-called Baseline Estimation And Denoising with Sparsity (BEADS) method. An illustration of result obtained by BEADS on real chromatogram data is displayed in Figure 1.6.

<figure><img src="image_7.png" /><figcaption>Figure 1.6 Baseline correction using BEADS. Baselines estimated in left and estimated peaks (obtained by subtraction of estimated baseline from data) in right.</figcaption></figure>

EXAMPLE 1.14 Let us consider the resolution of Problem (1.12) where

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \frac{1}{2} \|\mathbf{H}\mathbf{u} - \mathbf{z}\|^2, \quad (1.28) $$

$\mathbf{H} \in \mathbb{R}^{M \times N}$, $\mathbf{z} \in \mathbb{R}^M$. According to Property 1.10(ii), for every $\mathbf{v} \in \mathbb{R}^N$, a separable majorant function of $f$ at $\mathbf{v}$ is:

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h(\mathbf{u}|\mathbf{v}) = f(\mathbf{v}) + \nabla f(\mathbf{v})^\mathrm{T}(\mathbf{u} - \mathbf{v}) + \frac{1}{2} \sum_{n=1}^{N} a_n (u_n - v_n)^2, $$

with, for every $n \in \{1, \dots, N\}$,

$$ a_n = [|\mathbf{H}|^\mathrm{T}|\mathbf{H}| \mathbf{1}]_n, $$

where $|\mathbf{H}|$ denotes a matrix with same size than $\mathbf{H}$, whose entries are equal to the absolute value of those of $\mathbf{H}$. In that case, Algorithm (1.30) becomes equivalent to the simultaneous algebraic reconstruction technique (SART) [41] for tomography reconstruction:

$$ \begin{array}{l} \mathbf{u}_0 \in \mathbb{R}^N \\ \text{For } k = 0, 1, \dots \\ \quad \left[ \begin{array}{l} \text{For } n = 1, 2, \dots, N \\ \quad \left[ u_{n,k+1} = u_{n,k} - \frac{[\mathbf{H}^\mathrm{T}(\mathbf{H}\mathbf{u}_k - \mathbf{z})]_n}{a_n} \right]. \end{array} \right. \end{array} \quad (1.29) $$
Note that the separable form of the quadratic majorant function in that case leads to a fully parallelizable algorithm, in the sense that each component of $\mathbf{u}$ is updated in a parallel manner. Generalizations of the SART method to the minimization of a smoothed regularized least squares criterion can be found for instance in [42, 43].

**Remark 1.1** In the context of large scale optimization, the exact inversion of the majorant matrix involved at each iteration of Algorithm (1.26) can be computationally intensive. In practice, one can instead solve approximately the associated linear system using a linear solver, such as the conjugate gradient algorithm. A nice feature of Algorithm (1.26) is that its convergence properties are not modified if we replace the update by

$$ \text{For } k = 0, 1, \dots \left[ \begin{array}{l} \theta_k \in ]0, 2[ \\ \mathbf{u}_{k+1} = \mathbf{u}_k - \theta_k \mathbf{d}_k. \end{array} \right. \qquad (1.30) $$

where $\mathbf{d}_k$ results from $J_k \ge 1$ iterations of a conjugate gradient algorithm (possibly preconditioned) applied to the system $A(\mathbf{u}_k)\mathbf{d} = \nabla f(\mathbf{u}_k)$ [44]. Surprisingly, the best performance in practice for the modified Algorithm (1.30) are achieved for small values of $J_k$ (typically, $J_k = 5$).

**Remark 1.2** Let us consider the minimization of a twice differentiable, strictly convex function $f: \mathbb{R}^N \to \mathbb{R}$, by using the Newton algorithm:

$$ \text{For } k = 0, 1, \dots \left[ \begin{array}{l} \mathbf{u}_{k+1} = \mathbf{u}_k - (\nabla^2 f(\mathbf{u}_k))^{-1} \nabla f(\mathbf{u}_k). \end{array} \right. \qquad (1.31) $$

The strict convexity of $f$ is however not sufficient to guarantee the convergence of Algorithm (1.31), nor to ensure the monotonicity of the sequence $(f(\mathbf{u}_k))_{k \in \mathbb{N}}$, which makes the Newton method with unit stepsize inadequate for general use [10]. The convergence of the algorithm can actually be ensured, according to Theorem 1.12, if, for every $\mathbf{u} \in \mathbb{R}^N$, the following majorization holds:

$$ f(\mathbf{u}) \le f(\mathbf{v}) + (\mathbf{u} - \mathbf{v})^{\mathrm{T}} \nabla f(\mathbf{v}) + \frac{1}{2} \|\mathbf{u} - \mathbf{v}\|_{\nabla^2 f(\mathbf{v})}^2. \qquad (1.32) $$

However, this assumption is rarely satisfied in practice. Three practical strategies to secure the convergence of Algorithm (1.31) are discussed in [14, Chap.6],[39, Sec.1.4]: (i) the introduction of a backtracking linesearch ensuring that the Newton step satisfies some sufficient decrease condition, (ii) the addition of a definite positive matrix, possibly varying along the iterations, to the Hessian matrix in order to satisfy a majorizing condition similar to (1.32), or (iii) the combination of the previous technique with a trust region approach that check the majorization inequality only within a neighborhood of the current iterate.
#### 1.2.3.2 Half-quadratic algorithms

Half-quadratic methods [26, 27, 28, 29, 30, 31, 32] can be interpreted as MM quadratic algorithms for minimizing on $\mathbb{R}^N$ the following class of functions:

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \frac{1}{2} \|\mathbf{H}\mathbf{u} - \mathbf{z}\|^2 + \sum_{s=1}^{S} \psi(\|\mathbf{L}_s \mathbf{u} - \mathbf{c}_s\|) \quad (1.33) $$

with $\mathbf{H} \in \mathbb{R}^{M \times N}$, $\mathbf{z} \in \mathbb{R}^M$, $\psi: \mathbb{R} \to \mathbb{R}$, and $(\forall s \in \{1, \dots, S\}) \mathbf{L}_s \in \mathbb{R}^{N \times P_s}$, $\mathbf{c}_s \in \mathbb{R}^{P_s}$. Two particular choices for matrix $\mathbf{A}(\cdot)$ are employed in this context, depending on the assumptions on $\psi$.

The first strategy [30], relies on the assumption that $\psi$ has a $\beta$-Lipschitzian gradient on $\mathbb{R}$. Let us define

$$ (\forall \mathbf{v} \in \mathbb{R}^N) \quad \mathbf{A}(\mathbf{v}) = \mathbf{H}^\mathrm{T}\mathbf{H} + \mu \mathbf{L}^\mathrm{T}\mathbf{L}, \quad \mu \ge \beta, \quad (1.34) $$

with $\mathbf{L} = [\mathbf{L}_1^\mathrm{T} | \dots | \mathbf{L}_S^\mathrm{T}]^\mathrm{T} \in \mathbb{R}^{P \times N}$, $P = \sum_{s=1}^S P_s$. Then, according to Property 1.3(i), (1.25) is a majorant function of (1.33).

The second strategy [29] assumes that $\psi$ is differentiable on $]0, +\infty[$, $\psi(\sqrt{\cdot})$ is concave, $(\forall x \in [0, +\infty[$), $\dot{\psi}(x) \ge 0$, and $\lim_{x \to 0^+} (\omega(x) = \frac{\dot{\psi}(x)}{x}) \in \mathbb{R}$. Then, Property 1.7 implies that (1.25) is a majorant function of (1.33), if

$$ (\forall \mathbf{v} \in \mathbb{R}^N) \quad \mathbf{A}(\mathbf{v}) = \mathbf{H}^\mathrm{T}\mathbf{H} + \mathbf{L}^\mathrm{T} \mathrm{Diag} \left( (\omega(\|\mathbf{L}_s \mathbf{v} - \mathbf{c}_s\|)\mathbf{1})_{1 \le s \le S} \right) \mathbf{L}. \quad (1.35) $$

The later approach leads to a better convergence rate in theory [26, 32], but its computational cost may be higher. Indeed, in this case, the majorant matrix (1.35), as well as its inverse, must be recomputed at each iteration, which can be prohibitive in the context of large scale problems. Figure ?? illustrates on a simple scalar example the two different approaches for half-quadratic majorant construction.

The term *half-quadratic* stems from the fact that the application of the MM quadratic algorithm to the minimization of (1.33), with majorant matrix (1.34) or (1.35) can be written as an alternating minimization algorithm:

$$ \text{For } k = 0, 1, \dots \quad \left| \begin{array}{l} \mathbf{b}_k \in \underset{\mathbf{b} \in \mathcal{D}}{\operatorname{Argmin}} \Phi(\mathbf{u}, \mathbf{b}) \\ \bar{\mathbf{u}}_k \in \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{Argmin}} \Phi(\mathbf{u}, \mathbf{b}) \\ \mathbf{u}_{k+1} = \theta_k \bar{\mathbf{u}}_k + (1 - \theta_k) \mathbf{u}_k, \quad \theta_k \in ]0, 2[, \end{array} \right. \quad (1.36) $$

where $\Phi: \mathbb{R}^N \times \mathcal{D} \to (\mathbb{R} \cup \{+\infty\}) \times (\mathbb{R} \cup \{+\infty\})$ is such that for every $\mathbf{b} \in \mathcal{D}$, $\mathbf{u} \mapsto \Phi(\mathbf{u}, \mathbf{b})$ is a *quadratic* function.

#### 1.2.3.3 Subspace acceleration strategy

In the context of large scale optimization, the minimization of a quadratic majorant function at each iteration may become untractable. The main idea of subspace acceleration is to restrict the minimization space to a subspace spanned by a small
<figure><img src="image_8.png" /><figcaption>**Figure 1.7** Half-quadratic majorant function. Function $f(u) = \sqrt{1 + u^2}$ and quadratic majorizing functions $h(\cdot|v)$ for different values of the tangency point $v$, built according to the strategy (1.34) (left) and (1.35) (right). All the majorizing functions obtained from (1.34) share the same curvature, corresponding to the maximum curvature of $f$ (here, obtained at $u = 0$). In contrast, the majorizing functions resulting from (1.35) are even functions, with curvature depending on the tangency point $v$.</figcaption></figure>

number of vectors, instead of minimizing the majorant over the whole space $\mathbb{R}^N$.
The general form of an MM subspace algorithm is thus

$$ (\forall k \in \mathbb{N}) \quad \mathbf{u}_{k+1} \in \underset{\mathbf{u} \in \operatorname{span}\{\mathbf{d}_k^1, \mathbf{d}_k^2, \dots, \mathbf{d}_k^{M_k}\}}{\operatorname{Argmin}} h(\mathbf{u}|\mathbf{u}_k), $$

where, for every $k \in \mathbb{N}$, $M_k \ge 1$. Equivalently, we obtain:

$$ \begin{array}{l} \mathbf{u}_0 \in \mathbb{R}^N \\ \text{For } k = 0, 1, \dots \\ \quad \left[ \begin{array}{l} \text{Choose } \mathbf{D}_k \in \mathbb{R}^{N \times M_k}, \\ \quad s_k \in \underset{s \in \mathbb{R}^{M_k}}{\operatorname{Argmin}} h(\mathbf{u}_k + \mathbf{D}_k s | \mathbf{u}_k), \\ \quad \mathbf{u}_{k+1} = \mathbf{u}_k + \mathbf{D}_k s_k, \end{array} \right. \tag{1.37} \end{array} $$

where, for every $k \in \mathbb{N}$, $\mathbf{D}_k = [\mathbf{d}_k^1 | \mathbf{d}_k^2 | \dots | \mathbf{d}_k^{M_k}] \in \mathbb{R}^{N \times M_k}$ gathers some search directions (see [45, Tab.1] for examples) and $s_k \in \mathbb{R}^{M_k}$ is a multi-dimensional stepsize allowing to optimally combined these directions. By adopting a subspace spanned by a small number of directions, one can expect to reduce the computational cost of the algorithm.

Assume that, for every $v \in \mathbb{R}^N$, there exists a strongly positive self-adjoint operator $\mathbf{A}(v)$ such that the quadratic function (1.25) is a majorant function of $f$ at $v$. The MM subspace quadratic algorithm reads as follows, where $\dagger$ denotes the
pseudo-inverse operation:

$$
\begin{array}{ll}
u_0 \in \mathbb{R}^N & \\
\text{For } k = 0, 1, \dots & \\
\quad \text{Choose } D_k \in \mathbb{R}^{M_k}, & (1.38) \\
\quad \left\{
\begin{array}{l}
s_k = -(D_k^\top A(u_k) D_k)^\dagger (D_k^\top \nabla f(u_k)), \\
u_{k+1} = u_k + D_k s_k.
\end{array}
\right. &
\end{array}
$$

A convergence result of Algorithm (1.38) similar to Theorem 1.12 has been es-
tablished in [45], under the assumption that, for every $k \in \mathbb{N}$, the gradient direction
$\nabla f(u_k)$ belongs to the vector subspace spanned by the columns of $D_k$. This re-
sult has recently been extended to the case of non-convex functions [46], where the
convergence of the iterates of Algorithm (1.38) to a critical point of the objective
function is proved. Note that Algorithm (1.38) has recently been extended to the on-
line case when only a stochastic approximation of the criterion is employed at each
iteration [47].

The MM subspace algorithm shows very good performance in practice, when
compared with nonlinear conjugate gradient algorithms [48], low memory BFGS
strategy [49], and also with graph-cut based discrete optimization methods, as well
as proximal algorithms [45, 50, 51]. All the related works illustrate the fact that the
choice of the subspace has a major impact on the practical convergence speed of
Algorithm (1.38) (see, for instance [45, Section 5], [46, Section 5.1]). In particular,
it seems that the worse performance is obtained in the case of a gradient-like algo-
rithm, i.e., when, for every $k \in \mathbb{N}$, $M_k = 1$ and $D_k = -\nabla f(u_k)$. On the opposite,
when the search subspace is the full space, Algorithm (1.38) becomes equivalent to
the MM quadratic method mentioned in Section 1.2.3.1. This methods can converge
in a very small number of iterations, but each of them requires the inversion of an
$N \times N$ matrix, which may have a high computational cost. A good compromise
is obtained for the memory gradient subspace [52], spanned by the current gradient
and the previous direction, leading to the so-called MM Memory Gradient (3MG)
algorithm. These practical observations are confirmed by the theoretical analysis of
the rate of convergence of Algorithm (1.38) [53].

1.2.4

Variable metric forward-backward algorithm

Let us now focus on the minimization of $f = f_1 + f_2$ with $f_1$ differentiable and $f_2$
convex non necessarily differentiable (for instance, a regularization term enforcing
sparsity). Assume that, for every $v \in \mathbb{R}^N$, there exists a positive definite matrix
$A(v) \in \mathbb{R}^{N \times N}$ such that the quadratic function

$$
(\forall u \in \mathbb{R}^N) \quad h(u|v) = f_1(v) + \nabla f_1(v)^\top (u - v) + \frac{1}{2} \|u - v\|_{A(v)}^2
$$
is a majorant function of $f_1$ at $v$. The variable metric forward-backward (VMFB) algorithm [54, 55] reads

$$ 
\begin{array}{l}
u_0 \in \mathbb{R}^N \\
\text{For } k = 0, 1, \dots \\
\quad \theta_k \in [0, 2[, \\
\quad \mathbf{u}_{k+1} = \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{prox}} \theta_k^{-1} \mathbf{A}(\mathbf{u}_k) (\mathbf{u}_k - \theta_k \mathbf{A}(\mathbf{u}_k)^{-1} \nabla f_1(\mathbf{u}_k)),
\end{array}
\tag{1.39} $$

where the proximity operator has been defined in (1.15). Algorithm (1.39) can actually be viewed as a relaxed form of an MM algorithm. Let $\theta_k \equiv 1$. According to the definition of the proximity operator,

$$ 
\begin{align*}
\mathbf{u}_{k+1} &= \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \|\mathbf{u} - \mathbf{u}_k + \mathbf{A}(\mathbf{u}_k)^{-1} \nabla f_1(\mathbf{u}_k)\|_{\mathbf{A}(\mathbf{u}_k)}^2 + f_2(\mathbf{u}) \\
&= \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{argmin}} (\mathbf{u} - \mathbf{u}_k)^{\top} \nabla f_1(\mathbf{u}_k) + \frac{1}{2} \|\mathbf{u} - \mathbf{u}_k\|_{\mathbf{A}(\mathbf{u}_k)}^2 + f_2(\mathbf{u}) \\
&= \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{argmin}} h(\mathbf{u}|\mathbf{u}_k) + f_2(\mathbf{u}),
\end{align*}
$$

which identifies with a MM algorithm to minimize $f$.

A convergence result similar to Theorem 1.12 can be stated for the VMFB algorithm, under the assumption that $f_2$ is convex, with non empty domain, and continuous on its domain. An extension of the convergence properties of VMFB in the context of non convex optimization can be found in [55].

**Remark 1.3** Assume that $f_1$ is $\beta$-Lipschitz differentiable. According to Property 1.3(i), a valid choice is $\mathbf{A}(\mathbf{u}_k) \equiv \beta^{-1}\mathbf{I}$. In that case, VMFB algorithm becomes equivalent to the forward-backward algorithm [56] given below

$$ 
\begin{array}{l}
u_0 \in \mathbb{R}^N \\
\text{For } k = 0, 1, \dots \\
\quad \theta_k \in [0, 2\beta^{-1}], \\
\quad \mathbf{u}_{k+1} = \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{prox}}_{\theta_k f_2} (\mathbf{u}_k - \theta_k \nabla f_1(\mathbf{u}_k)).
\end{array}
\tag{1.40} $$

The iterative soft thresholding algorithm (ISTA) [57] is a special case of this algorithm in the context of sparse signal restoration when $f_2$ is a (possibly weighted) $\ell_1$ norm.

**Remark 1.4** There exist plenty of properties for the proximity operator, so that a number of functions have a closed form expression for their prox [15]$^4$. Moreover, if $f$ is a separable function, i.e.

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \sum_{n=1}^{N} f_n(u_n), \tag{1.41} $$

with $(f_n)_{1 \le n \le N}$ functions defined on $\mathbb{R}$, then

$$ \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{prox}} \operatorname{Diag}\{\mathbf{a}\}(\mathbf{u}) = \left( \underset{\mathbf{u} \in \mathbb{R}^N}{\operatorname{prox}}_{a_n^{-1} f_n}(u_n) \right)_{1 \le n \le N}. \tag{1.42} $$

<footnote n="4">4) See also the web repository http://proximity-operator.net/.</footnote>
Note that the later property is the cornerstone of the so-called separable surrogate functionals algorithm proposed in [42] in the context of sparse signal reconstruction, which identifies with a VMFB algorithm with a diagonal metric.

When the form of the function and/or the metric is more involved, one must resort to an iterative strategy in order to compute the proximity step. Attention must be paid to this problem since the overall computation cost of the optimization method becomes then strongly dependent on the efficiency of the subiterations implemented for computing the proximity operator. A very efficient algorithm for performing the inner loop is the dual forward-backward algorithm [58] (sometimes also known as the dual ascent algorithm) which simply consists in applying the above forward-backward algorithm to the dual of the proximal problem. Note that several recent works [59, 60, 61, 62] have proposed block-alternating, parallel, or distributed versions of this method, in order to reduce its computational and memory cost in the context of large scale optimization problems.

**EXAMPLE 1.15** Assume that

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f_2(\mathbf{u}) = \iota_C(\mathbf{u}), \qquad (1.43) $$

where $C$ is a closed non empty convex subset of $\mathbb{R}^N$. Then, Algorithm (1.39) becomes equivalent to the scaled projected gradient algorithm

$$ \begin{array}{l} \mathbf{u}_0 \in \mathbb{R}^N \\ \text{For } k = 0, 1, \dots \\ \quad \left[ \begin{array}{l} \theta_k \in ]0, 2[, \\ \tilde{\mathbf{u}}_k = \mathbf{u}_k - \theta_k \mathbf{A}(\mathbf{u}_k)^{-1} \nabla f_1(\mathbf{u}_k), \\ \mathbf{u}_{k+1} = \underset{v \in C}{\operatorname{argmin}} \|v - \tilde{\mathbf{u}}_k\|_{\mathbf{A}(\mathbf{u}_k)}. \end{array} \right. \end{array} \qquad (1.44) $$

The convergence of Algorithm (1.44) has been studied, for instance, in [63, 64] under various assumptions on the metric matrices, and the stepsize values. Note that the standard projected gradient [65] is recovered if the preconditioning matrices are scaled versions of the identity matrix.

**EXAMPLE 1.16** The recent work [66, 67] presents an application of the variable metric forward-backward approach in the context of mass spectrometry (MS). MS is a fundamental technology of analytical chemistry. In the context of protein analysis, the MS spectra $\mathbf{y} \in \mathbb{R}^M$ can be written as $\mathbf{y} = \mathbf{K}\mathbf{u} + \mathbf{w}$ with $\mathbf{K} \in \mathbb{R}^{M \times N}$ the MS average dictionary, $\mathbf{u} \in \mathbb{R}^N$ a sparse positive-valued signal to be estimated and $\mathbf{w}$ a noise term. The dictionary size $N$ is simply equal to $N = MZ$, with $M$ the grid size for the mass and $Z$ the number of charges to be explored. The estimation of $\mathbf{u}$ can be performed by minimizing $f = f_1 + f_2$ with $f_1$ a smooth term favorizing the sparsity of $\mathbf{u}$ and $f_2(\mathbf{u}) = \iota_C(\mathbf{u})$ with

$$ C = \{ \mathbf{u} \in \mathbb{R}^N \mid \| \mathbf{K} \mathbf{u} - \mathbf{y} \| \le \xi \text{ and } \mathbf{u} \ge 0 \} \qquad (1.45) $$
where $\xi > 0$ depends of the noise level. The authors in [66, 67] introduced the so-called SPOQ ((Smoothed p-Over-q) penalty:

$$f_1(\mathbf{u}) = \log \left( \frac{(\ell_{p,\alpha}^{p}(\mathbf{u}) + \beta^p)^{1/p}}{\ell_{q,\eta}(\mathbf{u})} \right), \qquad (1.46)$$

where $\beta > 0$, $\ell_{p,\alpha}$ and $\ell_{q,\eta}$ are smoothed versions of $\ell_p$ and $\ell_q$ norms parametrized by $(\alpha, \eta) \in ]0, +\infty[$. SPOQ penalty consists of a non-convex Lipschitz-differentiable surrogate for $\ell_p$-over-$\ell_q$ quasi-norm/norm ratios, with the ability to enforce sparsity property while preserving the signal scale. An accelerated version of Algorithm (1.39), is then proposed to minimize efficiently the resulting non-convex non-smooth objective function $f$, where the majorization condition is relaxed, using a trust-region based strategy. Figure 1.8 illustrates the results obtained when using SPOQ procedure to restore synthetic MS spectra.

<figure><img src="image_9.png" /><figcaption>Figure 1.8 Original sparse signal (top), observed measurements (middle) and restored signal (bottom) using SPOQ with $p = 0.75$ and $q = 2$, for a synthetic dataset with $M = 1000$ and mono-charged profile, i.e. $Z = 1$.</figcaption></figure>

### 1.2.5

#### Block-coordinate MM algorithms

##### 1.2.5.1 General principle

As emphasized in Section 1.1, most of the objective functions encountered in the context of source separation involve blocks of variables (very often, two blocks) with different significance from the application standpoint, which makes the resulting optimization problems intrinsically adapted to block alternating minimization methods. Such methods can also be very appealing in the context of large scale problems when the whole variable cannot be stored entirely in the local memory of
the machine, so that only update subsets of it can be loaded and updated at each iteration of the algorithm. The aim of this section is to show how MM algorithms can be combined with block alternating schemes.

For every $\mathbf{v} \in \mathbb{R}^N$, for every $j \in \{1, \dots, J\}$, assume that there exists a majorant function $h_j(\cdot|\mathbf{v}) : \mathbb{R}^{N_j} \to \mathbb{R}$, which majorizes the restriction of $f$ to its $j$-th block at $\mathbf{v}^{(j)}$, i.e. $h_j(\mathbf{v}^{(j)}|\mathbf{v}) = f(\mathbf{v})$ and,

$$ (\forall \mathbf{s} \in \mathbb{R}^{N_j}) \quad h_j(\mathbf{s}|\mathbf{v}) \le f(\mathbf{v}^{(1)}, \dots, \mathbf{v}^{(j-1)}, \mathbf{s}, \mathbf{v}^{(j+1)}, \dots, \mathbf{v}^{(J)}). $$

Then, a block-coordinate MM algorithm for minimizing $f$ reads:

$$ \begin{array}{ll} \mathbf{u}_0 \in \mathbb{R}^N & \\ \text{For } k = 0, 1, \dots & \\ \quad \left[ \begin{array}{ll} \text{Select } j_k \in \{1, \dots, J\}, & \\ \mathbf{u}_{k+1}^{(j_k)} = \underset{\mathbf{s} \in \mathbb{R}^{N_{j_k}}}{\operatorname{argmin}} h_j(\mathbf{s}|\mathbf{u}_k), & (1.47) \\ \mathbf{u}_{k+1}^{(\bar{j}_k)} = \mathbf{u}_k^{(j_k)}. & \end{array} \right. & \end{array} $$

where $\bar{j}_k = \{1, \dots, J\} \setminus \{j_k\}$. When $J = 1$, one recovers the standard MM algorithm (1.13). For $J > 1$, at each iteration $k \in \mathbb{N}$, $j_k \in \{1, \dots, J\}$ can be chosen according to:

* the cyclic rule:

$$ (\forall k \in \mathbb{N}) \quad j_k - 1 = k \mod (J), $$

* a quasi-cyclic rule:

There exists a constant $K \ge J$ such that, for every $k \in \mathbb{N}$,

$$ \{1, \dots, J\} \subset \{j_k, \dots, j_{k+K-1}\}, $$

* a random rule:

For every $k \in \mathbb{N}$, $j_k$ is a realization of a random variable, for example uniformly distributed on $\{1, \dots, J\}$.

The convergence properties of the general scheme (1.47), also known as *block successive upper-bound minimization* [12], have been studied in [68, 69] under various assumptions on the function $f$, on the majorant function sequence, and on the block selection rule. Discussions regarding the practical implementation of Algorithm (1.47) can be found for instance in [70] in the context of image reconstruction and in [11] in the context of machine learning.

### 1.2.5.2 Block-coordinate quadratic MM algorithm

Let us now describe a block-coordinate version of the quadratic MM algorithm (1.26). For every $\mathbf{v} \in \mathbb{R}^N$, for every $j \in \{1, \dots, J\}$, assume that there exists a symmetric positive definite matrix $A_j(\mathbf{v}) \in \mathbb{R}^{N_j \times N_j}$ such that the quadratic function

$$ (\forall \mathbf{s} \in \mathbb{R}^{N_j}) \quad h_j(\mathbf{s}|\mathbf{v}) = f(\mathbf{v}) + (\mathbf{s} - \mathbf{v}^{(j)})^\mathrm{T} \nabla_j f(\mathbf{v}) + \frac{1}{2} \|\mathbf{s} - \mathbf{v}^{(j)}\|_{A_j(\mathbf{v})}^2 $$
is a majorant function at $v^{(j)}$ of the restriction of $f$ to its $j$-th block. Then, the block-coordinate MM quadratic algorithm reads:

$$ 
\begin{array}{l}
u_0 \in \mathbb{R}^N \\
\text{For } k = 0, 1, \dots \\
\quad \left\{
\begin{array}{l}
\text{Select } j_k \in \{1, \dots, J\}, \\
\theta_k \in ]0, 2[, \\
\mathbf{u}_{k+1}^{(j_k)} = \mathbf{u}_k^{(j_k)} - \theta_k \mathbf{A}_{j_k}(\mathbf{u}_k)^{-1} \nabla_{j_k} f(\mathbf{u}_k), \\
\mathbf{u}_{k+1}^{(\bar{j}_k)} = \mathbf{u}_k^{(\bar{j}_k)}.
\end{array}
\right.
\end{array}
\tag{1.48} $$

For $J = 1$, we recover Algorithm (1.26). For $J > 1$, the convergence properties of Algorithm (1.48), also known as coordinate-ascent paraboloidal surrogate method [70], depend on the block selection rule. For example, in the case of a deterministic quasi-cyclic rule, a convergence result similar to Theorem 1.12 holds, under the assumption that there exists $(\underline{\nu}, \bar{\nu}) \in ]0, +\infty[^{2}$ such that $(\forall k \in \mathbb{N})$ $\underline{\nu}I \preceq A_{j_k}(\mathbf{u}_k) \preceq \bar{\nu}I$.

**EXAMPLE 1.17** Let us consider the problem of minimizing the following function:

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \frac{1}{2} \|\mathbf{H}\mathbf{u} - \mathbf{z}\|^2 + \sum_{n=1}^{N} \sqrt{1 + u_n^2}, \tag{1.49} $$

with $\mathbf{H} \in \mathbb{R}^{M \times N}$ and $\mathbf{z} \in \mathbb{R}^M$. According to Property (1.7), for every $\mathbf{v} \in \mathbb{R}^N$,

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) \le f(\mathbf{v}) + \nabla f(\mathbf{v})^\top (\mathbf{u} - \mathbf{v}) + \frac{1}{2} (\mathbf{u} - \mathbf{v})^\top \mathbf{A}(\mathbf{v}) (\mathbf{u} - \mathbf{v}), \tag{1.50} $$

with

$$ (\forall \mathbf{v} \in \mathbb{R}^N) \quad \mathbf{A}(\mathbf{v}) = \mathbf{H}^\top \mathbf{H} + \text{Diag}\{(1 + v_n^2)_{1 \le n \le N}\}^{-1/2}. \tag{1.51} $$

Consequently, for every $\mathbf{v} \in \mathbb{R}^N$, a majorant approximation of the restriction of $f$ to the $n$-th component of $\mathbf{u}$, i.e. $u \mapsto f([v_1, \dots, v_{n-1}, u, v_n, \dots, v_N]^\top)$, is:

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad h_n(u_n|\mathbf{v}) = f(\mathbf{v}) + \nabla_n f(\mathbf{v})(u_n - v_n) + \frac{a_n(\mathbf{v})}{2}(u_n - v_n)^2, \tag{1.52} $$

with

$$ (\forall \mathbf{v} \in \mathbb{R}^N)(\forall n \in \{1, \dots, N\}) \quad a_n(\mathbf{v}) = [\mathbf{H}^\top \mathbf{H}]_{n,n} + (1 + v_n^2)^{-1/2}, \tag{1.53} $$

and $\nabla_n f(\mathbf{v})$ the $n$-th component of $\nabla f(\mathbf{u})$. We can deduce the following alternating algorithm for minimizing $f$, referred to as the coordinate ascent with parabola
surrogate in [70]:

$$
\begin{array}{ll}
\mathbf{u}_0 \in \mathbb{R}^N & \\
\text{For } k = 0, 1, \dots & \\
\quad \left\{
\begin{array}{l}
\quad \text{Select } n \in \{1, \dots, N\}, \\
\quad \theta_k \in ]0, 2[, \\
\quad u_{n,k+1} = u_{n,k} - \theta_k a_n(\mathbf{u}_k)^{-1} \nabla_n f(\mathbf{u}_k), \\
\quad (\forall n' \in \{1, \dots, N\} \setminus \{n\}) \quad u_{n',k+1} = u_{n',k}.
\end{array}
\right. & (1.54)
\end{array}
$$

1.2.5.3 Block-coordinate VMFB algorithm

The VMFB algorithm (1.39) can also be extended to a block alternating version of
it. Consider the minimization of $f = f_1 + f_2$ with $f_1$ differentiable and $f_2$ convex
non necessarily differentiable. Assume additionally that $f_2$ is block-separable, i.e.

$$
(\forall \mathbf{u} = (\mathbf{u}^{(j)})_{1\leq j \leq J}) \quad f_2(\mathbf{u}) = \sum_{j=1}^{J} f_{2,j}(\mathbf{u}^{(j)}).
$$

For every $\mathbf{v} \in \mathbb{R}^N$, for every $j \in \{1, \dots, J\}$, assume that there exists a symmetric positive definite matrix $\mathbf{A}_j(\mathbf{v}) \in \mathbb{R}^{N_j \times N_j}$ such that the quadratic function

$$
(\forall \mathbf{u}^{(j)} \in \mathbb{R}^{N_j}) h_j(\mathbf{u}^{(j)} | \mathbf{v}) = f_1(\mathbf{v}) + (\mathbf{u}^{(j)} - \mathbf{v}^{(j)})^\mathrm{T} \nabla_j f_1(\mathbf{v}) + \frac{1}{2} \|\mathbf{u}^{(j)} - \mathbf{v}^{(j)}\|_{\mathbf{A}_j(\mathbf{v})}^2
$$

is a majorant function at $\mathbf{v}^{(j)}$ of the restriction of $f_1$ to its $j$-th block. Then, the
block-coordinate VMFB (BC-VMFB) algorithm reads:

$$
\begin{array}{ll}
\mathbf{u}_0 \in \mathbb{R}^N & \\
\text{For } k = 0, 1, \dots & \\
\quad \left\{
\begin{array}{l}
\text{Select } j_k \in \{1, \dots, J\}, \\
\theta_k \in ]0, 2[, \\
\mathbf{u}_{k+1}^{(j_k)} = \text{prox}_{f_{2,j_k}}^{\theta_k^{-1} \mathbf{A}_{j_k}(\mathbf{u}_k)} (\mathbf{u}_k^{(j_k)} - \theta_k \mathbf{A}_{j_k}(\mathbf{u}_k)^{-1} \nabla_{j_k} f_1(\mathbf{u}_k)), \\
\mathbf{u}_{k+1}^{(\bar{j}_k)} = \mathbf{u}_k^{(\bar{j}_k)}.
\end{array}
\right. & (1.55)
\end{array}
$$

Like its non alternating version, this algorithm can be understood as a relaxed MM
algorithm. Its convergence properties have been analyzed in the non convex case in
[71]. In the convex case, the convergence properties are the same as for the VMFB
algorithm, assuming that the update rule is quasi-cyclic.

EXAMPLE 1.18 *Let us consider the minimization of*

$$
\underset{\mathbf{A} \in \mathbb{R}^{P \times R}, \mathbf{S} \in \mathbb{R}^{R \times T}}{\operatorname{minimize}} F(\mathbf{A}, \mathbf{S}) = \frac{1}{2} \|\mathbf{A}\mathbf{S} - \mathbf{X}\|_{F}^{2} + \iota_{\mathbf{A}}(\mathbf{A}) + \iota_{\mathbf{S}}(\mathbf{S}), \quad (1.56)
$$

where $\mathcal{A} = [a_{\min}, a_{\max}]^{P \times R}$, with $0 < a_{\min} \le a_{\max}$, and $\mathcal{S} = [s_{\min}, s_{\max}]^{R \times T}$, with $0 < s_{\min} \le s_{\max}$. In order to apply Algorithm (1.55), we need to define, for a
given $(\mathbf{A}', \mathbf{S}') \in \mathcal{A} \times \mathcal{S}$, quadratic majorants of $F(\cdot, \mathbf{S}')$ (resp. $F(\mathbf{A}', \cdot)$). Since matrices $\mathbf{A}' \in \mathcal{A}$ and $\mathbf{S}' \in \mathcal{S}$ have positive elements, we can derive the following quadratic majorant functions (consequence of Lemma 1.10)

$$ (\forall \mathbf{A} \in \mathcal{A}) \quad F(\mathbf{A}, \mathbf{S}') \le F(\mathbf{A}', \mathbf{S}') + \operatorname{tr} \left( (\mathbf{A} - \mathbf{A}') \nabla_1 F(\mathbf{A}', \mathbf{S}')^\mathrm{T} \right) \\ \qquad + \frac{1}{2} \operatorname{tr} \left( (((\mathbf{A} - \mathbf{A}') \odot (\mathbf{A}'(\mathbf{S}'\mathbf{S}'^\mathrm{T}))) \otimes \mathbf{A}') (\mathbf{A} - \mathbf{A}')^\mathrm{T} \right), \quad (1.57) $$

$$ (\forall \mathbf{S} \in \mathcal{S}) \quad F(\mathbf{A}', \mathbf{S}) \le F(\mathbf{A}', \mathbf{S}') + \operatorname{tr} \left( (\mathbf{S} - \mathbf{S}') \nabla_2 F(\mathbf{A}', \mathbf{S}')^\mathrm{T} \right) \\ \qquad + \frac{1}{2} \operatorname{tr} \left( (((\mathbf{S} - \mathbf{S}') \odot (\mathbf{A}'^\mathrm{T} \mathbf{A}' \mathbf{S}') \otimes \mathbf{S}') (\mathbf{S} - \mathbf{S}')^\mathrm{T} \right), \quad (1.58) $$

where $\operatorname{tr}(\cdot)$ denotes the trace operator, and $\odot$ (resp. $\otimes$) the component-wise product (resp. division) between matrices of the same size. Then, the BC-VMFB algorithm, using a cyclic rule, reads:

$$ \mathbf{A}_0 \in \mathcal{A}, \mathbf{S}_0 \in \mathcal{S} $$

For $k = 0, 1, \dots$

$$ \left[ \begin{array}{l} \tilde{\mathbf{A}}_k = \mathbf{A}_k - (\mathbf{A}_k \oslash (\mathbf{A}_k \mathbf{S}_k \mathbf{S}_k^\mathrm{T})) \odot \nabla_1 F(\mathbf{A}_k, \mathbf{S}_k), \\ \mathbf{A}_{k+1} = \mathrm{P}_{\mathcal{A}}(\tilde{\mathbf{A}}_k) = \min\left(\max(\tilde{\mathbf{A}}_k, \mathrm{a}_{\min}), \mathrm{a}_{\max}\right), \\ \tilde{\mathbf{S}}_k = \mathbf{S}_k - (\mathbf{S}_k \oslash (\mathbf{A}_{k+1}^\mathrm{T} \mathbf{A}_{k+1} \mathbf{S}_k)) \odot \nabla_2 F(\mathbf{A}_{k+1}, \mathbf{S}_k), \\ \mathbf{S}_{k+1} = \mathrm{P}_{\mathcal{S}}(\tilde{\mathbf{S}}_k) = \min\left(\max(\tilde{\mathbf{S}}_k, \mathrm{s}_{\min}), \mathrm{s}_{\max}\right). \end{array} \right. \quad (1.59) $$

Let us formulate a few remarks about the above algorithm:

* The projections onto sets $\mathcal{A}$ and $\mathcal{S}$ are here computed in the standard Euclidean metric. This is a direct consequence of the separability of the majorant functions, and of the simple form of the constraint sets. More complex constraints on the sought solution (e.g., sum-to-one, sparsity) can be included at the price of modifying the underlying metric in the projection steps (see for example [72]).

* Algorithm (1.59) is reminiscent from the alternating minimization algorithm proposed in [2] in the context of NMF with Euclidean distance, the main difference being the strict positivity constraint on the components of the involved matrices, introduced in the projection step. This constraint allows us to guarantee the convergence of the iterates of the algorithm to a local minimizer of $F$, according to the convergence theorem from [71].

* A cyclic rule has been used in Algorithm (1.59) for updating matrices $\mathbf{A}$ and $\mathbf{S}$. In fact, the convergence of the algorithm still holds when using a quasi-cyclic rule i.e., when each block is updated at least once every $K \ge 0$ iteration. A particularly interesting strategy is to adopt an inner looping approach, where $\mathbf{A}$ (resp. $\mathbf{S}$) is updated several times before updating $\mathbf{S}$ (resp. $\mathbf{A}$). Such a technique has been applied for instance in [73, 74] in the context of blind deconvolution.
## 1.3 Primal-dual methods

In this section, we present algorithms making use of primal-dual formulations for dealing with the minimization of a composite cost function of the form

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = f_1(\mathbf{u}) + f_2(C\mathbf{u}) \quad (1.60) $$

where $f_1: \mathbb{R}^N \to \mathbb{R} \cup \{+\infty\}$, $f_2: \mathbb{R}^M \to \mathbb{R} \cup \{+\infty\}$, and $C \in \mathbb{R}^{M \times N}$. Subsequently, functions $f_1$ and $f_2$ will be assumed to be convex so as to facilitate our presentation. The presented algorithms can however be applied to the non convex case, even if convergence guarantees become hazardous.

### 1.3.1 Lagrange duality

By introducing an auxiliary variable $v$, the minimization of $f$ can be reformulated as the minimization of a separable sum of functions subject to the linear constraint $v = Cu$, that is

$$ \underset{\mathbf{u} \in \mathbb{R}^N, \mathbf{v} \in \mathbb{R}^M}{\text{minimize}} \quad f_1(\mathbf{u}) + f_2(\mathbf{v}) + \iota_{\{0\}}(Cu - \mathbf{v}), \quad (1.61) $$

with $\iota_{\{0\}}$ is the indicator function of the singleton $\{0\}$. A well-known strategy in optimization to address constrained problems is to resort to the Lagrange multiplier method. This classical technique makes use of the Lagrange function defined as

$$ (\forall \mathbf{u} \in \mathbb{R}^N)(\forall (\mathbf{v}, \mathbf{w}) \in (\mathbb{R}^M)^2) \quad \mathcal{L}(\mathbf{u}, \mathbf{v}, \mathbf{w}) = f_1(\mathbf{u}) + f_2(\mathbf{v}) + \mathbf{w}^\mathrm{T}(Cu - \mathbf{v}), \quad (1.62) $$

where the parameter $w$ is the so-called Lagrange multiplier vector. More precisely, it can be shown that under some mild qualification conditions (for example, the intersection of the interior of the domain of $f_2$ and the image of the domain of $f_1$ by $C$ is nonempty), if $(\hat{\mathbf{u}}, \hat{\mathbf{v}}, \hat{\mathbf{w}})$ is a saddle point of $\mathcal{L}$, then $\hat{\mathbf{u}}$ is a minimizer of $f$. Recall that such a saddle point is defined as follows:

$$ (\forall \mathbf{u} \in \mathbb{R}^N)(\forall (\mathbf{v}, \mathbf{w}) \in (\mathbb{R}^M)^2) \quad \mathcal{L}(\hat{\mathbf{u}}, \hat{\mathbf{v}}, \mathbf{w}) \le \mathcal{L}(\hat{\mathbf{u}}, \hat{\mathbf{v}}, \hat{\mathbf{w}}) \le \mathcal{L}(\mathbf{u}, \mathbf{v}, \hat{\mathbf{w}}). \quad (1.63) $$

It is characterized by the Karush-Kuhn-Tucker (KKT) optimality conditions:

$$ \left\{ \begin{array}{l} -C^\mathrm{T}\hat{\mathbf{w}} \in \partial f_1(\hat{\mathbf{u}}) \\ \hat{\mathbf{w}} \in \partial f_2(\hat{\mathbf{v}}) \\ C\hat{\mathbf{u}} = \hat{\mathbf{v}}. \end{array} \right. \quad (1.64) $$
1.3.2

Alternating direction method of multipliers

1.3.2.1 Basic form

The saddle point property of the Lagrangian suggests to find a minimizer of the cost function $f$ by alternating between a minimization step with respect to the primal variables $u$ and $v$, and a maximization step with respect to the dual variable $w$. However, in order to obtain good convergence properties of the resulting algorithm, the Lagrange function must be modified as follows:

$$ (\forall u \in \mathbb{R}^N)(\forall (v, z) \in (\mathbb{R}^M)^2) \quad \tilde{\mathcal{L}}(u, v, z) = \mathcal{L}(u, v, \gamma z) + \frac{\gamma}{2} \|Cu - v\|^2, \quad (1.65) $$

where, for convenience, we have performed a variable change for the Lagrange multiplier by setting $w = \gamma z$ with $\gamma \in ]0, +\infty[$. Due to the additional quadratic term, $\tilde{\mathcal{L}}$ is called an augmented Lagrangian. It can be noticed that $\tilde{\mathcal{L}}$ is a majorant function of $\mathcal{L}$, which coincides with it when the constraint $v = Cu$ is satisfied. The proposed optimization algorithm then generates a sequence $(u_k, v_k, z_k)_{k \ge 1}$ as follows:

$v_0 \in \mathbb{R}^M, z_0 \in \mathbb{R}^M$

For $k = 0, 1, \dots$

$$ \left. \begin{array}{l} u_{k+1} = \underset{u \in \mathbb{R}^N}{\operatorname{argmin}} \tilde{\mathcal{L}}(u, v_k, z_k) \\ v_{k+1} = \underset{v \in \mathbb{R}^N}{\operatorname{argmin}} \tilde{\mathcal{L}}(u_{k+1}, v, z_k) \\ z_{k+1} \text{ such that } \tilde{\mathcal{L}}(u_{k+1}, v_{k+1}, z_{k+1}) \ge \tilde{\mathcal{L}}(u_{k+1}, v_{k+1}, z_k). \end{array} \right\} (1.66) $$

If the last step is performed through a gradient ascent step with stepsize $1/\gamma$, we obtain the Alternating Direction Method of Multipliers (ADMM) [75, 76, 77, 78], which reads

$v_0 \in \mathbb{R}^M, z_0 \in \mathbb{R}^M$

For $k = 0, 1, \dots$

$$ \left. \begin{array}{l} u_{k+1} = \underset{u \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \|Cu - v_k + z_k\|^2 + \frac{1}{\gamma} f_1(u) \\ s_k = Cu_{k+1} \\ v_{k+1} = \underset{v \in \mathbb{R}^N}{\operatorname{prox}} \frac{f_2}{\gamma} (z_k + s_k) \\ z_{k+1} = z_k + s_k - v_{k+1}. \end{array} \right\} (1.67) $$

The convergence of the sequence $(u_k)_{k \in \mathbb{N}}$ to a minimizer of the cost function $f$ is then secured provided that $C$ has full column rank, that is $C^T C$ is invertible. Note that, by duality arguments, ADMM is strongly related to another famous algorithm in convex optimization, the Douglas-Rachford algorithm [79, 80].
### 1.3.2.2 Minimizing a sum of more than two functions

In practice, one may be interested in more involved cost functions of the form:

$$ (\forall \mathbf{u} \in \mathbb{R}^N) \quad f(\mathbf{u}) = \sum_{j=1}^{J} f_j(C_j \mathbf{u}), \qquad (1.68) $$

where, for every $j = 1, \dots, J$, $f_j$ is a convex function from $\mathbb{R}^{M_j}$ to $\mathbb{R} \cup \{+\infty\}$, and $C_j \in \mathbb{R}^{M_j \times N}$. One might think of a direct extension of ADMM to this context by introducing auxiliary variables $v_j = C_j \mathbf{u}$, but one has to be very cautious in doing so, since convergence guarantees may be lost [81]. A better approach may consist in resorting to consensus-like techniques leading to parallel forms of ADMM [82]. Such an efficient algorithm is the Parallel ProXimal Algorithm (PPXA) which was designed in [83]. An extended version of this algorithm (PPXA+) [84] is described next:

$$ \begin{array}{l} (\mathbf{y}_{0,j})_{1 \le j \le J} \in \mathbb{R}^N, \mathbf{u}_0 = (\sum_{j=1}^J C_j^\top C_j)^{-1} \sum_{j=1}^J C_j^\top \mathbf{y}_{0,j} \\ \text{For } k = 0, 1, \dots \\ \quad \mathbf{v}_{k,j} = \text{prox}_{\gamma f_j}(\mathbf{y}_{k,j}), \quad j = 1, \dots, J \qquad (1.69) \\ \quad \mathbf{c}_k = (\sum_{j=1}^J C_j^\top C_j)^{-1} \sum_{j=1}^J C_j^\top \mathbf{v}_{k,j} \\ \quad \mathbf{y}_{k+1,j} = \mathbf{y}_{k,j} + \lambda(C_j(2\mathbf{c}_k - \mathbf{u}_k) - \mathbf{v}_{k,j}), \quad j = 1, \dots, J \\ \quad \mathbf{u}_{k+1} = \mathbf{u}_k + \lambda(\mathbf{c}_k - \mathbf{u}_k), \end{array} $$

where $\gamma \in [0, +\infty[$ and $\lambda \in [0, 2[$ are two parameters of the algorithm. The use of this algorithm requires that $\sum_{j=1}^J C_j^\top C_j$ be invertible. One of the key additional advantages of PPXA+ is that many operations can be performed in parallel on $J$ processors.

**EXAMPLE 1.19** An illustration of the great performance of Algorithm (1.69) can be found in [85], in the context of signal restoration from DOSY NMR measurements. The applicative context for this modality is explained in detail in Section 1.4, but we introduce it briefly here for the sake of clarity. We are giving measurements $\mathbf{y} \in \mathbb{R}^M$, related to the sought DOSY spectra $\bar{\mathbf{u}} \in \mathbb{R}^N$ through the linear relation $\mathbf{y} = \mathbf{K}\bar{\mathbf{u}} + \mathbf{w}$, with $\mathbf{K} \in \mathbb{R}^{M \times N}$ the observation matrix associated to the DOSY physical model, and $\mathbf{w} \in \mathbb{R}^M$ a noise vector. The authors propose to solve the inverse problem of retrieving an estimate of $\bar{\mathbf{u}}$ from $\mathbf{y}$ and $\mathbf{K}$, by minimizing (1.68) where $J = 2$, $C_1 = \mathbf{K}$, $f_1 = \iota_C$ with

$$ C = \{ \mathbf{z} \in \mathbb{R}^M \mid \|\mathbf{z} - \mathbf{y}\| \le \xi \}, \qquad (1.70) $$

$\xi > 0$, $C_2 = \mathbf{I}_N$ and $f_2$ is an hybrid penalty combining an entropy term and an $\ell_1$ regularization [86]. The application of Algorithm (1.69) leads to the so-called PALMA algorithm, standing for “Proximal Algorithm for L1 combined with MAxent prior”. An example of result obtained by PALMA on the DOSY NMR analysis of crude ethanolic plant extract obtained from brown algae is provided in Figure 1.9.
<figure><img src="image_10.png" /><figcaption>**Figure 1.9** Result of PALMA method on an NMR DOSY experiment performed on a brown algae methanol/water extract.</figcaption></figure>

### 1.3.3

#### Primal-dual proximal algorithms

It can be observed that the computation of $u_{k+1}$ at the $k$-th iteration of ADMM is generally non explicit. Even when $J = 2$, PPXA+ does not present this short-coming, but it requires a matrix inversion for computing each variable $c_k$. When the involved matrices do not have a simple structure and they are of large-size, such inversion has a high computational cost. A number of primal-dual proximal algorithms have been proposed to circumvent this difficulty [87, 88, 89], which are applicable to possibly non smooth optimization problems. A simple way for deriving one of the most popular primal-dual proximal algorithms consists in starting from (1.67) and to replace the update of $u_{k+1}$ by the following semi-implicit subgradient step:

$$ u_{k+1} = u_k - \tau \gamma (C^T(Cu_k - v_k + z_k) + \frac{1}{\gamma}t_k) \quad (1.71) $$

where $t_k$ is a subgradient of $f_1$ at $u_{k+1}$ and $\tau$ is a positive stepsize. This is also equivalent to

$$ y_k = \gamma(Cu_k - v_k + z_k) \quad (1.72) $$

$$ u_{k+1} = \text{prox}_{\tau f_1} (u_k - \tau C^T y_k). \quad (1.73) $$

Now reverting the order of the updates of $v_k$ and $z_k$ with respect to the original form of ADMM and performing some algebra lead to

$$ u_0 \in \mathbb{R}^N, y_0 \in \mathbb{R}^M $$

For $k = 0, 1, \dots$

$$ \left\{ \begin{array}{l} u_{k+1} = \text{prox}_{\tau f_1} (u_k - \tau C^T y_k) \\ d_k = y_k + \gamma C(2u_{k+1} - u_k) \\ y_{k+1} = d_k - \gamma \text{prox}_{\frac{f_2}{\gamma}} \left( \frac{d_k}{\gamma} \right). \end{array} \right. \quad (1.74) $$

The convergence of $(u_k)_{k \in \mathbb{N}}$ to a minimizer of $f$ can then be shown if $\tau \gamma |||C|||^2 \le 1$, where we recall that $||| \cdot |||$ denotes the spectral norm. Let us emphasize that no
matrix inversion is required in the previous algorithm. Various extensions of this framework are possible. In particular, other forms of primal-dual proximal algorithms can be obtained [89]. It is also possible to add in the original criterion a Lipschitz-continuous term, which is addressed in the algorithm through its gradient [90, 91, 92, 93]. Variants of this algorithm have been proposed so as to minimize composite function of the form (1.68) [89]. Block-alternating implementations of primal-dual proximal algorithms, along with their convergence properties, can be found for instance in [94, 95].

**EXAMPLE 1.20** The primal-dual proximal algorithm is applied in [96] to mass spectrometry. Following a framework similar to the one of Example 1.16, the authors consider now the $l_1$ norm penalty function $f_1(\mathbf{u}) = \sum_{n=1}^{N} |u_n|$. The resulting optimization problem is then convex non-differentiable, and thus Algorithm 1.74 is particularly well suited. Example of results obtained for the analysis of a real MS dataset measured on a Brucker Solarix 15 T, FT-ICR instrument with an ESI source, is provided in Figure 1.10. The considered sample was constituted of 3 µM of the peptide EVEALEKKVAALESKVKQALEKKVEALEHG-NH2 (C140H240N38O45) in its trimer form within 50 mM of NH4OAc, acquired in native conditions. Despite the very large size of the problem, the processing time was of about 108 min on a standard laptop.

<figure><img src="image_11.png" /><figcaption>**Figure 1.10** Analysis of the real FT-ICR-MS spectrum of a peptide in trimer form: (top) zoom on the acquired data; (bottom) recovered spectrum for the charge $z = 5$, using Algorithm 1.74 for a mass grid size of $M = 8130981$ and total charge number of $Z = 5$.</figcaption></figure>

### 1.3.4 Primal-dual interior point algorithm

Let us consider the optimization problem

$$ \min_{\mathbf{u} \in C} f_1(\mathbf{u}), \qquad (1.75) $$
<figure><img src="image_12.png" /><figcaption>**Figure 1.11** Logarithmic barrier function $b(u) = -\log(u) - \log(1 - u)$ associated to the constraint $u \in [0, 1]$.</figcaption></figure>

subject to the linear equality constraints

$$C = \{ \mathbf{u} \in \mathbb{R}^N \mid C\mathbf{u} + \rho \in ]0, +\infty[ \} \quad (1.76)$$

$f_1$ is a convex twice differentiable function, $C \in \mathbb{R}^{M \times N}$, and $\rho = (\rho_m)_{1 \le m \le M} \in \mathbb{R}^M$. The existence of a minimizer $\hat{\mathbf{u}}$ will be subsequently assumed. Interior point methods, introduced in [97], solve the constrained optimization problem (1.75) by introducing a sequence of unconstrained optimization subproblems:

$$(\forall k \in \mathbb{N}) \quad \underset{\mathbf{u} \in \mathbb{R}^N}{\text{minimize}} \quad f_{\mu_k}(\mathbf{u}) = f_1(\mathbf{u}) + \mu_k b(\mathbf{u}) \quad (1.77)$$

for positive barrier parameter values $(\mu_k)_{k \in \mathbb{N}}$ decaying to 0. The auxiliary function $b$, called *barrier function*, penalizes the closeness to the constraint boundaries and maintains the iterates inside the strict interior of the constrained domain. The most widely used auxiliary function is the logarithmic barrier

$$(\forall \mathbf{u} \in \mathbb{R}^N) \quad b(\mathbf{u}) = f_2(C\mathbf{u}), \quad (1.78)$$

where

$$(\forall \mathbf{v} = (v_m)_{1 \le m \le M} \in \mathbb{R}^M) \qquad f_2(\mathbf{v}) = \begin{cases} -\sum_{m=1}^M \log(v_m + \rho_m) & \text{if } \mathbf{v} \in ]0, +\infty[ \\ +\infty & \text{otherwise.} \end{cases} \quad (1.79)$$

For every $\mu \in ]0, +\infty[$, this choice makes the penalized criterion $f_\mu$ unbounded at the boundary of the feasible region so that its minimizers $\hat{\mu}_\mu$ fulfill strictly the constraints. An example is presented in Figure ??.

Although various classical interior-point methods can be envisaged to solve Problem (1.77), a primal-dual approach can be followed by reformulating the problem under the form (1.61). The resulting primal-dual interior-point methods simultaneously estimate the primal variables $\mathbf{u}$ and the dual Lagrange multiplier vector
$\lambda = -w$ associated with the constraints [98]. The joint estimation of primal and dual variables is performed through KKT optimality conditions (1.64):

$$ egin{cases} 
\nabla f_1(\hat{u}_\mu) - C^T \hat{\lambda}_\mu = \mathbf{0} \\
\hat{\Lambda}_\mu(C\hat{u}_\mu + \rho) = \mu \\
C\hat{u}_\mu + \rho > \mathbf{0}. 
\end{cases} 
\quad (1.80) $$

where $\hat{\Lambda}_\mu = \text{Diag}\{\hat{\lambda}_\mu\}$ and $\mu = \mu[1, \dots, 1]^T \in ]0, +\infty[^{M}$.

Numerically, a sequence of optimization problems is solved for a sequence of penalization (also called perturbation) parameters $(\mu_k)_{k \in \mathbb{N}}$ converging to 0. At each iteration $k \in \mathbb{N}$ of the algorithm, a pair of primal-dual variables $(u_{k+1}, \lambda_{k+1})$ is firstly computed from an approximate solution to (1.80) through a Newton algorithm step on the equality conditions, in association with a linesearch strategy on a merit function incorporating some barrier terms, allowing to ensure the inequality condition [9, Chap.11]. The update strategy is then given by

$$ (u_{k+1}, \lambda_{k+1}) = (u_k - \alpha_k c_k, \lambda_k - \alpha_k d_k), \quad (1.81) $$

where $\alpha_k$ is the stepsize and $(c_k, d_k)$ are the primal and dual Newton directions. Finally, the perturbation parameter $\mu_{k+1}$ is updated in order to ensure the algorithm convergence. Based on the iterative scheme (1.81), several primal-dual interior point methods have been proposed in the literature, each of them calling for its own strategy for the computation of the primal-dual directions, the derivation of a suitable stepsize, and the update of the perturbation parameter (see [98, 99] for a review). We focus here on the iterative scheme introduced in [100, 101]. Note that the algorithm described above is reminiscent from the one from [102], with the introduction of novel tools based on the MM principle aiming at accelerating the practical convergence, as well as reducing the computational cost per iteration.

#### 1.3.4.1 Primal-dual directions

The Newton directions $(c_k, d_k)$ are obtained by solving the linear system

$$ \begin{bmatrix} \nabla^2 f_1(u_k) & -C^T \\ \Lambda_k C & \text{Diag}\{C u_k + \rho\} \end{bmatrix} \begin{bmatrix} c_k \\ d_k \end{bmatrix} = r(u_k, \lambda_k, \mu_k), \quad (1.82) $$

where $\Lambda_k = \text{Diag}\{\lambda_k\}$ and $r(u, \lambda, \mu)$ is the primal-dual residual defined as

$$ r(u, \lambda, \mu) = \begin{pmatrix} r^{\text{prim}}(u, \lambda) \\ r^{\text{dual}}(u, \lambda, \mu) \end{pmatrix} = \begin{pmatrix} \nabla f_1(u) - C^T \lambda \\ \Lambda(Cu + \rho) - \mu \end{pmatrix}. \quad (1.83) $$

As pointed out in [103, 104], the primal-dual matrix in the left side of (1.82) suffers from ill-conditioning as soon as some $[Cu_k + \rho]_m$ or some $\lambda_m$ gets closer to zero. Moreover, this matrix is not guaranteed to be symmetric nor definite positive [99, 105], so that the linear system (1.82) is difficult to solve. Therefore, rather than solving directly (1.82), variable substitution is used [102, 106]. From the second
equation of (1.82), one calculates the dual direction according to

$$ 
\begin{aligned}
d_k &= \operatorname{Diag}\{C u_k + \rho\}^{-1} (\Lambda_k (C u_k + \rho) - \mu_k - \Lambda_k C c_k) \\
    &= \lambda_k - \operatorname{Diag}\{C u_k + \rho\}^{-1} (\mu_k + \Lambda_k C c_k).
\end{aligned} 
\quad (1.84) 
$$

Then, the primal direction is obtained by solving the linear system

$$ H_k c_k = g_k \quad (1.85) $$

with

$$ H_k = \nabla^2 f_1(u_k) + C^\top \operatorname{Diag}\{C u_k + \rho\}^{-1} \Lambda_k C \quad (1.86) $$

$$ g_k = \nabla f_1(u_k) - C^\top \operatorname{Diag}\{C u_k + \rho\}^{-1} \mu_k. \quad (1.87) $$

The computation cost of the primal-dual interior point algorithm is almost exclusively due to the resolution of the primal system (1.85). It is shown in [107] that an approximate solution of this system is sufficient to ensure the convergence of the method. Several solutions exist to calculate such approximation (see [108] and references therein). For instance, [100] proposes to perform an approximate solution using a preconditioned bi-conjugate gradient algorithm based on an incomplete LU factorization of matrix $H_k$. Moreover, [101] emphasizes that, when $H_k$ has a block-diagonal matrix, the resolution of (1.85) reduces to the resolution of a family of linear systems of small-size and such structure is well suited for parallel computing as the blocks can be processed independently. When $H_k$ is not easily invertible, a separable quadratic majorization of it can be employed, and an MM algorithm can be applied to solve (1.85). More precisely, let us remark that solving (1.85) is equivalent to the resolution of a quadratic minimization problem of the form

$$ \underset{c \in \mathbb{R}^M}{\operatorname{minimize}} \quad \frac{1}{2} c^\top H_k c - g_k^\top c. \quad (1.88) $$

Let $B_k$ be a symmetric positive definite matrix such that $H_k \succeq B_k$, and whose inverse is simple to compute (for instance, when $B_k$ is a block diagonal matrix). Then, the solution of (1.88) is computed thanks to the following iterations:

$$ 
\begin{aligned}
& c_{k,0} = 0, \\
& \text{For } j = 0, 1, \dots \\
& \quad \left[ c_{k,j+1} = c_{k,j} - B_k^{-1} (H_k c_{k,j} - g_k) \right],
\end{aligned} 
\quad (1.89) 
$$

until the fulfillment of the following stopping criterion on the primal residual: $||H_k d_{k,j} - g_k|| \le \mu_k ||g_k||$ where $\mu_k$ is the barrier parameter. The benefit of this method has been illustrated in [109] through an example of spectral unmixing. Note that this MM-like approach is not only faster, it is also better suited for parallel implementation [110].
#### 1.3.4.2 Linesearch

At the k-th iteration, the stepsize value $\alpha_k$ must be chosen so as to ensure the convergence of the algorithm and the fulfillment of the inequalities of the pertubed KKT system (1.80). The convergence study of the primal-dual algorithm presented in [102] requests that $\alpha_k$ ensures a sufficient decrease of the primal-dual merit function $\Psi_{\mu_k}$, defined, for every $u \in \mathbb{R}^N$ and $\lambda \in \mathbb{R}^M$, as

$$ \Psi_{\mu_k}(u, \lambda) = f_{\mu_k}(u) + \lambda^T(Cu + \rho) - \mu_k \sum_{m=1}^{M} \log(\lambda_m[Cu + \rho]_m), \quad (1.90) $$

$$ = \mathcal{L}(u, Cu, -\lambda) - \mu_k \sum_{m=1}^{M} \log \lambda_m. \quad (1.91) $$

The sufficient decrease is assessed using the Armijo condition,

$$ \psi_{\mu_k}(\alpha_k) - \psi_{\mu_k}(0) \le \sigma \alpha_k \nabla \psi_{\mu_k}(0) \quad \text{with} \quad \sigma \in [0, 1/2[, \quad (1.92) $$

where, for every $\alpha \in [0, +\infty[$, $\psi_{\mu_k}(\alpha) \triangleq \Psi_{\mu_k}(u_k - \alpha c_k, \lambda_k - \alpha d_k)$. One can note that (1.91) contains two logarithmic barrier functions enforcing the fulfillment of the KKT inequalities, the positivity of the Lagrange multipliers $(\lambda_m)_{1 \le m \le M}$ being a straightforward consequence of (1.80). The presence of the barrier function may cause the inefficiency of general purpose backtracking line search methods for finding a stepsize satisfying (1.92) and, thus, the slowdown of the algorithm convergence. Several strategies have been proposed to override this issue [111, 112, 113]. In particular, the majorization-minimization strategy from [112], relying on the construction of a majorizing function made of a quadratic term and a logarithmic barrier term, was shown to lead to good practical performance through numerical applications to 2D nuclear magnetic resonance signal reconstruction under positivity constraints [114], and to sparse signal deconvolution [113].

#### 1.3.4.3 Perturbation parameter update

The convergence analysis of interior point methods requires that the sequence $(\mu_k)_{k \in \mathbb{N}}$ tends to 0 as $k$ tends to infinity. An efficient update strategy for the choice of the barrier parameter is the $\mu$-criticity rule introduced in [115]:

$$ (\forall k \in \mathbb{N}) \quad \mu_k = \theta \frac{\delta_k}{M}, \quad (1.93) $$

where $\delta_k = (Cu + \rho)^T \lambda_k$ is the duality gap and $\theta \in [0, 1[$. The barrier parameter is updated as soon as the primal and dual directions fulfill the inner stopping rule ([106, 116]):

$$ \|r^{\text{prim}}(u_k, \lambda_k)\|_{\infty} \le \eta^{\text{prim}} \mu_k \text{ and } \frac{\delta_k}{M} \le \eta^{\text{dual}} \mu_k, \quad (1.94) $$

with $\eta^{\text{prim}} > 0$ and $\eta^{\text{dual}} \in [1, 1/\theta[$.
#### 1.3.4.4 Resulting algorithm

The main steps of the proposed optimization method are summarized below:

$$ 
\begin{array}{l}
u_0 \in \mathbb{R}^N \quad \text{s.t.} \quad Cu_0 + \rho \in [0, +\infty)^M, \lambda_0 \in [0, +\infty)^M, \mu_0 \in [0, +\infty[. \\
\text{For } k = 0, 1, \dots \\
\quad \text{If (1.94) holds} \\
\quad \quad \lfloor \text{Update } \mu_k \text{ according to (1.93)} \\
\quad \quad \text{Calculate } c_k \text{ by solving (approximately) (1.85)} \\
\quad \quad \text{Deduce } d_k \text{ from (1.84)} \\
\quad \quad \text{Find } \alpha_k > 0 \text{ satisfying (1.92)} \\
\quad \quad \text{Update } (u_{k+1}, \lambda_{k+1}) \text{ according to (1.81).} \\
\end{array}
\tag{1.95}
$$

The convergence properties of Algorithm (1.95) have been studied in [102, 107, 101], under the assumption that the linear constraint set is nonempty and bounded. In particular, it is shown that if $f_1$ is strictly convex, then the sequence $(u_k)_{k \in \mathbb{N}}$ converges to the unique solution of (1.75).

<figure><img src="image_13.png" /><figcaption>**Figure 1.12** Estimated $T_1-T_2$ NMR distributions from real data (apple) using Algorithm (1.95).</figcaption></figure>

**EXAMPLE 1.21** The primal-dual interior point algorithm (1.95) has been applied to two dimensional nuclear magnetic resonance in [114]. Classical NMR experiments analyze the spin relaxation process independently, either in terms of longitudinal or transverse relaxation, leading to the so-called $T_1$ and $T_2$ relaxation spectra. Joint measurements with respect to both relaxation parameters allow to build two-dimensional distribution which is of high interest for chemical structure determination [117]. Experimental data consist of a series of discrete noisy samples $Y \in \mathbb{R}^{m_1 \times m_2}$ such that $Y = K_1 U K_2^\top + W$ with $U \in \mathbb{R}^{N_1 \times N_2}$ the sought distribution, $K_1, K_2$ matrices associated to the acquisition model, and $W$ a noise term. $T_1-T_2$ NMR reconstruction aims at estimating $U$ given $Y, K_1, K_2$. This problem can be solved by minimizing a penalized least-squares term:

$$ f_1(U) = \frac{1}{2} \|K_1 U K_2^\top - Y\|_F^2 + \lambda \|U\|_F^2 \tag{1.96} $$
with $\lambda > 0$, subject to the positivity constraint $U \ge 0$. Figure 1.12 presents numerical results of 2D spectra obtained from the analysis of an apple matter sample ($m_1 = 50, m_2 = 10000, N_1 = 300, N_2 = 300$), using Algorithm (1.95) [114].

## 1.4 Illustration in the context of DOSY NMR signal restoration

This section illustrates the applicability of the presented optimization tools for signal restoration in nuclear magnetic resonance (NMR). The measurement of diffusion by NMR is used in various application fields (agroalimentary sector, pharmaceutical industry, ecology) to analyze the properties of complex chemical mixtures in order to determine their molecular structure and dynamics [85, 117, 118]. After the immersion of the matter in a strong magnetic field, all the nuclear spins align to an equilibrium state along the field orientation. The application of a short magnetic pulse, i.e. the pulsed field gradient, in resonance with the spin motion disturbs the spin orientation. NMR aims at analyzing the process which corresponds to the re-establishment of the spin into its equilibrium state. During the DOSY (Diffusion Order Spectroscopy) experiment [119], a series of measurements is acquired for different pulsed field gradient strengths. The data are then analyzed with the aim to separate different species according to their diffusion coefficient.

The DOSY NMR data $\mathbf{y} = (y^{(m)})_{1\le m\le M} \in \mathbb{R}^M$ gathers the results of $M$ experiments characterized by a vector $\mathbf{t} = (t^{(m)})_{1\le m\le M} \in \mathbb{R}^M$ related to the pulsed field gradient strength and to the acquisition time. The relation between $\mathbf{y}$ and $\mathbf{t}$ can be expressed as the following Laplace transform:

$$ (\forall m \in \{1, \dots, M\}) \quad y^{(m)} = \int \chi(T) \exp(-t^{(m)}T) dT, \qquad (1.97) $$

where $\chi(T)$ is the unknown diffusion distribution. The problem is then to reconstruct $\chi(T)$ on the sampled grid $\mathbf{T} = (T^{(n)})_{1\le n\le N}$, from the measurements $\mathbf{y}$. After discretization and appropriate normalization, the observation model reads

$$ \mathbf{y} = K\bar{\mathbf{u}} + \mathbf{w}, \qquad (1.98) $$

where $K \in \mathbb{R}^{M \times N}$ is given by

$$ (\forall m \in \{1, \dots, M\})(\forall n \in \{1, \dots, N\}) \quad K^{(m,n)} = \exp(-T^{(n)}t^{(m)}), \qquad (1.99) $$

$\bar{\mathbf{u}} \in \mathbb{R}^N$ is the sought signal related to $(\chi(T^{(n)}))_{1\le n\le N}$ (up to a scaling factor depending on the discretization grid), and $\mathbf{w} \in \mathbb{R}^M$ is a perturbation noise.
<figure><img src="image_14.png" /><figcaption>Figure 1.13 Original signal $\bar{u}$ (left) and noisy acquired measure $y$ (right).</figcaption></figure>

We propose here to find an estimate $\hat{u} \in \mathbb{R}^N$ of $\bar{u}$ by solving the following minimization problem, requiring the knowledge of $K$ and $y$:

$$ \underset{u \in \mathbb{R}^N}{\text{minimize}} \quad \frac{1}{2} \|Ku - y\|^2 + \beta g(u) \qquad (1.100) $$

where $g \in \Gamma_0(\mathbb{R}^N)$ denotes a regularization term and $\beta > 0$ is a regularization parameter. Note that, in practice, a very large number of DOSY NMR acquisitions (typically, $10^4$) are conducted for various settings of the pulsed gradient field, so that Problem (1.100) must be solved many times, which motivates the search for a fast minimization algorithm [85]. In the following, we present optimization solutions, for various choices of function $g$, in the line of the work published in [86]. For illustration purpose, we will consider a test example with $\bar{u}$ of size $N = 256$, and $y$ of size $M = 50$, both represented in Figure 1.13. The noise $w$ is a realization of a zero-mean white Gaussian noise, with standard deviation equals to $10^{-2}y_0 \approx 0.34$. The grid $T = (T^{(n)})_{1 \le n \le N}$ is sampled following a logarithmic rule, within the interval $[T_{\text{min}}, T_{\text{max}}] = [1, 10^3]$ $\mu m^2 s^{-1}$, while $t = (t^{(m)})_{1 \le m \le M}$ is regularly sampled on $[t_{\text{min}}, t_{\text{max}}] = [0, 1.5]$ seconds. The quality of the results obtained by the different tested restoration approaches will be assessed quantitatively by means of the normalized mean square error (NMSE):

$$ \text{NMSE} = \frac{\|\bar{u} - \hat{u}\|^2}{\|\bar{u}\|^2}, \qquad (1.101) $$

with $\hat{u}$ the result of each method.

### 1.4.1 Quadratic penalization

Let us start with the following regularization term serving to promote the reconstruction of smooth signals:

$$ (\forall u \in \mathbb{R}^N) \quad g(u) = \frac{1}{2} \|Du\|^2 \qquad (1.102) $$

where $D \in \mathbb{R}^{N \times N}$ is the discrete gradient operator such that,

$$ (\forall n \in \{1, \dots, N\}) \quad [Du]^{(n)} = u^{(n)} - u^{(n-1)} \qquad (1.103) $$
<figure><img src="image_15.png" /><figcaption>Figure 1.14 Restored signal (1.104) for $\beta = 1$, NMSE = 0.29.</figcaption></figure>

with the circular convention $u^{(0)} = u^{(N)}$. In this case, the cost function involved in Problem (1.100) is convex and quadratic. Assuming that $K^\top K + \beta D^\top D$ is invertible, the solution to Problem (1.100) is unique and reads:

$$ \hat{u} = (K^\top K + \beta D^\top D)^{-1} K^\top y. \quad (1.104) $$

An example of such restored signal is displayed in Figure 1.14.

Now, assume that we want to impose some value range constraints on the sought signal. Then a possible regularization function is:

$$ (\forall u \in \mathbb{R}^N) \quad g(u) = \frac{1}{2} \|Du\|^2 + \iota_{[u_{\min}, u_{\max}]^N}(u) \quad (1.105) $$

with $0 \le u_{\min} < u_{\max}$ the minimum and maximum values of the original signal $(\bar{u}^{(n)})_{1 \le n \le N}$. The quadratic function $u \mapsto 1/2 \|Ku - y\|^2 + \beta/2 \|Du\|^2$ is Lipschitz differentiable, with constant $\|K^\top K + \beta D^\top D\|$. Problem (1.100) can thus be solved using the projected gradient algorithm introduced in Example 1.15, which reads in this case:

$$ \begin{array}{l} u_0 \in \mathbb{R}^N \\ \text{For } k = 0, 1, \dots \\ \quad \left[ \begin{array}{l} \alpha_k \in [0, 2/(\|K^\top K + \beta D^\top D\|)], \\ \tilde{u}_k = u_k - \alpha_k (K^\top (K u_k - y) + \beta D^\top D u_k), \\ u_{k+1} = P_{[u_{\min}, u_{\max}]^N} (\tilde{u}_k) = \min (\max (\tilde{u}_k, u_{\min}), u_{\max}) \end{array} \right. \end{array} \quad (1.106) $$

where the min and max operations are performed componentwise. An example of solution obtained using the above algorithm is displayed in Figure 1.15.

## 1.4.2

### Sparsity prior in the signal domain

Another strategy for regularization is to enforce the sparsity of the sought signal, in addition to the positivity constraint. In that respect, one can use

$$ (\forall u = (u_n)_{1 \le n \le N} \in \mathbb{R}^N) \quad g(u) = \|u\|_1 + \iota_{[0, +\infty]^N}(u). \quad (1.107) $$
<figure><img src="image_16.png" /><figcaption>Figure 1.15 Restored signal using Algorithm (1.106) for $\beta = 0.0625$, NMSE = 0.14.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 1.16 Restored signal using Algorithm (1.108) for $\beta = 1$, NMSE = 0.11.</figcaption></figure>

Function $g$ is convex, but it is not differentiable. Function $u \mapsto 1/2 \|Ku - y\|^2$ is Lipschitz differentiable with constant $||K||^2$. Problem (1.100) can be solved using the forward-backward algorithm presented in Remark 1.3:

$$ \begin{array}{ll} u_0 \in \mathbb{R}^N & \\ \text{For } k = 0, 1, \dots & \\ \quad \alpha_k \in [0, 2/||K||^2[, & (1.108) \\ \quad \tilde{u}_k = u_k - \theta_k K^\top (K u_k - y), & \\ \quad u_{k+1} = \text{prox}_{\alpha_k \beta (\|\cdot\|_1 + \iota_{[0, +\infty]^N})} (\tilde{u}_k) = \max (|\tilde{u}_k| - \alpha_k \beta, 0). & \end{array} $$

The convergence rate of the above method can be improved, using a variable metric approach ([55], Section 1.2.4) or a Nesterov-based scheme (see for instance, [120]). An example of solution obtained using Algorithm (1.108) is displayed in Figure 1.16.

### 1.4.3

#### Sparsity prior in a transformed domain

It can finally be useful to impose the sparsity of the signal in a transformed domain, in order to impose some regularity properties. For instance, it can be assumed that $W\bar{u}$ is sparse, where $W$ is a (possibly overcomplete) wavelet analysis operator [121]. In this case, we may use

$$ (\forall u \in \mathbb{R}^N) \quad g(u) = \|Wu\|_1. \qquad (1.109) $$

Because of the presence of matrix $W$, the proximal step in the forward-backward algorithm is not explicit anymore. In such context, an efficient strategy is to resort
<figure><img src="image_18.png" /><figcaption>Figure 1.17 Restored signal using Algorithm (1.110) for $\beta = 27$, NMSE = 0.39.</figcaption></figure>

to a primal-dual optimization approach as described in Section 1.3. For instance, the primal-dual proximal algorithm (1.74) [89] would read in this case:

$u_0 \in \mathbb{R}^N, v_0 \in \mathbb{R}^M$ \\
$(\tau, \gamma)$ such that $\tau\gamma \|W\|^2 \le 1$

For $k = 0, 1, \dots$

$$
\left\{
\begin{aligned}
u_{k+1} &= \operatorname{prox}_{\tau/2 \|K \cdot y\|^2} (u_k - \tau W^\top v_k) \\
&= (I + \tau K^\top K)^{-1} (u_k - \tau W^\top v_k + \tau K^\top y) \\
d_k &= v_k + \gamma W (2u_{k+1} - u_k) \\
v_{k+1} &= d_k - \gamma \operatorname{prox}_{\beta\gamma^{-1}\|\cdot\|_1} \left(\frac{d_k}{\gamma}\right) \\
&= d_k - \gamma \operatorname{sign}(d_k) \max\left(|\frac{d_k}{\gamma}| - \beta\gamma^{-1}, 0\right).
\end{aligned}
\right.
\quad (1.110)
$$

An example of solution obtained using Algorithm (1.110), where $W$ represents the Symlet orthonormal wavelet analysis operator, with level 3, is displayed in Figure 1.17.

1.4.4

Sparsity prior and range constraints

Let us end this section by focusing on the case when the regularization term in-
cludes both sparsity constraints (possibly in a transformed domain), and value range
constraints on the sought signal. The problem then becomes:

$$
\underset{\mathbf{u} \in \mathbb{R}^N}{\text{minimize}} \frac{1}{2} \|\mathbf{K}\mathbf{u} - \mathbf{y}\|^2 + \beta \|\mathbf{W}\mathbf{u}\|_1 + \iota_{[u_{\min}, u_{\max}]^N}(\mathbf{u}). \quad (1.111)
$$

The cost function reads as the sum of three terms, namely a quadratic term, and two
non-differentiable terms. An efficient solution to Problem (1.111) can be obtained
<figure><img src="image_19.png" /><figcaption>Figure 1.18 Restored signal using Algorithm (1.112) for $\beta = 20$, NMSE = 0.15.</figcaption></figure>

by using the PPXA+ algorithm (1.69) [84], which reads in this case:

$$
\begin{array}{l}
(z_{0,j})_{1 \le j \le 3} \in \mathbb{R}^N \\
\gamma \in ]0, +\infty[, \quad \lambda \in ]0, 2[ \\
\mathbf{u}_0 = (\mathbf{I} + \mathbf{K}^\top \mathbf{K} + \mathbf{W}^\top \mathbf{W})^{-1} (\mathbf{K}^\top \mathbf{z}_{0,1} + \mathbf{W}^\top \mathbf{z}_{0,2} + \mathbf{z}_{0,3}) \\
\text{For } k = 0, 1, \dots \\
\quad \mathbf{v}_{k,1} = \text{prox}_{\frac{\gamma}{2} \parallel\cdot\parallel} (-\mathbf{y})^2 (\mathbf{z}_{k,1}) = \mathbf{y} + (\mathbf{z}_{k,1} - \mathbf{y})/(\gamma + 1) \\
\quad \mathbf{v}_{k,2} = \text{prox}_{\gamma\beta\parallel\cdot\parallel_1} (\mathbf{z}_{k,2}) = \text{sign}(\mathbf{z}_{k,2}) \odot \max(|\mathbf{z}_{k,2}| - \beta\gamma, 0). \\
\quad \mathbf{v}_{k,3} = \text{P}_{[u_{\min}, u_{\max}]^N} (\mathbf{z}_{k,3}) = \min(\max(\mathbf{z}_{k,3}, u_{\min}), u_{\max}) \\
\quad \mathbf{c}_k = (\mathbf{I} + \mathbf{K}^\top \mathbf{K} + \mathbf{W}^\top \mathbf{W})^{-1} (\mathbf{K}^\top \mathbf{v}_{k,1} + \mathbf{W}^\top \mathbf{v}_{k,2} + \mathbf{v}_{k,3}) \\
\quad \mathbf{z}_{k+1,1} = \mathbf{z}_{k,1} + \lambda(\mathbf{K}(2\mathbf{c}_k - \mathbf{u}_k) - \mathbf{v}_{k,1}) \\
\quad \mathbf{z}_{k+1,2} = \mathbf{z}_{k,2} + \lambda(\mathbf{W}(2\mathbf{c}_k - \mathbf{u}_k) - \mathbf{v}_{k,2}) \\
\quad \mathbf{z}_{k+1,3} = \mathbf{z}_{k,3} + \lambda(2\mathbf{c}_k - \mathbf{u}_k - \mathbf{v}_{k,3}) \\
\quad \mathbf{u}_{k+1} = \mathbf{u}_k + \lambda(\mathbf{c}_k - \mathbf{u}_k).
\end{array}
\tag{1.112}
$$

An example of solution obtained using Algorithm (1.112) is displayed in Fig-
ure 1.18.

1.4.5

**Concluding remarks**

We have illustrated the applicability of the introduced optimization methods, in the
context of the restoration of a signal from NMR DOSY measurements. For various
choices of the regularization function, we have provided a possible minimization
scheme. It is worth noting that there is rarely a single technique available for the
resolution of an optimization problem. For a given application, it is always recom-
mended to test and compare different strategies, as we did in the presented example,
in order to reach the best possible tradeoff in terms of computational complexity and
convergence rate.
1.5

Conclusion

Table 1.2 presents a comprehensive list of the algorithms that have been presented in this chapter and the types of problem they can address. Available acceleration strategies to improve their computational efficiency are also mentioned. Since optimization constitutes a very rich and continuously evolving scientific area, it was not possible to provide an exhaustive description of all the existing methods and their various, sometimes empirical, extensions. We think however that the presented approaches provide a good basis for sound applicative developments grounded on strong theoretical foundations.

<table><caption>Table 1.2 Summary of the algorithms presented in this chapter.</caption><thead><tr><td>Algorithm name</td><td>Cost function</td><td>Particular cases</td><td>Acceleration techniques</td></tr></thead><tbody><tr><td>Generic MM algorithm (1.13)</td><td>Any proper lower-continuous function</td><td>Expectation-minimization, IRL1, IRLS, Richardson-Lucy, DC programming.</td><td>Linesearch strategy [122], block-coordinate approach [12], fixed-point schemes [123].</td></tr><tr><td>Quadratic MM algorithm (1.26)</td><td>Lipschitz differentiable function.</td><td>Gradient descent, SART, half-quadratic (in the smooth case), modified Newton.</td><td>Inexact version (1.30), subspace strategy (1.37), block coordinate approach (1.54).</td></tr><tr><td>Variable metric forward-backward algorithm (1.39)</td><td>Sum of a Lipschitz differentiable and of a convex non differentiable function.</td><td>Forward-backward, ISTA, (scaled) projected gradient, proximal point.</td><td>Inexact version [55], block-coordinate approach (1.55).</td></tr><tr><td>Primal-dual proximal algorithms [89]</td><td>Sum of convex functions composed with linear operators</td><td>ADMM (1.67), PPXA+ (1.69), primal-dual algorithm (1.74).</td><td>Preconditioning [124], block-coordinate approach [125], distributed strategy [95].</td></tr><tr><td>Primal-dual interior point algorithm (1.95)</td><td>Twice differentiable function subject to linear constraints</td><td></td><td>Inexact MM resolution (1.89), linesearch strategy [112].</td></tr></tbody></table>
References

1. Rockafellar, R. and Wets, R. (1998) *Variational Analysis, Grundlehren der Mathematischen Wissenschaften*, vol. 317, Springer, Berlin, 1st edn.

2. Lee, D.D. and Seung, H.S. (2001) Algorithms for non-negative matrix factorization, in *Advances in Neural and Information Processing Systems*, vol. 13, vol. 13, pp. 556–562.

3. Févotte, C., Bertin, N., and Durrieu, J.L. (2009) Nonnegative matrix factorization with the itakura-saito divergence. with application to music analysis. *Neural Computation*, 21 (3), 793–830.

4. Févotte, C. and Idier, J. (2011) Algorithms for nonnegative matrix factorization with the beta-divergence. *Neural Computation*, 23 (9), 2421–2456.

5. Pham, D.T. and Garat, P. (1997) Blind separation of mixtures of independent sources through a quasi maximum likelihood approach. *IEEE Transactions on Signal Processing*, 45 (7), 1712–1725.

6. Cichocki, A., Mandic, D., Lathauwer, L.D., Zhou, G., Zhao, Q., Caiafa, C., and PHAN, H.A. (2015) Tensor decompositions for signal processing applications: From two-way to multiway component analysis. *IEEE Signal Processing Magazine*, 32 (2), 145–163.

7. Hunter, D.R. and Lange, K. (2004) A tutorial on MM algorithms. *Amer. Stat.*, 58 (1), 30–37.

8. Lange, K., Hunter, D.R., and Yang, I. (2000) Optimization transfer using surrogate objective functions with discussion. *Journal of Computational and Graphical Statistics*, 9 (1), 1–59.

9. Boyd, S. and Vandenberghe, L. (2004)

*Convex Optimization*, Cambridge University Press, New York, 1st edn.

10. Bohning, D. and Lindsay, B.G. (1988) Monotonicity of quadratic-approximation algorithms. *Annals of the Institute of Statistical Mathematics*, 40 (4), 641–663.

11. Zhang, Z., Kwok, J.T., and Yeung, D.Y. (2007) Surrogate maximization/minimization algorithms and extensions. *Mach. Learn.*, 69, 1–33.

12. Hong, M., Razaviyayn, M., Luo, Z.Q., and Pang, J.S. (2016) A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing. *IEEE Signal Processing Magazine*, 33 (1), 57–77.

13. Horst, R. and Thoai, N.V. (1999) De programming: Overview. *Journal of Optimization Theory and Applications*, 103 (1), 1–43.

14. Nocedal, J. and Wright, S.J. (1999) *Numerical Optimization*, Springer-Verlag, New York, NY.

15. Combettes, P.L. and Pesquet, J.C. (2010) Proximal splitting methods in signal processing, in *Fixed-Point Algorithms for Inverse Problems in Science and Engineering* (eds H.H. Bauschke, R. Burachik, P.L. Combettes, V. Elser, D.R. Luke, and H. Wolkowicz), Springer-Verlag, New York, pp. 185–212.

16. Shor, N.Z. (1985) *Minimization Methods for Non-Differentiable Functions*, Springer Series in Computational Mathematics, vol. 3, Springer Berlin Heidelberg.

17. Ahmad, R. and Schniter, P. (2015)
Iteratively reweighted $ℓ_1$ approaches to sparse composite regularization, *Tech. Rep.*. http://arxiv.org/pdf/1504.05110.pdf.

18 Wipf, D. and Nagarajan, S. (2010) Iterative reweighted $ℓ_1$ and $ℓ_2$ methods for finding sparse solutions. *IEEE Journal on Selected Topics in Signal Processing*, **4** (2), 317–329.

19 Carillo, R.E., McEwen, J.D., Van de Ville, D., Thiran, J.P., and Wiaux, Y. (2013) Sparsity averaging for compressive imaging. *IEEE Signal Processing Letters*, **20** (6), 591–594.

20 Candes, E.J., Wakin, M.B., and Boyd, S.P. (2008) Enhancing sparsity by reweighted $ℓ_1$ minimization. *Journal of Fourier Analysis and Applications*, **14**, 877–905.

21 Fuchs, J.J. (2007) Convergence of a sparse representations algorithm applicable to real or complex data. *IEEE Transactions on Signal Processing*, **1** (4), 598–605.

22 Gorodnitsky, I.F. and Rao, B.D. (1997) Sparse signal reconstruction from limited data using focuss : a re-weighted minimum norm algorithm. *IEEE Transactions on Signal Processing*, **45** (3), 600–616.

23 Weiszfeld, E. and Plastria, F. (2009) On the point for which the sum of the distances to n given points is minimum. *Ann Oper Res*, **167**, 7–41.

24 Byrd, R.H. and Payne, D.A. (1979) Convergence of the iteratively reweighted least squares algorithm for robust regression, *Tech. Rep. 131*, The Johns Hopkins University, Baltimore, MD.

25 Rao, B.D., Engan, K., Cotter, S.F., Palmer, J., and Kreutz-Delgao, K. (2003) Subset selection in noise based on diversity measure minimization. *IEEE Transactions on Signal Processing*, **51** (3), 760–770.

26 Allain, M., Idier, J., and Goussard, Y. (2006) On global and local convergence of half-quadratic algorithms. *IEEE Transactions on Image Processing*, **15** (5), 1130–1142.

27 Charbonnier, P., Blanc-Féraud, L., Aubert, G., and Barlaud, M. (1997) Deterministic edge-preserving regularization in computed imaging. *IEEE Transactions on Image Processing*, **6**, 298–311.

28 Chan, T.F. and Mulet, P. (1999) On the convergence of the lagged diffusivity

fixed point method in total variation image restoration. *SIAM Journal on Numerical Analysis*, **36** (2), 354–367.

29 Geman, D. and Reynolds, G. (1992) Constrained restoration and the recovery of discontinuities. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, **14** (3), 367–383.

30 Geman, D. and Yang, C. (1995) Nonlinear image recovery with half-quadratic regularization. *IEEE Transactions on Image Processing*, **4** (7), 932–946.

31 Idier, J. (2001) Convex half-quadratic criteria and interacting auxiliary variables for image restoration. *IEEE Transactions on Image Processing*, **10** (7), 1001–1009.

32 Nikolova, M. and Ng, M.K. (2005) Analysis of half-quadratic minimization methods for signal and image recovery. *SIAM Journal on Scientific Computing*, **27**, 937–966.

33 Beck, A., Teboulle, M., and Chikishev, Z. (2008) Iterative minimization schemes for solving the single source localization problem. *SIAM Journal on Optimization*, **19** (3), 1397–1416.

34 Bissantz, N., Dumbgen, L., Munk, A., and Stratmann, B. (2009) Convergence analysis of generalized iteratively reweighted least squares algorithms on convex function spaces. *SIAM Journal on Optimization*, **19** (4), 1828–1845.

35 Lefkimmiatis, S., Bourquard, A., and Unser, M. (2012) Hessian-based norm regularization for image restoration with biomedical applications. *IEEE Transactions on Image Processing*, **21** (3), 983–995.

36 Fish, D.A., Walker, J.G., Brinicombe, A.M., and Pike, E.R. (1995) Blind deconvolution by means of the richardson-lucy algorithm. *J. Opt. Soc. Am. A*, **12**, 58–65.

37 Vardi, Y., Shepp, L.A., and Kaufman, L. (1985) A statistical model for positron emission tomography (with discussion). *Journal of American Statistical Association*, **80**, 8–38.

38 Lanteri, H., Roche, M., and Aime, C. (2002) Penalized maximum likelihood image restoration with positivity constraints : multiplicative algorithms. *Inverse Problems*, **18**, 1397–1419.
39 Bertsekas, D.P. (1999) *Nonlinear Programming*, Athena Scientific, Belmont, MA, 2nd edn.

40 Ning, X., Selesnick, I.W., and Duval, L. (2014) Chromatogram baseline estimation and denoising using sparsity (beads). *Chemometrics and Intelligent Laboratory Systems*, **139**, 156–167.

41 Jiang, M. and Wang, W. (2003) Convergence of the simultaneous algebraic reconstruction technique (SART). *IEEE Transactions on Image Processing*, **12** (8), 957–961.

42 Zibulevsky, M. and Elad, M. (2010) $l_2 - l_1$ optimization in signal and image processing. *IEEE Signal Processing Magazine*, **27**, 76–88.

43 De Pierro, A.R. (1995) A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography. *IEEE Trans. Med. Imag.*, **14** (1), 132–137.

44 Labat, C. and Idier, J. (2007) Convergence of truncated half-quadratic and newton algorithms, with application to image restoration, *Tech. Rep.*, IRCCyN, Nantes, France.

45 Chouzenoux, E., Idier, J., and Moussaoui, S. (2011) A Majorize-Minimize strategy for subspace optimization applied to image restoration. *IEEE Transactions Image Processing*, **20** (18), 1517–1528.

46 Chouzenoux, E., Jezierska, A., Pesquet, J.C., and Talbot, H. (2013) A majorize-minimize subspace approach for $l_2-l_0$ image regularization. *SIAM Journal Imaging Science*, **6** (1), 563–591.

47 Chouzenoux, E. and Pesquet, J.C. (2017) A stochastic majorize-minimize subspace algorithm for online penalized least squares estimation. *IEEE Transactions on Signal Processing*, **65** (18), 4770–4783.

48 Hager, W.W. and Zhang, H. (2006) A survey of nonlinear conjugate gradient methods. *Pacific Journal on Optimization*, **2** (1), 35–58.

49 Liu, D.C. and Nocedal, J. (1989) On the limited memory BFGS method for large scale optimization. *Mathematical Programming*, **45** (3), 503–528.

50 Chouzenoux, E., Pesquet, J.C., Talbot, H., and Jezierska, A. (2011) A memory gradient algorithm for $l_2-l_0$ regularization with applications to image restoration, in *18th IEEE International Conference on Image Processing (ICIP 2011)*, Brussels, Belgium, pp. 2717–2720.

51 Florescu, A., Chouzenoux, E., Pesquet, J.C., Ciuciu, P., and Ciochina, S. (2014) A majorize-minimize memory gradient method for complex-valued inverse problem. *Signal Processing*, **103**, 285–295. Special issue on Image Restoration and Enhancement: Recent Advances and Applications.

52 Miele, A. and Cantrell, J.W. (1969) Study on a memory gradient method for the minimization of functions. *Journal on Optimization Theory and Applications*, **3** (6), 459–470.

53 Chouzenoux, E. and Pesquet, J.C. (2016) Convergence rate analysis of the majorize-minimize subspace algorithm. *IEEE Signal Processing Letters*, **23** (9), 1284–1288.

54 Combettes, P.L. and Vū, B.C. (2014) Variable metric forward-backward splitting with applications to monotone inclusions in duality. *Optimization: A Journal of Mathematical Programming and Operations Research*, **63** (9), 1289–1318.

55 Chouzenoux, E., Pesquet, J.C., and Repetti, A. (2014) Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function. *Journal of Optimization Theory and Applications*, **162** (1), 107–132.

56 Combettes, P.L. and Wajs, V.R. (2005) Signal recovery by proximal forward-backward splitting. *Multiscale Modeling and Simulation*, **4** (4), 1168–1200.

57 Daubechies, I., Defrise, M., and De Mol, D. (2004) An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. *Communications on Pure and Applied Mathematics*, **57**, 1413–1457.

58 Combettes, P., Dung, D., and Vu, B. (2011) Proximity for sums of composite functions. *J. Math. Anal. Appl.*, **380** (2), 680–688.

59 Abboud, F., Chouzenoux, E., Pesquet, J.C., Chenot, J.H., and Laborelli, L.
(2015) A dual block coordinate proximal algorithm with application to deconvolution of interlaced video sequences, in *Proceedings of the 22nd IEEE International Conference on Image Processing (ICIP 2015)*, Quebec City, Canada, pp. 4917–4921.

60 Chambolle, A. and Pock, T. (2015) A remark on accelerated block coordinate descent for computing the proximity operators of a sum of convex functions. *SMAI-Journal of computational mathematics*, **1**, 29–54.

61 Jaggi, M., Smith, V., Takac, M., Terhorst, J., Krishnan, S., Hofmann, T., and Jordan, M.I. (2014) Communication-efficient distributed dual coordinate ascent, in *Advances in Neural Information Processing Systems 27* (eds Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Weinberger), Curran Associates, Inc., pp. 3068–3076.

62 Abboud, F., Chouzenoux, E., Pesquet, J.C., Chenot, J.H., and Laborelli, L. (2015) A distributed strategy for computing proximity operators, in *Proceedings of the 49th Asilomar Conference on Signals, Systems and Computers (ASILOMAR 2015)*, pp. 396–400.

63 Bertsekas, D.P. (1981) Projected Newton methods for optimization problems with simple constraints. *SIAM Journal Control and Optimization*, **20**, 762–767.

64 Bonettini, S., Zanella, R., and Zanni, L. (2009) A scaled gradient projection method for constrained image deblurring. *Inverse Problems*, **25** (1), 015 002+.

65 Iusem, A.N. (2003) On the convergence properties of the projected gradient method for convex optimization. *Computational Applied Mathematics*, **22** (1), 37–52.

66 Cherni, A., Chouzenoux, E., Duval, L., and Pesquet, J.C. (2019) A novel smoothed norm ratio for sparse signal restoration application to mass spectrometry, in *Proceedings of Signal Processing with Adaptive Sparse Structured Representations (SPARS 2019)*, Toulouse, France.

67 Cherni, A., Chouzenoux, E., Duval, L., and Pesquet, J.C. (2020) SPOQ Ip-Over-lq

regularization for sparse signal recovery applied to mass spectrometry, *Tech. Rep.*. Https://arxiv.org/abs/2001.08496.

68 Jacobson, M.W. and Fessler, J.A. (2007) An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms. *IEEE Transactions on Image Processing*, **16** (10), 2411–2422.

69 Razaviyayn, M., Hong, M., and Luo, Z. (2013) A unified convergence analysis of block successive minimization methods for nonsmooth optimization. *SIAM Journal on Optimization*, **23** (2), 1126–1153.

70 Sothivirat, S. and Fessler, J.A. (2002) Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms. *IEEE Transactions on Signal Processing*, **11** (3), 306–317.

71 Chouzenoux, E., Pesquet, J.C., and Repetti, A. (2016) A block coordinate variable metric forward-backward algorithm. *Journal on Global Optimization*, **66** (3), 457–485.

72 Repetti, A., Chouzenoux, E., and Pesquet, J.C. (2014) A preconditioned forward-backward approach with application to large-scale nonconvex spectral unmixing problems, in *Proceedings of the 39th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2014)*, Firenze, Italy, pp. 1498–1502.

73 Repetti, A., Pham, M.Q., Duval, L., Chouzenoux, E., and Pesquet, J.C. (2015) Euclid in a taxicab: Sparse blind deconvolution with smoothed 11/12 regularization. *IEEE Signal Processing Letters*, **22** (5), 539–543.

74 Abboud, F., Chouzenoux, E., Pesquet, J., Chenot, J., and Laborelli, L. (2019) An alternating proximal approach for blind video deconvolution. *Signal Processing: Image Communication*, **70**, 21–36.

75 Gabay, D. and Mercier, B. (1976) A dual algorithm for the solution of nonlinear variational problems via finite elements approximations. *Computers and Mathematics with Applications*, **2**, 17–40.

76 Fortin, M. and Glowinski, R. (eds) (1983) *Augmented Lagrangian Methods: Applications to the Numerical Solution of*
*Boundary-Value Problems*, Elsevier
Science Ltd, Amsterdam: North-Holland.

77 Figueiredo, M.A.T. and Bioucas-Dias, J.M. (2010) Restoration of Poissonian images using alternating direction optimization. *IEEE Transactions on Image Processing*, **19** (12), 3133–3145.

78 Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2011) Distributed optimization and statistical learning via the alternating direction method of multipliers. *Foundations and Trends in Machine Learning*, **3** (1), 1–222.

79 Lions, P.L. and Mercier, B. (1979) Splitting algorithms for the sum of two nonlinear operators. *SIAM Journal on Numerical Analysis*, **16**, 964–979.

80 Combettes, P.L. and Pesquet, J.C. (2007) A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery. *IEEE Journal of Selected Topics in Signal Processing*, **1** (4), 564–574.

81 Chen, C., He, B., Ye, Y., and Yuan, X. (2014) The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent. *Mathematical Programming*, pp. 1–23.

82 Setzer, S., Steidl, G., and Teuber, T. (2010) Deblurring poissonian images by split bregman techniques. *Journal on Visual Communication and Image Representation*, **21** (3), 193–199.

83 Combettes, P.L. and Pesquet, J.C. (2008) A proximal decomposition method for solving convex variational inverse problems. *Inverse Problems*, **24** (6).

84 Pesquet, J.C. and Pustelnik, N. (2012) A parallel inertial proximal optimization method. *Pacific Journal of Optimization*, **8** (2), 273–305.

85 Cherni, A., Chouzenoux, E., and Delsuc, M.A. (2016) PALMA, an improved algorithm for DOSY signal processing. *Analyst*, **142** (5), 772–779.

86 Cherni, A., Chouzenoux, E., and Delsuc, M.A. (2016) Proximity operators for a class of hybrid sparsity+entropy priors. application to DOSY NMR signal reconstruction, in *Proceedings of the International Symposium on Signal, Image, Video and Communications (ISIVC 2016)*, Tunis, Tunisia.

87 Chambolle, A. and Pock, T. (2011) A

first-order primal-dual algorithm for
convex problems with applications to
imaging. *Journal of Mathematical
Imaging and Vision*, **40** (1), 120–145.

88 Esser, E., Zhang, X., and Chan, T. (2010) A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science. *SIAM Journal on Imaging Sciences*, **3** (4), 1015–1046.

89 Komodakis, N. and Pesquet, J.C. (2014) Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems. *IEEE Signal Processing Magazine*, **32** (6), 31–54.

90 Combettes, P.L. and Pesquet, J.C. (2012) Primal-dual splitting algorithm for solving inclusions with mixtures of composite, lipschitzian, and parallel-sum type monotone operators. *Set-Valued and Variational Analysis*, **20** (2), 307–330.

91 Condat, L. (2013) A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms. *Journal of Optimization, Theory and Applications*, **158** (2), 460–479.

92 Vũ, B.C. (2013) A splitting algorithm for dual monotone inclusions involving cocoercive operators. *Advances in Computational Mathematics*, **38** (3), 667–681.

93 Raguet, H., Fadili, J., and Peyré, G. (2013) A generalized forward-backward splitting. *SIAM Journal on Imaging Sciences*, **6** (3), 1199–1226.

94 Bricenos-Arias, L., Chierchia, G., Chouzenoux, E., and Pesquet, J.C. (2019) A random block-coordinate douglas-rachford splitting method with low computational complexity for binary logistic regression. *Computational Optimization and Applications*, **72** (3), 707–726.

95 Pesquet, J.C. and Repetti, A. (2015) A class of randomized primal-dual algorithms for distributed optimization. *Journal of nonlinear and convex analysis*, **16** (12), 2453–2490.

96 Cherni, A., Chouzenoux, E., and Delsuc, M.A. (2018) Fast dictionary-based approach for mass spectrometry data
analysis, in *Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2018)*, Calgary, Canada.

97 Fiacco, A.V. and McCormick, G.P. (1967) The sequential unconstrained minimization technique (SUMT) without parameters. *Operations Research*, **15** (5), 820–827.

98 Wright, S.J. (1997) *Primal-dual interior-point methods*, SIAM, Philadelphia, PA, 1st edn.

99 Forsgren, A., Gill, P.E., and Wright, M.H. (2002) Interior methods for nonlinear optimization. *SIAM Review*, **44** (4), 525–597.

100 Moussaoui, S., Chouzenoux, E., and Idier, J. (2012) Primal-dual interior point optimization for penalized least squares estimation of abundance maps in hyperspectral imaging, in *Proceedings of the 4th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS 2012)*, Shangai, China.

101 Chouzenoux, E., Legendre, M., Moussaoui, S., and Idier, J. (2012) Fast constrained least squares spectral unmixing using primal-dual interior-point optimization. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, **7** (1), 59–69.

102 Armand, P., Gilbert, J.C., and Jan-Jégou, S. (2000) A feasible BFGS interior point algorithm for solving strongly convex minimization problems. *SIAM Journal on Optimization*, **11**, 199–222.

103 Wright, M.H. (1994) Some properties of the Hessian of the logarithmic barrier function. *Mathematical Programming*, **67** (2), 265–295.

104 Wright, M.H. (1998) Ill-conditioning and computational error in interior methods for nonlinear programming. *SIAM Journal on Optimization*, **9** (1), 84–111.

105 Friedlander, M.P. and Orban, D. (2012) A primal-dual regularized interior-point method for convex quadratic programs. *Mathematical Programming Computation*, **4** (1), 71–107.

106 Conn, A., Gould, N., and Toint, P.L. (1996) A primal-dual algorithm for minimizing a nonconvex function subject

to bounds and nonlinear constraints, in *Nonlinear Optimization and Applications* (eds G. Di Pillo and F. Giannessi), Kluwer Academic Publishers, 2nd edn.

107 Armand, P., Benoist, J., and Dussault, J. (2012) Local path-following property of inexact interior methods in nonlinear programming. *Computational Optimization and Applications*, **52** (1), 209–238.

108 Bonettini, S., Galligani, E., and Ruggiero, V. (2007) Inner solvers for interior point methods for large scale nonlinear programming. *Computational Optimization and Applications*, **37** (1), 1–34.

109 Legendre, M., Moussaoui, S., Chouzenoux, E., and Idier, J. (2014) Primal-dual interior-point optimization based on majorization-minimization for edge preserving spectral unmixing, in *Proceedings of the 21st IEEE International Conference on Image Processing (ICIP 2014)*, Shangai, China, pp. 4161–4165.

110 Legendre, M., Moussaoui, S., Schmidt, F., and Idier, J. (2013) Parallel implementation of a primal-dual interior-point optimization method for fast abundance maps estimation, in *Proceedings of the 5th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS 2013)*, Gainesville, FL, USA.

111 Murray, W. and Wright, M.H. (1994) Line search procedures for the logarithmic barrier function. *SIAM Journal on Optimization*, **4** (2), 229–246.

112 Chouzenoux, E., Moussaoui, S., and Idier, J. (2012) Majorize-minimize linesearch for inversion methods involving barrier function optimization. *Inverse Problems*, **28** (6), 065011.

113 Chouzenoux, E., Moussaoui, S., and Idier, J. (2011) Efficiency of line search strategies in interior point methods for linearly constrained optimization, in *Proceedings of the IEEE Workshop on Statistical Signal Processing (SSP 2011)*, Nice, France, pp. 101–104.

114 Chouzenoux, E., Moussaoui, S., Idier, J., and Mariette, F. (2013) Primal-dual
interior point optimization for a
regularized reconstruction of NMR
relaxation time distributions, in
*Proceedings of the 38th IEEE
International Conference on Acoustics,
Speech, and Signal Processing (ICASSP
2013)*, Vancouver, Canada, pp.
8747–8750.

115 El-Bakry, A.S., Tapia, R.A., Tsuchiya, T., and Zhang, Y. (1996) On the formulation and theory of the Newton interior-point method for nonlinear programming. *Journal of Optimization Theory and Applications, **89***, 507–541.

116 Johnson, C.A., Seidel, J., and Sofer, A. (2000) Interior-point methodology for 3-D PET reconstruction. *IEEE Transactions on Medical Imaging, **19*** (4).

117 Chouzenoux, E., Moussaoui, S., Idier, I., and Mariette, F. (2010) Efficient maximum entropy reconstruction of nuclear magnetic resonance t1-t2 spectra. *IEEE Transactions on Signal Processing, **58*** (12), 6040–6051.

118 Kaipin, X. and Shanmin, Z. (2013) Trust-region algorithm for the inversion of molecular diffusion nmr data. *Analytical chemistry, **86*** (1), 592–599.

119 Johnson, C.S. (1999) Diffusion ordered nuclear magnetic resonance spectroscopy: principles and applications. **34**, 203–256.

120 Beck, A. and Teboulle, M. (2009) A fast

iterative shrinkage-thresholding algorithm
for linear inverse problems. *SIAM Journal
on Imaging Sciences, **2*** (1), 183–202.

121 Pustelnik, N., Benazza-Benhayia, A., Zheng, Y., and Pesquet, J.C. (2016) Wavelet-based image deconvolution and reconstruction. *Wiley Encyclopedia of Electrical and Electronics Engineering*.

122 Lange, K. (1995) A gradient algorithm locally equivalent to the em algorithm. *Journal of the Royal Statistical Society. Series B (Methodological)*, pp. 425–437.

123 Varadhan, R. and Roland, C. (2008) Simple and globally convergent methods for accelerating the convergence of any em algorithm. *Scandinavian Journal of Statistics, **35*** (2), 335–353.

124 Pock, T. and Chambolle, A. (2011) Diagonal preconditioning for first order primal-dual algorithms in convex optimization, in *IEEE International Conference on Computer Vision (ICCV)*, Barcelona, Spain, pp. 1762–1769.

125 Repetti, A., Chouzenoux, E., and Pesquet, J.C. (2015) A random block-coordinate primal-dual proximal algorithm with application to 3d mesh denoising, in *Proceedings of the 40th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP 2015)*, Brisbane, Australia, pp. 3561–3565.