An entropic Landweber method for linear ill-posed problems

TL;DR

This paper introduces an entropic Landweber iterative method for regularizing linear ill-posed problems, providing a closed-form solution, convergence analysis, and numerical validation for reconstructing nonnegative unknowns and probability distributions.

Contribution

It presents a novel entropic projection approach with explicit iterates and analyzes its convergence, extending the Landweber method for specific applications.

Findings

Effective in reconstructing nonnegative unknowns

Converges reliably in numerical experiments

Relates to and extends existing methods

Abstract

The aim of this paper is to investigate the use of an entropic projection method for the iterative regularization of linear ill-posed problems. We derive a closed form solution for the iterates and analyze their convergence behaviour both in a case of reconstructing general nonnegative unknowns as well as for the sake of recovering probability distributions. Moreover, we discuss several variants of the algorithm and relations to other methods in the literature. The effectiveness of the approach is studied numerically in several examples.

Equations

$$u_{k+1}\in \text{mbox}arg\underset{u}{\min }\left\{\frac{1}{2}\Vert Au-y{\Vert }^2+cd(u,u_k)-\frac{1}{2}\Vert Au-A{u}_k{\Vert }^2\right\},k\in \mathbb{N},$$

$$u_{k+1}\in \text{mbox}arg\underset{u}{\min }\left\{⟨Au-y,A{u}_k-y⟩+cd(u,u_k)\right\},$$

$$f(u)={\int }_{\Omega }u(t)\ln u(t)dt.$$

$$\tilde{f}(u)={\int }_{\Omega }u(t)\ln u(t)-u(t)+1dt,$$

$${\int }_{\Omega }u(t)dt=1,$$

$$u_{k+1}\in \text{mbox}arg\underset{u}{\min }\left\{⟨Au-y,A{u}_k-y⟩+cd(u,u_k)+{\chi }_m(u)\right\},m\in \{0,1\},$$

$$\chi_{1}(u)=\left\{\begin{array}[]{ll}0&\mbox{if~{}}\int_{\Omega}u(t)~{}dt=1,\\ +\infty&\mbox{else},\end{array}\right.$$

$$f(u)=\left\{\begin{array}[]{ll}{\int_{\Omega}u(t)\ln u(t)\,dt}&\mbox{if $u\geq 0$ a.e. and $u\ln u\in{{L}}^{1}(\Omega)$},\\ &\\ {+\infty}&\mbox{otherwise.}\end{array}\right.$$

$$d(v,u)=f(v)-f(u)-f^{\prime }(u,v-u)$$

$$d(v,u)={\int }_{\Omega }\left[v(t)\ln \frac{v(t)}{u(t)}-v(t)+u(t)\right]dt$$

$$f^{\prime }(u,v)={\int }_{\Omega }v(t)[1+\ln u(t)]dt,$$

$$\Vert u-v{\Vert }_1^2\le \left(\frac{2}{3}\Vert v{\Vert }_1+\frac{4}{3}\Vert u{\Vert }_1\right)d(v,u).$$

$$\text{mbox}dom\partial f=\{u\in L^1(\Omega ):u\text{mbox}boundedandboundedawayfromzeroa.e.\}.$$

$$\{x\in L^1(\Omega ):d(x,u)\le C\}.$$

$$⟨u,v⟩={\int }_{\Omega }u(t)v(t)dt,$$

$$⟨\ell ,u⟩+d(u,v)+{\chi }_m(u)\to \underset{u\in \text{mbox}domf}{\min }$$

$$u_{m}=c_{m}ve^{-\ell},\qquad c_{j}=\left\{\begin{array}[]{ll}1&\mbox{if~{}}m=0,\\ \frac{1}{\int_{\Omega}ve^{-\ell}~{}dt}&\mbox{if~{}}m=1,\end{array}\right.$$

$$⟨\ell ,u⟩+d(u,v)+{\chi }_m(u)$$

$$={\int }_{\Omega }\left[u(t)\ln \frac{u(t)}{u_m(t)}-u(t)+v(t)+u(t)\ln c_m\right]dt+{\chi }_m(u)$$

$$=d(u,u_m)+\ln c_m\left({\int }_{\Omega }u(t)dt-1\right)+{\chi }_m(u)+C_m,$$

$$\ln c_m\left({\int }_{\Omega }u(t)dt-1\right)+{\chi }_m(u)={\chi }_m(u).$$

$$\Vert Au-Av\Vert \le \gamma \sqrt{d(u,v)}$$

$$\Vert Au-Av{\Vert }^2\le \Vert A{\Vert }^2\Vert u-v{\Vert }_{L^1(\Omega )}^2.$$

$$\Vert Au-Av{\Vert }^2\le \Vert A{\Vert }^{2}2d(u,v),$$

$$D(u,v)=cd(u,v)-\frac{1}{2}\Vert Au-Av{\Vert }^2,$$

$$u_{k+1}=u_{k}c_{k}^{m}e^{\lambda A^{*}(y-Au_{k})},\qquad c_{k}^{m}=\left\{\begin{array}[]{ll}1&\mbox{if~{}}m=0,\\ \frac{1}{\int_{\Omega}u_{k}e^{\lambda A^{*}(y-Au_{k})}~{}dt}&\mbox{if~{}}m=1,\end{array}\right.$$

$$\ln u_{k+1}=\ln u_k+\ln c_k^m+\lambda A^{\ast }(y-A{u}_k),$$

$$u_{k+1}\in \text{mbox}arg\underset{u}{\min }\left\{\frac{1}{2}\Vert Au-y{\Vert }^2+{\chi }_m(u)+D(u,u_k)\right\}.$$

$$v_k=w_0+\sum _{j=0}^{k-1}(y-A{u}_j)$$

$${\xi }_k={\xi }_0+\sum _{j=0}^{k-1}\ln c_j^m+\lambda A^{\ast }\left(\sum _{j=0}^{k-1}(y-A{u}_j)\right)=\sum _{j=0}^{k-1}\ln c_j^m+\lambda A^{\ast }{v}_k.$$

$$\frac{1}{2}\Vert A{u}_{k+1}-y{\Vert }^2+D(u_{k+1},u_k)\le \frac{1}{2}\Vert A{u}_k-y{\Vert }^2,$$

Full text

1 Introduction

This work deals with linear ill-posed equations $Au=y$ with $A:X\rightarrow Y$ acting between a Banach space $X$ and a Hilbert space $Y$, for which solutions with specific properties (such as positivity) are sought. In this respect, we consider iterative regularization methods of the following type

[TABLE]

where $d=D_{f}$ denotes the Bregman distance [8] associated with a convex functional $f:X\rightarrow\mathbb{R}\cup\{+\infty\}$ which is nonnegative and $c$ is some positive number. The term $d(u,u_{k})$ acts as a penalty enforcing the desired features for the solutions.

Note that we can rewrite (1) as

[TABLE]

which shows that the scheme can be obtained also by linearizing the quadratic data-fitting term $\displaystyle{\frac{1}{2}\|Au-y\|^{2}}$ at the current iterate $u_{k}$. This class of methods incorporates several procedures that have been proposed so far in the literature. For instance, the classical case when $f$ is quadratic in Hilbert spaces reduces to the Landweber method, as emphasized by [14] and as investigated for nonlinear operator equations by means of surrogate functionals in [25], see also the discussion in [28]. The case of quadratic $f$ in reflexive Banach spaces has been studied by [29]. The setting when $f$ is the total variation functional smoothed by a quadratic has been analyzed by [4], requiring fine analysis tools due to the bounded variation function space context. The case of $\ell^{1}$-penalties has been treated in [30, 10], resulting in the so-called linearized Bregman algorithm. In all those cases however, some quadratic term had to be part of $f$ to guarantee even well-definedness of the iterates and subsequently convergence.

We are interested here in the Shannon entropy setting without any quadratic term, i.e.

[TABLE]

We mention that one can alternatively consider a linear shift to

[TABLE]

which is a nonnegative functional inducing the same Bregman distance.

This raises challenges to analyze the problem in the $L^{1}$ setting without quadratic terms in the functional, but provides a simple closed iterative method with preserving the sign of the starting point (function) along the iterations. The latter formulation involves entropic projections, as one can see in the following section.

Moreover, we shall also be interested in the solution of inverse problems with unknowns being probability densities, i.e. we minimize on the domain of $f$ subject to the constraint

[TABLE]

which again results in a simple closed iterative form.

The advantages of using Bregman projections for solving variational problems with unknown probability densities have been exploited by several authors before, e.g. in optimal transport (cf. [6, 24]).

In order to write both problems in a closed form, we will use the equivalent formulation

[TABLE]

where $\chi_{0}\equiv 0$ denotes the original problem without integral constraint, and

[TABLE]

is employed for enforcing probability densities. The minimization is taken here over the domain of the entropy functional.

One finds the above entropy based algorithm in the finite dimensional optimization literature, as well as in the machine learning one, under quite different names. One could mention the mirror descent type algorithms for function minimization introduced in [22] and the Bregman-distance version with emphasis on entropy in [5], and the exponentiated gradient descent method for linear predictions - see [21]. The work [19] investigated three versions of the so-called approximate (linearized) proximal point methods for optimization in combination with line search strategies. The reader is referred to Sections 6.6 - 6.9 in [12] for other iterative optimization methods employing the Shannon entropy.

The main contribution of our work is the convergence of the iterates (4) to a solution of the equation $Au=y$ even in an infinite dimensional setting of such a nonquadratic penalty version, by stating also error estimates in the sense of a distance between the solution and the iterates, as opposed to the classical situation encountered in optimization, where the error for the objective function values is highlighted.

This manuscript is organized as follows. Section 2 provides the necessary background on entropy functionals, as well as on well-definedness of the proposed iterative procedure. Section 3 analyzes (weak) convergence of the method when both a priori and a posteriori stopping rules are considered, while Section 4 deals with error estimates only for the former rule. Section 5 explores a version of the entropic Landweber method for nonquadratic data fidelity terms. The theoretical results are tested in Section 6 on several integral equation examples, in comparison with the Expectation-Maximization algorithm and the projected Landweber method - see [13] for an overview on regularization methods for nonnegative solutions of ill-posed equations.

2 Preliminaries

In the following we collect some basic results and assumptions needed for the analysis below. We start with properties of the entropy and then proceed to the operator $A$.

2.1 Entropy and Entropic Projection

Let $\Omega$ be an open and bounded subset of ${\mathbb{R}}^{d}$. The negative of the Boltzmann-Shannon entropy is the function $f:L^{1}(\Omega)\rightarrow(-\infty,+\infty]$, given by111We use the convention $0\ln 0=0$.

[TABLE]

Here and in what follows ${{L}}_{+}^{p}(\Omega)$, $p\in[1,\infty]$, stands for the set $\{u\in{{L}}^{p}(\Omega):u(t)\geq 0\,\,\,\mbox{a.e.}\}$, while $\|\cdot\|_{p}$ denotes, as usual, the norm of the space $L^{p}(\Omega)$.

The Kullback-Leibler functional or the Bregman distance with respect to the Bolzmann-Shannon entropy can be defined as $d:{\mbox{dom}\,f}\times{\mbox{dom}\,f}\rightarrow[0,+\infty]$ by

[TABLE]

where $f^{\prime}(u,\cdot)$ is the directional derivative at $u$. Here $\mbox{dom}\,f=\{u\in{{L}}^{1}(\Omega):f(u)<\infty\}$ denotes the domain of $f$. One can also write

[TABLE]

if $d(v,u)$ is finite, as one can see below.

Some properties of the entropy functionals $f$ and $d$ are recalled below (see, e.g., [26, 27]).

Lemma 2.1

The function defined by (5) has the following properties:

1. (i)

The domain of the function $f$ is strictly included in ${{L}}_{+}^{1}(\Omega)$. 2. (ii)

The interior of the domain of the function $f$ is empty. 3. (iii)

The set $\partial f(u)$ is nonempty if and only if $u$ belongs to $L_{+}^{\infty}(\Omega)$ and is bounded away from zero. Moreover, $\partial f(u)=\{1+\ln u\}$. 4. (iv)

The directional derivative of the function $f$ is given by

[TABLE]

whenever it is finite. 5. (v)

For any $u,v\in{\mbox{dom}\,f}$, one has

[TABLE]

Based on Lemma 2.1 (iii) we define in the following

[TABLE]

Lemma 2.2

The statements below hold true:

1. (i)

The function $(v,u)\mapsto d(v,u)$ is convex; 2. (ii)

The function $d(\cdot,u^{*})$ is lower semicontinuous with respect to the weak topology of $L^{1}(\Omega)$, whenever $u^{*}\in{\mbox{dom}\,f}$; 3. (iii)

For any $C>0$ and any nonnegative $u\in L^{1}(\Omega)$, the following sets are weakly compact in $L^{1}(\Omega)$:

[TABLE] 4. (iv)

The set $\partial d(\cdot,u^{*})(u)$ is nonempty for $u^{*}\in{\mbox{dom}\,f}$ if and only if $u$ belongs to $L_{+}^{\infty}(\Omega)$ and is bounded away from zero. Moreover, $\partial d(\cdot,u^{*})(u)=\{\ln u-\ln u^{*}\}$.

Denote

[TABLE]

for $u,v\in L^{1}(\Omega)$, when the integral exists.

A key observation for obtaining well-definedness of the iterative scheme as well as an explicit form for the iterates is the following result on the entropic projection.

Proposition 2.3

Let $\ell\in L^{\infty}(\Omega)$ and $v\in$ dom $\partial f$. Then the problem

[TABLE]

has a unique solution in the cases $m=0$ and $m=1$, respectively, given by

[TABLE]

which satisfies $u_{m}\in$ dom $\partial f$.

Proof: We simply rewrite the functional as

[TABLE]

where $C_{m}=\ln c_{m}-\int_{\Omega}u_{m}(t)~{}dt$ is a constant independent of $u$. It is straightforward to notice that

[TABLE]

Hence, the problem is equivalent to minimizing $d(u,u_{m})+\chi_{m}(u)$. Since both terms are nonnegative and vanish for $u=u_{m}$, we see that $u_{m}$ is indeed a minimizer in ${\mbox{dom}\,f}$. Strict convexity of $d$ implies the uniqueness and since $u_{m}$ is the product of $v$ with a function strictly bounded away from zero it also satisfies $u_{m}\in$ dom $\partial f$.

2.2 Forward operators and entropy

In this paper we always assume that $A:L^{1}(\Omega)\rightarrow Y$ is a linear and bounded operator with $Y$ being a Hilbert space. In addition to the norm boundedness of $A$, we assume a continuity property in terms of the Bregman distance. More precisely we assume that

[TABLE]

holds on dom$(f+\chi_{m})$ in the respective cases $m=0$ or $m=1$ for some positive number $\gamma$. It is easy to see that the latter is already implied by the boundedness of $A$ in case $m=1$:

Lemma 2.4

Let $A$ be as above with $\|A\|$ denoting its operator norm, let $u,v\in$ dom$(f+\chi_{1})$, and $v\in$dom $\partial f$. Then (12) is satisfied with $\gamma=\sqrt{2}\|A\|$.

Proof: By the boundedness of $A$ we have

[TABLE]

Lemma 2.1 (v) further implies

[TABLE]

which yields the assertion.

We define the nonlinear functional

[TABLE]

which will be useful for the further analysis. Note that $D(u,v)\geq 0$ for $u,v\in$ dom$(f+\chi_{1})$, and $v\in$dom $\partial f$, whenever $c\geq\frac{\gamma^{2}}{2}$ (cf. Lemma 2.4). In case $m=0,$ we restrict the analysis to the class of operators $A$ for which $D(u,v)\geq 0$ for $u\in$ dom$f$ and $v\in$dom $\partial f$.

3 Convergence of the Entropic Landweber Method

In the following we consider the iterative method defined by (4), where $d$ is the Kullback-Leibler divergence given by (7).

Noticing that $A^{*}$ maps to $L^{\infty}(\Omega)$, we can equivalently rewrite the minimization in (4) in the form of Proposition 2.3, which implies the following result.

Proposition 3.1

Let $u_{0}\in$ dom $\partial f$. Then there exists a unique minimizer in (4) for any $k\geq 0$, given by ($\lambda=\frac{1}{c}$)

[TABLE]

which further satisfies $u_{k+1}\in$ dom $\partial f$.

Note that, from pointwise manipulation of (14) we rigorously obtain the first-order optimality condition for the variational problem in each step, i.e.,

[TABLE]

where $\ln c_{k}^{0}=0$ and $\ln c_{k}^{1}$ is to be interpreted as a Lagrange multiplier for the integral constraint. In the latter case, this constant term is orthogonal to all functions of the form $w-v$, where $\chi_{1}(w)=\chi_{1}(v)=0$. Since most estimates below for iterates will be based on taking duality products of (15) with such functions, they can be carried out in the same way for $m=0$ and $m=1$.

The analysis of the above method ressembles the one for proximal point methods, which is apparent from rewriting (4) as

[TABLE]

However, the quantity $D$ is neither a metric distance nor necessarily a Bregman distance of a convex function, rather a weighted difference of Bregman distances. This and the involved Kullback-Leibler divergence in an infinite dimensional setting require thus a careful investigation.

We choose $u_{0}\in\mbox{dom}\,\partial f$ such that $\xi_{0}:=1+\ln u_{0}\in\mathcal{R}(A^{*})$, that is $\xi_{0}=\lambda A^{*}w_{0}$ for some $w_{0}\in Y$, and denote

[TABLE]

for $k\geq 1$. Then (15) can be expressed as $\displaystyle{\xi_{k}=\xi_{k-1}+\ln c_{k-1}^{m}+\lambda A^{*}(y-Au_{k-1})}$ which implies

[TABLE]

We show next that the entropic Landweber method converges in the exact data case.

Proposition 3.2

Let $A:L^{1}(\Omega)\rightarrow Y$ be a bounded linear operator which satisfies (12) and such that the operator equation $Au=y$ has a positive solution $z$ verifying $\chi_{m}(z)=0$ if $m=1$. Let $u_{0}\in\mbox{dom}\,\partial f$ be an arbitrary starting element such that $1+\ln u_{0}\in\mathcal{R}(A^{*})$. Moreover, let $\chi_{m}(u_{0})=0$ if $m=1$. Then the following statements are true:

1. (i)

The residual $\|Au_{k}-y\|$ decreases monotonically. 2. (ii)

The term $D(z,u_{k})$ decreases monotonically. 3. (iii)

The sequences $\{u_{k}\}_{k\in\mathbb{N}}$ generated by the iterative method (14) converge weakly on subsequences in $L^{1}(\Omega)$ to solutions of the equation $Au=y$, with $\chi_{m}(u)=0$ if $m=1$.

Proof: We will use the proximal point method techniques in order to prove the statements, by taking care of the fact that $D$ is a nonnegative functional satisfying $D(u,u)=0$ for any $u$ in this function’s domain.

(i) We have for all $k\in\mathbb{N}$,

[TABLE]

which implies that the sequence $\{\|Au_{k}-y\|\}$ is nonincreasing, since $D(u_{k+1},u_{k})\geq 0$.

(ii) Consider first the case $m=0$. Let $z$ verify $Az=y$ and denote

[TABLE]

By using (15), one has for all $k\in\mathbb{N}$:

[TABLE]

This implies the typical inequality for a proximal-like method:

[TABLE]

which yields the conclusion.

(iii) Let $m=0$. Inequality (19) leads to

[TABLE]

which yields

[TABLE]

since the sequence $\{\|Au_{k}-y\|^{2}\}_{k\in\mathbb{N}}$ is monotone. We show now that $\{f(u_{k})\}_{k\in\mathbb{N}}$ is bounded. To this end, due to nonnegativity of $D(z,u_{k})$, $k\in\mathbb{N}$, and to (18), one has

[TABLE]

The right hand side is bounded by (20) and (21), thus ensuring boundedness of $\{f(u_{k})\}_{k\in\mathbb{N}}$. Consequently, there exists a subsequence $\{u_{l}\}_{l\in\mathbb{N}}$ in $\mbox{dom}\,\partial f$ which is $L^{1}$-weakly convergent to some $u\in{\mbox{dom}\,f}$, cf. Lemma 2.2. Then one has $Au_{l}\rightarrow Au$ weakly in $Y$ and moreover $Au_{l}\rightarrow y$ in the $Y$-norm since

[TABLE]

due to inequality (19) and to monotonicity of $\{D(z,u_{k})\}_{k\in\mathbb{N}}$. Hence, $u$ satisfies $Au=y$.

The proof of the statements above for the case $m=1$ is similar, the main difference being the optimality condition (15) with $c_{k}^{m}$ satisfying $c(\xi_{k}-\xi_{k+1})=A^{*}(Au_{k}-y)+c\ln c_{k}^{m}$. In more detail, the term $\langle c\ln c_{k}^{m},z-u_{k+1}\rangle$ vanishes when evaluating $a$ and does not influence further calculations, while other terms containing $c_{k}^{m}$ behave similarly in the remaining argumentation.

Let us consider now the iterative method based on the noisy data, that is

[TABLE]

We propose first a discrepancy principle for stopping the algorithm in this case. Before detailing how it works, denote

[TABLE]

for $k\geq 1$. Then the optimality condition for (22) yields

[TABLE]

Proposition 3.3

Assume that $A:L^{1}(\Omega)\rightarrow Y$ is a bounded linear operator which satisfies (12) and such that the operator equation $Au=y$ has a positive solution $z$ verifying $\chi_{m}(z)=0$ if $m=1$. Let $y^{\delta}\in Y$ be noisy data satisfying $\|y-y^{\delta}\|\leq\delta$, for some noise level $\delta$. Let $u_{0}\in\mbox{dom}\,\partial f$ be an arbitrary starting element with the properties $1+\ln u_{0}\in\mathcal{R}(A^{*})$ and $\chi_{m}(u_{0})=0$ if $m=1$. Then

1. (i)

The residual $\|Au_{k}-y^{\delta}\|$ decreasesmonotonically and the following inequalities hold

[TABLE]

[TABLE] 2. (ii)

The term $D(z,u_{k})$ decreases as long as $\|y^{\delta}-Au_{k}\|^{2}>\delta^{2}.$ 3. (iii)

The index $k_{*}(\delta)$ defined by

[TABLE]

is finite. 4. (iv)

There exists a weakly convergent subsequence of $\{u_{k_{*}(\delta)}\}_{\delta}$ in $L^{1}(\Omega)$. If $\{{k_{*}(\delta})\}_{\delta}$ is unbounded, then each limit point is a solution of $Au=f$.

Proof: We consider only the case $m=0$, since for $m=1$ one can use similar arguments, as explained in the previous proof.

First part of (i) follows by the definition of the iterative procedure. For proving the remaining inequalities in (i), and (ii), we consider as in the previous proof

[TABLE]

Inequality (26) can be obtained by writing (25) for $k=0,...,n-1$ und calculating the telescope sum:

[TABLE]

Moreover, (ii) follows from (25) by neglecting $D(u_{k+1},u_{k})$.

(iii) follows from (26) and the definition of $k_{*}(\delta)$:

[TABLE]

which implies

[TABLE]

(iv) can be shown similarly to Proposition 3.2 (iii). Due to nonnegativity of $D(z,u_{k})$ for any $k\in\mathbb{N}$ and to (24), one has

[TABLE]

The right hand side written for $k=k_{*}(\delta)$ is bounded by (28), (27), the monotonicity of the residual and by (29), thus ensuring boundedness of $\{f(u_{k_{*}(\delta)}\}_{\delta>0}$ for $\delta$ small enough. The conclusion follows then as in the proof of Proposition 3.2 (iii).

A convergence result can be established also in case of an a priori stopping rule with $k_{*}(\delta)\sim\frac{1}{\delta}$ by following the lines of Proposition 3.3 (iv).

Proposition 3.4

Assume that $A:L^{1}(\Omega)\rightarrow Y$ is a bounded linear operator which satisfies (12) and such that the operator equation $Au=y$ has a positive solution $z$ verifying $\chi_{m}(z)=0$ if $m=1$. Let $y^{\delta}\in Y$ be noisy data satisfying $\|y-y^{\delta}\|\leq\delta$, for some noise level $\delta$. Let $u_{0}\in\mbox{dom}\,\partial f$ be an arbitrary starting element with the properties $1+\ln u_{0}\in\mathcal{R}(A^{*})$ and $\chi_{m}(u_{0})=0$ if $m=1$. Let the stopping index $k_{*}(\delta)$ be chosen of order $1/{\delta}$. Then $\{f(u_{k_{*}(\delta)})\}_{\delta}$ is bounded and hence, as $\delta\rightarrow 0$, there exists a weakly convergent subsequence $\{u_{k(\delta_{n})}\}_{n}$ in $L^{1}(\Omega)$ whose limit is a solution of $Au=y$. Moreover, if the solution of the equation is unique, then $\{u_{k_{*}(\delta)}\}_{\delta>0}$ converges weakly to the solution as ${\delta\rightarrow 0}$.

4 Error estimates

In this section we derive error estimates under a specific source condition (on a solution) for the entropy type penalty. We proceed first with the case of exact data on the right-hand side of the operator equation and then with the noisy data case, by employing an a priori rule for stopping the algorithm.

4.1 Exact data case

Proposition 4.1

Assume that $A:L^{1}(\Omega)\rightarrow Y$ is a bounded linear operator which satisfies (12) and such that the operator equation $Au=y$ has a positive solution $z$ verifying $\chi_{m}(z)=0$ if $m=1$. Let $u_{0}\in\mbox{dom}\,\partial f$ be an arbitrary starting element with the properties $1+\ln u_{0}\in\mathcal{R}(A^{*})$ and $\chi_{m}(u_{0})=0$ if $m=1$. Additionally, let the following source condition hold:

[TABLE]

Then one has

[TABLE]

Moreover, $\|u_{k}-z\|_{1}=O(1/\sqrt{k})$ if $m=1$.

Proof: We consider only the case $m=0$ (similar arguments for the other case).

First, we symmetrize $D$ by considering $D^{s}(x,y)=D(x,y)+D(y,x)$. Let $\xi=1+\ln z=\lambda A^{*}v$ for some $v\in Y$.

One can use similar techniques as in [9] for deriving the announced error estimates, by carefully dealing with the setting of the $D$ distance penalty. Based on (17), one has

[TABLE]

By writing the last inequality also for $k-1,k-2,...,1$, by summing up and by combining with monotonicity of $\{D(z,u_{k})\}$, one obtains

[TABLE]

and thus, due to (21),

[TABLE]

holds. The announced convergence rate in the $L^{1}$-norm holds in case $m=1$ by Lemma 2.1 (v).

4.2 Noisy data case

Proposition 4.2

Assume that $A:L^{1}(\Omega)\rightarrow Y$ is a bounded linear operator which satisfies (12) and such that the operator equation $Au=y$ has a positive solution $z$ verifying $\chi_{m}(z)=0$ if $m=1$. Let $u_{0}\in\mbox{dom}\,\partial f$ be an arbitrary starting element with the properties $1+\ln u_{0}\in\mathcal{R}(A^{*})$ and $\chi_{m}(u_{0})=0$ if $m=1$. Let $y^{\delta}\in Y$ be noisy data satisfying $\|y-y^{\delta}\|\leq\delta$, for some noise level $\delta$. Let the stopping index $k_{*}(\delta)$ be chosen of order $1/\delta$. and let the source condition (30) hold. Then one has

[TABLE]

Moreover, $\|u_{k_{*}(\delta)}-z\|_{1}=O(\sqrt{\delta})$ if $m=1$.

Proof: Note that $c\xi-A^{*}Az=A^{*}q$ for $q=v-A^{*}z$. With this notation, one can show the following estimate as in Theorem 4.3 in [9]:

[TABLE]

Then one has

[TABLE]

which yields (32) when written for $k=k_{*}(\delta)$, due to (26).

Establishing convergence rates by means of a discrepancy rule remains an open issue.

5 General data fidelities

Before we conclude with numerical examples, we want to emphasize that Problem (1) can easily be generalised to

[TABLE]

Here $F_{y^{\delta}}:Y\rightarrow[0,+\infty)$ is a more general data fidelity term that is assumed to be convex and Fréchet-differentiable and $g:L^{1}(\Omega)\times L^{1}(\Omega)\rightarrow[0,+\infty)$ is the Bregman distance with respect to the function $F_{y^{\delta}}$, i.e.

[TABLE]

Note that (33) is an instance of the Bregman proximal method [11, 17]. The update for (33) can be written, in analogy to (14), as

[TABLE]

for $\lambda=1/c$ and $m\in\{0,1\}$. We want to emphasise that a more general data fidelity term that satisfies the assumptions mentioned above together with

[TABLE]

for all $u,v\in dom\,f$, is no restriction in terms of Fejér-monotonicity. In analogy to [7, Lemma 6.11] we can conclude

[TABLE]

for all $k<k_{\ast}(\delta)$, with $k_{\ast}(\delta)$ chosen according to a modified version of (27) that reads as

[TABLE]

However, we can also derive a monotonicity result for $d$ directly. First of all we observe that (35) implies

[TABLE]

Inserting (34) into the inequality above then yields

[TABLE]

As mentioned earlier in Section 3, we either have $\ln(c_{k}^{m})=0$ for $m=0$, or orthogonality of $\ln(c_{k}^{m})$ to all functions of the form $w-v$ with $\chi_{m}(w)=\chi_{m}(v)=0$ for $m=1$. Since $d(u_{k+1},u_{k})+d(u_{k},u_{k+1})=\langle\ln u_{k+1}-\ln u_{k},u_{k+1}-u_{k}\rangle$, we therefore estimate

[TABLE]

With the three-point identity we then observe

[TABLE]

Together with (36) we can then conclude

[TABLE]

for $k<k^{\ast}(\delta)$.

6 Example Problems

We finally discuss several types of problems that satisfy the conditions used in the analysis and present numerical illustrations for some of these situations.

6.1 Integral Equations

Let $\Omega\subset I\!\!R^{d}$ and $\tilde{\Omega}\subset I\!\!R^{\tilde{d}}$ be open and bounded sets and let $k\in L^{\infty}(\tilde{\Omega}\times\Omega)$. Then the integral operator

[TABLE]

is a well-defined and bounded linear operator. Thus, the convergence analysis is applicable due to Lemma 2.4.

We mention that in the case of $k$ being a nonnegative function, and hence $A$ and $A^{*}$ preserving nonnegativity, standard schemes preserving nonnegativity are available. In particular for $k$ including negative entries, the entropic Landweber scheme offers a straightforward alternative, since it does not depend on the positivity preservation of $A$ respectively its adjoint. For comparison we consider the EM-Algorithm

[TABLE]

and the projected Landweber iteration

[TABLE]

We implement the forward operator by discretization of $u$ on a uniform grid and a trapezoidal rule for integration. We use the following examples of kernels and initial values, all on $\Omega=(0,1)$, the first two being standard test examples used in the literature on maximum entropy methods (cf. [1])

Kernel $k_{1}(x,y)=e^{xy}$, exact solution $z_{1}(x)=e^{-\frac{x^{2}}{2\sigma^{2}}}$ 2. 2.

Kernel $k_{2}(x,y)=3e^{-\frac{(x-y)^{2}}{0.04}}$, exact solution

[TABLE] 3. 3.

Kernel $k_{3}(x,y)=1$ if $x\geq y$ and $k_{3}(x,y)=0$ else, exact solution $z_{3}(x)=e^{-\frac{x^{2}}{2\sigma^{2}}}$ .

In all examples, we have chosen $\sigma^{2}=0.01$ and a constant intial value $u^{0}$. In order to illustrate the behaviour of the iteration methods we plot the the error $\|u^{k}-z\|_{L^{1}}$ vs. the iteration number $k$ in Figure 1.

We observe that the entropic projection is at least competitive to the other schemes in all examples, it outperforms the EM and projection method in the first example, which is a combination of severe ill-posedness with an exact solution having many entries close to zero (which is a particularly difficult case for the EM algorithm).

In the second case, again severely ill-posed, the projected Landweber iteration performs better, mainly due to the strong initial decrease when the solution is positive and no projection is applied.

The third case corresponding to numerical differentiation, i.e. a very mildly ill-posed problem, is characterized by fast convergence of the schemes, but again the projection method converges significantly slower. For comparison we also include the stochastic version of the entropic projection method, with only one equation used in each iteration step, hence a highly efficient computation. That is, the operator $A$ is divided in $M$ blocks $A=(A_{1},A_{2},\ldots,A_{M})^{T}$, and the data $y$ are partitioned in the same way: $y=(y_{1},y_{2},\ldots,y_{M})$. With $J(k)$ a discrete uniform random variable in $\{1,\ldots,M\}$, we compute the iterates

[TABLE]

The initial convergence curve is similar to the other method, with much lower computational effort, then the asymptotic convergence close to the exact solution becomes significantly slower. Hence, it might be very attractive to use the stochastic version at least for the first phase of the reconstruction.

6.2 Discrete sampling of continuous probability densities

Suppose that our forward operator is the Fourier integral of a real-valued function evaluated at discrete samples $\xi_{1},\ldots,\xi_{n}$ on a compact domain $\Omega\subset I\!\!R^{d}$, i.e.

[TABLE]

Then the adjoint operator $A^{\ast}$ that satisfies $\sum_{j=1}^{n}(Au)_{j}\overline{v_{j}}=\int_{\Omega}u(x)\overline{(A^{\ast}v)(x)}dx=\int_{\Omega}u(x)(A^{\ast}v)(x)dx$ is given as

[TABLE]

where Re denotes the real part of a complex function. For this choice of $A$ the iterates of (4) read

[TABLE]

where $c_{k}^{m}$ is defined as in (14). Note that this update can also be written as

[TABLE]

Due to $\tilde{y}_{j(k+1)}=\tilde{y}_{jk}+(2\pi)^{-\frac{d}{2}}\int_{I\!\!R^{d}}u_{k}(t)e^{-it\cdot\xi_{j}}\,dt$, this formulation has the advantage that the numerical costs for evaluating the integrals remains constant.

In the following we consider a one-dimensional setting ($d=1$) with $\Omega=[-a,a]$ for $a=10$, where we measure $n=16$ samples $\{y_{j}\}_{j=1}^{n}$ of the Fourier integral for coordinates $\xi_{j}=(2\pi(j-1))/n$, $j=1,\ldots,n$. We assume that these measurements are of the form

[TABLE]

for a function $z\in L^{1}([-10,10])$ and where $n_{j}\in\mathcal{N}(0,\sigma^{2})$ are normal-distributed random variables with mean zero and variance $\sigma^{2}$, for all $j\in\{1,\ldots,n\}$. We consider numerical experiments for two choices of $z$. The first choice is the following Gaußian-mixture model,

[TABLE]

that is constructed as a linear combination of three normalised Gaußians, i.e.

[TABLE]

Note that $\tilde{z}$ does not satisfy the source condition (30), which is why we design a second function

[TABLE]

where $c$ is chosen to ensure $\int_{-10}^{10}z(x)\,dx=1$, which by construction satisfies (30). We design two functions $z_{1}$ and $z_{2}$; $z_{1}$ is defined as $z_{1}:=z$ for $\tilde{z}$ with means $\mu_{1}=0,\mu_{2}=-1$ and $\mu_{3}=1/2$, standard deviations $\sigma_{1}=1,\sigma_{2}=1/10$ and $\sigma_{3}=1/4$, and coefficients $c_{1}=1/10,c_{2}=3/5$ and $c_{3}=3/10$. The function $z_{2}:=\tilde{z}$ has the same means and standard deviations as $\tilde{z}$ in the previous example, but coefficients $c_{1}=1/10,c_{2}=2/5$ and $c_{3}=1/2$ instead. Both functions are visualised in Figure 2. Subsequently we create data samples via (38) with noise levels $\sigma=0$ and $\sigma=1/500$ for $z_{1}$, respectively $\sigma=0$ and $\sigma=1/50$ for $z_{2}$.

In the following we run the entropic projection method (14) for $j=1$, $\tau=9/(10\sqrt{2\pi})$ and with the initial function

[TABLE]

either until the discrepancy principle (27) is violated (for $\tau=1$) or until we reach a certain maximum number of iterations. We first investigate the algorithm for the function $z_{1}$ as seen in Figure 2(a), for perfect data ($\sigma=0$) and for noisy data $\sigma=1/500$. For perfect data we run the algorithm for 201 iterations and observe that we are converging towards $z_{1}$ as can be seen in Figure 3(a) as well as in Figure 5(a), which is a numerical confirmation of Proposition 4.1. For the non-trivial noise-level $\sigma=1/500$ the algorithm stops after 65 iterations according to the discrepancy principle (Figure 3(c)).

To conclude, we run the same numerical experiments for $z_{2}$ as shown in Figure 2(b). As we mentioned earlier, (30) is violated and even for perfect data (i.e. $\sigma=0$) we cannot expect the results of Proposition 4.1 to hold true. It can be seen in Figure 4(a) that $u^{k}$ does not seem to converge towards $z_{2}$ despite a decrease of the objective to values in the order of $10^{-5}$. In fact, if we compare the $L^{1}$-norm of the difference $u^{k}-z_{2}$, we also see in Figure 5(b) that $u^{k}$ does not seem to converge towards $z_{2}$. For noisy data with $\sigma=1/500$ the discrepancy principle is violated after 46 iterations, with its result being visualised in Figure 4(c).

6.3 Initial Densities for Stochastic Differential Equations

An interesting problem in several applications, e.g. in data assimilation scenarios (cf. e.g [18]), is the reconstruction of the initial density for a system evolving via stochastic differential equations with drift $b$ and volatility $a$. The density evolves via the Fokker-Planck equation (cf. [16])

[TABLE]

in $\Omega\times(0,T)$ with no-flux boundary conditions. Under appropriate smoothness conditions on $a$ and $b$ as well as positivity of $a$ is is well-known that the Fokker-Planck equation has a unique nonnegative solution $\rho\in C(0,T;L^{1}(\Omega))$ for nonnegative initial values $u\in L^{1}(\Omega)$ such that $\int_{\Omega}u\ln u~{}dx<\infty$. In problems related to reconstructing $u$ it is hence rather natural to use methods penalizing its entropy.

The forward operator $A$ maps the initial density to indirect measurements of the density $\rho$ over time, e.g. moments or local integrals. Parametrizing the measurements by values $\sigma$ in a bounded set $\sigma$ we obtain

[TABLE]

It is well-known that Fokker-Planck equations satisfy an $L^{1}$-contractivity property on the domain of the entropy functional (cf. [20]) i.e. for $\rho_{i}$ denoting the solution with initial value

[TABLE]

for almost all $t\in(0,T)$. Thus, the map $u\mapsto\rho$ is Lipschitz continuous with unit modulus when considered as a map into $L^{\infty}(0,T;L^{1}(\Omega))$ on the domain of the entropy. Hence, if $k\in L^{\infty}(\Sigma\times\Omega)$ we can easily verify that the operator $A$ satisfies (12).

We finally mention that in the case of stationary coefficients $a$ and $b$, the Fokker-Planck equation has a unique stationary solution $\rho_{\infty}$ among nonnegative functions with unit mass (cf. [15]), to which it converges with exponential speed in the relative entropy (cf. [2, 3]), i.e. the Bregman distance related to the entropy functional. Hence, it is natural to use $\rho_{\infty}$ as an initial value for the reconstruction of $u$, since we may expect them to be close in particular in the relative entropy.

7 Conclusions and remarks

In this study we have investigated a multiplicative entropic type method for ill-posed equations, which preserves nonnegativity of the iterates. Historically, the underlying strategy has spreading roots in the inverse problems literature: Landweber iterates, surrogate functionals and linearized Bregman, to quote a few approaches. In parallel, this has been treated in different contexts in finite dimensional optimization, e.g., as a mirror descent or as a steepest descent (linearized proximal) algorithm with generalized distances, or in machine learning - as an exponentiated gradient descent algorithm for online prediction via linear models.

The closed form algorithm is shown to converge weakly in $L^{1}$ to a solution of the ill-posed problem and convergence rates are obtained by means of the Kullback-Leibler (KL) distance. All the results are quite naturally established when imposing ”mean one” restriction to the unknown, while the case without restrictions relies on a norm combined with KL distance based Lipschitz condition, in which case operators satisfying it remain to be found.

Methods of this type involving other interesting fidelity terms, nonlinear operators and eventually stochastic versions and line search strategies might be considered in more detail for future research.

8 Acknowledgements

Martin Burger acknowledges support from European Union’s Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 777826 (NoMADS). Martin Benning acknowledges support from the Leverhulme Trust Early Career Fellowship ECF-2016-611 ’Learning from mistakes: a supervised feedback-loop for imaging applications’.