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

HAL
open science

# Kullback-Leibler residual and regularization for inverse problems with noisy data and noisy operator

Bruno Sixou, Cyril Mory

► To cite this version:

Bruno Sixou, Cyril Mory. Kullback-Leibler residual and regularization for inverse problems with noisy data and noisy operator. Inverse Problems and Imaging, 2019, 13 (5), pp.1113-1137. 10.3934/ipi.2019050. hal-04177985

HAL Id: hal-04177985

https://hal.science/hal-04177985

Submitted on 7 Aug 2023

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
KULLBACK-LEIBLER RESIDUAL AND REGULARIZATION FOR
INVERSE PROBLEMS WITH NOISY DATA AND NOISY
OPERATOR

BRUNO SIXOU AND CYRIL MORY

CREATIS, CNRS UMR 5220; INSERM U1044; INSA de Lyon; Université de Lyon 1;
Université de Lyon, 69621, Villeurbanne Cedex France

ABSTRACT. We study the properties of a regularization method for inverse problems with joint Kullback-Leibler data term and regularization when the data and the operator are corrupted by some noise. We show the convergence of the method and we obtain convergence rates for the approximate solution of the inverse problem and for the operator when it is characterized by some kernel, under the assumption that some source conditions are satisfied. Numerical results showing the effect of the noise levels on the reconstructed solution are provided for Spectral Computerized Tomography.

1. **Introduction.** In this paper, we are interested in the inversion of a nonlinear operator $F$. Our aim is thus to find an approximation of a solution $x$ of the inverse problem formulated as:

$$F(x) = y. \tag{1}$$

We assume that the data are corrupted by Poisson noise. For this type of noise, the relevant distance is the Kullback-Leibler (KL) distance defined in the following. Therefore, we assume that the distance between the non noisy data $y$ and the noisy data $y^\delta$ can be estimated with the Kullback-Leibler (KL) distance:

$$KL(y^\delta, y) \le \delta^2 \tag{2}$$

where $\delta$ is a positive constant. This assumption replaces the classical one, $\|y - y^\delta\|^2 \le \delta^2$ used for Tikhonov regularization. In order to use the Kullback-Leibler (KL) distance, we assume that $F: D \to D$ is defined on the domain $D = \{v \in L_1(\Omega), d_1 \le v \le d_2\}$, where $\Omega$ is a bounded open set of $\mathbb{R}^n$ and where $d_1$ and $d_2$ are two strictly positive constants, $0 < d_1 < d_2 < \infty$. For the sake of simplicity, we consider the case where the domain and the codomain of the operator are the same, but different strictly positive upper and lower bounds could be considered for these spaces. In the classical inverse problem theory, it is assumed that the operator $F$ is known exactly. We assume here that only an approximation $F_\delta: D \to D$ of the exact operator is known. An example of such problems is the spectral computerized tomography (SPCT) inverse problem when the detector response is unknown. We are interested to estimate an approximate solution of the inverse problem from the noisy data and also to study some methods to recover at the same time this solution and the unknown operator. In the following, we investigate variational

<footnote n="2000">2000 Mathematics Subject Classification. Primary: 65J22, 65J20, 65K10; Secondary: 52A41.
Key words and phrases. Inverse Problems, Kullback-Leibler divergence.</footnote>
approaches with regularization functionals based on the Kullback-Leibler distance as discrepancy term and regularization term. There is no detailed study of this type of inverse problems with Poisson noise on the data and an inexact operator. The joint additive Kullback-Leibler (KL) residual minimization and regularization for linear inverse problems has been investigated in [21] but the operator is assumed to be well determined. Our study generalizes this work and several proofs will be similar to the ones detailed in this work. In order to investigate linear problems with inexact operators, Regularized Total Least Squares methods have been proposed with cost functionals based on the Frobenius norm and the $L_2$ norm [12]. In the framework of the Regularized Total Least Squares, some estimates of $(x, F, y)$ are determined by solving the constrained minimization problem:

$$ \text{Minimize } \|F - F_{\delta}\|^2 + \|y - y^{\delta}\|^2 \text{ subject to } Fx = y, \|Bx\| \le R \quad (3) $$

where $B$ is some unbounded densely defined self-adjoint strictly positive definite operator used to restrict the set of admissible solutions and $R$ a positive constant. Dual Regularized Total Least Squares are studied in [17], [26]. In these cases, the inverse problem is formulated as the following constrained minimization problem:

$$ \text{Minimize } \|Bx\| \text{ subject to } Fx = y, \|F - F_{\delta}\| \le h, \|y - y_{\delta}\| \le \delta \quad (4) $$

where $h$ is some bound on the noise on the operator. A double regularization approach is considered in [6] for inverse problems with bilinear operators of the unknown function $x$ and of a kernel $k$. The kernel function determines the behaviour of the operator and a double penality term is used to stabilize the reconstruction of the unknown solution and of the characteristic function governing the operator. A convergence rate analysis of Tikhonov regularization for nonlinear inverse problems with noisy operators is detailed in [18] for uniform and non uniform noise bounds. In both cases, the discrepancy term is the square of the $L_2$ norm. In this work, we extend these studies and we investigate the properties of the regularization methods minimizing functionals based on the KL distance. We estimate error bounds for the solution of the inverse problem and convergence rates depending on noise bounds on the data and the operator. These convergence rates are obtained if source conditions are satisfied, as usual for linear and nonlinear inverse problems [11, 22, 8, 24]. In order to determine also the direct operator, we consider the case where it depends linearly on a kernel. In this case, the noise of the operator is due to this unknown characteristic function. A stabilizing term is included in the regularization functional. Assuming that source conditions holds, we estimate the convergence rate for both the solution of the inverse problem and the operator.

This article is organized as follows. In the second section, we present some preliminaries about the Kullback-Leibler functional and the operator noise and in the third section we present the mathematical framework with two different regularization functionals. In the fourth section, we analyze the well-posedness and convergence properties of the model, and we detail the convergence rates for the different cases investigated. Finally, in the last section, we present some simulations results for Spectral Computerized Tomography inverse problem illustrating the effect of the noise on the data or on the operator.

## 2. Notation and preliminary results on Kullback-Leibler divergence.
The Kullback-Leibler divergence is the Bregman divergence associated with the Poisson noise distribution [4]. The Kullback-Leibler divergence between two functions $u$ and
1 $v$ in its domain is given by:

$$KL(u, v) = \int_{\Omega} (v - u + u\ln(\frac{u}{v}))d\mathbf{x}. \qquad (5)$$

2 To avoid divergencies, it is possible to define a regularized distance:

$$KL_{\epsilon}(u, v) = \int_{\Omega} (v - u + u\ln(\frac{u+\epsilon}{v+\epsilon}))d\mathbf{x} \qquad (6)$$

3 where $\epsilon$ is a small parameter [29, 30]. In the following, we will consider restricted Kullback-Leibler distances. We consider functions in $D$ for which the KL divergence leads to well-defined regularization methods[20]. The assumption that the regularized solutions belongs to the set $D$ is quite natural. It has been shown in [21] that under the assumptions that the solution $\underline{x}$ is bounded and bounded away from zero almost everywhere and that the forward operator is a linear integral operator with a non negative kernel, then the regularized solutions obtained from the minimization of functionals based on the KL distance are also bounded and bounded away from zero. In this case, the set $D$ is well-defined.

12 From the definition of the Kullback-Leibler divergence, we have the following equality for any function $u$, $v$ and $w$ in the domain of Kullback-Leibler divergence:

$$KL(u, v) + KL(w, u) - KL(w, v) = \int_{\Omega} (\ln(u) - \ln(v))(u - w)d\mathbf{x}. \qquad (7)$$

14 The following properties will be useful in the following [22, 10]:

15 **Proposition 2.1.** For any $u, v$ in the domain of the KL divergence, $\|v - u\|_1^2 \le (\frac{2}{3}\|v\|_1 + \frac{4}{3}\|u\|_1)KL(u, v)$, where $\|\cdot\|_1$ denotes the $L_1(\Omega)$ norm.

16 **Proposition 2.2.** (i) The function $(u, v) \to KL(u, v)$ is convex.

17 (ii) For any fixed $v$, the function $KL(\cdot, v)$ is lower semicontinuous with respect to the weak topology of $L^1(\Omega)$. For any fixed $v$, the function $KL(v, \cdot)$ is lower semicontinuous with respect to the weak topology of $L^1(\Omega)$. $KL$ is also lower semicontinuous with respect to the product topology of $L^1(\Omega) \times L^1(\Omega)$.

18 (iii) Let $A: L^1(\Omega) \to L^1(\Omega)$ a compact, linear operator, with a positive range. For any $C > 0$ and any strictly positive $u \in L^1(\Omega)$, the following sets are weakly compact in $L^1(\Omega): \{x \in L^1(\Omega) : KL(Ax, u) \le C\}$.

25 **Proposition 2.3.** Let $u, v \in D \subset L_1(\Omega)$, then there exists a positive constant $C$ such that:

$$KL(u, v) \le C\|u - v\|_1. \qquad (8)$$

28 Proof. The functions $u$ and $v$ are uniformly bounded away from zero and

$$KL(u, v) = \int_{\Omega} (v - u + u\ln(1 + \frac{u-v}{v}))d\mathbf{x}. \qquad (9)$$

29 With $\ln(1 + t) \le t$ for all $t > -1$, we obtain:

$$KL(u, v) \le C\|u - v\|_1 \qquad (10)$$

30 for a positive constant $C$.

31 □
1. 3. **Problem formulation.** We assume that the distance between the noisy data and the non noisy data can be evaluated with the Kullback-Leibler distance with:

3. Moreover, we assume we have a uniform operator noise and that there exists a positive constant $\delta_F$ such that:

5. We consider two different regularization functionals depending on whether we want to estimate the solution of the inverse problem from noisy data and from a noisy operator or at the same time this solution and the forward operator.

8. The first regularization functional considered $J_1$ is:

9. where $\alpha$ is a regularization parameter and $x^* \in D$ an initial guess.

10. In order to estimate at the same time the solution of the nonlinear inverse problem and the operator, we will investigate the case where the nonlinear operators $F$ and $F_\delta$ can be characterized by functions $k$ and $k^\delta \in L_2(\Omega)$. We assume also that they depend linearly on them [6, 7]. We will consider the operator $\tilde{F}$

14. such that $\tilde{F}(k, x) = F(x)$ and $\tilde{F}(k^\delta, x) = F_\delta(x)$ where $\tilde{F}$ is linear with respect to $k$ and nonlinear with respect to $x$. In this case, the error on the operators in Eq.12 can be replaced by a bound on the error on the kernel $k$.

18. We will use the following regularization functional:

19. where $\alpha$ and $\beta$ are two regularization parameters. The kernel $k^*$ and the function $x^*$ are used as initial guesses. The two terms stabilize the inversion with respect to $x$ and $k$.

22. We will rewrite the regularization functional as:

23. where $\eta$ is a positive constant. In the following, we will denote $\Phi$ the functional including the regularization terms for $x$ and $k$

25. 3.1. **Well-posedness.** In this section, we analyze the properties of the regularization functionals $J_1$ and $J_2$. We first show the well-posedness of the method and that the minimizers of the functionals exist for every parameters $\alpha$ and $\beta$.

29. **Theorem 3.1.** Assume that $\alpha > 0$, and $y^\delta \in D$. Assume that the operator $F_\delta$ is weak-to-norm continuous with respect to the weak topology of $L_1(\Omega)$. Then there exists a global minimizer of the functional $J_1$ over $D$.

$$KL(y^\delta, y) \leq \delta^2. \quad (11)$$

$$\sup_{x \in D} \|F(x) - F_\delta(x)\|_1 \leq \delta_F \quad (12)$$

$$J_1(x) = KL(y^\delta, F_\delta(x)) + \alpha KL(x, x^*) \quad (13)$$

$$\tilde{F}: L_2(\Omega) \times D \to D \\ (k, x) \mapsto \tilde{F}(k, x)$$

$$J_2(x, k) = KL(y^\delta, F_\delta(x)) + \alpha KL(x, x^*) + \beta \|k - k^*\|_2^2 \quad (14)$$

$$\Phi(x, k) = KL(x, x^*) + \eta \|k - k^*\|_2^2. \quad (16)$$
1. **Proof.** The proof is along the line of Theorem 3.1 in [14], Theorem 4.2 [6]. We restrict the KL distance to functions in $D$ [20]. The regularization functional is weak lower semicontinuous, positive (Proposition 2.2). Since $x \in D$, we can not consider that $\|x\|_1 \to \infty$ and show that the functional is coercive. Yet, since $\Omega$ is bounded and with Eq.6, we can see that $J_1$ and $KL(x, x^*)$ are bounded from above and below for $x \in D$. It is possible to find a minimizing sequence of $J_1$ and with the weak compactness of the level sets of KL intersected with $D$, it is also possible to extract a subsequence converging weakly in the $L_1(\Omega)$ topology. The weak lower semicontinuity of $J_1$ concludes the proof. □

10. **Theorem 3.2.** Assume that $\alpha > 0$, and $y^\delta \in D$. Assume that the operator $\tilde{F}$ is weak-to-norm continuous with respect to the topologies of $L_2(\Omega) \times L_1(\Omega)$ and $L_1(\Omega)$. Then there exists a global minimizer of the functional $J_2$ over $L_2(\Omega) \times D$.

11. **Proof.** The functional $J_2$ is positive, proper and coercive since it follows with Proposition 2.1 that:

$$J_2(k, x) \geq \alpha KL(x, x^*) + \beta \|k - k^*\|_2^2 \to \infty \quad (17)$$
as $\|(k, x)\| \to \infty$.

Let $m = \inf\{J_2(k, x), (k, x) \in \text{dom}(J_2)\}$ where $\text{dom}(J_2)$ is the domain of the functional $J_2$. There exists a sequence $(k_j, x_j)$ such that $J_2(k_j, x_j) \to m$. Thus the sequence $(k_j, x_j)$ is bounded. The kernels $k_j$ belongs to $L_2(\Omega)$ which is reflexive Hilbert space and from Proposition 2.2 (iv), there exist subsequences also denoted as $(k_j, x_j)$ such that $k_j \rightharpoonup \underline{k}$ and $x_j \rightharpoonup \underline{x}$ in the $L_2(\Omega)$ and $L_1(\Omega)$ topologies respectively. With the weak lower semicontinuity of the functional $J_2$ with respect to the product topology, we obtain:

$$m \leq J_2(\underline{k}, \underline{x}) \leq \liminf J_2(k_j, x_j) = \lim J_2(k_j, x_j) = m. \quad (18)$$
Thus $(\underline{k}, \underline{x})$ is a global minimizer. □

23. **3.2. KL-minimal and $\Phi$-minimal solutions.** We will consider KL-minimal solutions of Eq.1, in the following sense: an element $\underline{x} \in D$, is called a KL-minimizing solution of (1) if $F(\underline{x}) = y$ and $KL(\underline{x}, x^*) = \min\{KL(x, x^*) : x \in D, F(x) = y\}$.

24. **Proposition 3.3.** Assume that there exists a solution of (1), and that the operator $F$ is weak-to-norm continuous. Then there exists a KL-minimal solution.

25. **Proof.** There exists a sequence $(x_k)$ of solutions of (1) in $D$, such that
$$KL(x_k, x^*) \to c = \inf\{KL(x, x^*) : x \in D, F(x) = y\} \quad (19)$$

With Proposition 2.2, it is possible to extract a weakly convergent subsequence which is also denoted by $(x_k)$, with a weak limit denoted by $\underline{x}$. From the weak lower semicontinuity of KL (Proposition 2.2 (ii)), it follows that $KL(\underline{x}, x^*) \leq c$. Moreover, for all $k$, $F(x_k) = y$ and $F$ is weak-to-norm continuous, thus it follows that $F(\underline{x}) = y$. □

26. Similarly, we can define a $\Phi$-minimal solution of $\tilde{F}(h, x) = y$ with $\Phi$ given in Eq.14. The $L_2$ norm and the KL distance are weakly lower semicontinuous with respect to the weak topologies of $L_2(\Omega)$ and $L_1(\Omega)$ respectively and therefore it follows:

27. **Proposition 3.4.** Assume that there exists a solution $(\underline{k}, \underline{x})$ of $\tilde{F}(\underline{k}, \underline{x}) = y$, and that the operator $\tilde{F}$ is weak-to-norm continuous. Then there exists a $\Phi$-minimal solution.
4. **Convergence properties and convergence rates.** In this section, we study the convergence properties and the convergence rates for linear and nonlinear inverse problems solved with the approximate solutions $x_{\alpha}^{\delta}$ or $(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})$ obtained with the regularization functionals $J_1$ or $J_2$. We show that the solutions converge to some solution of the inverse problems $F(x) = y$ or $\tilde{F}(\underline{k}, x) = y$ as the noise levels tend to zero, provided the regularization parameter $\alpha$ is chosen appropriately. We will use the following proposition:

**Proposition 4.1.** Let us assume that $F$ is Fréchet differentiable in a ball $B$ around $\underline{x}$ and that there is a positive constant $L$ such that for all $x$ in $B$:

$$ \|F'(x) - F'(\underline{x})\| \leq L\|x - \underline{x}\|_1 \quad (20) $$

where $\|\cdot\|$ is the norm of the linear operators $F'(x) : L_1(\Omega) \to L_1(\Omega)$ and $F'(\underline{x}) : L_1(\Omega) \to L_1(\Omega)$. Then we have for all $x$ in $D$:

$$ \|F(x) - F(\underline{x}) - F'(\underline{x})(x - \underline{x})\|_1 \leq L\|x - \underline{x}\|_1^2/2. \quad (21) $$

This result is obtained with the mean value theorem.

4.1. **Convergence properties for the regularization functional $J_1$ and non-linear inverse problems.** We first show that under an appropriate regularization parameter choice rule, the minimizer $x_{\alpha}^{\delta}$ of the functional $J_1$ converges to a KL-minimal exact solution as the noise levels $\delta$ and $\delta_F$ converge towards zero.

**Proposition 4.2.** Let $(y_j^{\delta})$ and $(F_j^{\delta})$ two sequences of noisy data and weak-to-norm continuous operators $F_j^{\delta} : D \to D$ such that $KL(y_j^{\delta}, y) \leq \delta_j^2$ and $\sup_{x \in D} \|F_j^{\delta}(x) - F(x)\|_1 \leq \delta_{F,j}$ with $\delta_j \to 0$ and $\delta_{F,j} \to 0$. Assume that the regularization parameter $\alpha_j$ satisfies $\alpha_j \to 0$ and

$$ \lim_{j \to \infty} \frac{\delta_j + \delta_{F,j}}{\alpha_j} = 0. \quad (22) $$

Let $x_{\alpha_j}^{\delta_j}$ be the minimizer of $J_1$ obtained with the noisy data $y_j^{\delta}$, the noisy operator $F_j^{\delta}$, and the regularization parameter $\alpha_j$. There exists a convergent subsequence of $(x_{\alpha_j}^{\delta_j})$. The limit of every convergent subsequence of $(x_{\alpha_j}^{\delta_j})$ is a KL-minimal solution of $F(x) = y$.

**Proof.** The minimizing property of $x_{\alpha_j}^{\delta_j}$ guarantees that:

$$ 0 \leq J_1(x_{\alpha_j}^{\delta_j}) = KL(y_j^{\delta}, F_j^{\delta}(x_{\alpha_j}^{\delta_j})) + \alpha_j KL(x_{\alpha_j}^{\delta_j}, x^*) \leq KL(y_j^{\delta}, F_j^{\delta}(\underline{x})) + \alpha_j KL(\underline{x}, x^*). \quad (23) $$

where $\underline{x}$ is a KL-minimal solution of $F(x) = y$.

With Eq. 7, we get:

$$ KL(y_j^{\delta}, F_j^{\delta}(\underline{x})) - KL(y_j, F_j^{\delta}(\underline{x})) - KL(y_j^{\delta}, y) = - \int_{\Omega} (\ln(y) - \ln(F_j^{\delta}(\underline{x}))) (y - y_j^{\delta}) d\underline{x} \quad (24) $$

With Propositions 2.2 and 2.3, there exist positive constants $C, C_1$ such that:

$$ KL(y, F_j^{\delta}(\underline{x})) \leq C\|F(\underline{x}) - F_j^{\delta}(\underline{x})\|_1 \leq C\delta_{F,j} \quad (25) $$

$$ \left| \int_{\Omega} (\ln(y) - \ln(F_j^{\delta}(\underline{x}))) (y - y_j^{\delta}) d\underline{x} \right| \leq \|(\ln(y) - \ln(F_j^{\delta}(\underline{x})))\|_{\infty} \|y - y_j^{\delta}\|_1 \leq C_1 \delta_j. \quad (26) $$
1 where we have used Proposition 2.1 and $KL(y_j^\delta, y) \le \delta_j^2$. With Eq.22, 23 and 24

2 and with $KL(y_j^\delta, y) \le \delta_j^2$, we obtain that, as $\delta_j \to 0$, there exists a positive constant

3 $C_2$ such that:

$$KL(y_j^\delta, F_j^\delta(\bar{x})) \le C_2(\delta_j + \delta_{F,j}). \quad (27)$$

4 The sequence $KL(x_{\alpha_j}^{\delta_j}, x^*)$ is bounded and therefore there exists a weakly con-

5 vergent subsequence also denoted by $x_{\alpha_j}^{\delta_j}$ converging weakly towards $\bar{x}$. We will

6 prove that $\bar{x}$ is a KL-minimal solution. The nokse on the operator is bounded and

7 thus we have:

$$\|F(\bar{x}) - F_j^\delta(\bar{x})\|_1 \le \sup_{x \in D} \|F(x) - F_j^\delta(x)\|_1 = \delta_{F,j}. \quad (28)$$

8 $F_j^\delta$ is a weak-to-norm continuous operator and thus $\|F_j^\delta(\bar{x}) - F_j^\delta(x_{\alpha_j}^{\delta_j})\|_1 \to 0$, and

9 thus we obtain that $F_j^\delta(x_{\alpha_j}^{\delta_j}) \to F(\bar{x})$ in $L_1(\Omega)$.

10 With the weak lower semicontinuity of KL, it follows:

$$KL(y, F(\bar{x})) \le \liminf_{j \to \infty} KL(y_j^\delta, F_j^\delta(x_{\alpha_j}^{\delta_j})). \quad (29)$$

11 Therefore, we obtain:

$$0 \le KL(y, F(\bar{x})) \le \liminf_{j \to \infty} KL(y_j^\delta, F_j^\delta(x_{\alpha_j}^{\delta_j})) \quad (30)$$

$$\le \liminf_{j \to \infty} KL(y_j^\delta, F_j^\delta(x_{\alpha_j}^{\delta_j})) + \alpha_j KL(x_{\alpha_j}^{\delta_j}, x^*) \quad (31)$$

$$\le \liminf_{j \to \infty} C_2(\delta_{F,j} + \delta_j) + \alpha_j KL(\underline{x}, x^*) = 0 \quad (32)$$

12 and thus $F(\bar{x}) = y$.

13 With Eq.21 and with the weak lower semicontinuity of KL, we obtain:

$$KL(\bar{x}, x^*) \le \liminf_{j \to \infty} KL(x_{\alpha_j}^{\delta_j}, x^*) \quad (33)$$

$$\le \liminf_{j \to \infty} \frac{C_2(\delta_j + \delta_{F,j})}{\alpha_j} + KL(\underline{x}, x^*) \quad (34)$$

$$\le KL(\underline{x}, x^*) \quad (35)$$

14 and thus $\bar{x}$ is a KL-minimal solution and $KL(x_{\alpha_j}^{\delta_j}, x^*) \to KL(\underline{x}, x^*)$. With Propo- sition 2.1, it follows that $x_{\alpha_j}^{\delta_j} \to \bar{x}$ in the strong topology of $L_1(\Omega)$. $\square$

16 The former result can be applied to a sequence $(F_j^\delta)$ of bounded linear operators.

17 We now detail convergence rates based on some source condition.

18 **Theorem 4.3.** Let $F$ and $F_\delta$ two nonlinear operators satisfying the uniform noise

19 bound of Eq.12. We assume that $F$ is Fréchet differentiable and that there exists

20 a Lipschitz constant $L$ such that:

$$\|F'(x_1) - F'(x_2)\| \le L\|x_1 - x_2\|_1 \quad (36)$$

21 for all $x_1, x_2 \in D$, where $\|F'(x)\|$ denotes the norm of the linear operator $F'(x)$.

22 Moreover, we assume that the following source condition holds:

$$F'(\underline{x})^* w = \ln(\frac{\underline{x}}{x^*}) \quad (37)$$

23 for some $w \in L_\infty(\Omega)$ with $2L\|w\|_\infty \sup_{x \in D} \|x\|_1 \le 1$. In this formula, $\underline{x} \in D$ is a

24 KL-minimal solution of (1), $x^* \in D$ and $F'(\underline{x})^* : L_\infty(\Omega) \to L_\infty(\Omega)$ is the adjoint

25 of $F'(\underline{x})$. Let $y$ and $y^\delta$ such that $KL(y^\delta, y) \le \delta^2$, and $x_\alpha^\delta$ the minimizer of the
1 regularization functional $J_1$. Then with the parameter choice $\alpha \sim (\delta + \delta_F)^{1/2}$, we

2 obtain the convergence rate $KL(x_{\alpha}^{\delta}, \underline{x}) = O((\delta + \delta_F)^{1/2})$.

3 *Proof.* The proof is similar to the one of Proposition 4.2. With Eq.7, we get:

$$KL(y^{\delta}, F_{\delta}(\underline{x})) - KL(y, F_{\delta}(\underline{x})) - KL(y^{\delta}, y) = - \int_{\Omega} (\ln(y) - \ln(F_{\delta}(\underline{x}))(y - y^{\delta})dx \quad (38)$$

4 With Proposition 2.2 and 2.3, there exist positive constants $C'$ and $C_1$ such that:

$$KL(y, F_{\delta}(\underline{x})) \leq C' \|F(\underline{x}) - F_{\delta}(\underline{x})\|_{1} \leq C' \delta_{F} \quad (39)$$

$$|\int_{\Omega} (\ln(y) - \ln(F_{\delta}(\underline{x}))(y - y^{\delta})dx| \leq \|\ln(y) - \ln(F_{\delta}(\underline{x}))\|_{\infty} \|y - y^{\delta}\|_{1} \leq C_{1}\delta. \quad (40)$$

7 With $KL(y^{\delta}, y) \leq \delta^{2}$, we obtain that, as $\delta \to 0$, there exists a positive constant $C_2$ such that:

$$KL(y^{\delta}, F_{\delta}(\underline{x})) \leq C_2(\delta + \delta_F). \quad (41)$$

9 The Lipschitz continuity of the Fréchet derivative $F'$ and Proposition 4.1 imply that:

$$F(x_{\alpha}^{\delta}) = F(\underline{x}) + F'(\underline{x})(x_{\alpha}^{\delta} - \underline{x}) + r_{\alpha}^{\delta} \quad (42)$$

11 with

$$\|r_{\alpha}^{\delta}\|_{1} \leq \frac{L}{2}\|x_{\alpha}^{\delta} - \underline{x}\|_{1}^{2}. \quad (43)$$

12 The minimizers $x_{\alpha}^{\delta}$ satisfies:

$$KL(y^{\delta}, F_{\delta}(x_{\alpha}^{\delta})) + \alpha KL(x_{\alpha}^{\delta}, x^{*}) - \alpha KL(\underline{x}, x^{*}) \leq KL(y^{\delta}, F_{\delta}(\underline{x})) \quad (44)$$

13 where $\underline{x}$ is a KL-minimal solution of $F(x) = y$. With Eq.7, we have:

$$KL(\underline{x}, x^{*}) + KL(x_{\alpha}^{\delta}, \underline{x}) - KL(x_{\alpha}^{\delta}, x^{*}) = \int_{\Omega} (\ln(\underline{x}) - \ln(x^{*}))(\underline{x} - x_{\alpha}^{\delta})dx. \quad (45)$$

14 Thus with Eq.37 it follows that:

$$KL(y^{\delta}, F_{\delta}(x_{\alpha}^{\delta})) + \alpha KL(x_{\alpha}^{\delta}, \underline{x}) + \alpha \langle w, F'(\underline{x})(x_{\alpha}^{\delta} - \underline{x}) \rangle \leq KL(y^{\delta}, F_{\delta}(\underline{x})) \quad (46)$$

15 where $\langle \cdot, \cdot \rangle$ denotes the dual pairing of two functions. Moreover, we have

$$\alpha \langle w, F'(\underline{x})(x_{\alpha}^{\delta} - \underline{x}) \rangle \leq \alpha L \|w\|_{\infty} \|x_{\alpha}^{\delta} - \underline{x}\|_{1}^{2} / 2 + \alpha \|w\|_{\infty} \|F_{\delta}(x_{\alpha}^{\delta}) - F(x_{\alpha}^{\delta})\|_{1} + \alpha \|w\|_{\infty} \|y^{\delta} - F_{\delta}(x_{\alpha}^{\delta})\|_{1} + \alpha \|w\|_{\infty} \|F(\underline{x}) - y^{\delta}\|_{1} \quad (47)$$

16 and there is positive constant $C_3$ such that:

$$\|F_{\delta}(x_{\alpha}^{\delta}) - y^{\delta}\|_{1} \leq C_{3}KL(y^{\delta}, F_{\delta}(x_{\alpha}^{\delta}))^{1/2} \quad (48)$$

17 and

$$\|F(\underline{x}) - y^{\delta}\|_{1} \leq C_{3}KL(y^{\delta}, F(\underline{x}))^{1/2} \quad (49)$$

18 We obtain thus for positive constants $C_4 = 2 \sup_{x \in D} \|x\|_1$ and $C_5$:

$$|\alpha \langle w, F'(\underline{x})(x_{\alpha}^{\delta} - \underline{x}) \rangle| \leq \alpha C_4 L \|w\|_{\infty} KL(x_{\alpha}^{\delta}, \underline{x}) + \alpha \|w\|_{\infty} C_5 \delta_F + \alpha C_3 \|w\|_{\infty} \delta + \alpha C_3 \|w\|_{\infty} KL(y^{\delta}, F_{\delta}(x_{\alpha}^{\delta}))^{1/2} \quad (50)$$
1 and we can reformulate Eq.46 as:

$$
\begin{align}
& (KL(y^{\delta}, F_{\delta}(x_{\alpha}^{\delta}))^{1/2} - \alpha C_3 \|w\|_{\infty}/2)^2 + \alpha(1 - C_4 L \|w\|_{\infty}) KL(x_{\alpha}^{\delta}, \underline{x}) \le \alpha \|w\|_{\infty} C_5 \delta_F \nonumber \\
& \qquad + \alpha C_3 \|w\|_{\infty} \delta + C_2(\delta + \delta_F) + (\alpha C_3 \|w\|_{\infty})^2/4 \tag{51} \\
& \qquad \le C_6(\delta + \delta_F + (\alpha \|w\|_{\infty})^2) \tag{52}
\end{align}
$$

2 for a positive constant $C_6$. We obtain thus:

$$
KL(x_{\alpha}^{\delta}, \underline{x}) \le \frac{C_6}{1 - C_4 L \|w\|_{\infty}} \left( \frac{(\delta + \delta_F)}{\alpha} + \alpha \|w\|_{\infty}^2 \right). \quad (53)
$$

3 With the choice $\alpha \sim (\delta + \delta_F)^{1/2}$, we obtain that $KL(x_\alpha^\delta, \underline{x}) = O((\delta + \delta_F)^{1/2})$. $\square$

4 The former proposition gives a convergence rate for an a priori choice of the regularization parameter. Methods of choice of regularization parameters based on the
5 Morozov principle have been studied in detail for the minimization of regularization
6 functional based on least squares [1, 11]. A Generalized Discrepancy Principle was
7 proposed in [27] to take into account the discretization errors and the incompatibility of the data with the equation to be solved. Some rules for the choice of
8 the regularization parameter for Poisson noise have been proposed in [23, 3, 5, 15].
9 Generalizing these rules in the case where the operator is inexact will be the subject
10 of future work.
11
12

13 **4.2. Convergence properties for the regularization functional J₂ and non-**
14 **linear operators depending linearly on a kernel.** In this section, we investi-**
15 gate the regularization functional J₂ and we want to reconstruct also the forward**
16 operator. The convergence rates are thus given for k_α^δ → k, and x_α^δ → x. For the**
17 case investigated here, we will use the following source condition:

$$
\mathcal{R}(\tilde{F}'(\underline{k}, \underline{x})^*) \cap \partial\Phi(\underline{k}, \underline{x}) \neq \emptyset \quad (54)
$$

18 where $\partial\Phi$ is the subdifferential of $\Phi$ defined in Eq.16 and $\mathcal{R}$ denotes the range
19 of an operator. The former condition implies that there exists, $(\phi_x, \phi_k)$ in the
20 subdifferential of $\Phi$ at $(\underline{k}, \underline{x})$ such that:

$$
(\phi_k, \phi_x) = (\eta(\underline{k} - k^*), \ln(\frac{\underline{x}}{x^*})) = \tilde{F}'(\underline{k}, \underline{x})^* w, w \in L_1(\Omega) \quad (55)
$$

We assume that $\tilde{F}$ is Fréchet differentiable, with a Lipschitz continuous derivative.
21 For $(u, v) \in L_2(\Omega) \times L_1(\Omega)$, the Fréchet derivative of $\tilde{F}$ at a point $(k, x)$ is given
22 by:

$$
\tilde{F}'(k, x)(u, v) = \tilde{F}'(u, x) + \tilde{F}'(k, v). \tag{56}
$$

For $(u, v) \in L_2(\Omega) \times L_1(\Omega)$, the remainder of the Taylor expansion can be estimated
24 by:

$$
\| \tilde{F}(k+u, x+v) - \tilde{F}(k, x) - \tilde{F}'(k, x)(u, v) \|_1 \le C \|(u, v)\|_1^2 \quad (57)
$$

26 for a positive constant C.

We first show the convergence of the regularization method.

28 **Proposition 4.4.** Let $y_j^\delta$ a sequence of data such that $KL(y_j^\delta, y) \le \delta_j^2$ with $\delta_j \to 0$.
29 We assume that the operator $\tilde{F}$ is weak-to-norm continuous with respect to the
30 topology of $L_2(\Omega) \times L_1(\Omega)$. Assume also that the regularization parameter $\alpha_j$ satisfies
31 $\alpha_j \to 0$ and

$$
\lim_{j \to \infty} \frac{\delta_j^2}{\alpha_j} = 0 \qquad (58)
$$
1. Let $(x_{\alpha_j}^{\delta_j}, k_{\alpha_j}^{\delta_j})$ be the minimizer of $J_2$ obtained from the noisy data $y_j^{\delta}$ and regularization parameters $\alpha_j$. There exists a convergent subsequence of $(x_{\alpha_j}^{\delta_j}, k_{\alpha_j}^{\delta_j})$. The limit of every convergent subsequence of $(x_{\alpha_j}^{\delta_j}, k_{\alpha_j}^{\delta_j})$ is a $\Phi$-minimal norm solution of $\tilde{F}(k, x) = y$.

5. **Proof.** The minimizing property of $(k_{\alpha_j}^{\delta_j}, x_{\alpha_j}^{\delta_j})$ guarantees that:

$$ 
\begin{aligned}
0 \le J_2(x_{\alpha_j}^{\delta_j}, k_{\alpha_j}^{\delta_j}) &= KL(y_j^{\delta}, \tilde{F}(k_{\alpha_j}^{\delta_j}, x_{\alpha_j}^{\delta_j})) + \alpha_j KL(x_{\alpha_j}^{\delta_j}, x^*) + \beta_j \|k_{\alpha_j}^{\delta_j} - k^*\|_1^2 \\
&\le KL(y_j^{\delta}, \tilde{F}(\underline{k}, \underline{x})) + \alpha_j KL(\underline{x}, x^*) + \eta \alpha_j \|\underline{k} - k^*\|_1^2 \\
&\le \delta_j^2 + \alpha_j KL(\underline{x}, x^*) + \eta \alpha_j \|\underline{k} - k^*\|_1^2.
\end{aligned} 
\quad (59) $$

6. where $(\underline{k}, \underline{x})$ is a $\Phi$-minimal norm solution of $\tilde{F}(k, x) = y$.

7. The sequence $KL(x_{\alpha_j}^{\delta_j}, x^*)$ is bounded and therefore there exists a weakly convergent subsequence $x_{\alpha_j}^{\delta_j}$ converging weakly towards $\bar{x}$. Similarly, $\|k_{\alpha_j}^{\delta_j} - k^*\|_2$ is bounded and thus there is a subsequence $k_{\alpha_j}^{\delta_j}$ converging weakly towards $\bar{k}$. We will prove that $(\bar{k}, \bar{x})$ is a $\Phi$-minimal solution and $\tilde{F}(\bar{k}, \bar{x}) = y$.

11. We have

$$ \tilde{F}(k_{\alpha_j}^{\delta_j}, x_{\alpha_j}^{\delta_j}) - \tilde{F}(\bar{k}, \bar{x}) = \tilde{F}(k_{\alpha_j}^{\delta_j}, x_{\alpha_j}^{\delta_j}) - \tilde{F}(k_{\alpha_j}^{\delta_j}, \bar{x}) + \tilde{F}(k_{\alpha_j}^{\delta_j} - \bar{k}, \bar{x}). \quad (60) $$

12. and thus $\tilde{F}(k_{\alpha_j}^{\delta_j}, x_{\alpha_j}^{\delta_j}) \to \tilde{F}(\bar{k}, \bar{x})$ in $L_1(\Omega)$.

13. With the lower semicontinuity of KL with respect to the product topology, it follows that:

$$ KL(y, \tilde{F}(\bar{k}, \bar{x})) \le \liminf_{j \to \infty} KL(y_j^{\delta}, \tilde{F}(k_{\alpha_j}^{\delta_j}, x_{\alpha_j}^{\delta_j})). \quad (61) $$

16. Therefore, we obtain:

$$ 0 \le KL(y, \tilde{F}(\bar{k}, \bar{x})) \quad (62) $$

$$ \le \liminf_{j \to \infty} KL(y_j^{\delta}, \tilde{F}(k_{\alpha_j}^{\delta_j}, x_{\alpha_j}^{\delta_j})) + \alpha_j (KL(x_{\alpha_j}^{\delta_j}, x^*) + \eta \|k^* - k_{\alpha_j}^{\delta_j}\|_2^2) \quad (63) $$

$$ \le \liminf_{j \to \infty} \delta_j^2 + \alpha_j (KL(\underline{x}, x^*) + \eta \|k^* - \underline{k}\|_2^2) = 0 \quad (64) $$

17. thus $F(\bar{k}, \bar{x}) = y$.

18. With Eq. 59 and with the weak lower semicontinuity of KL, we obtain:

$$ KL(\bar{x}, x^*) + \eta \|k^* - \bar{k}\|_1^2 \le \liminf_{j \to \infty} KL(x_{\alpha_j}^{\delta_j}, x^*) + \eta \|k^* - k_{\alpha_j}^{\delta_j}\|_2^2 \quad (65) $$

$$ \le \liminf_{j \to \infty} \frac{\delta_j^2}{\alpha_j} + KL(\underline{x}, x^*) + \eta \|k^* - \underline{k}\|_2^2 \quad (66) $$

$$ \le KL(\underline{x}, x^*) + \eta \|k^* - \underline{k}\|_2^2 \quad (67) $$

19. Therefore $(\bar{k}, \bar{x})$ is a $\Phi$-minimal solution. $\square$

20. The next result gives some convergence rates for the regularization method based on the functional $J_2$.

22. **Theorem 4.5.** Let $y$ and $y^\delta$ such that $KL(y^\delta, y) \le \delta^2$, $(\underline{k}, \underline{x})$ a $\Phi$-minimal solution of $\tilde{F}(k, x) = y$, $(k_\alpha^\delta, x_\alpha^\delta)$ the minimizer of the regularization functional $J_2$. We assume the source condition Eq. 54 holds for $\|w\|_\infty$ small enough. Then $KL(x_\alpha^\delta, \underline{x}) + \eta \|k_\alpha^\delta - \underline{k}\|_2^2 \sim \delta^{1/2}$.
1 *Proof.* With Eq. 7, we obtain:

$$
KL(y, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) + KL(y^{\delta}, y) - KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) = \int_{\Omega} (\ln(y) - \ln(\tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}))) (y - y^{\delta}) dx \quad (68)
$$

2 and thus with Proposition 2.1, there exist some positive constants C, C₁ such that:

$$
KL(y, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) \le KL(y^{\delta}, y) + KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) + C\|y - y^{\delta}\|_{1} \quad (69)
$$

$$
\le KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) + \delta^2 + C_1 \delta. \quad (70)
$$

3 Let $U_{\alpha}^{\delta} = (k_{\alpha}^{\delta}, x_{\alpha}^{\delta})$, and $\underline{U} = (\underline{k}, \underline{x})$, the Lipschitz continuity of the Fréchet

4 derivative $\tilde{F}'$ implies that:

$$
\tilde{F}(U_{\alpha}^{\delta}) = \tilde{F}(\underline{U}) + \tilde{F}'(\underline{U})(U_{\alpha}^{\delta} - \underline{U}) + R_{x,k} \quad (71)
$$

$$
\|R_{x,k}\|_1 \le \frac{L}{2} \|U_{\alpha}^{\delta} - \underline{U}\|^2 = \frac{L}{2} (\|x_{\alpha}^{\delta} - \underline{x}\|_1^2 + \|k_{\alpha}^{\delta} - \underline{k}\|_2^2). \quad (72)
$$

6 The minimizers $(x_{\alpha}^{\delta}, k_{\alpha}^{\delta})$ satisfies:

$$
\begin{align}
KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) + \alpha KL(x_{\alpha}^{\delta}, x^{*}) - \alpha KL(\underline{x}, x^{*}) + \beta \|k_{\alpha}^{\delta} - k^{*}\|_{2}^{2} \nonumber \\
\le KL(y^{\delta}, \tilde{F}(\underline{k}, \underline{x})) + \beta \|\underline{k} - k^{*}\|_{2}^{2}. \tag{73}
\end{align}
$$

7 Let $(\phi_x, \phi_k) = (\ln(\underline{x}^*), \eta(\underline{k} - k^*))$ the subdifferential of the regularization functional

8 at $(\underline{k}, \underline{x})$:

$$
\begin{align}
& KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) + \alpha KL(x_{\alpha}^{\delta}, \underline{x}) + \eta \alpha \|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} + \alpha (\langle \phi_x, x_{\alpha}^{\delta} - \underline{x} \rangle + \langle \phi_k, k_{\alpha}^{\delta} - \underline{k} \rangle) \nonumber \\
& \le KL(y^{\delta}, \tilde{F}(\underline{k}, \underline{x})) \tag{74}
\end{align}
$$

9 and with the source condition of Eq.54:

$$
\begin{align}
& KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) + \alpha KL(x_{\alpha}^{\delta}, \underline{x}) + \eta \alpha \|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} \nonumber \\
& \quad + \alpha \langle \tilde{F}'(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})^{*} w, (x_{\alpha}^{\delta} - \underline{x}, k_{\alpha}^{\delta} - \underline{k}) \rangle \le KL(y^{\delta}, \tilde{F}(\underline{k}, \underline{x})). \tag{75}
\end{align}
$$

10 With Eq.70, we get:

$$
\begin{align}
& KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) + \alpha KL(x_{\alpha}^{\delta}, \underline{x}) + \eta \alpha \|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} \nonumber \\
& \quad + \alpha \langle \tilde{F}'(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})^{*} w, (x_{\alpha}^{\delta} - \underline{x}, k_{\alpha}^{\delta} - \underline{k}) \rangle \le KL(y^{\delta}, \tilde{F}(\underline{k}, \underline{x})) + \delta^2 + C_1 \delta. \tag{76}
\end{align}
$$

$$
\langle -\tilde{F}'(\underline{k}, \underline{x})^* w, (x_{\alpha}^{\delta} - \underline{x}, k_{\alpha}^{\delta} - \underline{k}) \rangle \le \langle w, -\tilde{F}'(\underline{k}, \underline{x})(x_{\alpha}^{\delta} - \underline{x}, k_{\alpha}^{\delta} - \underline{k}) \rangle \quad (77)
$$

$$
\le \langle w, -\tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}) + \tilde{F}(\underline{k}, \underline{x}) - R_{x,k} \rangle \quad (78)
$$

$$
\le \|w\|_{\infty} \|R_{x,k}\|
$$

$$
+ \|w\|_{\infty} \|\tilde{F}(\underline{k}, \underline{x}) - y^{\delta}\|_{1} + \|w\|_{\infty} \|y^{\delta} - \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})\|_{1}. \quad (79)
$$

11 There exist positive constants $C'$, $C'_1$ and $C_2$ such that:

$$
\|R_{x,k}\| \le C'(\|x_{\alpha}^{\delta} - \underline{x}\|_{1}^{2} + \|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2}) \le C'_{1} KL(x_{\alpha}^{\delta}, \underline{x}) + C' \|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} \quad (80)
$$

$$
\| \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}) - y^{\delta} \|_{1} \le C_{2} KL(\tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}), y^{\delta})^{1/2} \quad (81)
$$

$$
\| \tilde{F}(\underline{k}, \underline{x}) - y^{\delta} \|_{1} \le C_2 KL(\tilde{F}(\underline{k}, \underline{x}), y^{\delta})^{1/2} \le C_2 \delta. \quad (82)
$$
1 We obtain thus:

$$
\alpha \langle -\tilde{F}'(\underline{k}, \underline{x})^* w, (x_{\alpha}^{\delta} - \underline{x}, k_{\alpha}^{\delta} - \underline{k}) \rangle \le C_1' \alpha \|w\|_{\infty} K L(x_{\alpha}^{\delta}, \underline{x}) + C' \alpha \|w\|_{\infty} \|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} \\
+ C_2 \|w\|_{\infty} \alpha \delta + C_2 \alpha \|w\|_{\infty} K L(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}))^{1/2}. \quad (83)
$$

2 Eq.76 can be rewritten:

$$
\begin{align}
& KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) - C_2 \alpha \|w\|_{\infty} KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}))^{1/2} \nonumber \\
& \quad + \alpha(1 - C_1' \|w\|_{\infty}) KL(x_{\alpha}^{\delta}, \underline{x}) + \alpha(\eta - C' \|w\|_{\infty}) \|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} \nonumber \\
& \le \delta^{2} + \delta^{2} + C' \delta + C_2 \alpha \|w\|_{\infty} \delta \tag{84}
\end{align}
$$

$$
\le A_1 \delta. \tag{85}
$$

3 where $A_1$ is a positive constant. Assuming that $\|w\|_{\infty} < \min(1/C'_1, \eta/C')$, we have:

$$
KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta})) - C_2 \alpha \|w\|_{\infty} KL(\tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}), y^{\delta})^{1/2} - A_1 \delta \le 0 \quad (86)
$$

4

$$
KL(y^{\delta}, \tilde{F}(k_{\alpha}^{\delta}, x_{\alpha}^{\delta}))^{1/2} \le \frac{C_2 \alpha \|w\|_{\infty}}{2} + \frac{\sqrt{(C_2 \alpha \|w\|_{\infty})^2 + 4A_1 \delta}}{2}. \quad (87)
$$

5 With this inequality, we obtain:

$$
KL(x_{\alpha}^{\delta}, \underline{x}) \le \frac{1}{1 - C_{1}^{\prime} \|w\|_{\infty}} \left( \frac{C_{2} \alpha \|w\|_{\infty}^{2}}{2} + \frac{\|w\|_{\infty} \sqrt{(C_{2} \alpha \|w\|_{\infty})^{2} + 4A_{1} \delta}}{2} \right. \\
\left. + \frac{A_{1} \delta}{\alpha} \right) \quad (88)
$$

6

$$
\|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} \le \frac{1}{\eta - C' \|w\|_{\infty}} \left( \frac{C_2 \alpha \|w\|_{\infty}^2}{2} + \frac{\|w\|_{\infty} \sqrt{(C_2 \alpha \|w\|_{\infty})^2 + 4A_1 \delta}}{2} \right. \\
\left. + \frac{A_1 \delta}{\alpha} \right). \quad (89)
$$

7 For the parameter choice $\alpha \sim \delta^{1/2}$, we obtain:

$$
KL(x_{\alpha}^{\delta}, \underline{x}) = O(\delta^{1/2}) \qquad (90)
$$

8

$$
\|k_{\alpha}^{\delta} - \underline{k}\|_{2}^{2} = O(\delta^{1/2}). \qquad (91)
$$

9

□

10 5. **A numerical experiment: spectral Computerized Tomography.** In this
11 section, we present some numerical experiments showing the effect of the noise levels
12 on the forward operator or on the data on the quality of the approximate solution for
13 a nonlinear inverse problem obtained with joint KL data term and regularization.
14 Multi-energy Computerized Tomography, or spectral CT, is an imaging modality
15 where the information received depends on the energy of each photon hitting the
detector. This information is provided by a photon counting detector which gives an
16 image for each energy bin. Assuming that the objects attenuation can be obtained
17 by linear combination of the attenuations of only a few materials, this energy-
18 resolved information allows to reconstruct several volumes, each one representing
19 a different material concentration map (e.g soft tissues and bone) [16, 25, 9]. The
20 measured projections are corrupted by Poisson noise. The spectral CT problem
21 sets an ill-posed inverse problem that has to be regularized to obtain stable results
22 for material decomposition. In the first subsection, we present the forward model
23 relating the projected mass of the material to the number of photons detected and
24 the regularization functional. In the second section, we detail the numerical results.
25
1 5.1. Forward model for spectral CT. For the sake of simplicity, we consider a
2 2D object to reconstruct. The forward model presented here can be easily extended
3 to 3D. We denote $\Omega$ the space that contains the 2-dimensional (2-D) object studied
4 and $\mathbf{x}$ the voxel position. This object is imaged on a 1-D detector for several
5 projection angles. The number of photons transmitted at energy $E$ for each pixel
6 $u$ of the detector and for each angle $\theta$ is denoted $n(E, u, \theta)$. The energy range
7 is denoted as $B$ and the set of the pixels values on the detector as $\Sigma$. With the
8 Beer-Lambert law, it follows:

$$
n(E, u, \theta) = n^0(E) \exp \left( - \int_{\mathcal{L}(u, \theta)} \mu(E, \mathbf{x}) \mathrm{d}\mathbf{x} \right) \quad (92)
$$

where $n^0(E)$ is the source spectrum, $\mathcal{L}(u, \theta)$ is the line defined by the X-ray beam and $\mu(E, \mathbf{x})$ is the linear attenuation coefficient at energy $E$ for the voxel $\mathbf{x}$. The number of photons detected by the detector for energy $\mathcal{E}$, pixel $u$ and angle $\theta$ can be written as :

$$
s(\mathcal{E}, u, \theta) = \int_{\mathbb{R}} d(\mathcal{E}, E) n(E, u, \theta) dE \quad (93)
$$

13 where $d(\mathcal{E}, E) : B \times B \rightarrow [0, 1]$ is the detector response function and thus:

$$
s(\mathcal{E}, u, \theta) = \int_{\mathbb{R}} d(\mathcal{E}, E) n^0(E) \exp \left( - \int_{\mathcal{L}(u, \theta)} \mu(E, \mathbf{x}) \mathrm{d}\mathbf{x} \right) dE. \quad (94)
$$

14 It describes the probability that a photon with the energy $E$ is detected at the
15 $\mathcal{E}$. We assume that we can write the attenuation coefficient as a sum of M
16 basis functions for each material that are separable in energy and space:

$$
\mu(E, \mathbf{x}) = \sum_{m=1}^{M} \rho_m(\mathbf{x}) \tau_m(E) \quad (95)
$$

17 where $\rho_m(\mathbf{x})$ is the concentration of the material $m$ at the voxel $\mathbf{x}$ and $\tau_m(E)$ is a
18 well-defined function describing the attenuation effects in the material $m$ at energy
19 $E$. The total number of materials is denoted as $M$. It is possible to reformulate the
20 direct problem with Radon projection operator which maps a function $f \in L^1(\Omega)$
21 to its line integrals. Let $\mathcal{L}(\theta, u)$ be the line defined by $\mathcal{L}(\theta, u) = \{\tau \bar{\theta}^* + u \bar{\theta} : \tau \in \mathbb{R}\}$,
22 with $\bar{\theta} = (\cos(\theta), \sin(\theta))$ and $\bar{\theta}^* = (-\sin(\theta), \cos(\theta))$, the Radon transform for
23 $f \in L^1(\Omega)$ is given by:

$$
Rf(\theta, u) = R_{\theta}f(u) = R_{u,\theta}f = \int_{\mathbf{x} \in \mathcal{L}(\theta, u) \cap \Omega} f(\mathbf{x})d\mathbf{x} \quad (96)
$$

24 and it follows:

$$
s(\mathcal{E}, u, \theta) = \int_{\mathbb{R}} d(\mathcal{E}, E) n^0(E) \exp(-R_{u,\theta}\mu) dE. \quad (97)
$$

The nonlinear forward problem can be formulated as $s = \mathcal{F}(\rho)$ with a nonlinear operator $\mathcal{F}: L_1(\Omega)^M \rightarrow L_1(B \times \Sigma \times [0, \pi])$. It can also be formulated as $s = \mathcal{F}(d, \rho)$ with a nonlinear operator $\mathcal{F}: L_2(B \times B) \times L_1(\Omega)^M \rightarrow L_1(B \times \Sigma \times [0, \pi])$. The framework presented above can thus be applied. The following proposition shows that the operator $\mathcal{F}$ satisfies the continuity assumption of the former sections.

**Proposition 5.1.** The operators $\mathcal{F}: L_1(\Omega)^M \rightarrow L_1(B \times \Sigma \times [0, \pi])$ and $\mathcal{F}: L_2(B \times B) \times L_1(\Omega)^M \rightarrow L_1(B \times \Sigma \times [0, \pi])$ are weak-to-norm continuous.
1. *Proof.* Let us assume for the sake of simplicity, that $M = 1$. Let us consider a sequence $(\rho_n)$ in $L_1(\Omega)$ converging weakly towards $\rho$ in the weak topology of $L_1(\Omega)$, then, for all $\phi \in L_\infty(\Omega)$:

$$ \int_{\Omega} \rho_n(\mathbf{x})\phi(\mathbf{x})d\mathbf{x} \to \int_{\Omega} \rho(\mathbf{x})\phi(\mathbf{x})d\mathbf{x} \text{ as } n \to \infty \quad (98) $$

4. For $u \in B$ and $\theta \in [0, \pi]$, taking $\phi$ as the indicator function of the line $\mathcal{L}(\theta, u)$, it follows that:

$$ \int_{\mathcal{L}(\theta, u)} \rho_n(\mathbf{x})d\mathbf{x} \to \int_{\mathcal{L}(\theta, u)} \rho(\mathbf{x})d\mathbf{x} \quad (99) $$

6. and thus $\mathcal{F}(\rho_n) \to \mathcal{F}(\rho)$ almost everywhere. With the dominated convergence theorem in $L_1(\Omega)$, it follows that $\mathcal{F}(\rho_n) \to \mathcal{F}(\rho)$ in $L_1(B \times \Sigma \times [0, \pi])$. Similarly, the weak convergence of a sequene $(d_n)$ towards $d$ in $L_2(B \times B)$ implies the strong convergence of $\mathcal{F}(d_n)$ towards $\mathcal{F}(d)$ in $L_1(B \times \Sigma \times [0, \pi])$. $\square$

10. The convergence properties presented in the former sections are obtained with an operator mapping a space $D \subset L_1(\Omega)$ with itself. They can be easily generalized to a mapping between different spaces of strictly positive functions included in $L_1(\Omega)$. Moreover, they may be extended to a product of spaces included in $(L_1(\Omega))^M$ by considering the regularization term $\sum_{1 \le m \le M} KL(\rho_m, \rho_m^*)$. After discretization, the inverse problem considered is a finite dimensional one. An integrated circuit counts the number of photons that are detected in the $i^{th}$ energy bin $[\mathcal{E}_i, \mathcal{E}_{i+1}]$. The number of photons detected in the $i^{th}$, for the angle $\theta$ is:

$$ s_i(u, \theta) = \mathcal{F}_i(\rho) = \int_B d_i(E)n(E, u)dE \quad (100) $$

18. where

$$ d_i(E) = \int_{\mathcal{E}_i}^{\mathcal{E}_{i+1}} d(\mathcal{E}, E)d\mathcal{E}. \quad (101) $$

19. The function $d_i(E) \in L_1(B, [0, 1])$ is the response function of the $i^{th}$ energy bin of the detector. For the numerical simulation, we consider that our detector has $P$ pixels, $I$ energy bins, and the projections are measured for $S$ angles. The detector response function $d$ is characterized by a set $(d_i)_{1 \le i \le I}$ of response functions. The projections data are thus the vector $s \in \mathbb{R}^{IPS} = (s_{ijl})_{1 \le i \le I, 1 \le j \le P, 1 \le l \le S}$. The domain to reconstruct is also discretized with $Q$ pixels. We want to recover the set of the masses of the materials $\rho \in \mathbb{R}^{MQ} = (\rho_{ij})_{1 \le i \le M, 1 \le j \le Q}$. The discretization levels are thus $n = MP$ and $m = IPS$. We consider a cost function similar to the ones in Eq. 13 and Eq. 14 in order to recover the maps of the materials:

$$ J_1(\rho) = KL(s, \mathcal{F}(\rho)) + \alpha \sum_{1 \le m \le M} KL(\rho_m, \rho_m^*) \quad (102) $$

28. and

$$ J_2(\rho, d) = KL(s, \mathcal{F}(\rho, d)) + \alpha \sum_{1 \le m \le M} KL(\rho_m, \rho_m^*) + \beta \sum_{1 \le i \le I} \|d_i - d_i^*\|_2^2. \quad (103) $$

29. **5.2. Numerical experiments.**
1. 5.2.1. *Minimization algorithms*. In this subsection, we present the minimization algorithms used for the functionals $J_1$ and $J_2$.

2. Regularization functional $J_1$

3. For the minimization of the regularization functional $J_1$, we have to optimize the concentrations map $\rho$. The objective functional is non convex with respect to the materials concentrations. During the minimization, the simulation may be trapped in local minima. Recently, several methods have been proposed which reconstruct material-specific volumes directly from the photon counts. They are commonly referred to as one-step inversion methods. All of these methods are iterative: there is currently no analytical inversion formula for the material decomposition problem. They assume a well-known detector response function [16, 2, 28, 19]. In this work, we have used the Mechlem approach to minimize the regularization functional. The principle of the method is to minimize the functional based on the KL divergence which corresponds to the Poisson negative log-likelihood of the data by Separable Quadratic Surrogates (SQS). The surrogates are derived sequentially. The minimization also integrates Ordered Subsets (OS) and Nesterov acceleration. The minimization method is described in detail in [19].

4. Regularization functional $J_2$

5. In [7], an efficient and convergent alternating minimization scheme for the minimization of doubly regularized Total Least Square was presented. In each step, the regularization functional is minimized over one variable while keeping the second one fixed. In the case investigated in [7], for each subproblem, the functional is convex and a global minimum is obtained. The sequence obtained is converging towards a critical point of the regularization functional. In this work, we have not the same convergence results since the functionals considered are not convex but we use the same optimization methodology. In the case of the functional $J_2$, numerical experiments have been performed to reconstruct both the function and the kernel. For the minimization of the regularization functional $J_2$, we alternate the minimization with respect to $\rho$ and with respect to the detector response. For the minimization with respect to $\rho$, we use the former Mechlem algorithm. For the minimization with respect to the detector response, we use a gradient descent algorithm.

6. The operator $\mathcal{F}_i$ is linear with respect to the detector response and thus the Fréchet derivative of $\mathcal{F}_i$ (Eq. 100) for $d_i \in L_2(B, [0, 1])$ with respect to the detector response is:

$$ \mathcal{F}'_i(\rho, d) : L_2(B, [0, 1]) \to L_1(\Sigma \times [0, \pi]) \\ h \mapsto \mathcal{F}_i(\rho, h) $$

7. Let $k \in L_\infty(\Sigma \times [0, \pi])$

$$ \begin{aligned} \langle k, \mathcal{F}_i(h, \rho) \rangle &= \int_{[0, \pi] \times \Sigma} k(u, \theta) \int_B d\epsilon h(\epsilon) \exp(-R_{u,\theta} \sum_m x_m(\mathbf{x})\mu_m(\epsilon)) dud\theta \\ &= \int_B d\epsilon h(\epsilon) \int_{[0, \pi] \times \Sigma} k(u, \theta) \exp(-R_{u,\theta} \sum_m x_m(\mathbf{x})\mu_m(\epsilon)) dud\theta. \end{aligned} \quad (104) $$

8. The adjoint of $\mathcal{F}'_i$ is thus the operator $\mathcal{F}^*_i : L_\infty(\Sigma \times [0, \pi]) \to L_2(B, [0, 1])$ defined by:

$$ k \mapsto \int_{[0, \pi] \times \Sigma} k(u, \theta) \exp(-R_{u, \theta} \sum_m x_m(\mathbf{x}) \mu_m(\epsilon)) dud\theta. \quad (105) $$
1 The gradient of the KL data term with respect to the detector function $d_i$ is thus
2 $\mathcal{F}_i^* Z(\mathcal{F}_i(h, \rho) - s_i^\delta)$ where $Z$ is a diagonal matrice $\text{diag}(1/(\mathcal{F}_i(\rho, d) + \epsilon))$.

3 5.2.2. *Simulations details.* We present here simulations to illustrate the former con-
4 vergence results. We can not show that the solution satisfies the source condition,
5 but we use different noise levels for the projection data $s$ and for the forward oper-
6 ator and show the effect on the quality of the reconstruction. We assume that the
7 detector response function is approximately known. The approximation of the de-
8 tector response function $d$ will be denoted $d^\delta$. Our aim is first to reconstruct $\rho$ from
9 the noisy data $s^\delta$ and from the noisy kernel $d^\delta$ using the regularization functional
10 $J_1$. Then, we give approximate solutions for the concentrations and the kernel using
11 the regularization functional $J_2$.

12 We have designed a simple two-dimensional 3-materials phantom, consisting of
13 a large square of water at 1 g/ml, a small square of iodine at 10 g/ml, and a small
14 square of gadolinium at 10 g/ml. The iodine and gadolinium squares are inside the
15 water square, but do not overlap. The phantom has 2562 voxels. This phantom is
16 displayed in Figure 1.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Ground truth maps (a) iodine (b) gadolinium (c) water.</figcaption></figure>

17 Through this phantom, 725 parallel projections were simulated, with 362 rays
18 per projection, using the AIR toolbox [13] to generate the sparse forward projection
19 matrix. The line integrals obtained were then converted to photon counts, following
the classical polychromatic Beer-Lambert attenuation law. We consider in this work 5 energy bins. The detector response functions can thus be decomposed into 5 energy response functions $d_i$. The detector response was simulated according to the model presented in appendix A.2 of [25]. The detector responses for the five energy bins are displayed in Figure 2. In the end, the photon counts were corrupted with Poisson noise. This model therefore neglects pile-up, scatter, charge sharing and probably many other complex effects.

In the simulations, the noise on the data $s$ is characterized by the value $||s^\delta -

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Detector response functions $d_i(E)$ for the five energy bins.</figcaption></figure>

$s\|_1/\|s\|_1 = \delta'$. Some Gaussian noise is added to the detector response function $d$ to obtain the noisy detector response $d^\delta$ such that the relative error on the kernel measured with the $L_2$ norm $\frac{\|d-d^\delta\|_2}{\|d\|_2}$ is 0.05 or 0.1. The noise on the operator is evaluated by an approximate value of $\sup_{\rho \in D} \|\mathcal{F}(\rho, d) - \mathcal{F}(\rho, d^\delta)\|_1 / \|\mathcal{F}(\rho, d)\|_1 = \delta'_F$. This quantity is estimated by sampling the concentration maps $\rho$ around the initial concentration $\rho^*$. The initial guess concentration $\rho^*$ used in the simulations consists of squares of gadolinium and iodine centered at the positions corresponding to the ground truth. The size of the borders is three times the ones of the ground truth squares. In order to sample the set $D$ and evaluate $\delta'_F$, we have tested considered several phantoms. The support of the materials are varied between the empty set and the support of the first guess concentration $\rho^*$ enlarged by a factor 1.5. The concentrations are chosen between 0.5 and 1.5 times the true gadolinium, water and iodine concentrations. The initial guess detector response $d^*$ is obtained from the true detector response function $d$ with a convolution with a Gaussian kernel. It is such that $\sup_{\rho \in D} \|\mathcal{F}(\rho, d) - \mathcal{F}(\rho, d^*)\|_1 / \|\mathcal{F}(\rho, d)\|_1 = 0.1$.

The iodine and gadolinium concentrations are not bounded away from 0 everywhere. In order to avoid divergencies, we have just taken into account the pixels in the images with values different from 0 for these materials. The lower box constraint is thus ensured even if the lower bound value $d_1$ is not precisely determined for these materials. Similar results are obtained with Eq.6 and a small constant $\epsilon$.

The regularization parameter $\alpha$ has been chosen in order to obtain to largest decrease of the regularization functional, with an extensive sweeping of the parameter values, multiplying or dividing it by powers of 2, and choosing the largest
possible regularization that did not cause significant amount of cross-talk. The pa-
rameter $\eta$ is chosen such that the two regularization terms in $\Phi$ have the same
order of magnitude at the end of the optimization. Several values have been
tested. The best simulation results have been obtained with $\eta = 10$. In the re-
construction of the simulated data, we have compared the concentrations of the
materials with the ground truth concentrations. In order to evaluate the quality
of the reconstruction for each material $m$, the $l_1$ relative error is calculated with
$$ \|\rho_m - \bar{\rho}_m\|_1 / \|\rho_m\|_1 = \sum_{p=1}^P |\rho_{mp} - \hat{\rho}_{mp}| / \sum_{p=1}^P |\rho_{mp}|, \text{ where } \bar{\rho}_m \text{ be the ground truth concentration map and } \rho_m \text{ the estimated concentration map} $$

10 5.3. **Results and discussion.** We have tested different noise levels for the data
11 and the detector response. We present successively the results for the functional $J_1$
12 and the functional $J_2$.
13 Functional $J_1$

14 The decrease with the iterations of the KL data term measuring the mismatch be-
15 tween the measured photons counts and those simulated through the reconstructed
16 volume are displayed in Fig.3 for the different noise levels investigated.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Evolution of the data term for different noise levels.
Bold line ($\delta' = 0.1, \delta'_F = 0$), dashed line ($\delta' = 0.1, \delta'_F = 0.1$), thin line ($\delta' = 0.1, \delta'_F = 0.2$).</figcaption></figure>

The decrease of the relative error measured with the $l_1$ norm for the concentra-
tions of iodine, water and gadolinium with the iterations are displayed on Figure 4,
5 and 6 for the different noise levels $\delta'$ and $\delta'_F$ studied.
17 18 19

Figure 7 and 8 show the last iterate for the iodine, gadolinium and water con-
centrations for the different noise levels investigated. The gray window is a very
20 broad one: for iodine and gadolinium, the values 0.0 are displayed in black and
21 values above 0.015 are displayed in white. For water, the values 0.0 are displayed
22 in black and values above 1.5 are displayed in white. Large decreases are obtained
23 for the KL data term and for the reconstruction errors for the three materials. As
24 expected, the quality of the reconstruction degrades when the noise level on the
25 forward operator increases. Some cross-talk artifacts are also visible when the noise
26 level increases.
27 28
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Evolution of the iodine relative reconstruction error. Bold line ($\delta' = 0.1, \delta'_F = 0$), dashed line ($\delta' = 0.1, \delta'_F = 0.1$), thin line ($\delta' = 0.1, \delta'_F = 0.2$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Evolution of the gadolinium relative reconstruction error. Bold line ($\delta' = 0.1, \delta'_F = 0$), dashed line ($\delta' = 0.1, \delta'_F = 0.1$), thin line ($\delta' = 0.1, \delta'_F = 0.2$).</figcaption></figure>

Functional $J_2$

2. For the functional $J_2$, we started the simulation with a noisy detector response $d^*$. A gradient descent step for the detector response is performed every 5 iterations. The decrease of the relative errors for the iodine, water and gadolinium concentrations obtained with the alternate minimization algorithm are displayed in Figure 9, 10 and 11. The final concentrations maps are displayed in Figure 12.

7. Large decreases of the errors for the three materials are observed with the minimization of the functional $J_2$. The increase of the noise level on the forward operator
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Evolution of the water relative reconstruction error. Bold line ($\delta' = 0.1, \delta'_F = 0$), dashed line ($\delta = 0.1, \delta'_F = 0.1$), thin line ($\delta' = 0.1, \delta'_F = 0.2$).</figcaption></figure>

1 degrades the accuracy of the reconstruction. The comparison of the thin and dashed lines on these figures shows the improvement achieved with the alternate minimization and the optimization of the detector response function. With these results, 2 it can be seen that better reconstruction results are obtained with a few gradient 3 steps for the detector response function. Yet, the solution seems to be trapped in a 4 local minima after a few iterations. 5

6. **Conclusion.** In this work, we have studied inverse problems with noisy data and unknown operator. We have considered regularization functionals based on the Kullback-Leibler divergence as data term and regularization term. We have first considered the case where we want only to recover the solution of linear or nonlinear inverse problems. Then we consider a more general functional when the aim is to estimate also the forward operator. We have studied the convergence properties of the regularization methods based on these functionals. We have derived convergence rates based on source conditions and restrictions on the linearity of the direct operator. Some numerical experiments illustrate the effect of the noise levels on the data and on the operator for Spectral Computerized Tomography. 7 8 9 10 11 12 13 14 15 16

17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

REFERENCES

[1] Anzengruber S.W. and Ramlau R. Discrepancy principle for Tikhonov-type functionals with nonlinear operators. Inverse Problems. 2010; 26:025001.

[2] Barber R.F., Sidky E.Y., Schmidt T.G. and Pan X. An algorithm for constrained one-step of spectral CT data. Physics in Medicine and Biology. 2016; 61: 3784.
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Reconstruction maps for the noise levels (a) iodine (b) gadolinium (c) water ($\delta' = 0.1, \delta_F = 0.1$).</figcaption></figure>

[3] Bardsley J.M. and Goldes J. Regularization parameter selection methods for ill-posed Poisson maximum likelihood estimation. Inverse Problems. 2009; 25:095005.

[4] Barndorff-Nielsen O., Information and Exponential Families in Statistical Theory. Wiley Publishers; 1978.

[5] Bertero M., Boccacci P., Talenti C., Zanella R. and Zanni L. A discrepancy principle for Poisson data. Inverse Problems. 2010; 25: 105004-105023.

[6] Bleyer J.R. and Ramlau R. A double regularization approach for inverse problems with noisy data and inexact operator. Inverse Problems. 2013; 29:025004.

[7] Bleyer I.R. and Ramlau R. An alternating iterative minimisation algorithm for the double-regularised total least square functional. Inverse Problems. 2015; 31:075004.

[8] Burger M. and Osher S. Convergence rates of convex variational regularization. Inverse Problems. 2004; 10: 1411-1422.

[9] Ducros N., Abascal J.F.P., Sixou B., Rit S. and Peyrin F. Regularization of nonlinear decomposition of spectral X-ray projection images. Med.Phys. 2017; 44: 174-187.

[10] Eggermont PPB. Maximum entropy regularization for Fredholm integral equations of the first kind. SIAM Journal on Mathematical Analysis. 1993; 24:1557-1576.

[11] Engl H.W., Hanke M. and Neubauer A. Regularization of Inverse Problems. Kluwer Academic Publishers:Dordrecht; 1996.

[12] Golub G.H, Hansen P.C and O'Leary D.P. Tikhonov regularization and total least squares. SIAM J.Numer.Anal.Appl.1999; 21:185-194.

[13] Hansen P.C. and Saxild-Hansen M. AIR Tools-A MATLAB package of algebraic iterative reconstruction methods. Journal of Computational and Applied Mathematics. 2012; 236:2067-2178.
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Reconstruction maps for the noise levels (a) iodine (b) gadolinium (c) water ($\delta' = 0.1, \delta_F = 0.2$).</figcaption></figure>

<footnote>[14] Hofmann B., Kaltenbacher B., Poschl C. and Scherzer O. A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators. Inverse Problems. 2007; 23:987-1010.</footnote>

<footnote>[15] Hohage T. and Werner F. Inverse problems with Poisson data: statistical regularization theory, applications and algorithms. Inverse Problems. 2016; 32:093001.</footnote>

<footnote>[16] Long Y. and Fessler J.A. Multi-material decomposition using statistical image reconstruction for spectral CT. IEEE Transactions on Medical Imaging. 2014; 33:1614-1626.</footnote>

<footnote>[17] Lu S., Pereverzev S.V. and Tautenhahn U. Regularized total least squares: computational aspects and error bounds. SIAM J.Numer.Anal.Appl. 2009; 31:918-41.</footnote>

<footnote>[18] Lu S. and Flemming J. Convergence rate analysis of Tikhonov regularization for nonlinear ill-posed problems with noisy operators. Inverse Problems. 2012; 28:104003.</footnote>

<footnote>[19] Mechlem K., Ehen S., Sellerer T., Braig E., Munzel D., Pfeibber F. and Noel P.B. Joint satistical iterative material image reconstruction for spectral computed tomography using a semi-empirical forward model. IEEE Transactions on Medical Imaging. 2018; 37: 68-80.</footnote>

<footnote>[20] Pöschl C. Tikhonov regularization with general residual term. PhD Thesis, Universtat Innsbruck; 2008.</footnote>

<footnote>[21] Resmerita E. and Anderssen R.S. Joint additive Kullback-Leibler residual minimization and regularization for linear inverse problem. Mathematical Methods in the Applied sciences. 2007; 30:1527-1544.</footnote>

<footnote>[22] Resmerita E. Regularization of ill-posed inverse problems in Banach spaces: convergence rates. Inverse Problems. 2005; 21:1301-1314.</footnote>
<figure><img src="image_10.png" /></figure>

FIGURE 9. Evolution of the water reconstruction error. Bold line ($\delta' = 0.1, \delta_F = 0$), thin line ($\delta' = 0.1, \delta_F = 0.1$), dashed line: solution obtained with alternate minimization.

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

FIGURE 10. Evolution of the iodine reconstruction error. Bold line ($\delta' = 0.1, \delta_F = 0$), thin line ($\delta' = 0.1, \delta_F = 0.1$), dashed line: solution obtained with alternate minimization.

<footnote>[23] Santos R.J. and De Pierro A.R. A new parameters choice method for ill-posed problem with Poisson data and its application to emission tomographic imaging. International Journal of tomography and Statistics. 2009; 11: 33–52.</footnote>

<footnote>[24] Scherzer O., Grassmair M., Grossauer H., Haltmaier M. and Lenzen F. Variational Methods in Imaging, Springer Verlag, New York; 2008.</footnote>

<footnote>[25] Schlomka J.P., Roessl E., Dorscheid R., Dill S., Martens G., Istel T., Bumer C., Herrmann C., Steadman R., Zeitler G., Livne A. and Proska R. Experimental feasiblitiy of multi-energy photon-counting k-edge imaging in pre-clinical computed tomography. Physics in Medicine and Biology. 2008; 53: 4031.</footnote>
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Evolution of the gadolinium relative reconstruction error. Bold line ($\delta' = 0.1, \delta_F = 0$), thin line ($\delta' = 0.05, \delta_F = 0.1$), dashed line: solution obtained with alternate minimization.</figcaption></figure>

1

2

3

4

5

6

7

8

9

10

11

12

13

[26] Tautenhahn U. Regularization of linear ill-posed problems with noisy right hand side and noisy operator. SIAM J.Inverse Ill-Posed Problems. 2008; 16:527-23.

[27] Wang Y., Yagola AC. and Yang C. Optimization and regularization for computational inverse problems and applications. Berlin: Springer; 2010.

[28] Weidinger T., Buzug T.M., Flohr T., Kappler S. and Stierstorfer K. Polychromatic iterative satistical material image reconstruction for photon-counting computed tomography. Interna-tional Journal of Biomedical Imaging. 2016; 2:1-15.

[29] Werner F. and Hohage T. Convergence rate in expectation ofr Tikhonov-type regularization of inverse problems with Poisson data. Inverse Problems. 2012; 28:104004.

[30] Werner F. Inverse problems with Poisson data: Tikhonov-type regularization and iteratively regularized Newton methods. PhD thesis, University of Göttingen; 2012.

Received xxxx 20xx; revised xxxx 20xx.

E-mail address: bruno.sixou@insa-lyon.fr
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Reconstruction maps for (a) iodine (b) gadolinium (c) water obtained with the iterative algorithm starting form ($\delta' = 0.1$, $\delta_F = 0.1$).</figcaption></figure>