[1,2]Marta Lazzaretti

1]Dipartimento di Matematica, Università di Genova, Via Dodecaneso 35, Genova,
16145, Italy

2]I3S Laboratory, CNRS-UniCA-Inria, 2000 Route des Lucioles, Sophia Antipolis,
06903, France

3]Laboratory of GTPases and Neurosignalling, Achucarro Basque Center for
Neuroscience, Parque Científico UPV/EHU, Edif. Sede, Planta 3, Leioa, E-48940,
Spain

# Off-the-grid regularisation for Poisson inverse problems

<EMAIL_ADDRESS>Claudio Estatico<EMAIL_ADDRESS>Alejandro
Melero Carrillo<EMAIL_ADDRESS>Luca Calatroni
<EMAIL_ADDRESS>[  [  [

###### Abstract

Off-the-grid regularisation has been extensively employed over the last decade
in the context of ill-posed inverse problems formulated in the continuous
setting of the space of Radon measures ${\mathcal{M}(\Omega)}$. These
approaches enjoy convexity and counteract the discretisation biases as well
the numerical instabilities typical of their discrete counterparts. In the
framework of sparse reconstruction of discrete point measures (sum of weighted
Diracs), a Total Variation regularisation norm in ${\mathcal{M}(\Omega)}$ is
typically combined with an $L^{2}$ data term modelling additive Gaussian
noise. To asses the framework of off-the-grid regularisation in the presence
of signal-dependent Poisson noise, we consider in this work a variational
model coupling the Total Variation regularisation with a Kullback-Leibler data
term under a non-negativity constraint. Analytically, we study the optimality
conditions of the composite functional and analyse its dual problem. Then, we
consider an homotopy strategy to select an optimal regularisation parameter
and use it within a Sliding Frank-Wolfe algorithm. Several numerical
experiments on both 1D/2D simulated and real 3D fluorescent microscopy data
are reported.

###### keywords:

Off-the-grid sparse regularisation, Poisson noise, Sliding Frank-Wolfe,
fluorescence microscopy imaging.

##  1 Introduction

Discrete (or _on-the-grid_) sparse optimisation approaches are nowadays
established techniques in the field of mathematical signal and image inverse
problems. In the context of linear ill-posed inverse problems, for instance,
they aim at retrieving a sparse approximation of a quantity of interest $x$
(e.g., an image) from blurred, noisy and potentially incomplete acquisitions
$y\in\mathbb{R}^{M}$ on a regular grid of size $N\geq M$, see Fig. 1(a). The
discretisation parameter $N$ determines the localisation precision, as only
the centres of the discretisation intervals are candidates for reconstruction.
As such, in case of rough discretisations, very approximate reconstructions
may be computed (see, e.g., Fig. 1(b) where small values $M=N$ are used). To
obtain higher precision, one typically chooses a grid-size parameter $N>M$ for
reconstruction, as in Fig. 1(c). Choosing a large value for $N$, however, may
cause instabilities in the reconstructions [1] due to the higher numerical
complexity.

Off-the-grid approaches aim at overcoming such difficulties. They can be
thought indeed as the natural framework to deal with the case $N\to+\infty$ of
on-the-grid formulations [2, 3]. In such framework, the spatial domain
$\Omega\subseteq\mathbb{R}^{d}$ is not discretised by a regular grid, but,
rather, the quantity of interest is modelled as an element of a suitable
functional space defined on $\Omega$. A natural framework for spike
reconstruction problems, is, for instance, the space of Radon measures
$\mu\in{\mathcal{M}(\Omega)}$ where the quantity of interest can be modelled
as the measure $\mu_{a,x}=\sum_{i=1}^{n}a_{i}\delta_{x_{i}}$ with $n$ being
the number of spikes and where the question is therefore how to estimate the
number of spikes $n$, and then how to retrieve both continuous positions
$x_{i}\in\Omega$ and amplitudes $a_{i}\geq 0$ for $i=1,\ldots,n$.

Inverse problems in the space of measures and off-the-grid optimisation
methods have been first proposed in [4, 5, 6, 7] and since then they have been
a topic of intense research activity for the mathematical community, both from
an analytical and numerical viewpoint, see, e.g., [8, 9, 10, 11, 12, 13, 14].
Off-the-grid methods have been proved to be particularly useful in
applications where fine-scale details need to be retrieved from noisy
acquisitions, such as spike detection in astronomy and microscopy [12], as
well as parameter estimation in in spectroscopy [15] and density mixture
estimation [16]. Standard approaches in this setting usually combine off-the-
grid regularisation with an additive (Gaussian) noise modelling on the
underlying signal, which is in general easier to work with, from both an
analytical and a computational point of view. Such modelling corresponds to
the well-studied variational formulation of the Beurling LASSO (BLASSO) model,
[5, 7].

In this work, we consider an off-the-grid modelling in ${\mathcal{M}(\Omega)}$
under the hypothesis of signal-dependent Poisson noise in the data. This
choice is motivated by some particular biological applications of interest,
such as fluorescence microscopy, where, due to the photon emission nature of
the light, Poisson noise is better suited than the Gaussian one to describe
the process of photon counts on acquired images [17]. While in a discrete
setting, a precise modelling of noise statistics if often not necessary due to
the inevitable biases introduced by the choice of the regularisation employed,
a natural question is whether whenever a refined off-the-grid regularisation
is used, a precise noise model could indeed be relevant. As a graphical
visualisation, we report in Figure 1 a visual comparison between
reconstructions obtained with on-the-grid approaches and the off-the-grid
approach proposed in this work to solve a spike-deconvolution problem under
the choice of the Poisson data term considered in this work.

(a) Acquisition $y\in\mathbb{R}^{M}$

(b) Discrete reconstruction $x\in\mathbb{R}^{M}$

(c) Discrete reconstruction $x\in\mathbb{R}^{N}$

(d) Off-the-grid rec. $\mu_{a,x}\in{\mathcal{M}(\Omega)}$

Figure 1: Comparison between discrete (on-the-grid) and off-the-grid sparse
reconstructions with a Poisson data term. In black: the ground-truth spikes to
retrieve. In Fig.1(a): the acquired blurred and noisy signal
$y\in\mathbb{R}^{M}$ lying on a low-resolution grid of size $M$. In Fig.1(b):
in red, a discrete reconstruction with support on a grid with $M$ pixels. In
Fig.1(c): in red, discrete reconstruction with support on a grid with $N>M$
pixels. In Fig.1(d): in green, off-the-grid reconstruction.

###  1.1 Inverse problems in the space of Radon measures
${\mathcal{M}(\Omega)}$

In this section we recall some important definitions and properties of the
space of Radon measures ${\mathcal{M}(\Omega)}$ following [18, 12, 3]. For
more details on Radon measures see [19, 20].

Let $\Omega\subseteq\mathbb{R}^{d}$, with $d\in\mathbb{N}$, $d\geq 1$, be a
compact subset of $\mathbb{R}^{d}$ with non-empty interior. We denote by
$\mathscr{C}_{0}(\Omega,\mathbb{R})$ the space of real continuous functions on
$\Omega$ vanishing at infinity, namely all the continuous maps
$\psi:\Omega\rightarrow\mathbb{R}$ such that:

$\forall\varepsilon>0,\quad\exists K\subset\Omega,\ K\text{ compact, s.t.
}\quad\sup_{x\in\Omega\backslash K}|\psi(x)|\leq\varepsilon.$

The space of Radon measures is defined through duality.

######  Definition 1.

The Banach space ${\mathcal{M}(\Omega)}$ of real signed Radon measures on
$\Omega$ is the topological dual of $\mathscr{C}_{0}(\Omega,\mathbb{R})$
endowed with the supremum norm $\|\cdot\|_{\infty,\Omega}$, defined by
$\|\psi\|_{\infty,\Omega}:=\sup_{x\in\Omega}|\psi(x)|$.

This definition allows to characterise any measure
$\mu\in{\mathcal{M}(\Omega)}$ as a continuous linear form evaluated on
functions $\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ in terms of the duality
pairing, so that

$\forall\mu\in{\mathcal{M}(\Omega)}\qquad\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}=\int_{\Omega}\psi\mathrm{d}\mu.$
(1)

A Radon measure $\mu\in{\mathcal{M}(\Omega)}$ is a positive measure if
$\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}$
is positive or null for any non-negative function
$\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})$. This specifies the meaning of the
term signed in the above definition, as the quantity
$\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}$
can be also negative. Note that the space of positive (or non-negative) Radon
measures ${\mathcal{M}^{+}(\Omega)}$ cannot be defined as the topological dual
of $\mathscr{C}_{0}(\Omega,\mathbb{R}^{+})$, since the latter is not a vector
space, since trivially if $\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R}^{+})$ then
$-\psi\not\in\mathscr{C}_{0}(\Omega,\mathbb{R}^{+})$.

It can be shown that ${\mathcal{M}(\Omega)}$ is a non-reflexive Banach space
endowed with the Total Variation (TV) norm, here defined
$\forall\mu\in{\mathcal{M}(\Omega)}$ as:

$|\mu|(\Omega)=\sup\Big{(}\int_{\Omega}\psi\mathrm{d}\mu\leavevmode\nobreak\
\Big{|}\leavevmode\nobreak\ \ \psi\in\mathscr{C}_{0}(\Omega,\mathbb{R}),\
\|\psi\|_{\infty,\Omega}\leq 1\Big{)}.$ (2)

The TV norm is convex and non-differentiable, hence it is possible to consider
its subdifferential [5, 8]:

$\partial|\mu|(\Omega)=\left\\{\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})|\
\|\psi\|_{\infty,\Omega}\leq 1\text{ and
}\int_{\Omega}\psi\mathrm{d}\mu=|\mu|(\Omega)\right\\}.$ (3)

Note that for sparse discrete measures, i.e. weighted sums of Diracs:

$\mu_{a,x}=\sum_{i=1}^{N}a_{i}\delta_{x_{i}}\text{ with }N\in\mathbb{N},\
a=(a_{1},\dots,a_{N})\in\mathbb{R}^{N},\ x=(x_{1},\dots,x_{N})\in\Omega^{N},$
(4)

the TV norm coincides with the $L^{1}$ norm of the amplitudes vector $a$, that
is $|\mu_{a,x}|(\Omega)=\|a\|_{1}$. This explains why the TV norm is
considered in this setting as a generalisation of the $L^{1}$ norm. Moreover,
in this special case
$\,\langle\psi,\mu\rangle=\sum_{i=1}^{N}a_{i}\psi(x_{i})\,$ and the
subdifferential has the following expression involving the sign function

$\partial|\mu_{a,x}|(\Omega)=\left\\{\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})|\
\|\psi\|_{\infty,\Omega}\leq 1,\ \forall i=1,\ldots,N\
\psi(x_{i})=\operatornamewithlimits{sign}(a_{i})\right\\}.$ (5)

We focus now on the formulation of linear inverse problems in
${\mathcal{M}(\Omega)}$. As an acquisition space we will consider a Hilbert
space $\mathcal{H}$. Let $\mu\in{\mathcal{M}(\Omega)}$ be the unknown source
measure. We consider an acquisition $\bar{y}\in\mathcal{H}$ being the result
of the action of the forward operator
$\Phi:{\mathcal{M}(\Omega)}\rightarrow\mathcal{H}$ evaluated on $\mu$, with a
fixed measurement kernel $\varphi:\Omega\rightarrow\mathcal{H}$, that is:

$\bar{y}:=\Phi\mu=\int_{\Omega}\varphi(x)\mathrm{d}\mu(x).$ (6)

Note that this integral should not be confused with the concept of duality
pairing (1), which is a scalar function defined as the integral over $\Omega$
of a continuous real function with respect to a measure
$\mu\in{\mathcal{M}(\Omega)}$. In (6), the integral is then a Böchner integral
[20] for vector-valued functions $\varphi(x)\in\mathcal{H}$, that is
$\varphi(x)$ is an element of $\mathcal{H}$, and it is well-defined if
$\varphi$ is continuous and bounded [21, 12].

The choice of the kernel $\varphi$ and of the acquisition space $\mathcal{H}$
depends on the specific physical acquisition process considered. In the
following, we consider a convolution kernel, which is of practical interest in
fluorescence microscopy (see, e.g., [12, 22]). In this setting, a natural
choice is therefore $\mathcal{H}=L^{2}(\Omega)$, with the convolution kernel
$\varphi:\Omega\rightarrow L^{2}(\Omega)$ defined in terms of a Point Spread
Function (PSF) $\tilde{\varphi}:\Omega\rightarrow L^{2}(\Omega)$ acting as:

$\varphi(x):=(s\mapsto\tilde{\varphi}(s-x))\in L^{2}(\Omega),\ \forall
x\in\mathbb{R}.$ (7)

Note that depending on the microscopy technique used one can have different
PSFs. For instance, the Gaussian PSF, centred in $c\in\Omega$ with radius
$\sigma>0$, is defined by

$s\mapsto\tilde{\varphi}(s-c):={(2\pi\sigma^{2})}^{-d/2}e^{-\|s-c\|_{2}^{2}/2\sigma^{2}}.$
(8)

For other possible choices of measurement kernels, we refer the reader to
[12].

We observe that the action of the forward operator $\Phi$ on finite linear
combination of Diracs (4) can be explicited as follows:

$\Phi\mu_{a,x}=\int_{\Omega}\varphi(x)\mathrm{d}\mu(x)=\sum_{i=1}^{N}a_{i}\varphi(x_{i}).$
(9)

For simplicity, the following notation will therefore be used in the sequel
$\Phi_{x}(a)=\sum_{i=1}^{N}a_{i}\varphi(x_{i})$ to denote $\Phi\mu_{a,x}$.

In fluorescence microscopy imaging applications, the objects of interest are
usually images of molecules, i.e. single point-sources emitting fluorescent
light. The unknown is therefore well-described by non-negative discrete
measures of the form (4). Other possible objects of interest in this field are
microtubules, 1-dimensional curve structures [23, 24], and cells
(2-dimensional), that can be modelled as piece-wise constant functions [25,
26, 27]. For simplicity, we will focus in the following only on inverse
problems aiming at recovering (0-dimensional) measures in the form (4). The
extension to more general regularisation models is left for future work.

###  1.2 The BLASSO problem: formulation, duality and optimality conditions

The standard sparse spike deconvolution problem, consists in recovering a
(small) finite linear combination of Diracs
$\mu_{a,x}=\sum_{i=1}^{N}a_{i}\delta_{x_{i}}$ from a blurred and noisy
acquisition

$y=\Phi\mu_{a,x}+\omega,$ (10)

where $\Phi:\Omega\rightarrow L^{2}(\Omega)$ is a convolution operator, as
defined in (7), and $\omega$ is an additive noise component, typically
describing white Gaussian noise. The variational formulation of the problem of
retrieving an estimate $\mu$ from $y$ is

$\operatornamewithlimits{argmin}_{\mu\in{\mathcal{M}(\Omega)}}\leavevmode\nobreak\
T_{\lambda}(\mu):=\frac{1}{2}\|\Phi\mu-y\|^{2}+\lambda|\mu|(\Omega),\quad\lambda>0$
($L^{2}-|\cdot|$)

which is the Beurling-LASSO (BLASSO) functional, named after the work of the
mathematician Beurling [28]. It is considered the generalisation of the
discrete LASSO variational problem

$\operatornamewithlimits{argmin}_{x\in\mathbb{R}^{N}}\frac{1}{2}\|\tilde{\Phi}x-y\|^{2}+\lambda\|x\|_{1},$
($L^{2}-L^{1}$)

where $\tilde{\Phi}\in\mathbb{R}^{M\times N}$ is a discretisation of
$\Phi:{\mathcal{M}(\Omega)}\rightarrow\mathcal{H}$ and $y\in\mathbb{R}^{M}$.
In [5], the functonal $T_{\lambda}:{\mathcal{M}(\Omega)}\rightarrow\mathbb{R}$
is proved to be proper, convex and coercive, which proves the existence of
solutions of ($L^{2}-|\cdot|$) and uniqueness under injectivity on $\Phi$.

By optimality, $\mu_{\lambda}$ minimises $T_{\lambda}$ if and only if
$0\in\partial T_{\lambda}(\mu_{\lambda})$, i.e.

$0\in\partial
T_{\lambda}(\mu_{\lambda})=\Phi^{*}(\Phi\mu_{\lambda}-y)+\lambda\partial|\mu_{\lambda}|(\Omega)\,\iff\,\frac{1}{\lambda}\Phi^{*}(y-\Phi\mu_{\lambda})\in\partial|\mu_{\lambda}|(\Omega),$
(11)

which can be written as

$\eta\in\partial|\mu_{\lambda}|(\Omega)\qquad\text{with}\qquad\eta:=\frac{1}{\lambda}\Phi^{*}(y-\Phi\mu_{\lambda}),$
(12)

where $\eta$ is the so-called dual certificate of ($L^{2}-|\cdot|$). The dual
certificate formally generalises the concept of Euler equation, playing a
crucial role in the characterisation of optimality conditions for
($L^{2}-|\cdot|$) and in devising optimisation algorithms in this setting, as
better specified in the following. If $\mu$ is a finite linear combination of
Dirac masses (4), recalling (5), we get that optimality conditions (12) become
simply

$\eta(x_{i})=\operatornamewithlimits{sign}(a_{i})\quad\wedge\quad\|\eta\|_{\infty}\leq
1.$ (13)

Note that the optimality conditions (13) can be also derived by studying the
dual problem of ($L^{2}-|\cdot|$):

$\operatornamewithlimits{argmax}_{\|\Phi^{*}p^{*}\|_{\infty,\Omega}\leq
1}\langle y,p^{*}\rangle-\frac{\lambda}{2}\|p^{*}\|^{2},$ (14)

which provides further insights to the meaning of the dual certificate $\eta$.
Given $\mu_{\lambda}\in{\mathcal{M}(\Omega)}$ solution of to ($L^{2}-|\cdot|$)
and denoting by $p^{*}_{\lambda}\in L^{2}(\Omega)$ the solution of the dual
problem (14), the following extremality conditions hold true

$\begin{cases}&\Phi^{*}p^{*}_{\lambda}\in\partial|\mu_{\lambda}|(\Omega),\\\
&-p^{*}_{\lambda}=\frac{1}{\lambda}\left(\Phi\mu_{\lambda}-y\right)\end{cases}.$
(15)

From (15), we retrieve in fact the optimality conditions (13) expressed in
terms of the dual certificate (12). In addition, when the dual certificate
satisfies (13), by duality we have $\eta=\Phi^{*}p^{*}_{\lambda}$ with
$p^{*}_{\lambda}$ being a solution of (14). Optimality conditions thus fully
characterise the solution(s) of the BLASSO problem ($L^{2}-|\cdot|$). They are
indeed crucial in devising algorithms for its minimisation and, in particular,
in the definition of good stopping criterions.

##  2 Off-the-grid Poisson inverse problems

We now present the signal-dependent Poisson modelling studied in this work,
which differs from (10) as it is not additive.

Let $\Phi:{\mathcal{M}(\Omega)}\rightarrow L^{2}(\Omega)$ be the forward
operator in (6). In the following, $\Phi$ is assumed to be a positive definite
operator, i.e.

$\mu\in{\mathcal{M}^{+}(\Omega)}\text{ positive measure
}\Rightarrow\Phi\mu(x)\geq 0\quad\forall x\in\Omega.$ (16)

Observe that whenever the measurement kernel $\varphi$ (6) is non-negative,
this always holds. This is the case, for instance, in standard image
deblurring problems. For simplicity, we introduce the notation
$L^{2}(\Omega)^{+}=\\{f\in L^{2}(\Omega)\text{ such that }f(x)>0\text{ a .e.
}x\in\Omega\\}$.

The Poisson sparse spike deconvolution problem in an off-the-grid setting
consists in finding a positive discrete measure
$\mu\in{\mathcal{M}^{+}(\Omega)}$ such that

$y=\mathcal{P}(\Phi\mu+b),$ (17)

where $b\in L^{2}(\Omega)^{+}$ is a strictly positive background term, and
$y\in L^{2}(\Omega)$ is a realisation of a Poisson distributed random process
$\mathsf{Y}\sim\mathcal{P}(\Phi\mu+b)$. Similarly as in the standard discrete
setting, we consider the Kullback-Leibler divergence as a suitable data term
[29, 30, 31, 32, 33].

######  Definition 2.

The Kullback-Leibler divergence $\mathcal{D}_{KL}:L^{2}(\Omega)^{+}\times
L^{2}(\Omega)^{+}\longrightarrow\mathbb{R}$ is defined by

$\mathcal{D}_{KL}(s,t):=\int_{\Omega}s(x)-t(x)+t(x)\Big{(}\log(t(x))-\log(s(x))\Big{)}\
\mathrm{d}x.$ (18)

To consider $\mathcal{D}_{KL}(\Phi\mu+b,y)$ as data fidelity term, two issues
have to be considered:

  * •

if $\mu$ is not a positive measure, then $\Phi\mu+b$ might not be positive;

  * •

the noisy acquisition $y$ might still vanish in a non negligible region of the
domain $\Omega$. This is due to the fact that a Poisson random variable with
mean $\alpha$ vanishes with positive probability equal to $e^{-\alpha}$. For
(18) to be well defined, it is thus required that $y>0$ almost everywhere.

To solve the first issue, we introduce the function
$\tilde{\mathcal{D}}_{KL}:L^{2}(\Omega)\times
L^{2}(\Omega)^{+}\longrightarrow\mathbb{R}\cup\\{+\infty\\}$ which extends
(18) as

$\tilde{\mathcal{D}}_{KL}(s,t)=\begin{cases}\mathcal{D}_{KL}(s,t)&s\in
L^{2}(\Omega)^{+}\\\ +\infty&s\not\in L^{2}(\Omega)^{+}\end{cases}.$ (19)

Moreover, we just restrict our study to positive acquisitions $y$, thus
requiring

$y\in L^{2}(\Omega)^{+}.$ (20)

Under assumptions (16) and (20) and using (19), the quantity
$\tilde{\mathcal{D}}_{KL}(\Phi\mu+b,y)$ is therefore well-defined for all
$\mu\in{\mathcal{M}(\Omega)}$.

We can thus introduce the following variational problem for solving (17)

$\boxed{\operatornamewithlimits{argmin}_{\mu\in{\mathcal{M}(\Omega)}}\leavevmode\nobreak\
\tilde{\mathcal{D}}_{KL}(\Phi\mu+b,y)+\lambda|\mu|(\Omega)+\mathds{1}_{\\{{\mathcal{M}^{+}(\Omega)}\\}}(\mu),\quad\lambda>0.}$
($\tilde{\mathcal{D}}_{KL}-|\cdot|$)

The Poisson fidelity (19) is there coupled with the TV norm, together with the
indicator function of the positive measures ${\mathcal{M}^{+}(\Omega)}$
defined by

$\mathds{1}_{\\{{\mathcal{M}^{+}(\Omega)}\\}}(\mu)=\begin{cases}0&\mu\geq 0\\\
+\infty&\text{otherwise}\end{cases},$ (21)

to ensure $\mu$ to be non-negative.

###  2.1 Dual problem and optimality conditions

We analyse in the following the dual problem of
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) and provide an analytical expression of
the convex conjugate of the involved functions. We refer to Appendix A for the
definition of convex conjugate [34] and classical concepts about duality [35].

The study of the dual problem of the problem
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) requires the computation of the convex
conjugate
$F^{*}:\mathscr{C}_{0}(\Omega,\mathbb{R})\rightarrow\mathbb{R}\cup\\{+\infty\\}$
of the penalty term

$F:{\mathcal{M}(\Omega)}\rightarrow\mathbb{R}\cup\\{+\infty\\},\quad
F(\cdot):=|\cdot|(\Omega)+\mathds{1}_{\\{{\mathcal{M}^{+}(\Omega)}\\}}(\cdot),$
(22)

and of the convex conjugate
$G^{*}:L^{2}(\Omega)\rightarrow\mathbb{R}\cup\\{+\infty\\}$ of the fidelity
term

$G:L^{2}(\Omega)\rightarrow\mathbb{R}\cup\\{+\infty\\},\quad
G(\cdot):=\frac{1}{\lambda}\tilde{\mathcal{D}}_{KL}(\cdot,y).$ (23)

To compute the convex conjugate of the Kullback-Leibler functional (23), we
start considering the one-dimensional Kullback-Leibler function, defined by

$g_{t}(s)=\frac{1}{\lambda}\big{(}s-t+t\log(t)-t\log(s)\big{)},\qquad
s,t>0\text{ and }\lambda>0.$ (24)

Applying the definition of convex conjugate (86) to $g_{t}$ yields

$\displaystyle g_{t}^{*}(s^{*})$ $\displaystyle=\sup_{s>0}\leavevmode\nobreak\
ss^{*}-g_{t}(s)=\sup_{s>0}\leavevmode\nobreak\
ss^{*}-\frac{1}{\lambda}\Big{(}s-t+t\log(t)-t\log(s)\Big{)}=$
$\displaystyle=\sup_{s>0}\leavevmode\nobreak\
\underbrace{s\Big{(}s^{*}-\frac{1}{\lambda}\Big{)}+\frac{t}{\lambda}\log(s)+\frac{t}{\lambda}-\frac{t}{\lambda}\log(t)}_{h(s)}.$

We have two cases:

  * (i)

If $s^{*}\geq\frac{1}{\lambda}$, then $\lim_{s\to+\infty}h(s)=+\infty$ implies
$\sup_{s>0}h(s)=+\infty\Rightarrow g_{t}\big{(}s^{*}\big{)}=+\infty$ for all
$t>0$.

  * (ii)

If $0<s^{*}<\frac{1}{\lambda}$, then $\lim_{s\to\pm\infty}h(s)=-\infty$. Thus,
being $h$ a convex and differentiable function its supremum is attained at
$\hat{s}$ such that $h^{\prime}(\hat{s})=0$, which can be computed by

$\displaystyle
h^{\prime}(\hat{s})=s^{*}-\frac{1}{\lambda}+\frac{t}{\lambda\hat{s}}=\frac{\lambda\hat{s}s^{*}-\hat{s}+t}{\lambda\hat{s}}=0$
$\displaystyle\iff\lambda\hat{s}s^{*}-\hat{s}+t=0\iff\hat{s}=\frac{t}{1-\lambda
s^{*}}.$

Thus,

$g_{t}^{*}(s^{*})=h\Big{(}\frac{t}{1-\lambda
s^{*}}\Big{)}=-\frac{t}{\lambda}\log(1-\lambda s^{*}).$ (25)

Observe that $g_{t}^{*}(s^{*})$ is well defined since
$s^{*}<\frac{1}{\lambda}$.

Hence, the convex conjugate $g_{t}^{*}$ of $g_{t}$ is

$g^{*}_{t}(s^{*})=\begin{cases}+\infty&s^{*}\geq\frac{1}{\lambda}\\\
-\frac{t}{\lambda}\log(1-\lambda s^{*})&s^{*}<\frac{1}{\lambda}\end{cases}\,.$
(26)

Since $\tilde{\mathcal{D}}_{KL}(\cdot,t)$ is defined also for non-positive
functions, its one-dimensional counterpart is given by
$\tilde{g}_{t}:\mathbb{R}\longrightarrow\mathbb{R}\cup\\{+\infty\\}$ where

$\tilde{g}_{t}(s)=\begin{cases}g_{t}(s)&s>0\\\ +\infty&s\leq 0\end{cases}\,,$
(27)

so that its convex conjugate coincides with (26). We thus have the following
lemma.

######  Lemma 1.

Let $G:L^{2}(\Omega)\rightarrow\mathbb{R}\cup\\{+\infty\\}$ be the function
defined by

$G(\cdot):=\frac{1}{\lambda}\tilde{\mathcal{D}}_{KL}(\cdot,y),$ (28)

where $\tilde{\mathcal{D}}_{KL}$ is given by (19). The convex conjugate of $G$
is given by $G^{*}:L^{2}(\Omega)\rightarrow\mathbb{R}\cup\\{+\infty\\}$
defined by

$\displaystyle G^{*}(p^{*})=\begin{cases}+\infty&p^{*}\geq\frac{1}{\lambda}\\\
-\frac{y}{\lambda}\log(1-\lambda
p^{*})&p^{*}<\frac{1}{\lambda}\end{cases}=-\frac{y}{\lambda}\log(1-\lambda
p^{*})\mathds{1}_{\\{z<1/\lambda\\}}(p^{*})\,,$ (29)

where $-\frac{y}{\lambda}\log(1-\lambda
p^{*})=\langle-\frac{y}{\lambda},\log(1-\lambda p^{*})\rangle\in\mathbb{R}$.

We compute now the convex conjugate of $F$ defined in (22). To simplify the
notation, we denote in the following by
$A,B:{\mathcal{M}(\Omega)}\rightarrow\mathbb{R}\cup{\\{+\infty\\}}$ the
quantities

$A(\cdot)=|\cdot|(\Omega),\qquad
B(\cdot)=\mathds{1}_{\\{{\mathcal{M}^{+}(\Omega)}\\}}(\cdot).$ (30)

The convex conjugate of $F$ can be obtained as infimal convolution of the
convex conjugate $A^{*}$ and $B^{*}$ [36, 37], that is

$F^{*}(\psi)=\min_{\psi_{1}+\psi_{2}=\psi}A^{*}(\psi_{1})+B^{*}(\psi_{2}),$
(31)

see Appendix A for more details and for the the computation of $A^{*}$, which
is a well-known result. We report here the computation of the convex conjugate
of $B(\cdot)$.

######  Lemma 2.

The convex conjugate of the indicator function of positive measures $B$ of
(30) is the function
$B^{*}:\mathscr{C}_{0}(\Omega,\mathbb{R})\rightarrow\mathbb{R}\cup\\{+\infty\\}$
given by

$B^{*}(\psi)=\begin{cases}0&\psi(x)\leq 0\quad\forall x\in\Omega\\\
+\infty&\exists x\in\Omega\text{ s.t. }\psi(x)>0\end{cases}.$ (32)

###### Proof.

By Definition 86, for any $\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ we get

$\displaystyle B^{*}(\psi)$
$\displaystyle=\sup_{\mu\in{\mathcal{M}(\Omega)}}\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}-B(\mu)=\sup_{\mu\in{\mathcal{M}^{+}(\Omega)}}\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}$
$\displaystyle\geq\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}\quad\forall\mu\in{\mathcal{M}^{+}(\Omega)}.$

If there exists $\bar{x}\in\Omega$ such that $\psi(\bar{x})>0$, consider the
positive measure $\bar{\mu}=\alpha\delta_{\bar{x}}$ with $\alpha>0$. Then,
$B^{*}(\psi)\geq\langle\psi,\bar{\mu}\rangle=\alpha\psi(\bar{x})$ and taking
the limit for $\alpha\rightarrow\infty$ of the latter inequality yields
$B^{*}(\psi)\geq+\infty$ if $\psi(x)>1$. On the other hand, if $\psi(x)\leq 0\
\forall x\in\Omega$, then
$\langle\psi,\mu\rangle=\int_{\Omega}\psi\mathrm{d}\mu\leq 0\
\forall\mu\in{\mathcal{M}^{+}(\Omega)}$. Moreover, $\langle\psi,0\rangle=0$.
Thus, $B^{*}(\psi)=0$ if $\psi(x)\leq 0\ \forall x\in\Omega$, which concludes
the proof. ∎

We thus get an expression for $F^{*}$ in the following lemma.

######  Lemma 3.

Let $F:{\mathcal{M}(\Omega)}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ be the
functional defined in (22). For $\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})$,
its convex conjugate
$F^{*}:\mathscr{C}_{0}(\Omega,\mathbb{R})\rightarrow\mathbb{R}\cup\\{+\infty\\}$
is defined by

$F^{*}(\psi)=\begin{cases}0&\text{if}\ \psi(x)\leq 1\quad\forall x\in\Omega\\\
+\infty&\text{otherwise}\end{cases}.$ (33)

###### Proof.

By (31), we have

$\displaystyle F^{*}(\psi)$
$\displaystyle=\min_{\psi_{1}+\psi_{2}=\psi}A^{*}(\psi_{1})+B^{*}(\psi_{2})$
$\displaystyle=\min_{\psi_{1}+\psi_{2}=\psi}\begin{cases}0&\|\psi_{1}\|_{\infty,\Omega}\leq
1\\\
+\infty&\|\psi_{1}\|_{\infty,\Omega}>1\end{cases}+\begin{cases}0&\psi_{2}(x)\leq
0\ \forall x\in\Omega\\\ +\infty&\exists x\in\Omega\text{ s.t.
}\psi_{2}(x)>0\end{cases}.$

Notice that $F^{*}(\psi)\geq 0\
\forall\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ and we have either
$F^{*}(\psi)=0$ or $F^{*}(\psi)=+\infty$. Thus, $F^{*}(\psi)=0$ if and only if
there exist $\psi_{1},\psi_{2}\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ such that
$\psi=\psi_{1}+\psi_{2}$ with $A^{*}(\psi_{1})=0$ and $B^{*}(\psi_{2})=0$. On
the contrary, $F^{*}(\psi)=+\infty$ if for all
$\psi_{1},\psi_{2}\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ such that
$\psi=\psi_{1}+\psi_{2}$ it holds $A^{*}(\psi_{1})=+\infty$ or
$B^{*}(\psi_{2})=+\infty$. Consider now the three following cases:

  * i)

If $\|\psi\|_{\infty,\Omega}\leq 1$, consider $\psi_{1}=\psi$. This implies
$A^{*}(\psi_{1})=0$, and $\psi_{2}=0$, which implies $B^{*}(\psi_{2})=0$. Thus
$F^{*}(\psi)=0$.

  * ii)

If $\|\psi\|_{\infty,\Omega}>1$ and there exists $\bar{x}$ such that
$\psi(\bar{x})>1$, we have $F^{*}(\psi)=+\infty$. Indeed, assume by
contradiction that $F^{*}(\psi)=0$, that is there exist
$\psi_{1},\psi_{2}\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ such that
$\psi_{1}+\psi_{2}=\psi$ with $\|\psi_{1}\|_{\infty,\Omega}\leq 1$ and
$\psi_{2}(x)\leq 0\ \forall x\in\Omega$. We would have

$\psi(\bar{x})=\underbrace{\psi_{1}(\bar{x})}_{\leq
1}+\underbrace{\psi_{2}(\bar{x})}_{\leq 0}\leq 1$

which is absurd, since $\psi(\bar{x})>1$. Then, $\forall\psi_{1},\psi_{2}$
such that $\psi_{1}+\psi_{2}=\psi$ either $\|\psi_{1}\|_{\infty,\Omega}>1$ or
$\psi_{2}(x)>0$ for some $x\in\Omega$. Hence, $F^{*}(\psi)=+\infty$.

  * iii)

If $\|\psi\|_{\infty,\Omega}>1$ and $\psi(x)\leq 1\ \forall x\in\Omega$,
consider $\psi_{1}=\psi^{+}$ and $\psi_{2}=\psi^{-}$, with

$\psi^{+}(x)=\begin{cases}\psi(x)&\psi(x)\geq 0\\\
0&\psi(x)<0\end{cases}\qquad\psi^{-}(x)=\begin{cases}0&\psi(x)>0\\\
\psi(x)&\psi(x)\leq 0\end{cases}\,.$

It is evident that $\psi=\psi^{+}+\psi^{-}$. Moreover,
$\|\psi^{+}\|_{\infty,\Omega}\leq 1$, thus $A^{*}(\psi^{+})=0$, and $\forall\
x\in\Omega\ \psi^{-}(x)\leq 0$, thus $B^{*}(\psi^{-})=0$, so that
$F^{*}(\psi)=0$.

Since $(i)$ and $(iii)$ are equivalent to $\psi(x)\leq 1,\;\forall
x\in\Omega$, and the three items cover all possible cases, the thesis is
proved. ∎

###  2.2 Dual problem

The study of dual problems of sparse-regularisation models with non-negativity
constraints has been carried out in [38, 39] in the discrete setting of LASSO,
and in [40] for the discrete counterpart of
($\tilde{\mathcal{D}}_{KL}-|\cdot|$), that is with the Kullback-Leibler
divergence as fidelity and the $L^{1}$ penalty. In the following we use these
results to compute the dual problem of ($\tilde{\mathcal{D}}_{KL}-|\cdot|$) in
${\mathcal{M}(\Omega)}$.

By Fenchel-Rockafellar duality (Lemma 94 in Appendix A), the dual problem of
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) can be obtained by plugging (29) and (33)
into (93). We thus have:

$\displaystyle\operatornamewithlimits{argmax}_{p^{*}\in
L^{2}(\Omega)}\leavevmode\nobreak\ -F^{*}(\Phi^{*}p^{*})-G^{*}(-p^{*})$ (34)
$\displaystyle=$ $\displaystyle\operatornamewithlimits{argmax}_{p^{*}\in
L^{2}(\Omega)}\leavevmode\nobreak\
-F^{*}(\Phi^{*}p^{*})+\begin{cases}-\infty&-p^{*}\geq\frac{1}{\lambda}\\\
-\frac{b}{\lambda}(1+\lambda p^{*})+\frac{y}{\lambda}\log(1+\lambda
p^{*})&-p^{*}<\frac{1}{\lambda}\end{cases}$ (35) $\displaystyle=$
$\displaystyle\operatornamewithlimits{argmax}_{p^{*}\in
L^{2}(\Omega)}\leavevmode\nobreak\
-F^{*}(\Phi^{*}p^{*})+\begin{cases}-\infty&p^{*}\leq-\frac{1}{\lambda}\\\
-\frac{b}{\lambda}(1+\lambda p^{*})+\frac{y}{\lambda}\log(1+\lambda
p^{*})&p^{*}>-\frac{1}{\lambda}\end{cases}$ (36) $\displaystyle=$
$\displaystyle\operatornamewithlimits{argmax}_{p^{*}\in L^{2}(\Omega)\text{
s.t.
}p^{*}>-\frac{1}{\lambda}}-F^{*}(\Phi^{*}p^{*})-\frac{b}{\lambda}(1+\lambda
p^{*})+\frac{y}{\lambda}\log(1+\lambda p^{*})$ (37) $\displaystyle=$
$\displaystyle\operatornamewithlimits{argmax}_{p^{*}\in L^{2}(\Omega)\text{
s.t. }p^{*}>-\frac{1}{\lambda}}\begin{cases}\;0&\forall
x\in\Omega,\;\,\Phi^{*}p^{*}(x)\leq 1\\\ -\infty&\exists
x\in\Omega,\;\,\Phi^{*}p^{*}(x)>1\end{cases}\quad-\frac{b}{\lambda}(1+\lambda
p^{*})+\frac{y}{\lambda}\log(1+\lambda p^{*})$ (38) $\displaystyle=$
$\displaystyle\operatornamewithlimits{argmax}_{p^{*}\in\mathcal{D}}\leavevmode\nobreak\
-\frac{b}{\lambda}(1+\lambda p^{*})+\frac{y}{\lambda}\log(1+\lambda p^{*})\,,$
(39)

where $\mathcal{D}=\\{p^{*}\in L^{2}(\Omega):p^{*}>-\frac{1}{\lambda}\text{
and }\forall x\in\Omega,\ \ \Phi^{*}p^{*}(x)\leq 1\\}$.

###  2.3 Extremality conditions

By Lemma 94, extremality conditions (94) can be obtained. Given
$\mu_{\lambda}\in{\mathcal{M}(\Omega)}$ solution of the primal problem
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) with regularisation parameter $\lambda>0$
and $p^{*}_{\lambda}\in L^{2}(\Omega)$ solution of the dual problem (39), they
read

$\begin{cases}\Phi^{*}p^{*}_{\lambda}\in{\partial
F(\mu_{\lambda})=\partial|\mu_{\lambda}|(\Omega)+\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu_{\lambda})}\\\
-p^{*}_{\lambda}\in\frac{1}{\lambda}\partial_{1}\tilde{\mathcal{D}}_{KL}(\Phi\mu_{\lambda}+b,y)=\frac{1}{\lambda}\Big{(}I-\frac{y}{\Phi\mu_{\lambda}+b}\Big{)}\end{cases},$
(40)

where $\partial_{1}\tilde{\mathcal{D}}_{KL}(\Phi\mu_{\lambda}+b,y)$ denotes
the subdifferential of $\tilde{\mathcal{D}}_{KL}(\cdot,\cdot)$ computed with
respect to the first variable and evaluated in $(\Phi\mu_{\lambda}+b,y)$ and
$I$ denotes the identity operator.

######  Remark 1.

If $\mu_{\lambda}\in{\mathcal{M}(\Omega)}$ is solution of the primal problem
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) and $p^{*}_{\lambda}\in L^{2}(\Omega)$ is
solution of the dual problem (39), then from (40) we have

$-p^{*}_{\lambda}=\frac{1}{\lambda}\Big{(}I-\frac{y}{\phi\mu_{\lambda}+b}\Big{)}\quad\Rightarrow\quad
p^{*}_{\lambda}=\frac{y-\Phi\mu_{\lambda}-b}{\lambda(\Phi\mu_{\lambda}+b)}.$
(41)

It follows that ${p^{*}_{\lambda}>-\frac{1}{\lambda}\iff y>0}$, which holds by
hypothesis (20).

In the extremality conditions (40), the set-sum
$\partial|\mu_{\lambda}|(\Omega)+\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu_{\lambda})$
requires attention. Recall that, in general, for two proper convex functions
$f_{1},f_{2}:\mathcal{X}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ there holds
$\partial(f_{1}+f_{2})\subseteq\partial f_{1}+\partial f_{2}$, with the
equality holding if and only if
$\text{int}(\text{dom}(f_{1}))\cap\text{int}(\text{dom}(f_{2}))\not=\emptyset$,
see, e.g., [41].

In the following we show that $\partial
F(\cdot)=\partial|\cdot|(\Omega)+\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\cdot)$.
Let $\mu\in{\mathcal{M}^{+}(\Omega)}$ and consider $B$ as in (30). By
definition, the subdifferential of $B$ at $\mu\in{\mathcal{M}^{+}(\Omega)}$ is

$\displaystyle\partial B(\mu)$
$\displaystyle=\\{\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R})|\
B(\tilde{\mu})\geq B(\mu)+\langle\eta,\tilde{\mu}-\mu\rangle\
\forall\tilde{\mu}\in{\mathcal{M}(\Omega)}\\}$ (42)
$\displaystyle=\\{\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R})|\
\langle\eta,\tilde{\mu}-\mu\rangle\leq 0\
\forall\tilde{\mu}\in{\mathcal{M}^{+}(\Omega)}\\},$ (43)

that is the normal cone to ${\mathcal{M}^{+}(\Omega)}$ at $\mu$. We now
consider (43) for some scenarios of interest:

  * •

If $\mu=0$, we have $\eta\in\partial B(0)$ if and only if
$\forall\tilde{\mu}\in{\mathcal{M}^{+}(\Omega)}$ there holds
$\langle\eta,\tilde{\mu}\rangle\leq 0$. Since $\tilde{\mu}$ is a positive
measure, we have that
$\langle\eta,\tilde{\mu}\rangle=\int_{\Omega}\eta\mathrm{d}\tilde{\mu}\leq 0$
for any $\eta$ such that $\eta(x)\leq 0$ for all $x\in\Omega$. Thus,

$\partial B(0)=\\{\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R})|\quad\eta(x)\leq
0\quad\forall x\in\Omega\\}.$

  * •

If $\mu=\bar{a}\delta_{\bar{x}}$ with $\bar{a}>0$, $\bar{x}\in\Omega$, since
$\eta\in\partial B(\mu)$ if and only if
$\langle\eta,\tilde{\mu}-\mu\rangle\leq
0\;\forall\tilde{\mu}\in{\mathcal{M}^{+}(\Omega)}$, then $\eta\in\partial
B(\mu)$ if and only if

$\displaystyle\int_{\Omega}\eta\mathrm{d}\tilde{\mu}=\langle\eta,\tilde{\mu}\rangle\leq\langle\eta,\mu\rangle=\int_{\Omega}\eta\mathrm{d}\mu=\bar{a}\eta(\bar{x}).$

Consider $\tilde{\mu}=2\mu$. Hence,
$\int_{\Omega}\eta\mathrm{d}\tilde{\mu}=\int_{\Omega}\eta
2\mathrm{d}\mu=2\bar{a}\eta(\bar{x})$ and
$2\bar{a}\eta(\bar{x})\leq\bar{a}\eta(\bar{x})$ if and only if
$\eta(\bar{x})=0$. Thus, $\eta\in\partial B(\mu)$ implies $\eta(\bar{x})=0$.

Let $\eta\in\partial B(\mu)$, which entails $\eta(\bar{x})=0$, and suppose
$\exists\tilde{x}\in\Omega$ s.t. $\eta(\tilde{x})>0$. Then,

$\forall\tilde{\mu}\in{\mathcal{M}^{+}(\Omega)}\qquad\int_{\Omega}\eta\mathrm{d}\tilde{\mu}\leq\int_{\Omega}\eta\mathrm{d}\mu=\bar{a}\eta(\bar{x})=0.$

The above inequality with $\tilde{\mu}=\tilde{a}\delta_{\tilde{x}}$,
$\tilde{a}>0$, yields
$0\geq\int_{\Omega}\eta\mathrm{d}\tilde{\mu}=\tilde{a}\eta(\tilde{x})>0$,
which is absurd. Thus, $\eta(x)\leq 0\;\forall x\in\Omega$ and

$\partial B(\mu)=\\{\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R})|\ \eta(x)\leq
0\;\forall x\in\Omega\ \text{ and }\eta(\bar{x})=0\\}.$

  * •

Similarly, if $\mu=\sum_{i=1}^{N}a_{i}\delta_{x_{i}}$ with $a_{i}>0$ and
$x_{i}\in\Omega$ for $i=1,\ldots,N$, then we have:

$\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu)=\partial
B(\mu)=\left\\{\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R})\Big{|}\ \eta(x)\leq
0\;\forall x\in\Omega\ \text{ and }\eta(x_{i})=0,\ i=1,\ldots,N\right\\}.$
(44)

We now have an analytical expression for both $\partial|\cdot|(\Omega)$ and
$\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\cdot)$, whenever $\mu$ is a
finite linear combination of Diracs. We thus state the following proposition.

######  Proposition 1.

If $\mu=\sum_{i=1}^{N}a_{i}\delta_{x_{i}}$ with $a_{i}>0$, $x_{i}\in\Omega$,
$i=1,\ldots,N$, we have

$\partial|\mu|(\Omega)+\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu)=\left\\{\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R}\Big{|}\
\eta(x)\leq 1\ \forall x\in\Omega\text{ and }\eta(x_{i})=1,\
i=1,\ldots,N\right\\}.$ (45)

###### Proof.

If $\mu=\sum_{i=1}^{N}a_{i}\delta_{x_{i}}$ with $a_{i}>0$, $x_{i}\in\Omega$,
by (5) and (44) we have

$\partial|\mu|(\Omega)=\\{\eta_{1}\in\mathscr{C}_{0}(\Omega,\mathbb{R}):\
\|\eta_{1}\|_{\infty}\leq 1\text{ and }\eta_{1}(x_{i})=1,\ i=1,\ldots,N\\},$
(46)
$\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu)=\\{\eta_{2}\in\mathscr{C}_{0}(\Omega,\mathbb{R}):\
\forall x\in\Omega\ \eta_{2}(x)\leq 0\text{ and }\eta_{2}(x_{i})=0,\
i=1,\ldots,N\\}.$ (47)

Thus, we have the following inclusion

$\partial|\mu|(\Omega)+\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu)\subseteq\\{\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R})\big{|}\
\eta(x)\leq 1\;\forall x\in\Omega\text{ and }\eta(x_{i})=1,\ i=1,\ldots,N\\}.$
(48)

To show that the latter is indeed an equality, we verify now the opposite
inclusion. Let $\eta\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ such that
$\eta(x)\leq 1\ x\in\Omega\text{ and }\eta(x_{i})=1,\ i=1,\ldots,N$. Our aim
is to prove that there exist $\eta_{1}\in\partial|\mu|(\Omega)$ and
$\eta_{2}\in\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu)$ such that
$\eta=\eta_{1}+\eta_{2}$. Let’s consider the two cases $\|\eta\|_{\infty}\leq
1$ and $\|\eta\|_{\infty}>1$.

  * i)

If $\|\eta\|_{\infty}\leq 1$, we can take $\eta_{1}=\eta$, $\eta_{2}=0$.

  * ii)

If $\|\eta\|_{\infty}>1$, we can take

$\eta_{1}(x)=\begin{cases}\eta(x)&-1\leq\eta(x)\leq 1\\\
-1&\eta(x)\leq-1\end{cases},\qquad\eta_{2}(x)=\begin{cases}0&-1\leq\eta(x)\leq
1\\\ \eta(x)+1&\eta(x)\leq-1\end{cases}.$ (49)

In both cases, we have found $\eta_{1}\in\partial|\mu|(\Omega)$ and
$\eta_{2}\in\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu)$ such that
$\eta=\eta_{1}+\eta_{2}$, so that the opposite inclusion also holds. This
concludes the proof. ∎

Thanks to (45) of the latter proposition, it is now possible to better
characterise the extremality conditions (40), under the assumption that
$\mu_{\lambda}$, solution of ($\tilde{\mathcal{D}}_{KL}-|\cdot|$), is a
discrete measure. Indeed, if
$\mu_{\lambda}=\sum_{i=1}^{N_{\lambda}}(a_{\lambda})_{i}\delta_{(x_{\lambda})_{i}}$
we have that (40) is equivalent to

$\Phi^{*}p^{*}_{\lambda}(x)\leq 1\quad\forall
x\in\Omega\quad\text{and}\quad\Phi^{*}p^{*}_{\lambda}\big{(}(x_{\lambda})_{i}\big{)}=1,\
i=1,\ldots,N_{\lambda},$ (50)

with $p^{*}_{\lambda}$ solution of the dual problem (39), explicitly given by
(41).

It is possible to obtain a result similar to (50) by analysing the optimality
conditions of ($\tilde{\mathcal{D}}_{KL}-|\cdot|$). The functional

$T_{\lambda}^{KL}(\mu):=\tilde{\mathcal{D}}_{KL}(\Phi\mu+b,y)+\lambda|\mu|(\Omega)+\mathds{1}_{\\{{\mathcal{M}^{+}(\Omega)}\\}}(\mu)$
(51)

is convex. Hence $\mu_{\lambda}$ minimises $T_{\lambda}^{KL}$ if and only if
$0\in\partial T_{\lambda}^{KL}(\mu_{\lambda})$, i.e.

$\displaystyle 0$
$\displaystyle\in\Phi^{*}\partial_{1}\tilde{\mathcal{D}}_{KL}(\Phi\mu_{\lambda},y)+\lambda\partial|\mu_{\lambda}|(\Omega)+\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu_{\lambda})$
(52)
$\displaystyle=\Phi^{*}\Big{(}I-\frac{y}{\Phi\mu_{\lambda}+b}\Big{)}+\lambda\partial|\mu_{\lambda}|(\Omega)+\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu_{\lambda})\,,$
(53)

or equivalently

$\frac{1}{\lambda}\Phi^{*}\Big{(}\frac{y}{\Phi\mu_{\lambda}+b}-I\Big{)}\in\partial|\mu_{\lambda}|(\Omega)+\frac{1}{\lambda}\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu_{\lambda})\,,$

which can be written as

$\eta\in\partial|\mu_{\lambda}|(\Omega)+\frac{1}{\lambda}\partial\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu_{\lambda})\,,$
(54)

by defining the dual certificate of ($\tilde{\mathcal{D}}_{KL}-|\cdot|$)

$\eta:=\frac{1}{\lambda}\Phi^{*}\Big{(}\frac{y}{\Phi\mu_{\lambda}+b}-I\Big{)}.$
(55)

If $\mu$ is a finite linear combination of Diracs, Proposition 1 yields the
explicit expression (45) for the sum of the two subdifferentials of (54).
Hence, inclusion (54) simply becomes

$\eta(x)\leq 1\;\;\forall
x\in\Omega\quad\text{and}\quad\eta(x_{i})=1,\;i=1,\ldots,N.$ (56)

Note the similarity between (56) and (50), and of the dual certificate (55)
and the quantity $\Phi^{*}p^{*}_{\lambda}$ in the extremality conditions (40).

##  3 The Sliding Frank-Wolfe algorithm

Both problems ($\tilde{\mathcal{D}}_{KL}-|\cdot|$) and ($L^{2}-|\cdot|$) are
defined over the space ${\mathcal{M}(\Omega)}$, an infinite dimensional non-
reflexive Banach space. Due to non-reflexivity, it is not simple to define
therein proximal-based algorithms, see [42] for some recent attempt. Any
solver for such problems shall take into account the infinite dimensional
nature of ${\mathcal{M}(\Omega)}$: popularly used algorithms in this setting
are semi-definite programming approaches (for Fourier measurements) [7],
conditional gradient algorithms [43, 44, 12] and particle gradient descent
[21, 45], which is an optimal-transport based algorithm. In this work, we will
only focus on conditional gradient strategies, namely the Frank-Wolfe and
Sliding Frank-Wolfe algorithms, see [46] for a survey.

Algorithm 1 Sliding Frank-Wolfe (SFW) algorithm [12]

Initialisation: $\mu^{0}=0$.

repeat for $k=1,2,\ldots,K_{\text{max}}$

  * $\mu^{k}=\sum_{i=1}^{N^{k}}a_{i}^{k}\delta_{x_{i}^{k}},a_{i}^{k}\in\mathbb{R},x_{i}^{k}\in\Omega$, find $x_{*}^{k}\in\Omega$ s.t.:

$x_{*}^{k}\in\operatornamewithlimits{argmax}_{x\in\Omega}\left|\eta(\lambda,\mu^{k})(x)\right|\quad\text{where}\quad\eta(\lambda,\mu^{k})\quad\text{is
defined by \eqref{eq_gen_dual_cert} with }\mu=\mu^{k}$ (57)

  * if $\left|\eta(\lambda,\mu^{k})\left(x_{*}^{k}\right)\right|\leq 1$

    * $\mu^{k}$ is a solution of ($\mathcal{P}(\lambda)$) $\Rightarrow$ stop

  * else

    * –

insertion step: add support for the new spike and amplitudes’ estimation

$\displaystyle x^{k+1/2}=\left(x_{1}^{k},\ldots,x_{Nk}^{k},x_{*}^{k}\right)$
$\displaystyle
a^{k+1/2}\in\operatornamewithlimits{argmin}_{a\in\mathbb{R}^{N^{k}+1}}f_{y^{\delta},b}\left(\Phi_{x^{k+1/2}}a\right)+\lambda\|a\|_{1}+\alpha\mathds{1}_{\geq
0}(a)$ (58)
$\displaystyle\text{Update:}\qquad\mu^{k+1/2}=\sum_{i=1}^{N^{k}}a_{i}^{k+1/2}\delta_{x_{i}^{k}}+a_{N^{k}+1}^{k+1/2}\delta_{x_{*}^{k}}$
(59)

    * –

sliding step: using a non-convex solver initialised with
$\left(a^{k+1/2},x^{k+1/2}\right)$

$\displaystyle\left(a^{k+1},x^{k+1}\right)\in\underset{(a,x)\in\mathbb{R}^{N^{k}+1}\times\Omega^{N^{k}+1}}{\arg\min}f_{y^{\delta},b}\left(\Phi_{x}a\right)+\lambda\|a\|_{1}+\alpha\mathds{1}_{\geq
0}(a)$ (60)
$\displaystyle\text{Update:}\qquad\mu^{k+1}=\sum_{i=1}^{N^{k}+1}a_{i}^{k+1}\delta_{x_{i}^{k+1}}$
(61)

    * –

pruning: eventually remove zero amplitudes Dirac masses from $\mu^{k+1}$

until convergence

In [12], the authors detail how the conditional gradient algorithm, also known
as Frank-Wolfe algorithm [43], can be used to minimise the BLASSO
($L^{2}-|\cdot|$) functional using an epigraphic lift, and then propose the
Sliding version, which significantly improves the reconstruction quality by
adding an extra step where both positions and amplitudes are re-optimised.
Similar strategies are here used to minimise
($\tilde{\mathcal{D}}_{KL}-|\cdot|$). We report here the statement of the
result without a proof as it follows from [12, Lemma 4].

######  Lemma 4.

The solution $\bar{\mu}\in{\mathcal{M}(\Omega)}$ to
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) is equivalent to the solution
$\bar{\mu}\in{\mathcal{M}(\Omega)}$ to the problem

$\operatornamewithlimits{argmin}_{(t,\mu)\in C}\leavevmode\nobreak\
\tilde{T}_{\lambda}^{KL}(\mu,t):=\tilde{\mathcal{D}}_{KL}(\Phi\mu+b,y)+\lambda
t+\mathds{1}_{\\{{\mathcal{M}^{+}(\Omega)}\\}}(\mu),$ (62)

where
$C:=\left\\{(t,\mu)\in\mathbb{R}^{+}\times{\mathcal{M}(\Omega)};|\mu|(\Omega)\leq
t\leq M\right\\}\quad\text{with
}M:=\frac{\tilde{\mathcal{D}}_{KL}(b,y)}{2\lambda}$.

Given a convex and smooth fidelity functional
$f_{y^{\delta},b}:{\mathcal{M}(\Omega)}\rightarrow\mathbb{R}$ defined in terms
of the observed data $y^{\delta}\in L^{2}(\Omega)$ and, potentially, a
background term $b^{\delta}\in L^{2}(\Omega)+$, we consider the following
general optimisation problem:

$\operatornamewithlimits{argmin}_{\mu\in{\mathcal{M}(\Omega)}}\leavevmode\nobreak\
f_{y^{\delta},b}(\Phi\mu)+\lambda|\mu|(\Omega)+\alpha\mathds{1}_{{\mathcal{M}^{+}(\Omega)}}(\mu),\qquad\lambda>0,\quad\alpha\in\\{0,1\\},$
($\mathcal{P}(\lambda)$)

which encompasses the BLASSO as well as problem ($\mathcal{P}(\lambda)$) and
where $\alpha\in\\{0,1\\}$ may enforce non-negativity constraints.

The pseudocode corresponding to the optimisation of the general problem
($\mathcal{P}(\lambda)$) is reported in Algorithm 1. Observe that the stopping
criterion is expressed in terms of the dual certificate $\eta(\lambda,\mu)\in
L^{2}(\Omega)$ of the general problem ($\mathcal{P}(\lambda)$), which is
defined in terms of the subgradient of the fidelity and reads

$\eta(\lambda,\mu)=\frac{1}{\lambda}\tilde{\eta}(\mu)\quad\text{with}\quad\tilde{\eta}(\mu):=\begin{cases}-\Phi^{*}\frac{\partial}{\partial\mu}f_{y^{\delta},b}(\Phi\mu)&\alpha=0\\\
\left(-\Phi^{*}\frac{\partial}{\partial\mu}f_{y^{\delta},b}(\Phi\mu)\right)_{+}&\alpha=1\end{cases},$
(63)

depending on the parameter $\alpha$ so that the dual certificate (12) of
($L^{2}-|\cdot|$) and the dual certificate (55) of
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) can be obtained by choosing $\alpha$ and
$f_{y^{\delta},b}$, respectively. The generalised optimality condition for
($\mathcal{P}(\lambda)$) reads

$\|\eta(\lambda,\mu)\|_{\infty}\leq 1,$ (64)

and, under the hypothesis that the solution is a finite linear combination of
Diracs as (4), the dual certificate of the solution satisfies
$\eta(\lambda,\mu)(x_{i})=\operatornamewithlimits{sign}(a_{i})$ where the
points $x_{i}\in\Omega$ are the support of $\mu$.

##  4 Parameter selection via algorithmic homotopy

(a) $\lambda\approx 0$

(b) $\lambda=1$

(c) $\lambda=3000$

Figure 2: Reconstructions obtained using SFW for a 1D sparse deconvolution
problem with Poisson noise. Ground truth spikes (black) and reconstructed ones
(green) using Algorithm 1 for some choices of $\lambda$ are shown. When
$\lambda\approx 0$, the number and intensities of spikes are overestimated,
while for $\lambda\gg 0$ they are underestimated.

The performance of Algorithm 1 for solving ($\mathcal{P}(\lambda)$) depends on
the choice of the regularisation parameter $\lambda>0$: it plays indeed a
fundamental role in the sparsity pattern of the solution and in the
enforcement of the stopping rule,
$\max_{x\in\Omega}\left|\eta(\lambda,\mu^{k})(x)\right|\leq 1$. Namely, a high
value of $\lambda$ forces only few iterations of the algorithm to be
performed. Such choice impacts also both the estimation and the sliding steps
of Algorithm 1, being it associated with the sparsity-promoting $L^{1}$
penalty used to compute the amplitude vector. On the other hand, smaller
values of $\lambda$ provide a better data fit, with more spikes with higher
intensities, see Figure 2.

In [47], the authors propose a method based on algorithmic homotopy [48, 49]
to choose an optimal regularisation parameter for the resolution of BLASSO
($L^{2}-|\cdot|$). The idea behind homotopy algorithms is to avoid the
exploration of the whole Pareto frontier by grid search, while providing an
iterative procedure computing only few parameters up to a target value. More
precisely, starting from an initial overestimated value $\lambda_{1}>0$, a
solution $\mu_{\lambda_{1}}$ to ($\mathcal{P}(\lambda)$) is computed by, e.g.,
Algorithm 1. At each homotopy iteration, if the solution does not fit well the
data up to some tolerance $\sigma_{\text{target}}$ depending on the noise
magnitude $\delta$, then $\lambda$ is decreased. A new solution
$\mu_{\lambda_{2}}$ to ($\mathcal{P}(\lambda)$) is thus computed in the next
homotopy step and so on. Homotopy algorithms thus explore the Pareto frontier
for a small set of values $\lambda$ and select its biggest value for which the
solution to ($\mathcal{P}(\lambda)$) meets a convergence criterion depending
on $\delta$. In Figure 3, one can see in red the discrete values produced by
the homotopy strategy we are going to describe, which stops when the fidelity
term goes under the value of $\sigma_{\text{target}}(\delta)$, in grey.

Figure 3: Homotopy algorithms explore the Pareto frontier iteratively for a
strictly decreasing sequence of regularisation parameters $\lambda$. In blue:
fine discretisation of the Pareto frontier with grid search. In red: homotopy
iterations, corresponding to different values of $\lambda$. In gray: target
value for the fidelity.

Input: $y\in L^{2}(\Omega)$, $b\in L^{2}(\Omega),\ b\geq 0$,
$\Phi\in\mathcal{L}({\mathcal{M}(\Omega)},L^{2}(\Omega))$

Parameters: $\gamma\in(0,1)$, $c>0$, $\sigma_{\text{target}}>0$

Output: optimal $\hat{\mu}\in{\mathcal{M}(\Omega)}$

Initialisation: $\hat{\mu}_{0}\in{\mathcal{M}(\Omega)}$ and
$\lambda_{1}=\gamma\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}$

repeat for $t=1,2,\ldots,T_{\text{max}}$

  * 1.

Compute $\hat{\mu}_{t}$ solution of ($\mathcal{P}(\lambda)$) with
$\lambda=\lambda_{t}$ with warm start $\mu_{t}^{[0]}=\hat{\mu}_{t-1}$.

  * 2.

Compute $\sigma_{t}$ from the residual:

$\sigma_{t}=f_{y^{\delta},b}(\Phi\hat{\mu}_{t})$ (65)

  * 3.

if $\sigma_{t}<\sigma_{\text{target}}$

    * $\hat{\mu}_{t}$ is a solution $\Rightarrow$ stop

  * 4.

else if $\sigma_{t}\geq\sigma_{\text{target}}$

$\text{Update}\quad\lambda_{t+1}=\frac{\lambda_{t}\left\|\eta(\lambda_{t},\hat{\mu}_{t})\right\|_{\infty}}{c+1}$
(66)

until $\sigma_{t}<\sigma_{\text{target}}$

Algorithm 2 Homotopy algorithm in ${\mathcal{M}(\Omega)}$

In th following we follow [47], where Sliding Frank-Wolfe with homotopy is
proposed for the BLASSO problem ($L^{2}-|\cdot|$). We extend the strategy in
Algorithm 2 for a general fidelity term within an off-the-grid setting. Note
that each homotopy iteration $t\in\mathbb{N}$ performs the following steps:

  * •

Compute $\mu_{\lambda_{t}}:=\hat{\mu}_{t}$, solution of
($\mathcal{P}(\lambda)$) with $\lambda=\lambda_{t}$.

  * •

Check if
$\sigma_{t}:=f_{y^{\delta},b}(\Phi\mu_{\lambda_{t}})\leq\sigma_{\text{target}}(\delta)$,
where the target value $\sigma_{\text{target}}(\delta)$ depends on the noise
level $\delta>0$.

  * •

If the condition above is not met, decrease $\lambda$:
$\lambda_{t+1}<\lambda_{t}$.

In the following we discuss the choice of the starting value for the sequence
of $\lambda$’s, its updating rule and the choice of a suitable value
$\sigma_{\text{target}}$ for the noise scenario considered. In Table 1, we
outline the different choices made for Alg.2 to both ($L^{2}-|\cdot|$) and
($\tilde{\mathcal{D}}_{KL}-|\cdot|$).

| ($L^{2}-|\cdot|$) | ($\tilde{\mathcal{D}}_{KL}-|\cdot|$)
---|---|---
$\lambda_{1}$ := | $\gamma\|\Phi^{*}y^{\delta}\|_{\infty}$ | $\gamma\left\|\left(\Phi^{\ast}\big{(}\frac{y-b}{b}\big{)}\right)_{+}\right\|_{\infty}$
$\sigma_{t}$ := | $\frac{1}{2}\|\Phi\hat{\mu}_{t}-y^{\delta}\|^{2}$ | $\tilde{\mathcal{D}}_{KL}(\Phi\hat{\mu}_{t},y^{\delta})$
$\sigma_{\text{target}}(\delta)$ := | $\frac{1}{2}\|y-y^{\delta}\|^{2}=\frac{\delta^{2}}{2}$ | $\tilde{\mathcal{D}}_{KL}(y,y^{\delta})$
Table 1: Homotopy algorithmic choices for problems ($L^{2}-|\cdot|$) and
($\tilde{\mathcal{D}}_{KL}-|\cdot|$).

#### Starting value

The element $\eta(\lambda,\mu)$ is crucial for the definition of Algorithm 2
as it allows to define a good starting value $\lambda_{1}$. We propose to
initialise as follows:

$\displaystyle\lambda_{1}$
$\displaystyle:=\gamma\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}=\gamma\left\|\tilde{\eta}(\hat{\mu}_{0})\right\|_{\infty},\qquad\gamma\in(0,1)$
(67)

where $\hat{\mu}_{0}\in{\mathcal{M}(\Omega)}$ is the initialisation of the
solution and $\eta$ is the dual certificate (63). The parameter
$\gamma\in(0,1)$ is a relaxation parameter usually chosen close to 1. The
choice (67) is motivated by optimality conditions (64). If one takes at the
first iteration
$\lambda_{1}\geq\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}$, then
$\hat{\mu}_{0}$ is an optimal solution for ($\mathcal{P}(\lambda)$) with
$\lambda=\lambda_{1}$, since

$\left\|\eta(\lambda_{1},\hat{\mu}_{0})\right\|_{\infty}=\left\|\frac{1}{\lambda_{1}}\tilde{\eta}(\hat{\mu}_{0})\right\|_{\infty}=\frac{1}{\lambda_{1}}\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}\leq
1\quad\iff\quad\lambda_{1}\geq\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}.$
(68)

In this case, the algorithm does not improve upon the initialisation
$\hat{\mu}_{0}$ since it does not perform any iteration. On the contrary,
choosing $\lambda_{1}<\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}$ ensures
that the initial measure $\hat{\mu}_{0}$ is updated since the dual certificate
computed with respect to the initialisation $\hat{\mu}_{0}$ and
$\lambda_{1}>0$ is such that (67) has supremum norm that satisfies

$\displaystyle\left\|\eta(\lambda_{1},\hat{\mu}_{0})\right\|_{\infty}$
$\displaystyle=\frac{1}{\lambda_{1}}\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}=\frac{1}{\gamma\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}}\left\|\eta(1,\hat{\mu}_{0})\right\|_{\infty}=\frac{1}{\gamma}>1.$
(69)

#### Updating rule

The updating rule (66) for $\lambda$ together with the choice of a strictly
positive parameter $c>0$ ensures that the measure $\hat{\mu}_{t}$, which is
used to initialise $(\mathcal{P}_{\lambda_{t+1}})$ as
$\mu_{t+1}^{[0]}=\hat{\mu}_{t}$, is not an optimal solution for the problem.
Indeed, the dual certificate computed in correspondence with $\lambda_{t+1}$
and $\hat{\mu}_{t}$ reads

$\eta(\lambda_{t+1},\hat{\mu}_{t})=\frac{\lambda_{t}}{\lambda_{t+1}}\eta(\lambda_{t},\hat{\mu}_{t})=\frac{1+c}{\|\eta(\lambda_{t},\hat{\mu}_{t})\|_{\infty}}\eta(\lambda_{t},\hat{\mu}_{t})\;\Longrightarrow\;\|\eta(\lambda_{t+1},\hat{\mu}_{t})\|_{\infty}=1+c>1.$
(70)

Thus, $\mu_{t+1}^{[0]}=\hat{\mu}_{t}$ does not satisfy (64) for
($\mathcal{P}(\lambda)$) with $\lambda=\lambda_{t+1}$, so the homotopy step
$t+1$ computes a new candidate solution $\hat{\mu}_{t+1}$. This is indeed
consistent with the choice of rejecting $\hat{\mu}_{t}$ at the previous step
$t$.

###  4.1 Descent property

In this section, we show that the homotopy algorithm 2 produces a strictly
decreasing sequence of residual errors $(\sigma_{t})_{t}$. This properties
gives insight on the good convergence of the algorithm. We remark that the
proof follows [47, Proposition 3.1] but is here adapted to more general
choices of the data term $f_{y^{\delta},b}$.

######  Proposition 2.

The homotopy algorithm (Alg.2) for the resolution of ($\mathcal{P}(\lambda)$)
produces a strictly decreasing sequence of residual distances
$(\sigma_{t})_{t\in\mathbb{N}}$ provided that $f_{y^{\delta},b}(\cdot)$ is
twice differentiable and strictly convex.

###### Proof.

We consider

$\mu^{\ast}(\lambda)\in\operatornamewithlimits{argmin}_{\mu}\leavevmode\nobreak\
C(\mu;y^{\delta},\lambda):=f_{y^{\delta},b}(\Phi\mu)+\lambda|\mu|(\Omega),$
(71)

where the $\operatornamewithlimits{argmin}$ is computed over
${\mathcal{M}(\Omega)}$ if $\alpha=0$ in ($\mathcal{P}(\lambda)$) and over
${\mathcal{M}^{+}(\Omega)}$ if $\alpha=1$. Denoting by $\partial/\partial\mu$
the Radon-Nikodym derivative w.r.t. $\mu$, we have that the minimiser
$\mu^{\ast}(\lambda)$ satisfies

$\frac{\partial
C(\mu;y^{\delta},\lambda)}{\partial\mu}\Big{|}_{\mu=\mu^{\ast}}=\Phi^{\ast}\partial
f_{y^{\delta},b}(\Phi\mu^{*})+\lambda=0,$ (72)

Deriving the expression above with respect to $\lambda$ and applying the chain
rule yields:

$\displaystyle 0$
$\displaystyle=\frac{\partial}{\partial\lambda}\Big{(}\frac{\partial
C(y,\lambda,\mu)}{\partial\mu}\Big{|}_{\mu=\mu^{*}}\Big{)}=\frac{\partial}{\partial\lambda}\Big{(}\Phi^{\ast}\partial
f_{y^{\delta},b}(\Phi\mu^{*})+\lambda\Big{)}$ (73) $\displaystyle 0$
$\displaystyle=\frac{\partial}{\partial\mu}\Big{(}\Phi^{\ast}\partial
f_{y^{\delta},b}(\Phi\mu)\Big{)}\Big{|}_{\mu=\mu^{\ast}}\cdot\frac{\partial\mu^{\ast}(\lambda)}{\partial\lambda}+1$
(74) $\displaystyle 0$
$\displaystyle=\Phi^{\ast}\Phi\cdot\partial^{2}f_{y^{\delta},b}(\Phi\mu^{*})\cdot\frac{\partial\mu^{\ast}(\lambda)}{\partial\lambda}+1.$
(75)

Since $\Phi$ is positive definite and $f_{y^{\delta},b}(\cdot)$ is stricly
convex, from the computations above it follows that
$\frac{\partial\mu^{\ast}(\lambda)}{\partial\lambda}<0$. Consider now the
following:

$\frac{\partial}{\partial\lambda}f_{y^{\delta},b}\big{(}\Phi\mu^{\ast}(\lambda)\big{)}=\Phi^{\ast}\partial
f_{y^{\delta},b}(\Phi\mu^{*})\cdot\frac{\partial\mu^{\ast}(\lambda)}{\partial\lambda}.$
(76)

Recalling (72), we notice that
$\frac{\partial}{\partial\mu}f_{y^{\delta},b}(\Phi\mu^{*})=-\lambda$.
Moreover, since $\frac{\partial\mu^{\ast}(\lambda)}{\partial\lambda}<0$, then
$\frac{\partial}{\partial\lambda}f_{y^{\delta},b}\big{(}\Phi\mu^{\ast}(\lambda)\big{)}>0$.
Since by homotopy $\lambda_{t+1}<\lambda_{t}$, and denoting by
$\hat{\mu}_{t}:=\mu^{\ast}(\lambda_{t})$, we have

$\sigma_{t}=f_{y^{\delta},b}(\Phi\hat{\mu}_{t})<\sigma_{t+1}=f_{y^{\delta},b}(\Phi\hat{\mu}_{t+1}),$

which concludes the proof. ∎

######  Remark 2.

For BLASSO ($L^{2}-|\cdot|$), the fidelity
$f_{y^{\delta},b}(\Phi\mu)=\frac{1}{2}\|\Phi\mu-y\|^{2}$ satisfies the
hypothesis required to prove Proposition 2. Indeed,

$\partial^{2}f_{y^{\delta},b}=\partial^{2}\left(\frac{1}{2}\|\cdot-y^{\delta}\|^{2}\right)=1,$
(77)

that was explicitly used in [47, Proposition 3.1].

For problem ($\tilde{\mathcal{D}}_{KL}-|\cdot|$) the fidelity
$f_{y^{\delta},b}(\Phi\mu)=\tilde{\mathcal{D}}_{KL}(\Phi\mu+b,y)$ is twice
differentiable with

$\partial^{2}f_{y^{\delta},b}=\partial\tilde{\mathcal{D}}_{KL}(\cdot+b,y)=\frac{y^{\delta}}{\|\cdot+b\|^{2}}.$
(78)

######  Remark 3.

The proof of Proposition 2 holds for any strictly decreasing updating rule.
However, using an updating rule different from (66) there is no guarantee that
at iteration $t+1$ of Algorithm 2 the measure $\mu_{t+1}^{[0]}=\hat{\mu}_{t}$
is not already optimal for $\mathcal{P}(\lambda_{t+1})$, thus requiring an
immediate update of $\lambda$.

##  5 Numerical tests

In this section, we report numerical results to
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) computed by means of the Sliding Frank-
Wolfe (SFW) Algorithm 1 with homotopy (Algorithm 2) for several off-the-grid
sparse deconvolution problems in simulated 1D/2D and real 3D fluorescence
microscopy data.

###  5.1 1D comparisons between Gaussian and Poisson modelling

Figure 4: 1D comparison between Gaussian (left) and Poisson (right) models. In
black: ground truth spikes. In green: reconstructed spikes. For both models,
$\lambda=8.82$.

The aim of this first set of experiments on simulated 1D blurred signals
corrupted with Poisson noise is to validate our model
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) and to compare its performance with
BLASSO ($L^{2}-|\cdot|$). In a discrete setting, several works (e.g., [17, 50,
51] and [52] for examples in microscopy) have proposed both analytical and
numerical approaches for precisely modelling signal-dependent Poisson noise
using discrete version of the KL divergence (18). Generally speaking, such
choices improves performance (in terms of, e.g., localisation/reconstruction
quality, especially in low-photon count regimes) in comparison to simpler
$L^{2}$ Gaussian models, but only slightly. This is due to the biases
introduced by the handcrafted regularisation employed, such as, e.g.,
$\ell_{1}$ or TV-type. Aiming at reconstructing weighted sums of Diracs, for
which the regulariser (2) is tailored, we wonder in the following experiments
whether in such infinite-dimensional setting the improvement is more evident.

For that, we test 10 different ground truth signals with 6 randomly located
spikes in the 1D domain $\Omega=[0,1]$. For each ground truth signal, the
position of each spike is sampled from a uniform distribution over $\Omega$,
as well as their amplitudes from a uniform distribution over $[1-d,1+d]$ with
$d=0.4$. The corresponding acquisitions are blurred by a Gaussian 1D PSF with
$\sigma=0.07$, a spatially constant background $b=0.01$ is considered, and
then several Poisson noise realisations are generated as acquired data. In
Figure 4, one simulated ground truth signal $\mu_{\text{gt}}$ is shown (black
spikes) together the corresponding Poisson noisy and blurred data (blue
signal). All 1D signals are then reconstructed by using both the Poisson
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) and the Gaussian noise model BLASSO
($L^{2}-|\cdot|$) with $\lambda\in(0,10]$ using Algorithm 1. Figure 4 shows an
example between the two reconstructions (BLASSO, left, Poisson model, right)
$\mu_{\text{rec}}$ (green spikes) for $\lambda=8.82$. For such illustrative
example, the Poisson model provides a better estimate (in terms both of
amplitudes and localisation) than the Gaussian model.

To assess such observation over different choices of regularisation parameters
(to which the models could be sensitive), we then performed a statistical test
comparing localisation/reconstruction performance for all the generated
signals. To evaluate the goodness of the reconstructions, we consider the
Jaccard index defined in terms of the number of True Positive (TP), False
Positive (FP) and False Negative (FN) spikes as follows

$\text{Jac}_{\delta}(\mu_{\text{gt}},\mu_{\text{rec}})=\frac{\\#\text{TP}}{\\#\text{TP}+\\#\text{FP}+\\#\text{FN}}\in[0,1]$
(79)

with tolerance radius $\delta>0$. TP are reconstructed spikes that are at a
distance less than $\delta$ from a ground truth spile, while reconstructed
spikes that are more than $\delta$ distant from each ground truth spike are
denoted by FP. FN are spikes in the ground truth which have not been
associated to any TP.

---
(a) TP, FN and TP
(b) Jaccard index and RMSE of amplitudes and positions
Figure 5: Mean values over 100 different randomly generated ground truth
signals with 6 spikes and their corresponding reconstructions. Shaded area
corresponds to standard deviation. Maximum number of iterations of SFW:
$2N_{\text{spikes}}$. Tolerance radius for computation of the Jaccard Index is
$\delta=0.05$.

Similarly as in [12], we used also the RMSE of the amplitudes $a$ and
positions $x$ of the TP spikes as a different quality metric:

$\displaystyle\text{RMSE}_{x}(\mu_{\text{gt}},\mu_{\text{rec}})=\sqrt{\frac{1}{\\#\text{TP}}\sum_{i\in\text{TP}}\big{(}(x_{\text{rec}})_{i}-(x_{gt})_{i}\big{)}^{2}}$
(80)
$\displaystyle\text{RMSE}_{a}(\mu_{\text{gt}},\mu_{\text{rec}})=\sqrt{\frac{1}{\\#\text{TP}}\sum_{i\in\text{TP}}\big{(}(a_{\text{rec}})_{i}-(a_{gt})_{i}\big{)}^{2}}.$
(81)

In Figure 5(a), we plot the number of TP, FN and FP reconstructed by models
($L^{2}-|\cdot|$) and ($\tilde{\mathcal{D}}_{KL}-|\cdot|$) for a large finely
discretised range of $\lambda$. In Figure 5(b) the Jaccard index (computed
with $\delta=0.05$) and the Root Mean Square Error (RMSE) of amplitudes and
positions are also reported. The proposed Poisson model
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) has a better performance in terms of TP
and FN and, overall, in terms of the Jaccard index and RMSE of amplitudes and
positions. Only for the number of FP, BLASSO ($L^{2}-|\cdot|$) results
slightly better than ($\tilde{\mathcal{D}}_{KL}-|\cdot|$) for small values of
$\lambda$. This is due to the fact that ($\tilde{\mathcal{D}}_{KL}-|\cdot|$)
usually requires more iterations of SFW before reaching convergence. This
results in a better estimation of the number of molecules, with TP being
closer to the actual number of spikes, which may cause an overestimation of
the number of spikes with a consequently higher value of FP.

Using the same dataset, we then compare the results obtained running the
homotopy Algorithm 2 for the automatic selection of the regularisation
parameter $\lambda$ for both problems ($L^{2}-|\cdot|$) and
($\tilde{\mathcal{D}}_{KL}-|\cdot|$) solved by Algorithm 1 as inner solver,
with the algorithmic parameters specified in Table 2. To reduce the
computational burdens, we observed that only one iteration of SFW was enough,
as the estimated measures is anyway updated in the next homotopy step. For the
same reason, we set the maximum number of homotopy outer iterations to be
equal to twice the number of peaks in the ground truth. As far as
$\sigma_{\text{target}}$ is concerned, being in a simulated environment, we
computed the exact value of the residual $\sigma_{\text{exact}}$ as
$f_{y^{\delta},b}(\Phi\mu_{\text{gt}})$, to be compared with the residual
$\sigma_{t}$ at current iteration $t$ of Algorithm 2. Since
$\sigma_{\text{exact}}$ is unknown in real situations, a possible strategy for
its estimation will be discussed in the next section.

Parameters | $L^{2}-|\mu|$ | $\tilde{\mathcal{D}}_{KL}-|\mu|$
---|---|---
Max. number of homotopy iterations | $2N_{\text{molecules}}$ | $2N_{\text{molecules}}$
Max. number of inner SFW iterations | 1 | 1
Homotopy parameter $c$ | 15 | 40
Homotopy parameter $\gamma$ | 0.9 | 0.9
Choice of $\sigma_{\text{target}}$ | $1.5\times\frac{1}{2}\|\Phi\mu_{\text{gt}}+b-y\|^{2}$ | $1.5\times\mathcal{D}_{KL}(\Phi\mu_{\text{gt}}+b,y)$
Table 2: Parameters used by the homotopy algorithm (Alg.2) in the 1D simulated comparison tests. | $L^{2}-|\mu|$ | $\tilde{\mathcal{D}}_{KL}-|\mu|$
---|---|---
Jaccard index | 0.74 | 0.76
Number of TP | 4.50 | 4.80
Number of FN | 1.50 | 1.20
Number of FP | 0.10 | 0.40
RMSE on amplitudes of TP | 0.41 | 0.44
RMSE on positions of TP | 0.014 | 0.015
Final estimated $\mathbf{\lambda}$ | 6.09 | 40.21
Number of homotopy iterations | 4.55 | 3.93
Value of $\sigma_{\text{target}}$ | 4.09 | 77.16
Table 3: Homotopy algorithm: comparison between BLASSO and the Poisson off-
the-grid modelling. Mean values over 100 different randomly generated ground
truths with 6 spikes.

In Table 3, we report the values of TP,FN,FP, Jaccard index and RMSE of the
reconstructed signals for both models. The final estimated $\lambda$ is also
reported in Table 3, together with the number of performed homotopy iterations
and of the value $\sigma_{\text{target}}$. Using homotopy, we observe that we
retrieve values which are comparable with the best ones obtained using SFW
with grid search. This shows the effectiveness of the homotopy strategy.
Overall, the algorithm applied to solve ($\tilde{\mathcal{D}}_{KL}-|\cdot|$)
yields better results than ($L^{2}-|\cdot|$) in the presence of Poisson noise,
with a reduction of the number of FN and an improvement of the accuracy in
terms of Jaccard index.

###  5.2 2D examples: choice of $\sigma_{\text{target}}$

To avoid the choice of $\sigma_{\text{target}}$ to depend on the ground truth
image, we propose in this section an heuristic strategy to estimate a
reasonable value $\sigma_{\text{target}}$ relying on the sole acquisition
$y^{\delta}$ and on the assumption that the signal is sparse. To better
illustrate the proposed strategy we consider a 2D example of simulated blurred
and noisy microscopy acquisitions on the domain $\Omega=[0,1]^{2}$ 111Blurring
operator: 2D Gaussian PSF with $\sigma=0.07$. Spatially constant background
$b=0.05$., see Fig. 6(a).

Under a suitable sparsity level of the ground truth image, it is safe to
assume that its corresponding noisy and blurred acquisition $y^{\delta}$
presents regions containing only background noise, which we denote by
$\Omega_{\text{bg}}\subset\Omega$. In Figure 6(a) we show $y^{\delta}$ and
highlight $y^{\delta}|_{\Omega_{\text{bg}}}$ in transparency, i.e. the
acquisition "masked" to the area of background noise in the external square-
ring. We propose to estimate the value of $\sigma_{\text{target}}$ as follows

$\sigma_{\text{target}}=f_{y^{\delta},b}(0)\big{|}_{\Omega_{\text{bg}}}\frac{|\Omega|}{|\Omega_{\text{bg}}|},$
(82)

where the restriction of the fidelity term to $\mu=0$ is due to the fact we
assume the desired image $\mu$ to be null in $\Omega_{\text{bg}}$, i.e.
$\mu|_{\Omega_{\text{bg}}}=0$. Note that considering
$\sigma_{\text{target}}=f_{y^{\delta},b}(0)\big{|}_{\Omega_{\text{bg}}}$ would
be equivalent to assume that the noise is null on
$\Omega\setminus\Omega_{\text{bg}}$, which is obviously not true. The formula
(82) is thus adjusted to account for noise not only on $\Omega_{\text{bg}}$
but on the whole domain $\Omega$. Note that it is also possible to
approximately estimate the (constant) background $b\in L^{2}(\Omega)$ in
$\Omega_{\text{bg}}$ by taking

$b=\frac{1}{|\Omega_{\text{bg}}|}\int_{\Omega_{\text{bg}}}y^{\delta}(t)\mathrm{d}t.$
(83)

Theoretical value (based on the ground truth) | $\tilde{\mathcal{D}}_{KL}(\Phi\mu_{gt}+b,y^{\delta})$ | 842.3
---|---|---
Poisson discrepancy principle (84) [30] | $\frac{|\Omega|}{2}$ | 8192
Estimation based only on $y^{\delta}$ and $\Omega_{\text{bg}}$ (82) | $\tilde{\mathcal{D}}_{KL}(b,y_{bg})\frac{|\Omega|}{|\Omega_{\text{bg}}|}$ | 846.4
Table 4: Different estimates of $\sigma_{\text{target}}$.

(a) Data

(b) Reconstruction

(c) Values of $\lambda_{t}$, $\sigma_{\text{target}}$ and cost functional
along the homotopy iterations.

Figure 6: 2D fluorescence microscopy simulated image. (a) 2D sparse image
(white crosses denote off-the-grid locations) and its corresponding noisy
blurred acquisition $y^{\delta}$ with $y^{\delta}|_{\Omega_{\text{bg}}}$
visualised in trasparency. (b) 2D reconstruction (green spikes) obtained with
homotopy algorithm (Alg. 2) compared with the ground truth spikes (black).

For the image shown in Figure 6(a), we compare in Table 4 different choices of
$\sigma_{\text{target}}$ for ($\tilde{\mathcal{D}}_{KL}-|\cdot|$) with
$f_{y^{\delta},b}$ being the Kullback-Leibler fidelity term. In particular, we
consider in the second row the estimate proposed in [30], where a discrepancy
principle for Poisson data is studied under the following approximation

$\tilde{\mathcal{D}}_{KL}(\Phi\mu_{gt}+b,y^{\delta})\approx\frac{|\Omega|}{2}.$
(84)

This value is obtained by computing the expected value for Kullback-Leibler
fidelity and by approximating it with a first order Taylor expansion. As
observed also in [30], the estimate (84) might not be optimal and, indeed, one
should consider

$\tilde{\mathcal{D}}_{KL}(\Phi\mu_{gt}+b,y^{\delta})\approx\frac{1-\epsilon}{2}|\Omega|,$
(85)

where $\epsilon\in\mathbb{R}$ is small. When $\Omega$ is big, (85) might lead
to bad estimates even if $\epsilon$ is very small. As showed in Table 4, the
approximation (84) does not lead indeed to an accurate estimation of
$\sigma_{\text{target}}$. On the contrary, the estimation (82) proposed (last
row), is close to the real value (first row) which is known given the
simulated setting. We remark that in recent work [53] similar masking
strategies are used to define parameter selection strategies for variational
noise in Poisson scenarios, with a detailed description of the modifications
arising to the underlying statistical laws when performing such choice. The
application of analogous strategies to the problem considered is an
interesting venue for future work.

By using the homotopy algorithm (Alg.2) with SFW (Alg. 1) 222Max. number of
(outer) homotopy it. 20, max. number of (inner) SFW it. 1, $c=30$,
$\gamma=0.9$. for the reconstruction of the data in Figure 6(a) with
background and $\sigma_{\text{target}}$ estimated by (83) and (82),
respectively, we obtain the results shown in Figure 6.

###  5.3 Real 3D wide-field dataset

Figure 7: ERES 3D real data (I and II columns: lateral views, III column: top
view).

We finally consider a real 3D fluorescence microscopy setting with blurred and
noisy volume acquisitions. The image was acquired in wide-field microscopy
with a 100x/1.49 objective, a Hamamatsu Flash 4 camera with 6.5um pixels, by
Alejandro Melero Carrillo at the MRC-LMB and was used in [54]. It is an
acquisition of yeasts expressing fluorescent proteins (*SEC16-sfGFP* and a
*SEC24-sfGFP*) localised at the Endoplasmic Reticulum exit sites (ERES). The
acquired volume $y^{\delta}$ is shown in Figure 7 with maximum intensities
projections over the $xz$, $yz$, $yx$ planes. The 3D volume blurred and noisy
acquisition has $190\times 190\times 17$ voxels with voxel size of 65nm in
$yz$ and 250nm in $z$. Signal dependency of the noise is observed.

To reconstruct a sparse volume from the 3D acquisition $y^{\delta}$, we use
the homotopy algorithm (Alg. 2) with an accelerated version of the SFW, called
Boosted SFW proposed in [47], as an inner solver to minimise
($\tilde{\mathcal{D}}_{KL}-|\cdot|$). Boosted SFW reduces computational costs
by limiting the number of sliding steps.

Parameters | $\eqref{poissonblasso}$
---|---
Max. number of homotopy iterations | 10
Max. number of inner SFW iterations | 50
Homotopy parameter $c$ | $0.5$
$\sigma_{\text{target}}$ given by (82) | 1102067.75
$\sigma_{\text{target}}$ given by (84) | 306850
Constant background estimate (83) | b = 337.77
Table 5: Parameters used in Algorithm 2 for the reconstruction of the 3D
volume

We consider a 3D convolution kernel (7) estimated as a 3D Gaussian PSF,

$\varphi\big{(}(x,y,z)\big{)}=\frac{1}{\sqrt{(2\pi)^{3}}\sigma_{x}\sigma_{y}\sigma_{z}}\exp\left[-\frac{x^{2}}{2\sigma_{x}^{2}}\right]\exp\left[-\frac{y^{2}}{2\sigma_{y}^{2}}\right]\exp\left[-\frac{z^{2}}{2\sigma_{z}^{2}}\right].$

The standard deviation $\sigma_{x},\sigma_{y}$ of the 3D PSF can be estimated
from the Full Width Half Maximum, which is given by
$\text{FWHM}=0.61\lambda_{\text{wavelength}}/\text{NA}$, where
$\lambda_{\text{wavelength}}$ is the emission wavelength of the fluorescent
proteins and NA is the numerical aperture of the microscope. Note that if the
FWHM is known, then it is possible to retrieve information about the variance
parameters of the PSF since $\text{FWHM}=2.355\cdot\sigma_{x}$ and
$\text{FWHM}=2.355\cdot\sigma_{y}$. For the standard deviation in the
$z$-axis, we assume $\sigma_{z}=2\cdot\sigma_{x}$. Since the value of the
numerical aperture is known, $\text{NA}=1.49$, and
$\lambda_{\text{wavelength}}=508$nm for the green fluorescent proteins under
test, we obtain a PSF estimation with $\sigma_{x}=\sigma_{y}=89$nm and
$\sigma_{z}=178$nm, which appears to be a good approximation of the underlying
blur model.

Figure 8: ERES 3D data. Values of $\lambda_{t}$, $\sigma_{\text{target}}$ and
cost functional along the homotopy iterations.

The parameters used to run the homotopy algorithm 2 with BSFW as inner solver
are reported in Table 5. Note that as far as the estimate of
$\sigma_{\text{target}}$ is concerned, the formulas (84) and (82) give very
different results, so that the estimate given by (82) was considered as more
accurate as showed in the previous section. The background is estimated by
(83) as $b=337.77$. We computed this reconstruction using 50 homotopy
iterations fixing to 10 the maximum number of inner loops of the BSFW
algorithm, using Google Colab CPUs for about 10 hours. The reconstructed
volume $\mu_{\text{rec}}$ counted 274 spikes. While probably underestimating
the exact result this first result is promising, since the use of Algorithm 2
yields very precise localisation of spikes automatically, with no need of
estimating the regularisation parameter, and no a-priori information about the
data. A better visualisation of the reconstruction under different views is
given in Figure 9 using the visualisation codes used in [12] and provided by
the authors at their GitHub page333https://github.com/qdenoyelle/sfw4blasso.

|  |
---|---|---
Figure 9: Sparse reconstruction of the 3D real ERES data

##  6 Conclusions

In this work, we considered sparse inverse problems in the off-the-grid
setting of the space of radon measures ${\mathcal{M}(\Omega)}$ under the
assumption of signal-dependent Poisson noise in the measurement. Such choice
is motivated by fluorescence microscopy applications, where the noise observed
is modelled as Poisson to account for photon-counting processes. We proposed a
variational approach where the Total Variation norm is coupled with the
Kullback-Leibler fidelity term and a non-negativity constraint. For the
resulting composite problem we derive analytically the optimality and
extremality conditions. Using the Sliding Frank-Wolfe algorithm as numerical
solver, we then discussed the choice of a good regularisation parameter by
means of an algorithmic homotopy strategy. Finally, we presented several
numerical experiments on simulated 1D/2D data to validate the proposed Poisson
model, to compare its performance with the analog Gaussian model and to
confirm the good performance of the homotopy strategy proposed. To conclude,
we tested the proposed framework on a 3D real fluorescence microscopy dataset.

Acknowledgements We thank B. Laville and L. Blanc-Féraud (Morpheme, Inria),
and R. Petit (MaLGa, University of Genoa) for thoughtful discussions and for
providing useful codes of the Sliding Frank-Wolfe algorithm for the resolution
of 1D and 2D BLASSO results. We thank J. Boulanger (MRC-LMB, Cambrisge) for
making this work possible by having put ML, LC and CE in touch with AM, who
acquired the 3D data and for thoughtful discussions. Funding LC and ML
acknowledge the support received by the projects ANR MICROBLIND
(ANR-21-CE48-0008) and ANR JCJC TASKABILE (ANR-22-CE48-0010). CE acknowledges
the support received by the MUR Excellence Department Project awarded to
Dipartimento di Matematica, Università di Genova, CUP D33C23001110001. CE and
ML acknowledge the support of the "PRIN 2022ANC8HL - Inverse Problems in the
Imaging Sciences (IPIS)" project, granted by the Italian Ministero
dell’Università e della Ricerca within the framework of the European Union -
Next Generation EU program and of the Italian INdAM group on scientific
calculus GNCS.

Data availability

The codes used for implementing the models and algorithms described in this
work are available at https://github.com/martalazzaretti/KL-TV-off-the-grid.

## Declarations

Conflict of interest The authors declare that they have no conflict of
interest.

##  Appendix A Technical Results

We report here the well-known definition of convex conjugate and some of its
properties.

######  Definition 3.

Given $f:\mathcal{X}\longrightarrow\mathbb{R}$, with $\mathcal{X}$ real Banach
space, the convex conjugate of $f$ is the function
$f^{*}:\mathcal{X}^{*}\longrightarrow\mathbb{R}\cup\\{+\infty\\}$ defined by

$f^{*}(x^{*}):=\sup_{x\in\mathcal{X}}\leavevmode\nobreak\ \langle
x^{*},x\rangle_{\mathcal{X}^{*}\times\mathcal{X}}-f(x),\quad
x^{*}\in\mathcal{X}^{*}.$ (86)

It is possible to obtain the convex conjugate of the sum of functions as
infimal convolution of the convex conjugate of each singular function [36,
37].

######  Proposition 3.

Let $f_{1},\ldots,f_{n}:\mathcal{X}\rightarrow\mathbb{R}$ with $\mathcal{X}$
real Banach space. Suppose there exists at least a point $x\in\mathcal{X}$
such that $f_{1},\ldots,f_{n}$ are continuous in $x$. Then,

$\left(f_{1}+\ldots+f_{n}\right)^{*}(x^{*})=\min_{x_{1}^{*}+\ldots+x_{n}^{*}=x^{*}}\
f_{1}^{*}(x_{1}^{*})+\ldots+f_{n}^{*}(x_{n}^{*}).$ (87)

The convex conjugate of any norm is a well-known result. We report here for
the sake of completeness, the convex conjugate of the TV norm.

######  Lemma 5.

Let $A:{\mathcal{M}(\Omega)}\rightarrow\mathbb{R}$,
$A(\cdot)=\left|\cdot\right|(\Omega)$ be the TV-norm. The convex conjugate of
$A$ is the function
$A^{*}:\mathscr{C}_{0}(\Omega,\mathbb{R})\rightarrow\mathbb{R}\cup\\{+\infty\\}$
given by

$A^{*}(\psi)=\begin{cases}0&\|\psi\|_{\infty,\Omega}\leq 1\\\
+\infty&\|\psi\|_{\infty,\Omega}>1\end{cases}\,.$ (88)

###### Proof.

By Definition 86 of convex conjugate, for any
$\psi\in\mathscr{C}_{0}(\Omega,\mathbb{R})$ we write

$\displaystyle A^{*}(\psi)$
$\displaystyle=\sup_{\mu\in{\mathcal{M}(\Omega)}}\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}-|\mu|(\Omega)\geq\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}-|\mu|(\Omega),\quad\forall\mu\in{\mathcal{M}(\Omega)}.$
(89)

Taking $\mu=\alpha\delta_{x}$ with $x\in\Omega$ and $\alpha>0$, we have

$\sup_{\mu\in{\mathcal{M}(\Omega)}}\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}-|\mu|(\Omega)\geq\alpha(\psi(x)-1),$

and the limit for $\alpha\rightarrow\infty$ of the latter inequality yields
$A^{*}(\psi)\geq+\infty$ if $\psi(x)>1$. A similar result for $\psi(x)<1$ is
obtained with the measure $\mu=-\alpha\delta_{x}$, with $\alpha>0$. Thus, we
have $A^{*}(\psi)=+\infty$ if $\|\psi\|_{\infty,\Omega}>1$.

Assume now that $\|\psi\|_{\infty,\Omega}\leq 1$. First notice that

$A^{*}(\psi)=\sup_{\mu\in{\mathcal{M}(\Omega)}}\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}-|\mu|(\Omega)\geq
0$ (90)

by considering $\mu=0$. Moreover,

$\displaystyle\langle\psi,\mu\rangle_{\mathscr{C}_{0}(\Omega,\mathbb{R})\times{\mathcal{M}(\Omega)}}-|\mu|(\Omega)$
$\displaystyle\leq\|\psi\|_{\infty,\Omega}|\mu|(\Omega)-|\mu|(\Omega)\leq|\mu|(\Omega)\left(\|\psi\|_{\infty,\Omega}-1\right)\leq
0$ (91)

since $\|\psi\|_{\infty,\Omega}\leq 1$. By evaluating the sup on both sides of
the last inequality, since $A(\psi)\geq 0$ one finally gets $A^{*}(\psi)=0$ if
$\|\psi\|_{\infty,\Omega}\leq 1$, which concludes the proof. ∎

######  Lemma 6.

_(Fenchel-Rockafellar duality.)_ Let $V,Y$ be two Banach spaces. For
$\Lambda:V\rightarrow Y$ linear operator, $F:V\rightarrow\mathbb{R}$ and
$G:Y\rightarrow\mathbb{R}$ convex functionals, we consider the following
primal problem:

$\operatornamewithlimits{argmin}_{u\in V}\leavevmode\nobreak\ F(u)+G(\Lambda
u).$ (92)

The corresponding dual problem reads

$\operatornamewithlimits{argmax}_{p^{*}\in
Y^{*}}-F^{*}(\Lambda^{*}p^{*})-G^{*}(-p^{*}),$ (93)

where $\Lambda^{*}:Y^{*}\rightarrow V^{*}$ is the adjoint operator of
$\Lambda$ and $F^{*}:V^{*}\rightarrow\mathbb{R}\cup\\{+\infty\\}$,
$G^{*}:Y^{*}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ are the convex conjugate
of $F$ and $G$.

Moreover, if $u\in V$ and $p^{*}\in Y^{*}$ are respectively solutions of the
primal (92) and dual (93) problems, the following extremality conditions hold:

$\begin{cases}&\Lambda^{*}p^{*}\in\partial F(u)\\\ &-p^{*}\in\partial
G(\Lambda u)\end{cases}.$ (94)

## References

  *   * Duval and Peyré [2017a] Duval, V., Peyré, G.: Sparse spikes super-resolution on thin grids II: the continuous basis pursuit. Inverse Probl. 33(9), 095008 (2017)
  * Duval and Peyré [2017b] Duval, V., Peyré, G.: Sparse regularization on thin grids I: the Lasso. Inverse Probl. 33(5), 055008 (2017)
  * Laville et al. [2021] Laville, B., Blanc-Féraud, L., Aubert, G.: Off-the-grid variational sparse spike recovery: Methods and algorithms. J. Imaging. 7(12) (2021)
  * Scherzer and Walch [2009] Scherzer, O., Walch, B.: Sparsity regularization for radon measures. In: Proc.-SSVM 2009, pp. 452–463. Springer, Berlin Heidelberg (2009)
  * Bredies and Pikkarainen [2013] Bredies, K., Pikkarainen, H.K.: Inverse problems in spaces of measures. ESAIM: COCV. (2013)
  * De Castro and Gamboa [2012] De Castro, Y., Gamboa, F.: Exact reconstruction using Beurling minimal extrapolation. J. Math. Anal. Appl. (2012)
  * Fernandez-Granda [2013] Fernandez-Granda, C.: Support detection in super-resolution. arXiv preprint 1302.3921 (2013)
  * Duval and Peyré [2013] Duval, V., Peyré, G.: Exact support recovery for sparse spikes deconvolution. Found. Comput. Math. 15, 1315–1355 (2013)
  * Denoyelle et al. [2017] Denoyelle, Q., Duval, V., Peyré, G.: Support recovery for sparse super-resolution of positive measures. J. Fourier Anal. Appl. (2017)
  * Boyer et al. [2017] Boyer, C., De Castro, Y., Salmon, J.: Adapting to unknown noise level in sparse deconvolution. Inf. Inference J. IMA (2017)
  * Poon and Peyré [2019] Poon, C., Peyré, G.: Multi-dimensional sparse super-resolution. SIAM J. Math. Anal. (2019)
  * Denoyelle et al. [2019] Denoyelle, Q., Duval, V., Peyré, G., Soubies, E.: The sliding frank–wolfe algorithm and its application to super-resolution microscopy. Inverse Probl. 36(1), 014001 (2019)
  * Flinth et al. [2021] Flinth, A., Gournay, F., Weiss, P.: On the linear convergence rates of exchange and continuous methods for total variation minimization. Mathematical Programming 190(1), 221–257 (2021)
  * Bénard et al. [2024] Bénard, P.-J., Traonmilin, Y., Aujol, J.-F., Soubies, E.: Estimation of off-the grid sparse spikes with over-parametrized projected gradient descent: theory and application. Inverse Problems (2024)
  * Dossal et al. [2017] Dossal, C., Duval, V., Poon, C.: Sampling the fourier transform along radial lines. SIAM J. Numer. Anal. 55(6), 2540–2564 (2017)
  * De Castro et al. [2021] De Castro, Y., Gadat, S., C., M., Maugis-Rabusseau, C.: SuperMix: Sparse regularization for mixtures. Ann. Statist. 49(3), 1779–1809 (2021)
  * Bertero et al. [2018] Bertero, M., Boccacci, P., Ruggiero, V.: Inverse Imaging with Poisson Data. 2053-2563. IOP Publishing, UK (2018)
  * Poon [2019] Poon, C.: An Introduction to Sparse Spikes Recovery via the BLASSO. Lecture notes (2019)
  * Rudin [1987] Rudin, W.: Real and Complex Analysis, 3rd Ed. McGraw-Hill, Inc., USA (1987)
  * Cohn [2013] Cohn, D.L.: Measure Theory. Springer, New York (2013)
  * Chizat and Bach [2018] Chizat, L., Bach, F.: On the global convergence of gradient descent for over-parameterized models using optimal transport. In: Proc. 32nd Int. Conf. Neural Inf. Process. Syst., NIPS’18, pp. 3040–3050 (2018)
  * Koulouri et al. [2021] Koulouri, A., Heins, P., Burger, M.: Adaptive superresolution in deconvolution of sparse peaks. IEEE Trans. Signal Process. 69, 165–178 (2021)
  * Laville et al. [2023a] Laville, B., Blanc-Féraud, L., Aubert, G.: Off-the-grid curve reconstruction through divergence regularization: An extreme point result. SIAM J. Imaging Sci. 16(2), 867–885 (2023)
  * Laville et al. [2023b] Laville, B., Blanc-Féraud, L., Aubert, G.: Off-the-grid charge algorithm for curve reconstruction in inverse problems. In: Proc.-SSVM 2023, pp. 393–405. Springer, Cham (2023)
  * De Castro et al. [2023] De Castro, Y., Duval, V., Petit, R.: Towards off-the-grid algorithms for total variation regularized inverse problems. J. Math. Imaging Vis. 65(1), 53–81 (2023)
  * De Castro et al. [2021] De Castro, Y., Duval, V., Petit, R.: Towards off-the-grid algorithms for total variation regularized inverse problems. In: Proceedings SSVM 2021, pp. 553–564. Springer, Cham (2021)
  * De Castro et al. [2023] De Castro, Y., Duval, V., Petit, R.: Exact recovery of the support of piecewise constant images via total variation regularization. arXiv preprint: 2307.03709 (2023)
  * Beurling [1938] Beurling, A.: Sur les intégrales de fourier absolument convergentes et leur application à une transformation fonctionnelle. In: Proc. Ninth Scand. Math. Congr., Helsinki, Finland, pp. 345–366 (1938)
  * Bertero et al. [2009] Bertero, M., Boccacci, P., Desidera, G., Vicidomini, G.: Image deblurring with Poisson data: from cells to galaxies. Inverse Probl. 25, 123006 (2009)
  * Bertero et al. [2010] Bertero, M., Boccacci, P., Talenti, G., Zanella, R., Zanni, L.: A discrepancy principle for Poisson data. Inverse Probl. 26, 105004 (2010)
  * Li et al. [2015] Li, J., Shen, Z., Yin, R., Zhang, X.: A reweighted $l^{2}$ method for image restoration with Poisson and mixed Poisson-Gaussian noise. Inverse Probl. Imaging. 9(3), 875–894 (2015)
  * Sawatzky et al. [2013] Sawatzky, A., Brune, C., Kösters, T., Wübbeling, F., Burger, M.: EM-TV Methods for Inverse Problems with Poisson Noise, pp. 71–142. Springer, Cham (2013)
  * Lanza et al. [2014] Lanza, A., Morigi, S., Sgallari, F., Wen, Y.-W.: Image restoration with Poisson–Gaussian mixed noise. "Comput. Methods Biomech. Biomed. Eng. Imaging Vis. 2(1), 12–24 (2014)
  * Rockafellar [1970] Rockafellar, R.T.: Convex Analysis. Princeton University Press, Princeton, NJ, USA (1970)
  * Ekeland and Témam [1999] Ekeland, I., Témam, R.: Convex Analysis and Variational Problems. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA (1999)
  * Hiriart-Urruty [2006] Hiriart-Urruty, J.-B.: A note on the legendre-fenchel transform of convex composite functions. In: Alart, P., Maisonneuve, O., Rockafellar, R.T. (eds.) Nonsmooth Mechanics and Analysis, pp. 35–46. Springer, Boston, MA (2006)
  * Fajardo et al. [2012] Fajardo, M.D., Vicente-Pérez, J., Rodríguez, M.M.L.: Infimal convolution, c-subdifferentiability, and Fenchel duality in evenly convex optimization. Inverse Probl. 20(2), 375–396 (2012)
  * Ndiaye et al. [2017] Ndiaye, E., Fercoq, O., Gramfort, A., Salmon, J.: Gap safe screening rules for sparsity enforcing penalties. J. Mach. Learn. Res. 18(1), 4671–4703 (2017)
  * Wang and Ye [2014] Wang, J., Ye, J.: Two-layer feature reduction for sparse-group lasso via decomposition of convex sets. In: Proc. 27th Int. Conf. Neural Inf. Process. Syst. - Vol. 2, NIPS’14, pp. 2132–2140. MIT Press, Cambridge, MA, USA (2014)
  * Dantas et al. [2021] Dantas, C.F., Soubies, E., Févotte, C.: Safe screening for sparse regression with the Kullback-Leibler divergence. In: Proc. - ICASSP IEEE Int. Conf. Acoust., Toronto, Canada (2021)
  * Bauschke and Combettes [2011] Bauschke, H.H., Combettes, P.L.: Convex Analysis and Monotone Operator Theory in Hilbert Spaces, (2011)
  * Valkonen [2023] Valkonen, T.: Proximal methods for point source localisation. J. Nonlinear Anal. Optim. Volume 4 (2023)
  * Frank and Wolfe [1956] Frank, M., Wolfe, P.: An algorithm for quadratic programming. Nav. Res. Logist. Q. 3(1-2), 95–110 (1956)
  * Bredies et al. [2009] Bredies, K., Lorenz, D.A., Maass, P.: A generalized conditional gradient method and its connection to an iterative shrinkage method. Computational Optimization and Applications 42(2), 173–193 (2009)
  * Chizat [2022] Chizat, L.: Sparse optimization on measures with over-parameterized gradient descent. Math. Program. 194(1–2), 487–532 (2022)
  * Pokutta [2024] Pokutta, S.: The Frank-Wolfe Algorithm: A Short Introduction. Jahresbericht der Deutschen Mathematiker-Vereinigung 126(1), 3–35 (2024)
  * Courbot and Colicchio [2021] Courbot, J.-B., Colicchio, B.: A fast homotopy algorithm for gridless sparse recovery. Inverse Probl. 37(2), 025002 (2021)
  * Osborne et al. [2000a] Osborne, M., Presnell, B., Turlach, B.: A new approach to variable selection in least squares problems. IMA J. Numer. Anal. 20(3), 389–403 (2000)
  * Osborne et al. [2000b] Osborne, M.R., Presnell, B., Turlach, B.A.: On the LASSO and Its Dual. J. Comput. Graph. Stat. 9(2), 319–337 (2000)
  * Harmany et al. [2012] Harmany, Z.T., Marcia, R.F., Willett, R.M.: This is SPIRAL-TAP: Sparse Poisson Intensity Reconstruction ALgorithms—Theory and Practice. IEEE Trans. Image Process. 21(3), 1084–1096 (2012)
  * Calatroni et al. [2017] Calatroni, L., De Los Reyes, J.C., Schönlieb, C.-B.: Infimal convolution of data discrepancies for mixed noise removal. SIAM J. Imaging Sci. 10(3), 1196–1233 (2017)
  * Reiser et al. [2023] Reiser, N., Guo, M., Shroff, H., Riviere, P.J.L.: Phase Diverse Phase Retrieval for Microscopy: Comparison of Gaussian and Poisson Approaches. arXiv preprint 2308.00734 (2023)
  * Bevilacqua et al. [2023] Bevilacqua, F., Lanza, A., Pragliola, M., Sgallari, F.: Masked unbiased principles for parameter selection in variational image restoration under poisson noise. Inverse Probl. 39(3), 034002 (2023)
  * Gomez-Navarro et al. [2020] Gomez-Navarro, N., Melero, A., Li, X.-H., Boulanger, J., Kukulski, W., Miller, E.A.: Cargo crowding contributes to sorting stringency in COPII vesicles. J. Cell. Biol. 219(7), 201806038 (2020)