# Reinterpreting EMML as Mirror Descent for Constrained Maximum Likelihood Estimation

Antonin Clerc, Ségolène Martin, Nicolas Papadakis, Gabriele Steidl

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

To cite this version:

Antonin Clerc, Ségolène Martin, Nicolas Papadakis, Gabriele Steidl. Reinterpreting EMML as Mirror Descent for Constrained Maximum Likelihood Estimation. 2026. hal-05534462

HAL Id: hal-05534462

https://hal.science/hal-05534462v1

Preprint submitted on 3 Mar 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Reinterpreting EMML as Mirror Descent for Constrained Maximum Likelihood Estimation

Antonin Clerc¹,², Ségolène Martin¹, Nicolas Papadakis², Gabriele Steidl¹

¹ Institute of Mathematics, Technische Universität Berlin, Germany

² Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France

antonin.clerc@math.u-bordeaux.fr, nicolas.papadakis@math.u-bordeaux.fr

martin@math.tu-berlin.de, steidl@math.tu-berlin.de

**Abstract**— The Expectation–Maximization Maximum Likelihood (EMML) algorithm belongs to the Expectation–Maximization family and is widely used for image reconstruction problems under Poisson noise. In this paper, we reinterpret EMML as a mirror descent method applied to a reparametrized objective function. This perspective allows us to incorporate convex constraints into the algorithm through appropriately chosen Bregman projections, while preserving the multiplicative structure of the EMML updates to ensure computational efficiency. We then establish the convergence of the resulting algorithm toward a solution of the constrained maximum-likelihood problem. Numerical experiments on hyperspectral unmixing problems demonstrate that the constrained EMML converges in fewer iterations than the classical EMML.

## 1 Introduction

Image reconstruction problems under Poisson noise naturally arise in scientific imaging domains, where measurements are photon-limited. Prominent examples include fluorescence microscopy [1,2], astronomical imaging [3,4], hyperspectral image restoration [5], and Positron Emission Tomography [6]. In these applications, photon counts are inherently low and follow discrete statistics.

Mathematically, the goal is to recover an unknown non-negative image $x \in \mathbb{R}_+^n$ from Poisson-degraded measurements $y \in \mathbb{R}_+^m$, modeled as

$$y \sim \text{Poisson}(Ax), \quad (1)$$

where $A \in \mathbb{R}_+^{m \times n}$ is a known degradation operator, typically ill-conditioned. Under Poisson noise, the negative log-likelihood is $F(x) := \text{KL}(y, Ax)$ where KL is the Kullback–Leibler divergence defined for all $u \in \mathbb{R}_{\ge 0}^m$ and $v \in \mathbb{R}_{> 0}^m$

$$\text{KL}(u, v) := \sum_{i=1}^{m} \left[ u_i \log \frac{u_i}{v_i} + v_i - u_i \right],$$

with the convention $0 \log 0 = 0$. To ensure $(Ax)_i > 0$ for all $i$, we assume that each row $a_i$ of A satisfies $a_i \in \mathbb{R}_{\ge 0}^n \setminus \{0\}$, and $x \in \mathbb{R}_{> 0}^n$. These conditions are generally satisfied by standard operators such as convolutions with smoothing filters.

Minimizing the Poisson log-likelihood is numerically challenging due to the lack of Lipschitz continuity of the gradient and the singular behavior near the boundary of $\mathbb{R}_{> 0}^n$,

which prevents the direct use of standard gradient descent methods. Early iterative reconstruction methods minimizing the KL divergence were introduced in several works [6–9] based on the Expectation–Maximization (EM) framework. In particular, the Richardson–Lucy algorithm [7,8] was originally proposed for linear image deconvolution under Poisson noise. Later, Byrne [9] and Shepp and Vardi [6] independently derived the same multiplicative updates within the EM formalism, introducing the Expectation–Maximization Maximum Likelihood (EMML) and Maximum Likelihood Expectation–Maximization (MLEM) algorithms, respectively. Although introduced in different communities under different names, these methods correspond to the same EM iterations for Poisson maximum-likelihood estimation. In the remainder of this work, we refer to all these algorithms as EMML.

The EMML iterations to minimize $F$ are given by

$$x^{(k+1)} = x^{(k)} \odot \frac{A^{\top} \frac{y}{Ax^{(k)}}}{A^{\top} \mathbf{1}}. \quad (2)$$

These iterations are simple, preserve non-negativity, monotonically increase the likelihood [9]. EMML typically provides high-quality reconstruction after only a few iterations, but the quality degrades afterward. Nevertheless, many works continue to rely on EMML [10, 11]. To improve reconstruction quality and mitigate the need for early stopping, several recent approaches extend the EMML iterations by incorporating regularization, see, e.g. [12–16], or use them to initialize a neural network [17].
The focus of this paper is on incorporating constraints into the EMML algorithm while preserving its fast multiplicative updates. Many imaging problems naturally involve the minimization of a convex objective under convex constraints [18, 19]. Yet, standard EMML iterations do not directly accommodate these constraints, motivating the developments proposed in this work. Existing constrained variants of EMML as [20] only take into account averages over selected pixels, corresponding to constraints of the form $\sum_{i \in \mathcal{I}} x_i = y$, and do not cover general constraints.

**Contributions** Building on the connection between EM algorithms for the exponential family and Mirror Descent (MD), we reinterpret EMML as a particular instance of MD, yielding a natural framework for incorporating general convex constraints while preserving both convergence guarantees and the characteristic multiplicative structure of the updates. Section 2 recalls MD results. Section 3 shows that EMML can be interpreted as MD on a reparametrized objective. Section 4 extends EMML to constrained problems using Bregman projections, preserving convergence guarantees. Section 5 presents numerical experiments on hyperspectral image restoration under Poisson noise, demonstrating improved PSNR and faster convergence with the constrained EMML. Conclusions are drawn in Section 6.

## 2 Background on Mirror Descent

Let $F: \mathbb{R}^d \to \mathbb{R} \cup \{+\infty\}$ be a convex and differentiable function on $\text{int}(\text{dom } F)$. Let $u: \mathbb{R}^d \to \mathbb{R} \cup \{+\infty\}$ be a strictly convex and differentiable function on $\text{int}(\text{dom } u)$. The mapping $\nabla u$ is referred to as a *mirror map*. For all $(x, y) \in \text{dom } u \times \text{int}(\text{dom } u)$, the Bregman divergence associated with the potential $u$ is defined by

$$D_u(x, y) = u(x) - u(y) - \langle \nabla u(y), x - y \rangle,$$

MD can be derived as a Bregman proximal gradient step. Given an iterate $x^{(k)} \in \text{int}(\text{dom } u)$ and a step size $\tau > 0$, the next iterate is defined as

$$x^{(k+1)} = \arg \min_{x \in \mathbb{R}^d} \left\{ \langle \nabla F(x^{(k)}), x - x^{(k)} \rangle + \frac{1}{\tau} D_u(x, x^{(k)}) \right\}. \quad (3)$$

The first-order optimality condition of (3) leads to the MD update

$$x^{(k+1)} = \nabla u^* \left( \nabla u(x^{(k)}) - \tau \nabla F(x^{(k)}) \right),$$

where $u^*$ denotes the convex conjugate of $u$ [21]. To ensure convergence of the iterates when the objective function is not Lipschitz continuous, we make the following assumptions:

**Assumption 1.**

1. $F$ is proper, lower semicontinuous, and convex, with $\text{dom } u \subseteq \text{dom } F$, and differentiable on $\text{int}(\text{dom } u)$.

2. $u$ is of Legendre type [22, Definition 1].

3. $F$ is relatively smooth with respect to $u$, i.e., there exists $L > 0$ such that $Lu - F$ is convex on $\text{int}(\text{dom } u)$.

4. For every $y \in \text{int}(\text{dom } u)$, the Bregman proximal gradient mapping

$$T_{\tau}(y) := \arg \min_{x \in \mathbb{R}^d} \left\{ \langle \nabla F(y), x - y \rangle + \frac{1}{\tau} D_u(x, y) \right\} \quad (4)$$

is nonempty, single-valued, and maps $\text{int}(\text{dom } u)$ into itself.

The Assumption is satisfied for $F(x) = \text{KL}(y, Ax)$ when using either Burg's entropy [17, 22] or Shannon entropy [23] as the mirror map. We recall below standard convergence results for MD following [22, 24, 25].

**Lemma 2** (Adapted from [22, Lemma 5]). Suppose that Assumption 1 holds. Let $\tau = 1/L$ and, for all $x \in \text{int}(\text{dom } u)$,

$$x^+ := T_{\tau}(x).$$

where $T_{\tau}$ is given by (4). Then, for all $s \in \text{dom } u$,

$$\tau(F(x^+) - F(s)) \leq D_u(s, x) - D_u(s, x^+).$$

**Theorem 3** ([22, Theorem 2]). Suppose that Assumption 1 holds. If $\tau \le 1/L$ and the solution set of

$$(\mathcal{P}) \quad \min\{F(x) : x \in \overline{\text{dom } u}\}$$

is nonempty and compact, then any limit point of the sequence $(x^{(k)})_{k \in \mathbb{N}}$ generated by the MD iteration is a solution of $(\mathcal{P})$.

## 3 The EMML Algorithm

Following the Poisson model (1), the EMML algorithm (2) computes the maximum-likelihood estimate of $x$ by minimizing the negative log-likelihood $F(x) := \text{KL}(y, Ax)$.

### 3.1 EMML as EM algorithm

To derive the EM algorithm, we introduce latent variables $z_{ij} \in \mathbb{N}$ representing the number of photons emitted from pixel $j$ and detected by detector $i$ as

$$z_{ij} \sim \text{Poisson}(A_{ij}x_j), \quad y_i = \sum_{j=1}^{n} z_{ij},$$

where the observed counts are obtained by marginalization. The joint complete-data likelihood can then be written as

$$p(z, y | x) = \frac{1}{\prod_{i,j} z_{ij}!} \exp \left( \sum_{i,j} z_{ij} \log(A_{ij}x_j) - A_{ij}x_j \right). \quad (5)$$

We now apply the EM algorithm.

**E-step:** We compute the expected complete-data negative log-likelihood:

$$Q_{x^{(r)}}(x) = -\mathbb{E}_{p(z|y, x^{(r)})} \left[ \log p(z, y | x) \right] \\ \propto - \sum_{i,j} \left( -A_{ij}x_j + y_i \frac{A_{ij}x_j^{(r)}}{(Ax^{(r)})_i} \log(A_{ij}x_j) \right).$$

Up to an additive constant independent of $x$, this expression coincides with $\text{KL}(y, Ax)$.
**M-step:** We minimize $Q_{x^{(r)}}$ with respect to $x$. Setting the gradient with respect to $x_j$ to zero yields

$$0 = \sum_{i} \left( A_{ij} - \frac{y_i A_{ij} x_j^{(r)}}{x_j (Ax^{(r)})_i} \right).$$

Solving for $x_j$ leads to the classical EMML iteration (2).

## 3.2 EMML as MD

EM applied to exponential families is known to admit a MD interpretation [23,26]. Since the Poisson distribution belongs to the exponential family [27], this framework applies to our setting. In particular, the joint complete-data log-likelihood (5) can be rewritten in the canonical exponential-family form as

$$p(z, y | x) = \frac{1}{h(z)} \exp \left( \langle S(z, y), \theta \rangle - \mathcal{A}(\theta) \right).$$

The corresponding components are:

* the natural parameter $\theta_j = \log x_j$, $j = 1, \dots, n$,

* the sufficient statistic $S(z, y)_j = \sum_i (1 + \log A_{ij}) z_{ij}$,

* the log-partition function $\mathcal{A}(\theta) = \sum_{i,j} A_{ij} e^{\theta_j} = \mathbf{1}^{\top} A e^{\theta}$,

* and the base measure

$$h(z) = \exp \left( \sum_{i,j} z_{ij} \log A_{ij} - \log(z_{ij}!) \right).$$

This formulation is central for interpreting EMML as a MD-type algorithm. By exploiting the mirror map and the reparametrization derived in [23], we obtain the following proposition for EMML.

**Proposition 4.** Define $\mathcal{L}: \mathbb{R}^n \to \mathbb{R}$ by $\mathcal{L}(\theta) = \text{KL}(y, Ae^{\theta})$ and assume that every row and every column of $A$ contains at least one strictly positive entry. Consider MD applied to $\mathcal{L}$ with mirror map

$$\mathcal{A}(\theta) = \sum_{i,j} A_{ij} e^{\theta_j}.$$

Then the corresponding update

$$\theta^{(r+1)} = \nabla \mathcal{A}^{-1} (\nabla \mathcal{A}(\theta^{(r)}) - \nabla \mathcal{L}(\theta^{(r)}))$$

is equivalent to the classical EMML iteration, up to the reparametrization $x^{(r)} = \exp(\theta^{(r)})$.

The convergence of the EMML algorithm was established earlier in [6, 9]. Convergence rates for EM algorithms associated with exponential family models were subsequently analyzed in [23, Proposition 2] and [26, Proposition 11]. The sublinear convergence of EMML is a direct consequence of Lemma 2 when viewed through the lens of MD theory.

**Proposition 5.** For any $N \in \mathbb{N}^*$ and any initial point $x^0 \in \mathbb{R}^n$, the EMML iterates satisfy

$$\text{KL}(y, Ax^N) - \text{KL}(y, Ax^*) \le \frac{D_A(\log(x^*), \log(x^0))}{N},$$

where $x^*$ is any minimizer of $x \mapsto \text{KL}(y, Ax)$.

EMML exhibits a worst-case sublinear $\mathcal{O}(1/N)$ decrease of the objective value. Faster local convergence rates have been established in the literature; for instance, [23] proved superlinear convergence in a neighborhood of a solution under certain assumptions. However, these improved rates rely on the evaluation of additional terms, such as the missing information matrix and its eigenvalues, which are generally inaccessible or costly to compute.

## 4 Constrained EMML

Interpreting the EMML algorithm as a MD allows incorporating constraints on the iterates while preserving convergence guarantees. Let $C \subset \mathbb{R}^m$ be a convex set. We make the following assumption:

**Assumption 6.** Each column and each row of $A$ contains at least one strictly positive entry and $C \cap \mathbb{R}^n_{>0}$ is nonempty.

We consider the constrained optimization problem

$$(\mathcal{P}_C) \quad \min_{x \in \mathbb{R}^m_{\ge 0}} \text{KL}(y, Ax) \quad \text{subject to} \quad x \in C.$$

As shown in Proposition 4, EMML is equivalent to MD applied to the variable $\theta = \log(x)$. A natural way to enforce constraints is thus to apply projected MD in the $\theta$-variable. Define the constraint set in the dual variables by

$$C_{\theta} := \{\theta \in \mathbb{R}^m : e^{\theta} \in C\}.$$

Denoting as $\iota_{C_{\theta}}$ the indicator function of $C_{\theta}$, we target

$$\min_{\theta} \text{KL}(y, Ae^{\theta}) + \iota_{C_{\theta}}(\theta).$$

The Bregman proximal gradient [22] update with respect to the Legendre function $\mathcal{A}$ reads

$$\theta^{(r+1)} = \text{prox}_{\iota_{C_{\theta}}}^{\mathcal{A}} \left( \nabla \mathcal{A}^{-1} \left( \nabla \mathcal{A}(\theta^{(r)}) - \tau \nabla \text{KL}(y, Ae^{\theta^{(r)}}) \right) \right).$$

Recall that the Bregman proximal operator of an indicator function reduces to a Bregman projection [28, 29]:

$$\text{prox}_{\iota_{C_{\theta}}}^{\mathcal{A}}(\theta^{(0)}) = \arg \min_{\theta \in C_{\theta}} D_{\mathcal{A}}(\theta, \theta^{(0)}), \quad (6)$$

where $D_{\mathcal{A}}$ is the Bregman divergence associated with $\mathcal{A}$. In practice, for many constrained sets defined for the primary variable $x$ (e.g., linear constraints), computing the Bregman projection in the $\theta$-variable may be difficult.

**Proposition 7.** Let $\theta^*$ denote the Bregman projection (6) and define $x = e^{\theta}$, $x^* = e^{\theta^*}$, $x^{(0)} = e^{\theta^{(0)}}$. Let

$$w_j := \sum_i A_{i,j} > 0, \quad \tilde{w}_j := w_j x_j^{(0)}, \quad u(x) = - \sum_j \tilde{w}_j \log x_j.$$

Then, solving (6) is equivalent to the Bregman projection

$$x^* = \arg \min_{x \in C} D_u(x, x^{(0)}),$$

where $D_u$ is the Bregman divergence associated with $u$.
Proof. The problem (6) can be rewritten as

$$x^* \in \arg\min_{x \in \mathbb{R}_{>0}^m} D_A(\log x, \log x^{(0)}) \quad \text{s.t.} \quad x \in C.$$

A direct computation shows that, up to a constant independent of $x$, the Bregman divergence $D_A(\log x, \log x^{(0)})$ writes

$$D_A(\log x, \log x^{(0)}) \propto \sum_j w_j x_j - \sum_j w_j x_j^{(0)} \log x_j,$$

which coincides, up to constants, with the Bregman divergence $D_u(x, x^{(0)})$ associated with $u(x) = -\sum_j \tilde{w}_j \log x_j$ and $\tilde{w}_j = w_j x_j^{(0)}$. The Bregman projection in the $\theta$-variable is thus equivalent to the Bregman projection of $x$ onto $C$ with divergence $D_u$. □

The resulting constrained EMML algorithm is summarized in Algorithm 1. The theoretical convergence rate of Algorithm 1 can also be established.

**Algorithm 1** Constrained EMML via Projected MD

**Input:** A, y, M, y satisfying Assumption 6; initial point $x^{(0)} \in \mathbb{R}_{>0}^m$; weights $w_j := \sum_i A_{i,j}$.

1: **for** $r = 0, 1, 2, \dots$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp; $\tilde{x}^{(r)} = x^{(r)} \odot \frac{A^{\top}(y/(Ax^{(r)}))}{A^{\top}1}$ ▷ EMML update

3: &nbsp;&nbsp;&nbsp;&nbsp; Define $\tilde{w} = w \odot \tilde{x}^{(r)}$ and $u(x) = -\tilde{w}^{\top} \log(x)$

4: &nbsp;&nbsp;&nbsp;&nbsp; $x^{(r+1)} = \text{Proj}_C^u(\tilde{x}^{(r)})$ ▷ Bregman projection

5: **end for**

6: **return** $x^{(r+1)}$

**Proposition 8.** Let $x^*$ be any minimizer of ($\mathcal{P}_C$). For any $N \in \mathbb{N}^*$ and any initial point $x^0 \in \mathbb{R}_{>0}^n$, the iterates of the constrained EMML algorithm 1 satisfy

$$\text{KL}(y, Ax^N) - \text{KL}(y, Ax^*) \le \frac{D_A(\log(x^*), \log(x^0))}{N}.$$

**Proof.** Let $\theta^{(r)} := \log x^{(r)}$, for $r \in \mathbb{N}$, with $(x^{(r)})_{r \in \mathbb{N}}$ the iterates of the constrained EMML algorithm 1. Define

$$\theta^* \in \arg\min_{\theta} \text{KL}(y, Ae^{\theta}) + \iota_{C_{\theta}}(\theta), \quad x^* = e^{\theta^*}.$$

By Lemma 2, the corresponding mirror-descent updates satisfy the three-point inequality

$$\text{KL}(y, Ax^{(r+1)}) - \text{KL}(y, Ax^*) \le D_A(\theta^*, \theta^{(r)}) - D_A(\theta^*, \theta^{(r+1)}).$$

Summing this inequality for $r = 0, \dots, N-1$ yields

$$\begin{aligned} & \min_{1 \le r \le N} (\text{KL}(y, Ax^{(r)}) - \text{KL}(y, Ax^*)) \\ & \le \frac{D_A(\log(x^*), \log(x^0)) - D_A(\log(x^*), \log(x^N))}{N}. \end{aligned} \quad (7)$$

Since Bregman divergences are nonnegative and $(\text{KL}(y, Ax^{(r)}))_{r \in \mathbb{N}}$ is nonincreasing [22, Theorem 1], the inequality (7) immediately yields the desired bound, proving the claim. □

# 5 Numerical Experiments

We consider hyperspectral unmixing [30] to compare the classical EMML algorithm with iterations given by (2) and its constrained version in Alg.1. Hyperspectral images record hundreds of contiguous spectral bands for each pixel of a scene, and due to the limited spatial resolution of the sensors, each pixel typically contains a mixture of spectral responses from several materials. Hyperspectral unmixing aims at decomposing each observed pixel spectrum into a set of pure spectral signatures, called endmembers, and estimating their associated abundance fractions [31]. The abundance vector $x = (x^q)^{\top}$, $q \in \{\text{Soil, Tree, Water}\}$, represents the spatial abundance maps. Estimating $x$ from $y$ and $A$ constitutes the hyperspectral unmixing problem [32]. We compare

Classical EMML: $\min_{x \ge 0} \text{KL}(y, Ax)$,

Constrained EMML: $\min_{x \ge 0} \text{KL}(y, Ax)$ s.t. $\sum_j x_j = 1$,

where the simplex constraint enforces that the abundances for each pixel sum to one [33–35]. Reconstruction quality is assessed using the Peak Signal-to-Noise Ratio (PSNR). Algorithms are run until $\|x^{(k+1)} - x^{(k)}\|_2^2 / \|x^{(k)}\|_2^2 < 10^{-5}$ or a maximum of $N = 1000$ iterations is reached.

We first empirically assess the convergence rates derived in Proposition 5. To this end, we approximate the optimal value $\text{KL}(y, Ax^*)$ by running EMML until convergence. At each iteration, we evaluate both sides of the inequality stated in Proposition 5 and verify that it is satisfied throughout the procedure. The results are reported in Figure 1. We perform the same analysis for the constrained EMML algorithm using the bound provided in Proposition 8. The results confirms that the theoretical inequality holds in this setting and that convergence to the optimum requires fewer. Consistent with the observations reported in [23], both EMML variants exhibit convergence rates that are significantly faster than the sublinear theoretical rate predicted by the bounds.

<figure><img src="image_3.png" /><figcaption>Figure 1: Comparison between the observed convergence behavior and the corresponding theoretical bounds for EMML and constrained EMML with $\sigma = 40$. In both cases, the empirical convergence is substantially faster than the predicted sublinear rate.</figcaption></figure>

We then assess the quality of the reconstructions. For all noise levels, the reconstructed abundance maps produced
by both algorithms are visually similar, as shown in Figure 3. Nevertheless, the reconstructions obtained with the unconstrained EMML appear noticeably noisier than those produced by the constrained EMML. This qualitative observation is supported by a quantitative evaluation: the constrained EMML consistently achieves more accurate reconstructions in fewer iterations, yielding higher PSNR values across all noise levels, as reported in Figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2: PSNR of the reconstructed abundance maps obtained with EMML and constrained EMML at different noise levels. The constrained algorithm consistently achieves higher PSNR values and outperforms EMML even at higher noise levels, while converging in nearly twice fewer iterations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Reconstructed abundance map for Soil with Poisson noise $\sigma = 40$. The reconstructions obtained with the EMML algorithm are visibly noisier than those obtained with the Constrained EMML algorithm.</figcaption></figure>

## 6 Conclusion and Perspectives

We revisited EMML as a MD method on a reparametrized objective, which allowed the introduction of convex constraints while preserving convergence guarantees. Numerical experiments show that these constraints accelerate convergence and improve reconstruction quality by acting as a regularization.

The method relies on Bregman projections, which can be challenging beyond simplex constraints [20, 36]. Extending it to more general convex or linear constraints would further broaden its applicability in inverse problems.

## Acknowledgment

Antonin Clerc acknowledges funding from the Berlin Mathematical School and the French Research Agency through the Holibrain project (ANR-23-CE45-0020). Ségolène Martin’s work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – The Berlin Mathematics Research Center MATH+ (EXC-2046/1, project ID: 390685689).

## References

[1] P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” *IEEE Signal Proces. Magazine*, vol. 23, no. 3, pp. 32–45, 2006.

[2] M. Bertero, P. Boccacci, G. Desidera, and G. Vicidomini, “Image deblurring with Poisson data: from cells to galaxies,” *Inverse Probl.*, 2009.

[3] F. Murtagh, *Astronomical image and data analysis*. Springer, 2006.

[4] F.-X. Dupé, J. Fadili, and J. L. Starck, “A proximal iteration for deconvolving Poisson noisy images using sparse representations,” *IEEE Trans. Image Process.*, vol. 18, no. 2, pp. 310–321, 2009.

[5] C. Zou and Y. Xia, “Restoration of hyperspectral image contaminated by Poisson noise using spectral unmixing,” *Neurocomputing*, vol. 275, pp. 430–437, 2018.

[6] L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” *IEEE Trans. Med. Imaging*, vol. 1, no. 2, pp. 113–122, 1982.

[7] L. B. Lucy, “An iterative technique for the rectification of observed distributions,” *The Astronomical Journal*, vol. 79, p. 745, 1974.

[8] W. H. Richardson, “Bayesian-based iterative method of image restoration,” *J. Opt. Soc. Am.*, vol. 62, no. 1, p. 55, 1972.

[9] C. L. Byrne, “Iterative image reconstruction algorithms based on cross-entropy minimization,” *IEEE Trans. Image Process.*, vol. 2, no. 1, pp. 96–103, 1993.

[10] A. Zunino, M. Castello, and G. Vicidomini, “Reconstructing the image scanning microscopy dataset: an inverse problem,” *Inverse Probl.*, vol. 39, no. 6, p. 064004, 2023.

[11] H. Zhang, J. Xu, P. Li, Z. Qin, D. Si, Y. Wang, Y. Wang, Q. Li, and Z. Xiao, “Probing the three-dimension emission source and neutron skin via $\pi - \pi$ correlations in heavy-ion collisions,” *arXiv preprint arXiv:2510.20554*, 2025.

[12] S. Setzer, G. Steidl, and T. Teuber, “Deblurring Poissonian images by split Bregman techniques,” *J. Visual Commun. Image Represent.*, vol. 21, no. 3, pp. 193–199, 2010.

[13] R. Kasai and H. Otsuka, “Iterative reconstruction with dynamic ElasticNet regularization for nuclear medicine imaging,” *J. Imaging*, vol. 11, p. 213, 2025.

[14] D. Deidda, N. A. Karakatsanis, P. M. Robson, Y.-J. Tsai, N. Efthimiou, K. Thielemans, et al., “Hybrid PET-MR list-mode kernelized expectation maximization reconstruction,” *Inverse Probl.*, vol. 35, no. 4, p. 044001, 2019.

[15] H. Marquis, K. P. Willowson, C. R. Schmidtlein, and D. L. Bailey, “Investigation and optimization of PET-guided SPECT reconstructions for improved radionuclide therapy dosimetry estimates,” *Front. Nucl. Med.*, vol. 3, 2023.

[16] T. Modrzyk, A. Etxebeste, E. Bretin, and V. Maxim, “Fast deconvolution using a combination of Richardson-Lucy iterations and diffusion regularisation,” in *Proc. 32nd Eur. Signal Process. Conf. (EUSIPCO)*, pp. 1901–1905, IEEE, 2024.

[17] C. Daniele, S. Villa, S. Vaiter, and L. Calatroni, “Deep equilibrium models for poisson imaging inverse problems via mirror descent,” *arXiv preprint arXiv:2507.11461*, 2025.
[18] M. V. Afonso, J. M. Bioucas-Dias, and M. A. Figueiredo, "An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems," *IEEE Trans. Imag. Proc.*, vol. 20, no. 3, pp. 681-695, 2010.

[19] W. Ha, E. Y. Sidky, R. F. Barber, T. G. Schmidt, and X. Pan, "Estimating the spectrum in computed tomography via Kullback-Leibler divergence constrained optimization," *Med. Phys.*, vol. 46, no. 1, pp. 81-92, 2019. arXiv:1805.00162 [physics].

[20] F. Beier, R. Beinert, and G. Steidl, "Multi-marginal Gro-mov-Wasserstein transport and barycentres," *Inf. Inference J. IMA*, vol. 12, pp. 2753-2781, 2023.

[21] R. T. Rockafellar, *Convex Analysis*. Princeton Univ. Press, 1997.

[22] H. H. Bauschke, J. Bolte, and M. Teboulle, "A descent lemma beyond Lipschitz gradient continuity: first-order methods revisited and applications," *Math. Oper. Res.*, vol. 42, no. 2, pp. 330-348, 2017.

[23] F. Kunstner, R. Kumar, and M. Schmidt, "Homeomorphic-invariance of EM: Non-asymptotic convergence in KL divergence for exponential families via mirror descent," in *Int. Conf. Artif. Intell. Stat. (AISTATS)*, pp. 3295-3303, PMLR, 2021.

[24] J. Bolte, S. Sabach, M. Teboulle, and Y. Vaisbourd, "First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems," *SIAM J. Optim.*, vol. 28, no. 3, pp. 2131-2151, 2018.

[25] H. Lu, R. M. Freund, and Y. Nesterov, "Relatively smooth convex optimization by first-order methods, and applications," *SIAM J. Optim.*, vol. 28, no. 1, pp. 333-354, 2018.

[26] P.-C. Aubin-Frankowski, A. Korba, and F. Léger, "Mirror descent with relative smoothness in measure spaces, with application to sinkhorn and em," in *Adv. Neural Inf. Process. Syst. 36th (NeurIPS 2022)*, pp. 17263-17275, 2022. Article No. 1255, Published 28 November 2022.

[27] C. Rossignol, F. Sureau, E. Chouzenoux, C. Comtat, and J.-C. Pesquet, "A Bregman majorization-minimization framework for PET image reconstruction," in *Proc. IEEE Int. Conf. Image Process. (ICIP)*, pp. 1736-1740, IEEE, 2022.

[28] L. Bregman, "The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming," *USSR Comput. Math. Math. Phys.*, vol. 7, no. 3, pp. 200-217, 1967.

[29] Y. Censor, *Parallel Optimization: Theory, Algorithms, and Applications*. New York, NY: Oxford Academic, online edn ed., 1998.

[30] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussoit, "Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches," *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*, vol. 5, no. 2, pp. 354-379, 2012.

[31] B. Zhang, "Advancement of hyperspectral image processing and information extraction," *J. of Remote Sens.*, vol. 20, no. 5, pp. 1062-1090, 2016.

[32] R. H. Chan, K. K. Kan, M. Nikolova, and R. J. Plemmons, "A two-stage method for spectral-spatial classification of hyperspectral images," *J. Math. Imaging Vis.*, vol. 62, pp. 790-807, 2020.

[33] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, "Total variation spatial regularization for sparse hyperspectral unmixing," *IEEE Trans. Geosci. Remote Sens.*, vol. 50, no. 11, pp. 4484-4502, 2012.

[34] E. Chouzenoux, M. Legendre, S. Moussaoui, and J. Idier, "Fast constrained least squares spectral unmixing using primal-dual interior-point optimization," *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*, vol. 7, no. 1, pp. 59-69, 2014.

[35] E. Chouzenoux, M.-C. Corbineau, and J.-C. Pesquet, "A proximal interior point algorithm with applications to image processing," *J. Math. Imaging Vision*, vol. 62, no. 6-7, pp. 919-940, 2020.

[36] V. Martin-Marquez, S. Reich, and S. Sabach, "Iterative methods for approximating fixed points of Bregman nonexpansive operators," *Discrete Contin. Dyn. Syst*, vol. 6, no. 4, pp. 1043-1063, 2013.