text:
\begin{document}

\newcounter{assum}

\begin{abstract}
We pass to continuum in optimisation problems associated to linear inverse problems $y = Ax$ with non-negativity constraint $x \geq 0$. We focus on the case where the noise model leads to maximum likelihood estimation through the so-called $\beta$-divergences, which cover several of the most common noise statistics such as Gaussian, Poisson and multiplicative Gamma. Considering~$x$ as a Radon measure over the domain on which the reconstruction is taking place,  we show a general sparsity result. In the high noise regime corresponding to $y \notin \setc{Ax}{x \geq 0}$, optimisers are typically sparse in the form of sums of Dirac measures. We hence provide an explanation as to why any possible algorithm successfully solving the optimisation problem will lead to undesirably spiky-looking images when the image resolution gets finer, a phenomenon well documented in the literature. We illustrate these results with several numerical examples inspired by medical imaging.
\end{abstract}

\maketitle
\section{Introduction}
Linear inverse problems $Ax = y$ often come with natural constraints, one of the most common being nonnegativity of the unknown, {i.e.}, $x \geq 0$. This happens in various applications, in particular in the imaging sciences. In this setting, $x \in \R^r$, $y \in \R^\ndet$ and $A \in \R^{\ndet \times r}$, where $\ndet$ is the number of data points, $r$ is the number of voxels.
One important example is that of Positron Emission Tomography (PET) where the sought-for image is the activity, which must be nonnegative~\cite{Shepp1982}. Deconvolution problems often also incorporate such constraints~\cite{Henrot2012}.

Depending on the noise model, the corresponding (negative) log-likelihood problem typically writes
\[ \min_{x \geq 0} \quad D(y|Ax),\]
where $D$ is some kind of divergence functional. If the noise model is Gaussian, for instance, then $D$ is simply the Euclidean distance, whereas if the noise model is Poisson, $D$ is the Kullback--Leibler divergence.

Analysing the effect of increased resolution leads to considering $x$ as a function (henceforth denoted $\mu$) in some functional space $X$, and $A$ now stands for some linear mapping from $X$ onto $\R^\ndet$, leading to the optimisation problem
\begin{equation}
\label{opt_int}
 \min_{\mu \geq 0} \quad D(y|A\mu).
 \end{equation}

The non-negativity constraint has been proved to provoke sparsity in various contexts in optimisation and optimal control~\cite{Clason2019}. As a result, the right functional space $X$ to be considered appears to be that of Radon measures, with discrete measures considered as sparse. Examples of this phenomenon can be seen in optimal control~\cite{Clason2017, Loheac2017}, and the same goes for the optimisation problem~\eqref{opt_int} when $D$ is the Kullback--Leibler divergence~\cite{Verdier2020}.

In the examples above, sparsity (which can arise in the form of Dirac masses) is undesirable as the sought-for image is expected to be at least piecewise smooth. In other contexts, sparsity of the signal must be enforced, as is the case for instance in sparse super-resolution~\cite{Denoyelle2017}. In the latter case, the situation is completely different as the unknown signal is known to be a sum of (nonnegative Dirac masses) and the aim is to recover their support.

The goal of the present work is to generalise the sparsity results of~\cite{Verdier2020} to the optimisation problem~\eqref{opt_int} when $D$ is \emph{any} $\beta$-divergence for $\beta \in [0,2]$. The $\beta$-divergences have attracted interest recently in non-negative matrix factorisation~\cite{Fevotte2011}, and are  now also advocated for in some medical imaging contexts, such as in PET~\cite{Cavalcanti2019}. They indeed have the appealing property of interpolating between three common divergences, namely the Itakura--Saito divergence ($\beta = 0$), the Kullback--Leibler divergence ($\beta = 1$) and the Euclidean distance ($\beta = 2$)~\cite{Cichocki2010}, which correspond to different noise models.

Some intermediate values of $\beta$ have shown to yield interesting results in instances where real data may have complex noise statistics, as is the case in PET~\cite{Cavalcanti2019}. The variable $\beta$ then adds flexibility to treating the inverse problem.

In this article, our main result may informally be stated as follows (see \autoref{sparse} for a precise statement), with mild assumptions on the operator $A$.
\begin{theorem}
\label{thm1}
If $y \notin \{A \mu, \; \mu \geq 0\}$, then any optimal solution $\mu^\star$ to~\eqref{opt_int} posed with a $\beta$-divergence is sparse.
\end{theorem}
In other words, if the data is not in the image of $\{\mu \geq 0\}$ under the operator $A$, any optimiser $\mu^\star$ will be sparse, typically a sum of Dirac masses. The condition $y \notin \{A \mu, \; \mu \geq 0\}$ should be interpreted as a condition on the level of noise: the more noise there is, the more likely it is that this condition be fulfilled. Sparsity will hence arise in the high-noise regime.

Our results show that the optimisation problem itself leads to sparse results. Thus, any algorithm successfully solving~\eqref{opt_int} will inevitably lead to undesirably spiky-looking images as one keeps iterating. In the context of medical imaging, this has been observed when using the Maximum-Likelihood-Expectation-Maximisation (ML-EM, also called the Richardson--Lucy algorithm) for solving~\eqref{opt_int}, and has been referred to in the literature as the ``night-sky'' or the ``draughtsboard'' effect~\cite{Vardi1985}.

Consequently, the same kind of artefacts will be observed for other likelihoods, hinting at the necessity of either early stopping when solving~\eqref{opt_int} (see \cite{Resmerita2007}), or adding appropriate regularisation terms in the form
\begin{equation}
\label{reg}
 \min_{\mu \geq 0} \quad D(y|A\mu) + \lambda R(\mu).
 \end{equation}
with $\lambda$ big enough in order to alleviate the issue~\cite{Cavalcanti2019}.

In this work, we also provide examples where, solving~\eqref{opt_int} for some values of $\beta \in [0,2]$, reconstructions exhibit the night-sky effet for inverse problems with sufficient noise. In agreement with our theoretical results, this is the case with enough iterations of either a convergent algorithm for solving~\eqref{opt_int}, or a convergent algorithm for solving~\eqref{reg} with $\lambda$ small.

\subsection*{Outline of the paper}
The paper is organised as follows.
In \autoref{sec2}, we define the inverse problem by setting the functional analytic framework as well as the noise model through $\beta$-divergences, leading to the corresponding maximum likelihood problem.
Then, in \autoref{sec3}, we analyse the resulting optimisation problem and prove~\autoref{thm1}.
Numerical simulations are presented in \autoref{sec4} in different contexts taken from medical imaging, confirming our results about sparsity.

\section{Setting the inverse problem}
 \label{sec2}

 \subsection{Linear inverse problem}
 \label{sec:lininv}
 We aim at reconstructing an image $\mu$ defined on a compact $K \subset \R^p$, $p \geq 1$. The sought-for image is taken to be an element of the space of Radon measures, denoted $\MM(K)$, which is the topological dual space of continuous fonctions over the compact, denoted $\mathcal{C}(K)$. We endow $\MM(K)$ with the weak-$\ast$ topology, making $\mathcal{C}(K)$ its dual space. The dual pairing between a function $\mu \in \MM(K)$ and a function $f \in \mathcal{C}(K)$ will be denoted
 $\scl{\mu}{f}$. Finally, $\MM_+(K)$ stands for the set of nonnegative Radon measures.

 The data is made of $\ndet$ measurements, stacked into a nonnegative vector $y \in \R_+^\ndet$.
 We shall use the notation
 \[I:= \setc{i \in \{1, \ldots, m\}}{ y_i>0}.\]
 This vector itself is the realisation of some random variable with mean $A \mu$, where $\mu \in \MM_+(K)$ is the image to be reconstructed, and $A$ is a linear mapping $A \colon \MM(K) \to \R^\ndet$.

We make the only regularity assumption that $A$ is continuous in the weak-$\ast$ topology.
From~\cite[Proposition 3.14]{Brezis2010}, this forces $A$ to be of the form
 \begin{equation}
\label{scalar} (A \mu)_i = \scl{\mu}{a_i}, \quad i = 1, \ldots, \ndet,
\end{equation} where the $a_i$ are elements of $\mathcal{C}(K)$. This covers the case of PET~\cite{Mair1996, Verdier2020} and more generally the setting of Hilbert--Schmidt operators: if the underlying operator in infinite dimension is of the form
\[ \mu \longmapsto  \int_K k(\cdot,y) \, d\mu(y),\] for some smooth kernel $k \in \mathcal{C}(K \times K)$, the operator $A$ typically is obtained from a sampling for $\ndet$ points $x_i \in K$ or integrating the kernel over some subdomains $\Omega_i \subset K$, namely
\[a_i = k(x_i, \cdot),  \text{ or } a_i = \int_{\Omega_i} k(x, \cdot) \, dx, \quad i = 1, \ldots, \ndet.\]

We will make the assumption that $A$ maps $\MM_+(K)$ into the set of (componentwise) nonnegative vectors denoted $\R_+^\ndet$, and that it is non-trivial, {i.e.}, $A \mu > 0$ componentwise for $\mu$ the uniform measure over $K$. These two assumptions may be written in the terms of the functions $a_i$ as
\begin{equation}
\label{positive}
a_i \geq 0, \; a_i \neq 0, \qquad  i = 1, \ldots, \ndet.
\end{equation}
 Another consequence of the simple continuity assumption on $A$ is that its adjoint $A^\ast \colon \R^\ndet \to \mathcal{C}(K)$ is simply defined as
 \[A^\ast \lambda =  \sum_{i=1}^\ndet \lambda_i a_i, \quad \lambda \in \R^\ndet.\]

 We shall sometimes need to know when $A^\ast$ is injective. This is of course equivalent to the linear independence of the family $(a_i)_{i = 1,\ldots, \ndet}$.

 Note that since the codomain of \(A\) is finite dimensional, we have
 \[A^\ast \text{ is injective}  \iff  A  \text{ is surjective}.\]

 In order to solve the inverse problem (written informally as $A \mu = y$), we aim at solving the optimisation problem
 \begin{equation}
\label{opt}
 \operatorname*{min}_{\mu \in \MM_+(K)}  \; D_\beta (y | A\mu),
 \end{equation}
 where $D_\beta$ is the $\beta$-divergence for $\beta  \in [0,2]$, see below for the definition.
 We will use the notation \[\loss(\mu) := D_\beta(y|A \mu).\]

 Finally, let us define the notion of support for the various relevant cases (all these cases can be covered in one single definition, but we prefer separating them for clarity):
 \begin{itemize}
 \item for a vector $w \in \R_+^\ndet$, \[\supp(w) = \setc*{i \in \{1, \ldots, \ndet\}}{w_i > 0}\]
 \item for a nonnegative function $f \in \mathcal{C}(K)$,  \[\supp(f) = \overline{\setc*{x \in K}{f(x) > 0}}\]
  \item for a nonnegative measure $\mu \in \MM_+(K)$,  \[\supp(\mu) := \setc[\big]{x \in K}{ \mu(N) > 0, \; \forall N \in N(x)}.
\]
where $N(x)$ is the set of all open neighbourhoods of $x$.

 \end{itemize}

  \subsection{$\beta$-divergences}

  \subsubsection{For scalar positive variables}
 For $u > 0$, $v>0$ scalar variables, $\beta \in [0,2]$, we define
\[ d_{\beta}(u | v) := \frac{1}{\beta(\beta - 1)}\left(u^\beta + (\beta - 1) v^\beta - \beta u v^{\beta - 1}\right),\]
which for $\beta = 2$ gives the Euclidean distance \[d_2(u|v) = \frac{1}{2}(u-v)^2,\]
and by continuity for $\beta = 1$ the Kullback--Leibler divergence
\[d_1(u|v) = u \log\Big(\frac{u}{v}\Big) - u + v,\]
and for $\beta = 0$ the Itakura--Saito divergence
\[d_0(u|v) = \frac{u}{v} - \log\Big(\frac{u}{v}\Big) - 1.\]

More precisely, using the convention $0/0 = 0$, $0 \log 0 = 0$, $d_\beta$ is defined for nonnegative scalars $u \geq 0$, $v\geq 0$ as follows.
\begin{description}
\item[Case $1< \beta\leq 2$] \[d_\beta(u|v) = \frac{1}{\beta\,(\beta-1)}\,\big(
u^{\beta} + (\beta-1)\,v^{\beta} - \beta\,u\,v^{\beta-1}
\big)\]

\item[Case $\beta = 1$]
  \[d_1(u|v) =

    \begin{cases}
+ \infty & \text{if \(v = 0,\, u > 0\)}
\\
v -u-u\,\log\big(\frac{v}{u}\big) & \text{otherwise}
\end{cases}
\]

\item[Case $0 < \beta < 1$]

  \[
    d_\beta(u|v) = \begin{cases}
+ \infty & \text{ if \(\; v = 0,\, u > 0\)}
\\
\frac{1}{\beta\,(\beta-1)}\,\big(
u^{\beta} + (\beta-1)\,v^{\beta} - \beta\,u\,v^{\beta-1}
\big) & \text{otherwise}
\end{cases}
\]

\item[Case $\beta = 0$]

  \[
    d_0(u|v) = \begin{cases}
+ \infty & \text{ if \(u = 0 \text{ or } v = 0\)}
\\
\frac{u}{v}- \log\big(\frac{u}{v}\big) - 1 & \text{otherwise}
\end{cases}
\]
\end{description}

 Let us review some of the most important properties of the $\beta$-divergences, which will be used throughout.
 \begin{description}
\item[Separation] for all $\beta \in [0,2]$,  \[d_\beta(u|v) = 0 \iff  u = v.\]
\item[Convexity] if (and only if) $\beta \in [1,2]$ and for any $u \in \R_+$, $v \mapsto d_\beta(u|v)$ is convex on $\R_+$.
\item[Lower semicontinuity] for all $\beta \in [0,2]$ and $u \in \R_+$, the mapping $v \mapsto d_\beta(u|v)$ is lower-semicontinuous on $\R_+$.
\end{description}
 Finally, it will be convenient to analyse when, for a given $u \in \R_+$, $v \mapsto d_\beta(u|v)$ is differentiable.
More precisely, we will insist on situations when there is a point~$v$, then necessarily $0$ in view of the definition, at which $v \mapsto d_\beta(u|v)$ is non-differentiable. This happens in the two following regimes:
\begin{itemize}
\item if $\beta \in (1,2)$, and $u>0$, the function $v \mapsto d_\beta(u|v)$ is not differentiable at $v=0$.
\item if $\beta \in (0,1)$, the function $v \mapsto d_\beta(u|v)$ is not differentiable at $v=0$.
\end{itemize}
Note that the function is not subdifferentiable either at such points.

 \subsubsection{For nonnegative vectors} The $\beta$-divergence for vectors $y, w \in \R_+^\ndet$ is then defined as
\[ D_\beta(y|w) := \sum_{i=1}^\ndet d_\beta(y_i | w_i),\]
and inherits the smoothness properties of $d_\beta$.

   \subsection{The noise model.}

Let us review the underlying statistical model, of which the minimisation problem~\eqref{opt} is (up to constants) the (negative-log) likelihood maximum problem, denoting $w = A \mu$, and $\phi$ a dispersion parameter. We write the statistical model for $y$ and $w$ as scalar variables, the full statistical model is defined component by component by independent draws of the same form.

A general way to write the noise model giving rise to $\beta$-divergences is to use the so-called Tweedie distributions. The $\beta$-divergences are a special case of such distributions, as the corresponding Tweedie distribution is given by
\[ y \mapsto h (x,\phi) \operatorname{exp}\left(-\frac{1}{\phi} d_\beta(y|w)\right),\]
so that minimising the negative log-likelihood problem indeed is equivalent to minimising $w \mapsto d_\beta(y|w)$. We refer to~\cite{Simsekli2013} for more details.

The underlying density is not always tractable (this is the case if $1< \beta < 2$), making the noise model unclear. In other cases, the noise model can be further identified as follows.

\begin{description}
\item[Case $\beta = 2$] the noise model is Gaussian:  \[y = \mathcal{N}(w, \phi).\]
\item[Case $1< \beta<  2$] no explicit model is known.
\item[Case $\beta = 1$] the noise model is Poisson:  \[y =\phi \, \mathcal{P}\Big(\frac{1}{\phi} w\Big).\]
\item[Case $0 < \beta < 1$] the noise model is compound Poisson, whose distribution is the sum of a singular measure at $0$ and an absolutely continuous measure on the positive reals. More precisely:  \[y = \sum_{i=1}^n g_i,\]
where $n = \mathcal{P}(\lambda)$, $g_i = \mathcal{G}(a,b)$ (here, $\mathcal{G}$ stands for the Gamma distribution, of density $x \mapsto \frac{b^a x^{a-1}e^{-b x}}{\Gamma(a)}$ over $(0, +\infty)$), and the parameters $\lambda$, $a$ and $b$ are given as functions of $\phi$ and $\beta$ through
\[\lambda = \frac{1}{\phi}\frac{w^\beta}{\beta}, \quad a = \frac{\beta}{1-\beta}, \quad b = \frac{1}{\phi}\frac{w^{\beta-1}}{1-\beta}.\]
\item[Case $\beta = 0$] the noise model is multiplicative Gamma:  \[y = w\, \mathcal{G}(a,b),\]
with \[a = \frac{1}{\phi}, \quad b = \frac{1}{\phi}.\]
\end{description}

A common property of these noise models is that,
the random variable $y$ tends to $w$
as the dispersion parameter $\phi$ goes to zero:
 \beq
 \label{estimator}
 y \to w \;\;  \text{a.s.}, \quad \text{ as } \, \phi \to 0.\eeq

 \section{Optimisation problem}
 \label{sec3}
 We now investigate the optimisation problem~\eqref{opt}, starting with a related optimisation problem and its dual.
 \subsection{Related optimisation problem and its dual}
 Let us define the cone
 \[
   C:= \imcone = \setc*{A \mu}{\mu \in \MM_+}
   ,
 \]
 which is closed and convex.
 The original optimisation problem~\eqref{opt} is related to the following optimisation problem:
 \begin{equation}
\label{cone_formulation}
\operatorname*{min}_{w \in \imcone}  \; D_\beta (y | w).
\end{equation}
More precisely, for any $w^\star$ optimal for the above problem, any measure $\mu^\star \in \MM_+(K)$  such that $A \mu^\star = w^\star$ is optimal for the original problem.

\begin{lemma}
\label{unique_min}
The minimum in~\eqref{opt} is attained.
\end{lemma}
\begin{proof}
To prove the claim, we prove that the optimal value in the minimisation problem~\eqref{cone_formulation} is attained. For $\beta >0$, the function $w \mapsto D_\beta(y|w)$ behaves at infinity like $w \mapsto \frac{1}{\beta} \norm{w}^\beta_\beta$ and hence is coercive. For $\beta = 0$, it is also coercive since it behaves live $w \mapsto - \sum_{i=1}^\ndet \log(w_i)$ at infinity. It is also lower semi-continuous, and since $\imcone$ is closed, there exists an optimal $w^\star \in \imcone$ for the problem~\eqref{cone_formulation}. Any $\mu^\star$ such that $A\mu^\star =w^\star$ then provides a minimiser for the original problem~\eqref{opt}.
\end{proof}

We now compute the dual problem to the original one. To be more precise, we compute the dual to~\eqref{cone_formulation}.
Finding the dual will be instrumental in proving our main result, \autoref{sparse}, via~\autoref{diff}.
These computations will also happen to be crucial in analysing results from numerical simulations, determining whether we should expect sparsity or not for some given $y \in \R_+^\ndet$.

We define the cone $\imcone^\ast := \setc*{\lambda \in \R^\ndet}{\scl{\lambda}{w} \geq 0, \;\, \forall w \in \imcone}$ dual to $\imcone$, which can be characterised as in~\cite{Georgiou2005} by \[\imcone^\ast = \setc*{\lambda \in \mathbb{R}^\ndet}{A^* \lambda \geq 0 \text{ on } K}.\]
The dual problem writes
\begin{equation}
\label{dual_problem}
\operatorname*{max}_{\lambda \in \imcone^\ast}  \; g(\lambda)
,
\end{equation}
where the function \(g \colon \R^\ndet \to \R\) is defined for $\lambda \in \imcone^\ast$ by
\begin{equation}
  \label{eq:defg}
  g(\lambda) := \min_{ w \in \R^\ndet} D_\beta(y|w)- \langle\lambda , w \rangle
  .
\end{equation}
The point of taking the dual in this form is that, since $D_\beta$ decomposes, so does $g$ and we find
\[ g(\lambda) = \sum_{i=1}^m \min_{w_i \in \R} \paren[\big]{d_\beta(y_i|w_i)- \lambda_i w_i}. \]
As a result, we focus on computing
\begin{equation}
  \label{eq:defh}
  h(y, \lambda) := \min_{w \in \R} \left(d_\beta(y|w)- \lambda w \right) =\min_{w \geq 0} \left(d_\beta(y|w)- \lambda w \right)
\end{equation}
for a scalar $y \geq 0$.
The resulting dual function will take the form
\(g(\lambda) = \sum_{i=1}^m h(y_i, \lambda_i).\)

As the explicit computation of the functions $h$ in \eqref{eq:defh} is quite involved, we postpone it to \autoref{appA}.
\begin{lemma}
\label{diff}
Let $w^\star$ be an optimal vector for the optimisation problem~\eqref{cone_formulation}.
Then we have \[ \operatorname{supp}(y) \subset  \operatorname{supp}(w^\star).\]
\end{lemma}

\begin{proof}
Let us argue depending on the value of $\beta$. For any $0 \leq \beta \leq 1$, the objective function would be infinite at any $w$ such that $w_i = 0$ with $y_i>0$, which obviously contradicts optimality.

When $1 < \beta \leq 2$, we shall make use of duality. The lemma is equivalent to proving that any $w^\star$ optimal for~\eqref{cone_formulation} will satisfy $w^\star_i > 0$ whenever~$i$ is such that $y_i>0$.
Since $1 < \beta \leq 2$, the problem at hand is convex. Slater's condition is satisfied since the cone $\imcone$ then has non-empty relative interior, as $A$ is non-trivial.
Thus, strong duality holds and the dual optimal value is attained.

Denoting $\lambda^\star$ a dual optimal variable, this implies that any primal optimal variable $w^\star$ must minimise $w \mapsto D_\beta(y|w) - \langle \lambda^\star, w\rangle$ over $\R^\ndet$.\footnote{This stems from the usual computations leading to the KKT conditions, see~\cite{Boyd2004}:
by strong duality, letting $\lambda^\star$ be dual optimal, we have for any $w^\star$ primal optimal,
$D_\beta(y|w^\star) = g(\lambda^\star) = \min_{w \in \R^\ndet} (D_\beta(y|w) - \langle \lambda^\star, w\rangle) \leq D_\beta(y|w^\star) - \langle \lambda^\star, w^\star\rangle \leq D_\beta(y|w^\star)$. All the inequalities must be equalities and in particular, $w^\star$ must minimise $w \mapsto D_\beta(y|w) - \langle \lambda^\star, w\rangle$ over $\R^\ndet$.}
 Thus, for each $i \in \{1, \ldots, m\}$,  $w_i^\star$ minimises $w \mapsto d_\beta(y_i|w) - \lambda_i^\star w$ over $\R$. Computations made to evaluate the dual problem impose that when $y_i >0$, we must have
\[(w^\star_i)^{\beta- 2}(w^\star_i - y_i) = \lambda_i^\star\] for such an $i$. Consequently, $w^\star_i > 0$ and the proof is finished.
\end{proof}

We also gather some results about optimal solutions to~\eqref{cone_formulation} in \autoref{appB}.

 \subsection{KKT conditions}
Let us compute the KKT conditions for problem~\eqref{opt}. We will do so at points at which $\loss$ is differentiable (endowing $\MM_+(K)$ with its strong topology), namely at measures $\mu \in \MM_+(K)$ such that
\begin{itemize}
\item when $\beta \in [1,2)$, $\operatorname{supp}(y) \subset  \operatorname{supp}(A\mu)$,
\item when $\beta \in [0,1)$, we have $A \mu> 0$.
\end{itemize}
We shall denote $S$ the set of such measures $\mu$.

For $\mu \in S$, we define
\begin{equation}
  \label{eq:deflambda}
  \lambda(\mu) := (A\bar \mu)^{\beta - 2} (A\bar \mu- y),
\end{equation}
where multiplications and divisions are to be understood componentwise. This formula is well-defined when $\mu \in S$, recalling the convention $0/0 = 0$.

\begin{proposition}
Let $\bar \mu \in S$.
Then the KKT conditions write
\beq
\label{kkt}
A^\ast \lambda(\bar \mu) \geq 0 \text { on } K, \qquad A^\ast \lambda(\bar \mu) = 0 \text{ on } \supp(\bar \mu).
\eeq
\end{proposition}
\begin{proof}
The KKT conditions here write
\[ \nabla \loss(\bar \mu) \in -N_{\MM_+}(\bar \mu),\]
where $N_{\MM_+(K)} (\mu)$  is the normal cone of $\MM_+(K)$ at $\bar \mu$, defined by \[N_{\MM_+(K)} (\mu):= \setc*{f \in \mathcal{C}(K)}{\forall \nu \in \MM_+(K), \; \scl{\mu-\nu}{f} \geq 0}.\]
On the one hand, the normal cone can be identified as \[ N_{\MM_+(K)}(\mu) = \setc*{ f \in \mathcal{C}(K)} {f\leq 0 \text{ on } K, \; f = 0 \text{ on } \supp(\mu)},\]
see~\cite{Verdier2020}, and on the other hand the gradient of $\loss$ at points of differentiability is readily computed as
\[\nabla l(\mu) = A^\ast \lambda(\mu).\]
The combination of these two results exactly leads to~\eqref{kkt}.
\end{proof}

Since $\loss$ is smooth at points of regularity and the constraint $\{\mu \geq 0 \}$ is convex, we have the following dichotomy for a point $\bar \mu$ satisfying the KKT conditions
\begin{itemize}
\item for any $\beta \in [0,2]$, \[\bar \mu\text{ is optimal }  \implies  \bar \mu \text{ satisfies }\eqref{kkt}.\]
\item for any $\beta \in [1,2]$, \[\bar \mu \text{ satisfies }~\eqref{kkt}\implies  \bar \mu\text{ is optimal}.  \]
\end{itemize}

 \subsection{Sparsity theorem}

 We are now in a position to prove the sparsity theorem.
 We will say that a measure $\mu$ is \emph{sparse} (with respect to an operator \(A \colon \MM \to \RR^{\ndet}\) as described in \autoref{sec:lininv}) if there exists $\varphi \in A^*{\RR^{\ndet}}$, $\varphi \geq 0$, $\varphi \neq 0$ such that
\[ \supp(\mu) \subset  \argmin(\varphi).\]

\begin{theorem}
\label{sparse}
For any $\mu^\star$ optimal for~\eqref{opt},

\beq
\operatorname{supp}(\mu^\star) \subset \argmin(\varphi^\star), \quad \varphi^\star = A^\ast \lambda(\mu^\star).
\eeq
In particular, if $A$ is surjective and $y \notin \imcone$, any optimal measure $\mu^\star$ is sparse.
\end{theorem}
\begin{proof}
Let $\mu^\star$ be optimal for~\eqref{opt}. \par
\begin{description}
\item[{First case, $\beta \in \{0\} \cup [1,2]$}]
we know from \autoref{diff} that the function $\loss$ is differentiable at $\mu^\star$. We may thus write the KKT conditions which are necessary (and sufficient for $\beta \in [1,2]$). Hence, we find that the support of $\mu^\star$ satisfies
\begin{equation*}
\operatorname{supp}(\mu^\star) \subset \argmin(\varphi^\star), \quad \varphi^\star = A^\ast \lambda(\mu^\star).
\end{equation*}

\item[Second case, $\beta \in (0,1)$]
let us denote \[I_{\mu^\star}:= \setc*{i \in \{1,\ldots, \ndet\}}{(A \mu^\star)_i > 0} =  \supp(A\mu^\star).\]
First, we notice that \autoref{diff} ensures that for any $i \notin I_{\mu^\star}$, we have $y_i = 0$.
Thus we have \[D_\beta(y|A\mu^\star) = \sum_{i \in I_{\mu^\star}} d_\beta(y_i|(A \mu^\star)_i).\]
We define a reduced compact
\[\widetilde K = \bigcup_{i \in I_{\mu^\star}} \, \supp(a_i)
  ,
\]
and, correspondingly, a reduced operator  \(\widetilde A \colon \MM_+(\widetilde K) \to \RR^{\#I_{\mu^\star}}\) by
\[
  \widetilde A \mu \coloneqq (A \mu)_{i \in I_{\mu^\star}}
  .
\]

Since $(A \mu^\star)_i = 0$ for all $i \notin I_{\mu^\star}$, $\supp(\mu^\star) \subset \widetilde K$.

Let us prove that $\mu^\star$, regarded as an element of $\MM_+(\widetilde K)$ is optimal for the problem
\begin{equation}
\label{reduced} \operatorname*{\min}_{\mu \in \MM_+(\widetilde K)} \tilde \loss (\mu), \qquad \tilde \loss(\mu):=  D_\beta\big(y |\widetilde A \mu\big).
\end{equation}
If not, and arguing by contradiction, we may find $\bar \mu \in \MM_+(\widetilde K)$ such that $\tilde \loss(\bar \mu) < \tilde \loss(\mu^\star)$. Let us extend $\bar \mu$ by $0$ outside of $\widetilde K$, still abusively denoting the corresponding measure of $\MM_+(K)$ as~$\bar \mu$. Given that $(A \bar \mu)_i = 0$, and $y_i = 0$ for all $i \notin I_{\mu^\star}$, we find
\[\loss(\mu^\star) = \tilde \loss(\mu^\star) > \tilde \loss(\bar \mu) = \loss(\bar \mu),\]
which contradicts the optimality of $\mu^\star$ for the original problem.

Considering that $\mu^\star$ is optimal for the optimisation problem~\eqref{reduced}, and that $\tilde \loss$ is differentiable at $\mu^\star$, we may this time use the necessity of KKT conditions for this problem, which leads to
\beq
\operatorname{supp}(\mu^\star) \subset \argmin\big(\widetilde A^\ast \lambda(\mu^\star)\big), \quad \lambda_i^\star = (A\mu_i^\star)^{\beta - 2} (A \mu_i^\star - y_i), \; i \in I_{\mu^\star},
\eeq
where $\mu^\star$ is regarded as a measure in $\MM_+(\widetilde K)$. With the convention $0/0 = 0$ and from the fact that the original measure $\mu^\star \in \MM_+(K)$ has support included in $\widetilde K$, the result above is equivalent to
\beq
\operatorname{supp}(\mu^\star) \subset  \argmin(A^\ast \lambda(\mu^\star)),
\eeq
as claimed.
\end{description}

Let us now assume that $A$ is surjective and $y \notin \imcone$ and we shall prove that sparsity holds. We already know that $\varphi^\star = A^\ast \lambda(\mu^\star) \geq 0$. It remains to show that $\varphi^\star \neq 0$. In order to do so, let us recall that the surjectivity of $A$ is equivalent to the injectivity of $A^\ast$ (itself equivalent to the linear independence of the $a_i$'s in $\mathcal{C}(K)$). Since $A^\ast$ is injective and $\varphi^\star = A^\ast \lambda(\mu^\star)$, all we need to prove is that $\lambda^\star \neq 0$. Given that $\lambda(\mu^\star) =  (A\mu^\star)^{\beta - 2} (A \mu^\star -y)$, this vector cannot vanish as otherwise we would get $A \mu^\star = y$, contradicting our initial assumption.
\end{proof}

For a discussion of why we may speak of a sparse measure when its support is included in $\argmin(\varphi)$ for some non trivial function $\varphi \in A^*\RR^{\ndet}$ satisfying $\varphi \geq 0$, we refer to~\cite{Verdier2020}.
For example, if $A$ is surjective and if the $a_i$ functions defining $A$ are analytic on $K$, any sparse measure in the sense defined here will be a sum of Dirac masses.

\begin{remark}
  \label{unique_fct}
  In the case \(\beta \in [1,2]\),
  since  $w \mapsto D_\beta (y | w)$ is strictly convex,
  we may say the following:
  \begin{itemize}
  \item
    if \(\mu_1\) and \(\mu_2\) are minimisers of \eqref{opt}, then  \(\lambda(\mu_1) = \lambda(\mu_2)\), where \(\lambda\) is given by~\eqref{eq:deflambda}.
    Let us denote this unique value by \(\lambda^{\star}\).
    In particular, this means that all the minimisers have a support included in the same set \(\argmin(A^*\lambda^{\star})\),
  \item
    the set $\setc*{A \mu^\star}{\mu^\star \text{ optimal for }\eqref{opt}}$ is reduced to a singleton,
  \item
    the optimisation problem~\eqref{cone_formulation} has a single minimum $w^\star$.
  \end{itemize}
  We do not know if these observations still hold

  for $\beta \in [0,1)$.
\end{remark}

 \subsection{Low-noise and high noise regimes}

For high values of the dispersion parameter $\phi$, {i.e.}, when the noise level is high, we expect that the random variable $y$ may be outside of the cone $\imcone$, leading to sparse measures as optimal solution to the optimisation problem~\eqref{opt}.
We may prove the surprising result that there will be a single optimal measure to~\eqref{opt} under quite general conditions, which require $y \notin \imcone$ and $\beta \in [1,2]$.

  Let us start with a formula for optimisers $\lambda(\mu^\star) =  (A\mu^\star)^{\beta - 2} (A \mu^\star - y)$ for an optimal measure $\mu^\star$.

\begin{lemma}
Let $\mu^\star$ be an optimal measure. Then
\beq
\label{form}
\sum_{i=1}^\ndet y_i (A \mu^\star)_i^{\beta-1} =  \sum_{i=1}^\ndet (A \mu^\star)_i^{\beta}.
\eeq
In particular, if $y \notin \imcone$, we find that $\lambda(\mu^\star)$ must have at least one positive and one negative component.
\end{lemma}
\begin{proof}
Proving the formula requires establishing that
\beq
\label{form}
\sum_{i=1}^\ndet y_i (w^\star)_i^{\beta-1} =  \sum_{i=1}^\ndet (w^\star)_i^{\beta},
\eeq
for any $w^\star$ optimal form for~\eqref{cone_formulation}, which is nothing but the result of~\autoref{relation_opt} in \autoref{appB}.

Now assume that $y \notin \imcone$ and by contradiction that all components of $\lambda^\star$ are nonnegative. Then $y_i \leq (A \mu^\star)_i$ for all $i \in \{1, \ldots, m\}$, and at least one of inequalities must be strict since $y \notin \imcone$,  leading to
\[\sum_{i =1}^{\ndet} (A \mu^\star)_i^{\beta} = \sum_{i =1}^{\ndet} y_i (A \mu^\star)_i^{\beta-1} < \sum_{i =1}^{\ndet} (A \mu^\star)_i^{\beta},\]
a contradiction with~\eqref{form}.

\begin{corollary}
\label{unique_opt}
Assume that $\beta \in [1,2]$  and $y \notin \imcone$.
Then, if the function $\varphi^\star \coloneqq A^*\lambda^{\star} \neq 0$, where \(\lambda^{\star}\) is defined in~\autoref{unique_fct}, has a single minimum $x^\star$, there is a single optimiser solving~\eqref{opt}, given by
\[\mu^\star = \xi \,\delta_{x^\star}, \qquad \xi := \frac{\sum_{i=1}^\ndet y_i a_i^{\beta-1}(x^\star)}{\sum_{i=1}^\ndet  a_i^{\beta}(x^\star)}, \]
with $\delta_{x^\star}$ the Dirac mass at $x^\star$.
\end{corollary}

\begin{proof}
Let $\bar \mu^\star$ be optimal. We know that $\supp(\mu^\star) \subset \argmin (\varphi^\star)$, with $\varphi^\star \neq 0$ since $y \notin \imcone$. The assumption imposes $\mu^\star = \xi \,\delta_{x^\star}$ for some $\xi> 0$. In virtue of the formula~\eqref{form}, we must have
\[\sum_{i=1}^\ndet y_i \xi^{\beta -1} (A \delta_{x^\star})_i^{\beta-1} =  \sum_{i=1}^\ndet \xi^\beta (A \delta_{x^\star})_i^{\beta},\]
equation with a unique solution for $\xi$, as given in the statement.

\end{proof}

On the other hand, in the low-noise regime and in view of~\eqref{estimator}, we expect $y \in \imcone$.
In this more favourable situation, a natural question is the regularity of optimisers, hoping that some measures solving of $A \mu = y$ will be absolutely continuous with respect to the Lebesgue measure on $K$ (in the sequel, \emph{absolutely continuous} will mean with respect to the Lebesgue measure).
This question is tackled in~\cite{Verdier2020}, based on the results of~\cite{Georgiou2005}, we end up with the following result where we use the notation
\[\widetilde K :=  \bigcup_{i \in I} \, \supp(a_i), \quad \widetilde A \colon \MM_+(\widetilde K) \to \RR^{\#I}, \;
  \widetilde A \mu \coloneqq (A \mu)_{i \in I}.\]

\begin{proposition}[{\cite[Proposition~3.13]{Verdier2020}}]
Assume that $y \in \imcone$, and further that
\beq
\label{smooth}
\tilde y := (y_i)_{i \in I} \in \mathrm{int}\Big(\widetilde \imcone\Big), \quad \widetilde \imcone:= \setc*{ \widetilde A \mu}{\mu \in \MM_+(\widetilde K)}.\eeq
Then there exists an absolutely continuous measure $\mu^\star$ solving~\eqref{opt}, {i.e.}, such that $A \mu^\star = y$.
\end{proposition}

Quantifying the probability to have $y \in \imcone$ is outside of the scope of this paper, we refer to~\cite{Verdier2020} for some results in that direction in the case of Poisson noise, namely for $\beta = 1$. Also outside the scope of this paper is to know whether a given algorithm will pick an absolutely continuous measure at the limit, rather than a sparse one, in the low-noise regime where $y \in \imcone$ with~\eqref{smooth} satisfied.

\end{proof}

 \section{Numerical experiments}
  \label{sec4}
We here present some simulations of algorithms solving~\eqref{opt} in different contexts, where the results exhibit sparsity as expected from the theoretical results. All simulations are run using Python.
 \subsection{Sparsity Certificates}
 As evidenced by our results, the relevant criterion for sparsity of optimisers for the optimisation problem~\eqref{opt} is independent of $\beta$, as the question reduces to:
 \[\text{do we have } y \in \imcone?\]
In practice, as one wants to solve~\eqref{opt} (or possibly a regularised version thereof) in the form of some iterative algorithm defined by iterates of the form
\[ \mu_{k+1} = G_k(\mu_k),\]
we are looking for methods allowing us to \emph{certify} that $y \notin \imcone$ along iterates, {i.e.}, a method to prove that $y \notin \imcone$ which writes as a function of $\mu_k$ and that gets better as $k \to +\infty$.

One approach towards this is to make use of duality: by weak duality (and recalling the definition of \(g\) in \eqref{eq:defg}) we always have
\[ \forall \mu \in \MM_+(K), \; \forall \lambda \in \imcone^\ast, \quad \loss(\mu) \geq g(\lambda).\]
Since $y \in \imcone \iff \loss(\mu) = 0$, this entails the following straightforward result:
\beq
\label{certificate}
(\exists \lambda \in \imcone^\ast, \; g(\lambda) > 0)  \quad \implies \quad y \notin \imcone.
\eeq
We will call a vector $\lambda \in \imcone^\ast$ such that $g(\lambda) > 0$ a \emph{dual certificate} of sparsity.

This provides a natural method when it comes to establishing sparsity: assume we have a convergent algorithm for solving~\eqref{opt}, in the sense that all subsequences of $(\mu_k)_{k \in \mathbb{N}}$ converge (in the weak-$\ast$ sense) to some minimiser $\mu^\star$ of~\eqref{opt}.

Then, a candidate of choice for a dual certificate is given by
\[\lambda_k := (A \mu_k)^{\beta - 2} (A \mu_k - y).\]
Indeed, at least for $\beta \in [1,2]$, strong duality holds and the (unique) optimal dual variable $\lambda^\star \in \imcone^\ast$ is related to the (unique) primal variable $w^\star$ by $\lambda^\star := (w^\star)^{\beta - 2} (w^\star - y)$.

In particular, $\lambda_k$ will converge to $\lambda^\star$ along subsequences, whence the convergence of $g(\lambda_k)$ to $\max_{\lambda \in \imcone^\ast} g(\lambda) = \min_{\mu \geq 0} \loss(\mu)$. Hence, if $y \notin \imcone$, we should have $\lim g(\lambda_k) > 0$ as $k \to + \infty$.

We will use this choice of $\lambda_k$ even for $\beta \notin [1,2]$, although no theoretical guarantee is of avail when it comes to the convergence to a certificate.

A caveat with our choice is that we should only expect $\lambda_k \in \imcone^\ast$ at the limit $k \to +\infty$, and not for a fixed iteration number $k$. In practice, if $\lambda_k \notin \imcone^\ast$, we set $\tilde \lambda_k := \lambda_k  + c$ where $c >0$ is a small constant, large enough to ensure that $A^\ast \tilde \lambda_k \geq 0$, {i.e.}, $\tilde \lambda_k \in \imcone^\ast$.

  \subsection{First example}
We first look at an example from PET, where the aim is to solve~\eqref{opt} with $y$.
A common way to do so is to use the following iterates, called \emph{multiplicative}~\cite{Lee2001, Fevotte2011}. Starting from some $\mu_0 \in \MM_+(K)$ (typically a positive constant over the domain), the iterates write
 \begin{equation}
 \label{multiplicative_updates}
 \mu_{k+1} = \mu_k \frac{A^\ast((A\mu_k)^{\beta - 2} y)}{A^\ast((A\mu_k)^{\beta - 1})}.
 \end{equation}
 These iterates have the ML-EM algorithm ($\beta = 1$) and the Iterative Image Space Reconstruction ($\beta = 2$) as particular cases~\cite{Depierro1993}, and proofs of convergence for these algorithms with any $\beta \in [0,2]$ can be found in~\cite{Yang2011}, in the finite-dimensional case.

One advantage of these algorithms is the decrease of the functional $\loss$ along iterates, a result we may generalise to the infinite-dimensional setting, when $\beta \in [1,2]$, see \autoref{appC}, adapting a proof of~\cite{Fevotte2011}. Up to our knowledge, such a result remains elusive even in the finite-dimensional case when $\beta \in [0,1)$, although simulations suggest it does hold true as well.

We now present the results of applying the algorithm in the case of a 2D PET operator $A$ with 90 views and 64 tangential positions (hence, $m = 5760$). Simulations are run using Python with the Operator Discretization Library~\cite{Adler2017}.
The sought-for image is taken to be the Derenzo phantom, denoted $\mu_r$. The data is obtained by (re-scaled) Poisson draws, with a time-variable (or dose-variable)~$t$ which accounts for the level of noise. In other words, $y \sim \frac{1}{t} \mathcal P(t A \mu_r)$, and the higher~$t$, the lower the noise. In order to approach the infinite-dimensional setting of our work, we take a fine image space discretisation of $512 \times 512$ pixels.

Finally, we take $\beta = 1.2$, on purpose not quite matching the noise statistics, as Poisson noise should lead one to take $\beta = 1$. We hence  mimic the situation of not knowing the exact noise statistics.

Figures~\ref{time0} for $t=1$ and~\ref{time-1} for $t=10^{-1}$ both display the evolution of the loss function $\loss$ along iterates, \textit{i.e}, $k \mapsto \loss(\mu_k)$, starting from $\mu_0 = 1$. As theoretically expected, the function decreases. We also plot the maximum attained for each reconstruction, namely $k \mapsto \max(\mu_k)$, which tends to increase. Finally, we show the reconstruction after $k=100$ and $k=1000$ iterates.

In the noisier case $t=0$, some pixels clearly take over as one keeps iterating. Moreover, we can certify that the we should indeed expect sparsity, as we may provide a dual certificates proving that $y \notin \imcone$. This is also suggested by the fact that $k \mapsto \loss(\mu_k)$ seems to converge to a positive value rather than to $0$.

\begin{figure}[h!]

\centering{
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=0.8\linewidth]{Figures/div_time0.png}
\caption{} \label{}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=0.8\linewidth]{Figures/max_time0.png}
\caption{} \label{}
\end{subfigure}
\\
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=1\linewidth]{Figures/100iter_time0_nobound.png}
\caption{} \label{}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=1\linewidth]{Figures/1000iter_time0_nobound.png}
\caption{} \label{}
\end{subfigure}
}
  \caption{Case $t=0$. (A) Divergence along iterates. (B) Maximum of reconstruction along iterates. (C) Reconstruction after $100$ iterates. (D) Reconstruction after $1000$ iterates.
  }
  \label{time0}
\end{figure}

In the less noisy case $t=10^{-1}$, it seems like the divergence $\loss(\mu_k)$ is not converging to $0$, which might be a hint that we should still expect sparsity. However, we are not able to certify it. Note that if sparsity were to hold true in this case as well, it may be that many more iterates are required to see sparsity arise more clearly.

\begin{figure}[h!]

\begin{subfigure}{0.49\textwidth}
\includegraphics[width=0.8\linewidth]{Figures/div_time-1.png}
\caption{} \label{}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=0.8\linewidth]{Figures/max_time-1.png}
\caption{} \label{}
\end{subfigure}
\\
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=1\linewidth]{Figures/100iter_time-1_nobound.png}
\caption{} \label{}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=1\linewidth]{Figures/1000iter_time-1_nobound.png}
\caption{} \label{}
\end{subfigure}
  \caption{Case $t=10^{-1}$. (A) Divergence along iterates. (B) Maximum of reconstruction along iterates. (C) Reconstruction after $100$ iterates. (D) Reconstruction after $1000$ iterates.}
  \label{time-1}
\end{figure}
  \subsection{Examples with \(\beta = 2\)}
\subsubsection{Toy Example}
We illustrate the sparsity with a toy example.
Here we choose \(\beta = 2\), \(K = [0,1]\), \(a_0 = 1\), \(a_1(x) = x\), with \(\ndet=2\).
In this case, one can compute the sparse solutions explicitly depending on the parameter \(y \in \R^{2}\), as shown in \autoref{fig:regions}.

We also look at the effect of regularisation.
In this case, we use total variation regularisation, that is, we solve
\[ \min_{\mu \geq 0}  \quad \loss(\mu) +\rho \operatorname{TV}(\mu),\] where \(\operatorname{TV}(\mu)\) is the total variation of the derivative of the measure \(\mu\) and \(\rho\) is a regularisation parameter.
In this discretised, one-dimensional setting, this is simply \(\operatorname{TV}(\mu) = \sum_i \abs{\mu^{i+1} - \mu^{i}}\), where \(\mu^i\) is the value of the discretised measure at pixel~\(i\).
We then compute the minimum using a primal-dual hybrid gradient method~\cite{ChaPo10}.
We plot the resulting minima for various values of the regularisation parameter \(\rho\) in \autoref{fig:toyexperiment}.

\begin{figure}
  \centering

  \begin{tikzpicture}[
    scale=1.5,
    vector/.style={draw,thick,->,},
    boundary/.style={draw,very thick,opacity=.6},
    zerob/.style={red},
    coneb/.style={green},
    caxis/.style={draw,->,very thick},
    isoval/.style={draw,thin,dashed},
    isolab/.style={sloped, fill=white,opacity=.9},
    reglab/.style={fill=white,opacity=.9,draw=black},
    ]
      \coordinate (O) at (0,0);
      \coordinate (A) at (2,0);
      \coordinate (M) at (2,2);
      \coordinate (P) at (-2,2);
      \coordinate (Q) at (0,-2);

      \draw[step=1cm,very thin, gray,opacity=.3] (-2,-2) grid(3,3);

      \draw[caxis] (-1,0) --  (3.5,0) node[below right]{$y_0$};
      \draw[caxis] (0,-1) -- (0,3.2) node[above]{$y_1$};

      \path[fill=green!40!white,draw=green,opacity=.3] (O) -- (3,0) -- (3,3) -- cycle;
      \path[fill=red!20!white,opacity=.3] (O) -- (P) -- (-2,-2) -- (Q) -- cycle;

      \path[boundary,coneb] (3,0) -- (O) -- (3,3);

      \path[boundary,zerob] (P) -- (O) -- (Q);

      \path[isoval] (1,1) -- node[isolab]{\(\xi=1\)} +(135:2.5);
      \path[isoval] (2,2) -- node[isolab]{\(\xi=2\)} +(135:1.5);

      \path[isoval] (0:1) -- node[isolab]{\(\xi=1\)} +(-90:2);
      \path[isoval] (0:2) -- node[isolab]{\(\xi=2\)} +(-90:2);

      \node[reglab] at (3,-2) {\(\mu=\xi \delta_0\)};
      \node[reglab] at (-1,2) {\(\mu=\xi \delta_1\)};
      \node[reglab] at (-1.2,.5) {\(\mu=0\)};
      \node[reglab] at (3,1) {\(A\mathcal{M}_+\)};

      \fill[blue,opacity=.75]  (0,1) circle (2pt);

    \end{tikzpicture}
    \caption{
      The sparse solutions of the problem \eqref{opt} with two detectors, \(a_0 = 1\), \(a_1(x)=x\) on the interval \(K=[0,1]\).
      There are three regions outside the cone \(A \MM_{+}\).
      In one of them, the optimal solution is zero, otherwise it takes the form \(\xi\delta_{x_0}\) with \(x_0\) being zero or one depending on the region.
    }
 \label{fig:regions}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=.8\textwidth]{Figures/regs-crop}
  \caption{
    The solutions for \(y_0=0\) and \(y_1=1\)  of the problem in \autoref{fig:regions} with total variation regularisation.
    Each curve corresponds to a different regularisation parameter, labelled by its base-10 logarithm.
    We see that when the regularisation parameter goes to zero, the computed solution converges to the expected exact solution, which,
    we see from \autoref{fig:regions},
    is \(\mu = .5 \delta_1\).
  }
  \label{fig:toyexperiment}
\end{figure}

\subsubsection{Tomography example}

The setting is similar to the previous section: in \autoref{fig:rhomax} we look at the effect of the regularisation parameter (with the same total variation regularisation), with the same optimisation algorithm (primal-dual hybrid gradient).
The image resolution is \(127 \times 127\), and there are \num{285} angles and \num{183} tangential coordinates.

\begin{figure}
  \centering
  \begin{subfigure}{.6\textwidth}
    \includegraphics[width=\textwidth]{Figures/rhomax}
    \caption{The maximum of the reconstruction for various regularisation parameters \(\rho\)}
  \end{subfigure}
  \\
  \begin{subfigure}{.4\textwidth}
    \includegraphics[width=\linewidth]{Figures/minregreco.png}
    \caption{Reconstruction with no regularisation (\(\rho=0\))}
  \end{subfigure}
  \begin{subfigure}{.4\textwidth}
    \includegraphics[width=\textwidth]{Figures/maxregreco.png}
    \caption{Reconstruction with high regularisation (\(\rho=6\))}
  \end{subfigure}
  \caption{Relation between regularisation parameter and the reconstruction of a Shepp--Logan phantom.}
  \label{fig:rhomax}
\end{figure}

\bibliography{biblio.bib}
\small
\bibliographystyle{acm}

\appendix
\normalsize
\section{Computation of the dual}
\label{appA}
Let us denote \[\psi_y (w) := d_\beta(y|w)- \lambda w.\]

\subsection{Case $\beta = 0$}
Recall that in this case, we always assume $y > 0$.
If $\lambda > 0$, we have $\psi_y(w) \sim - \lambda w$ as $w \to +\infty$, whence $h(y, \lambda) = - \infty$.
We now focus on the case $\lambda < 0$.
We have $\psi_y(w) \sim - \lambda w$ as $w \to +\infty$.  As $w \to 0$, we have $\psi_y(w) \sim \frac{y}{w}$, thus the function tends to $+\infty$ at both ends.
Since $\psi_y'(w) =  w^{- 2}(w - y) - \lambda$, $\psi_y^{(2)}(w) = w^{-3}( 2 y - w )$ for $w>0$, $w \mapsto \psi_y'(w)$ increases to the positive value $2 y$ and then decreases to $0^+$: there is a unique $w \geq 0$ minimising $\psi_y(w)$, which we denote  $w(y, \lambda)$, solving
\[w(y, \lambda)^{-2}(w(y, \lambda) - y) = \lambda.\]
A bit of algebra actually yields
$w(y, \lambda) = \frac{\sqrt{1-4 \lambda y} - 1}{-2 \lambda}.$

If $\lambda = 0$, $w \mapsto d_\beta(y|w)$ is minimised at $w = y$, and we may gather this case with $\lambda< 0$, since $w^{-2} (w - y) = \lambda$ indeed has solution $y$ for $\lambda = 0$, which amounts to setting $w(y, 0) = y$.

Summing up, we find

\[
  h(y, \lambda) = \begin{cases}
- \infty & \text{ if  \(\lambda > 0\)}
\\
d_\beta(y |w(y, \lambda)) - \lambda w(y, \lambda) & \text{ if \(\lambda \leq 0\)}
\end{cases}
\]
 In the last case, further computations lead to
 \[d_\beta(y |w(y, \lambda)) - \lambda w(y, \lambda) = \sqrt{1-4 \lambda y} - \ln\paren[\Big]{\frac{1}{2}\paren[\big]{\sqrt{1-4 \lambda y} + 1}}  - 1.\]

\subsection{Case $0 < \beta < 1$.}
If $\lambda > 0$, we have $\psi_y(w) \sim - \lambda w$ as $w \to +\infty$, whence $h(y, \lambda) = - \infty$ in this case.

We now focus on the case $\lambda < 0$.
We still have $\psi_y(w) \sim - \lambda w$ as $w \to +\infty$.
If $y>0$, $\psi_y(w) \sim - \frac{y}{\beta-1} w^{\beta - 1}$ as $w \to 0$, thus the function tends to $+\infty$ at both ends.
Since $\psi_y'(w) =  w^{\beta- 2}(w - y) - \lambda$, $\psi_y^{(2)}(w) = w^{\beta-3}((2 - \beta)y - (1-\beta)w )$ for $w>0$, $w \mapsto \psi_y'(w)$ increases to the positive value $\frac{2-\beta}{1-\beta} y$ and then decreases to $0^+$: there is a unique $w \geq 0$ minimising $\psi_y(w)$, which we denote  $w(y, \lambda)$, solving
\[w(y, \lambda)^{\beta- 2}(w(y, \lambda) - y) = \lambda.\]
If $y =0$, it is easily seen that the function $\psi_y(w) = \frac{1}{\beta}w^\beta - \lambda w$ is minimised at $w = 0$, with value $0$.

If $\lambda = 0$, $w \mapsto d_\beta(y|w)$ is minimised at $w = y$, and we may gather this case with $\lambda< 0$, since if $y=0$ we find $w = 0$ and if $y>0$, $w^{\beta- 2}(w - y) = 0$ leads to $w = y$, and we also here set $w(y, 0) = y$.

Summing up, we find

\[
  h(y, \lambda) = \begin{cases}
- \infty & \text{ if \(\lambda > 0\)}
\\
0 & \text{ if \(\lambda \leq 0, y = 0 \)}
\\
d_\beta(y |w(y, \lambda)) - \lambda w(y, \lambda)  & \text{ if \(\lambda \leq 0, y > 0 \)}
\end{cases}
\]

 \subsection{Case $ \beta = 1$.}
We have $\psi_y(w) \sim (1- \lambda) w$ as $w \to +\infty$, whence $h(y, \lambda) = - \infty$ if $\lambda >1$.
If $\lambda =  1$ and $y >0$, $\psi_y(w) \sim -y \ln(w)$ as $w \to 0$: the function tends to $-\infty$ as $w \to 0$ and $h(y, \lambda) = - \infty$. If $y = 0$, the function equals $0$ identically and its minimum is $0$.

We now focus on the case $\lambda < 1$.
We still have $\psi_y(w) \sim (1- \lambda) w$ as $w \rightarrow +\infty$.
If $y>0$, $\psi_y(w) \sim - y \ln(w)$ as $w \rightarrow 0$, thus the function tends to $+\infty$ at both ends.
Since $\psi_y$ is strictly convex in this case, it has a unique minimum for $w>0$, which we again denote $w(y, \lambda)$, solving
\[w(y, \lambda)^{-1}(w(y, \lambda) - y) = \lambda \iff w(y, \lambda) = \frac{y}{1- \lambda}. \]
If $y =0$, it is easily seen that the function $\psi_y(w) = \frac{1}{\beta}w^\beta - \lambda w$ is minimised at $w = 0$, with value $0$.

We may also gather the cases $y = 0$ and  $y>0$ whenever $\lambda \leq 1$, since the formula  for $w(y, \lambda)$ shows that it vanishes with $y$.

Summing up, we find

\[
  h(y, \lambda) = \begin{cases}
- \infty & \text{ if  \(\lambda > 1\)}
\\
d_\beta(y |w(y, \lambda)) - \lambda w(y, \lambda) & \text{ if \(\lambda \leq 1\)}
\end{cases}
\]
 In the last case, further computations lead to
 \[d_\beta(y |w(y, \lambda)) - \lambda w(y, \lambda) = y \ln(1-\lambda).\]

 \subsection{Case $ 1 < \beta < 2$.}
If $y>0$, $\psi_y(w) \sim \frac{w^\beta}{\beta}$ as $w \rightarrow +\infty$, and $\psi_y(w) \sim - \frac{y}{\beta-1} w^{\beta - 1}$ as $w \rightarrow 0$. The derivative $\psi'_y$ satisfies  $\psi'_y(w) \sim -y w^{\beta - 2}$ as $w \rightarrow 0$, and $\psi'_y(w) \sim w^{\beta - 1}$ as $w \rightarrow +\infty$. Thus the function $\psi'$ increases (as $\psi$ is convex) from $-\infty$ to $+\infty$. As a consequence, it has a unique minimum $w>0$, which we again denote $w(y, \lambda)$, solving
\[w(y, \lambda)^{\beta - 2}(w(y, \lambda) - y) = \lambda.\]
If $y =0$ and $\lambda \leq 0$, it is easily seen that the function $\psi_y(w) = \frac{1}{\beta}w^\beta - \lambda w$ is minimised at $w = 0$, with value $0$, whereas if $\lambda > 0$, it has a unique minimum (also defined by the equation for $w(y,\lambda)$).
Summing up, we find

\[
  h(y, \lambda) = \begin{cases}
0 & \text{ if \(y = 0, \lambda \leq 0\)}
\\
d_\beta(y |w(y, \lambda)) - \lambda w(y, \lambda) & \text{otherwise}.
\end{cases}
\]

  \subsection{Case $\beta = 2$.}
In this case, $\psi_y$ has a unique minimum $w$ given by $w = \lambda + y$, which gives the explicit formula
 \[ h(y, \lambda) = - \frac{1}{2}(\lambda  +y)^2 + \frac{1}{2} y^2.\]

\section{Further results regarding optimality}
\label{appB}
\begin{lemma}
\label{relation_opt}
Let $w^\star$ be optimal for the optimisation problem~\eqref{cone_formulation}.
Then there holds
\[ \sum_{i=1}^m (w_i^\star)^\beta = \sum_{i=1}^m y_i (w_i^\star)^{\beta -1}.\]
\end{lemma}

\begin{proof}
We consider small (multiplicative) variations around $w^\star$: in other words, we compute the first order expansion of the difference $\psi(\e) := D_\beta (y | w^\star) - D_\beta (y | (1-\e) w^\star)$ for $\e$ small (without any sign assumption), with $(1-\e) w^\star$ an admissible choice since $C$ is a cone.
We find
\begin{align*}\psi(\e) & = \frac{1}{\beta}  \sum_{i=1}^m(w_i^\star)^\beta \left(  1 - (1- \e)^{\beta} \right) + \frac{1}{1-\beta} \sum_{i=1}^m\,  y_i (w_i^\star)^{\beta - 1}\left(1 - (1- \e)^{\beta-1} \right) \\
 & = \e \left( \sum_{i=1}^m (w_i^\star)^\beta - \sum_{i=1}^m y_i (w_i^\star)^{\beta -1}\right) + o(\e).
\end{align*}
The first order term has to vanish for $(1-\varepsilon) w^\star$ not to violate the optimality of $w^\star$, whence the claim. Note that the final expansion obtained is also valid for $\beta = 0$ and $\beta = 1$.
\end{proof}

\begin{lemma}
\label{relation_opt}
Let $w^\star$ be optimal for the optimisation problem~\eqref{cone_formulation}.
Then for all $v \in C$ such that $\operatorname{supp}(v) \subset  \operatorname{supp}(w^\star)$, we have
\[ \sum_{i=1}^m (w_i^\star)^{\beta-2}(w_i^\star - y_i)(w_i^\star - v_i) \leq 0.\]
This inequality holds regardless of the assumption $\operatorname{supp}(v) \subset  \operatorname{supp}(w^\star)$ when $\beta \in \{0\} \cup [1,2]$.
\end{lemma}

\begin{proof}
Similarly to \autoref{relation_opt}, we perturb $w^\star$ but this time in the form $(1-\varepsilon)w^\star + \e v$ with $v \in C$ satisfying $\operatorname{supp}(v) \subset  \operatorname{supp}(w^\star)$, where $\e\in [0,1]$ is small. Denoting $\psi(\e) := D_\beta (y | w^\star) - D_\beta (y | ((1-\e)w^\star + \e v)$, we find $\psi(\e) = \sum_{i=1}^\ndet \psi_i(\e)$,
where
\begin{align*}\psi_i(\e):=\frac{1}{\beta} \big( (w_i^\star)^\beta - \left((1-\e)w_i^\star + \e v_i\right)^{\beta}\big) +\frac{1}{1-\beta} y_i\big( (w_i^\star)^{\beta-1} - \left((1-\e)w_i^\star + \e v_i\right)^{\beta-1}\big)
\end{align*}
For a given $i \in \{1,\ldots, m\}$, let us argue depending on whether $i \in \operatorname{supp}(w^\star)$ or not.

If $i \in \operatorname{supp}(w^\star)$, the functions $x \mapsto x^\beta$ and $x \mapsto x^{\beta -1}$ are  differentiable at $x = w_i^\star$ and we find by a Taylor expansion
\[\psi_i(\e) = \big( (w_i^\star)^{\beta-2}(w_i^\star - y_i)(w_i^\star - v_i) \big) \e + o(\e),\]

When $i \notin \operatorname{supp}(w^\star)$, \autoref{diff} ensures that $y_i = 0$, so the only remaining contribution is explicitly given by
\[\psi_i(\e) = - \frac{1}{\beta} v_i^\beta \e^\beta,\]

Note these two expansions also hold for $\beta = 0$ and $\beta = 1$ (the case $i \notin \operatorname{supp}(w^\star)$ is void when $\beta = 0$).

All in all, we find
\[\psi (\e) =  \left(\sum_{i \in\operatorname{supp}(w ^\star)} (w_i^\star)^{\beta-2}(w_i^\star - y_i)(w_i^\star - v_i) \right) \e + o(\e) - \frac{1}{\beta}\left(  \sum_{i \notin\operatorname{supp}(w ^\star)} v_i^\beta \right) \e^\beta  \]
For $\beta = 1$, the terms may be regrouped since they coincide: the leading term is $ \e\sum_{i=1}^m (w_i^\star)^{\beta-2}(w_i^\star - y_i)(w_i^\star - v_i) $, so the parenthesis has to be non-positive, as claimed.
For $\beta \in (1,2]$ and $\beta = 0$, the leading term is $ \e \sum_{i \in\operatorname{supp}(w ^\star)}  (w_i^\star)^{\beta-2}(w_i^\star - y_i)(w_i^\star - v_i) $ and has to be non-positive.

Finally, for $\beta \in (0,1)$, the leading term is  $- \frac{1}{\beta}\left(  \sum_{i \notin\operatorname{supp}(w ^\star)} v_i^\beta \right) \e^\beta$ but assumption $\operatorname{supp}(v) \subset  \operatorname{supp}(w^\star)$ shows that it vanishes.

Thus, in all cases, we obtain the expected inequality but with a sum ranging over $i \in\operatorname{supp}(w ^\star)$. That we may take the sum to range over $i \in \{1, \ldots, m\}$ is a consequence of \autoref{diff}: if $i \notin \operatorname{supp}(w ^\star)$, we have $y_i = 0$, whence $(w_i^\star)^{\beta -2} (y_i - w_i^\star) (v_i - w_i^\star) = 0$ and this term can be kept in the sum without affecting it.
\end{proof}

\section{Monotonicity along multiplicative updates}
\label{appC}

\begin{proposition}
Assume that $\beta\in [1,2]$. For any $\mu_0 \in \MM_+(K)$ with $\widetilde K \subset \supp(\mu_0)$, the iterates
 \[\mu_{k+1} = \mu_k \frac{A^\ast((A\mu_k)^{\beta - 2} y)}{A^\ast((A\mu_k)^{\beta - 1})},\]
 are well-defined and they make $\loss$ decrease, {i.e.},
 \[ \forall k \in \mathbb{N}, \quad D_\beta (y|A\mu_{k+1})  \leq D_\beta(y|A\mu_k).\]
\end{proposition}
\begin{proof}
Let us check that the assumption $\widetilde K \subset \supp(\mu_0)$ implies $\supp(\mu_k) = \widetilde K$ along iterates, making all divisions are well-defined (with the convention $0/0 = 0$).

Outside of $\supp(\mu_k)$, the value of the function $\frac{A^\ast((A\mu_k)^{\beta - 2} y)}{A^\ast((A\mu_k)^{\beta - 1})}$ does not matter (hence we do not care if there are divisions by $0$).
All we have to prove is that for all $k \in \mathbb{N}^*$, there holds
\begin{itemize}
\item $\supp(y) \subset \supp(A \mu_k)$,
\item $\supp(\mu_k) = \widetilde K$.
\end{itemize}
Then, the product $(A\mu_k)^{\beta - 2} y$ is indeed well-defined with the convention $0/0 = 0$, and the denominator $A^\ast((A\mu_k)^{\beta - 1})$ vanishes only on $\operatorname{int}(K \setminus \widetilde K)$, {i.e.}, outside of $\supp(\mu_k) = \widetilde K$. Hence $\mu_{k+1}$ is well-defined.

We proceed recursively. The second point is clear since we have $\supp(\mu_{k+1}) \subset \supp(\mu_k)$, and $\supp(A^\ast((A\mu_k)^{\beta - 2} y)) = \widetilde K$. For the first point, we fix $i \in I$ and prove that $(A\mu_{k+1})_i >0$.
Notice that the denominator$A^\ast((A\mu_k)^{\beta - 1})$ is bounded from above by some positive constant $M_k$.
We may write
\begin{align*} (A\mu_{k+1})_i &\geq y_i (A\mu_k)_i^{\beta-2} \int_K \frac{a_i^2}{A^\ast((A\mu_k)^{\beta - 1})} d\mu_k \geq \frac{y_i (A\mu_k)_i^{\beta-2}}{M_k} \int_K a_i^2  d\mu_k \\
 & \geq  \frac{y_i (A\mu_k)_i^{\beta-2}}{M_k}  \frac{ (A\mu_k)_i^2}{\mu_k(K)} = \frac{y_i (A\mu_k)_i^{\beta}}{M_k \mu_k(K)} > 0,
 \end{align*}

where we have used the Cauchy--Schwarz inequality $\int_K a_i^2 d\mu_k \geq \frac{(A\mu_k)_i^2}{ \mu_k(K)}$.

For $k \in \mathbb{N}$, and with the notation $w_k := A \mu_k$, consider the surrogate function, defined for $\mu \in \MM_+(K)$ absolutely continuous with respect to $\mu_k$ by
\[G_k(\mu) = \sum_{i=1}^\ndet \int_K d_\beta\bigg(y_i \bigg| (w_k)_i \frac{d \mu}{d \mu_k}\bigg) \frac{a_i}{(w_k)_i } \, d \mu_k.\]

\begin{enumerate}[label=\textit{Step \arabic*.}]
  \item
Let us prove that $G_k$ is above $\loss$ and coincides with it at $\mu_k$.
We compute
\begin{align*}G_k(\mu_k) & = \sum_{i=1}^\ndet \int_K d_\beta\Big(y_i \Big| (w_k)_i \Big) \frac{a_i}{(w_k)_i} \, d \mu_k \\
& = \sum_{i=1}^\ndet  d_\beta\Big(y_i \Big| (w_k)_i \Big)\underbrace{\int_K \frac{a_i}{(w_k)_i} \, d \mu_k}_{=1} = D_\beta (y |w_k) = \loss(\mu_k)
\end{align*}
Then, for a given $\mu \in \MM_+(K)$ absolutely continuous with respect to $\mu_k$ and by the convexity
of $v \mapsto d_\beta(u,v)$ with the probability measure  $ \nu_k := \frac{a_i}{(w_k)_i} \, \mu_k$,
\[ \int_K d_\beta\bigg(y_i \bigg| (w_k)_i \frac{d \mu}{d \mu_k}\bigg)\, d \nu_k \geq d_\beta\bigg(y_i \bigg| \int_K (w_k)_i \frac{d \mu}{d \mu_k} \, d\nu_k\bigg) = d_\beta(y_i |(A\mu)_i).\]
Summing on $i$ yields the result.

\item
We now prove
\[G_k(\mu_{k+1}) \leq G_k(\mu_{k}).\]
From the convexity of $v \mapsto d_\beta(u,v)$ and for $i \in I$,
\begin{align}
  d_\beta\bigg(y_i \bigg| (w_k)_i \frac{d \mu_{k}}{d \mu_k}\bigg) \geq & d_\beta\paren[\Big]{y_i \Big| (w_k)_i \frac{d \mu_{k+1}}{d \mu_k}} \\
  &+ \partial_v d_\beta\bigg(y_i \bigg|(w_k)_i \frac{d \mu_{k+1}}{d \mu_k}\bigg) (w_k)_i \frac{d(\mu_k - \mu_{k+1})}{d \mu_k},
  \end{align}
which yields
\begin{align}
G_k(\mu_k) & \geq G(\mu_{k+1}) + \sum_{i=1}^\ndet \int_K  \partial_v d_\beta\bigg(y_i \bigg|(w_k)_i \frac{d \mu_{k+1}}{d \mu_k}\bigg) (w_k)_i \frac{d(\mu_k - \mu_{k+1})}{d \mu_k} \frac{a_i}{(w_k)_i} \, d \mu_k  \\
                 & = G(\mu_{k+1}) +\int_K   \sum_{i=1}^\ndet a_i \, \partial_v d_\beta\bigg(y_i \bigg|(w_k)_i \frac{d \mu_{k+1}}{d \mu_k}\bigg) \, d(\mu_k - \mu_{k+1}).
\end{align}
We are done if we prove that \[\sum_{i=1}^\ndet a_i \, \partial_v d_\beta\bigg(y_i \bigg|(w_k)_i \frac{d \mu_{k+1}}{d \mu_k}\bigg) = A^\ast \partial_v d_\beta\bigg(y \bigg| w_k \frac{d \mu_{k+1}}{d \mu_k}\bigg)\] on $K$.
Since $\frac{d \mu_{k+1}}{d \mu_k} = \frac{A^\ast(w_k^{\beta - 2} y)}{A^\ast(w_k^{\beta - 1} y)}:= f_{k+1}$, and $\partial_v d(u|v) = v^{\beta-2}(v-u)$, we compute
\begin{align*} A^\ast \partial_v d_\beta\big(y \big| w_k f_{k+1} \big) & = f_{k+1}^{\beta - 2} A^\ast(w_k^{\beta-1} f_{k+1} - w_k^{\beta-2} y) \\
& =  f_{k+1}^{\beta - 2}\left(f_{k+1} A^\ast w_k^{\beta-1} - A^\ast (w_k^{\beta - 1} y)\right) = 0
\end{align*} by the definition of $f_{k+1}$. We may now conclude as the combination of the two results entails for all $k \in \mathbb{N}$
\[\loss (\mu_{k+1}) \leq G_k(\mu_{k+1}) \leq G_k(\mu_k) = \loss(\mu_k).\]
\end{enumerate}
\end{proof}
\end{document}

config: arxiv