Inverse Problems

Generally speaking, inverse problems typically consist in the reconstruction of causes for observed effects. In imaging applications the cause is usually a probe and the effect are observed data. The corresponding forward problems then consists in predicting experimental data given perfect knowledge of the probe. In some sense solving an inverse problems means “computing backwards”, which is usually more difficult then solving the forward problem. To model these kind of problems mathematically we describe the imaging system or experimental setup by a forward operator F : X → Y between Banach spaces X,Y, which maps a probe f ∈ X to the corresponding effect g ∈ Y. Then the inverse problem is, given data g ∈ Y, to find a solution f ∈ X to the equation F( f ) = g. (5.1)


What Is an Inverse Problem?
Generally speaking, inverse problems typically consist in the reconstruction of causes for observed effects. In imaging applications the cause is usually a probe and the effect are observed data. The corresponding forward problems then consists in predicting experimental data given perfect knowledge of the probe. In some sense solving an inverse problems means "computing backwards", which is usually more difficult then solving the forward problem.
To model these kind of problems mathematically we describe the imaging system or experimental setup by a forward operator F : X → Y between Banach spaces X, Y, which maps a probe f ∈ X to the corresponding effect g ∈ Y. Then the inverse problem is, given data g ∈ Y, to find a solution f ∈ X to the equation (5.1)

Ill-Posedness and Regularization
A first obvious question to ask is whether or not a probe f is uniquely determined by the data g, i.e. if the operator F is injective. Such questions are addressed e.g. in Sect. 13.6.6. Given uniqueness, one might try to find the solution f to (5.1) by just applying the inverse operator F −1 , which gives another reason to call the problem anss inverse problem. However in practice this can cause several problems, due to the forward operator being ill-posed. According to J. Hadamard a problem is called well-posed if the following conditions are satisfied: 1. There exists a solution.
2. The solution is unique. 3. The solution depends continuously on the data (stability).
Otherwise the problem is called ill-posed.
An inverse problem in the form of an operator equation (5.1) is well-posed if F is surjective (such that for all g there exists a solution), injective (such that the solution is unique) and if F −1 is continuous (guaranteeing stability). For many inverse problems in practice only the third condition is violated, and ill-posedness in the narrower sense often refers to this situation: The reconstruction of causes from observed effects is unstable since very different causes may have similar effects.
The remedy against ill-posedness is regularization: To obtain stable solutions of inverse problems one constructs a family of continuous operators R α : Y → X parameterized by a parameter α > 0 converging pointwise to the discontinuous inverse F −1 : We will discuss several generic constructions of such families of stable approximate inverses R α in §5.2.

Numerical Differentiation
In our first example we consider the forward operator given by integration. We fix the free integration constant by working in spaces of functions with mean zero, where c( f ) is such that F( f ) ∈ L 2 . The corresponding inverse problem described by the operator equation (5.1) is to compute the derivative g . If g = F( f ), then the existence of a solution to the inverse problem is guaranteed and this solution is unique. Now assume that instead of the exact data g we are given noisy data g obs that fulfill for some small noise level δ > 0. For example this could be one of the functions As the derivatives are given by (g obs This illustrates the typical ill-posedness of inverse problems: amplification of noise may be arbitrarily large for naive application of the inverse F −1 . This is the main difficulty one has to cope with. Our toy example illustrates another typical feature of inverse problems: Each function in the image of the operator F defined above is at least once differentiable. Also for many other inverse problems the forward operator is smoothing in the sense that the output function has higher smoothness than the input function, and this property causes the instability of the inverse problem.

Fluorescence Microscopy
In fluorescence microscopy one is interested in recovering the density f of fluorescent markers in some specimen in R 3 . The probe is sampled by a laser beam, and one detects fluorescent photons. In confocal microscopy, spatial resolution is achieved by focusing the laser beam by some lense and collecting fluorescent electrons by same lense such that out-of-focus fluorescent photons can be blocked by a pinhole. Let y ∈ R 3 be the focal point and assume that the probability (density) that for the focus point y ∈ R 3 a fluorescent photon emitted by a marker at point x ∈ R 3 is detected is k(x − y). k is called the point-spread function, and we assume here that it is spatially invariant. Then our problem is described by the operator equation i.e. the observation g is given by a convolution of the marker density f with point spread function k. As convolution will usually blur an image, the forward operator is smoothing. Smoother kernels will lead to stronger smoothing.

Choice of Regularization Parameters and Convergence Concepts
Due to the ill-posedness discussed above it is essential to take into account the effects of noise in the observed data. Let f † ∈ X denote the unknown exact solution, and first assume data g δ with deterministic errors such that As mentioned above, regularization of an ill-posed operator equation (5.1) with an injective operator F consists in approximating the discontinuous inverse operator F −1 by a pointwise convergent family of continuous operators R α : Y → X, α > 0. This immediately gives rise to the question which operator in the family should be chosen for the reconstruction, i.e how to choose the parameter α. Usually the starting point of deterministic error analysis in regularization theory is the following splitting of the reconstruction error: The first term on the right hand side is called propagated data noise error, and the second term is referred to as approximation error or bias. Due to pointwise convergence (see (5.2)), the bias tends to 0 as α → 0. Hence, to control this error term, we should choose α as small as possible. However, as R α converges pointwise to the discontinuous operator F −1 , the Lipschitz constant (or operator norm in the linear case) of R α will explode as α → 0, and hence also the propagated data noise error. Therefore, α must not be chosen too small. This indicates that the choice of the regularization parameter must be a crucial ingredient of a reliable regularization method. Probably the most well-known parameter choice rule in the deterministic setting is Morozov's discrepancy principle: with some parameter τ ≥ 1. In other words, among all estimators R α (g δ ) which can explain the data within the noise level (times τ ), we choose the most stable one.
Definition 5.1 A family of operators R α : Y → X parameterized by a parameter α > 0 together with some rule α : [0, ∞) × Y → (0, ∞) how to choose this parameter depending on the noise level δ and the data g δ is called a regularization method if the worst case error tends to 0 with the noise level in the sense that for all f † ∈ X.
The convergence (5.7) is a minimal requirement that one expects from a regularization method. However, it can be shown that for ill-posed problems this convergence may be arbitrarily slow depending on f † . This is of course not satisfactory. Fortunately, a-priori information on the solution f † , which is often available, may help. If we know a-priori that f † belongs to some set K ⊂ X, then it is often possible to derive explicit error bounds Let us now consider statistical noise models instead of the deterministic noise model (5.4). Often statistical data G t belong to a different space Y , e.g. a space of distributions. The distribution depends on some parameter t, and we assume that G t → F( f † ) in some sense as t → ∞. Let t e.g. denote the number of observations or in photonic imaging the expected number of photons. As the estimator R α (g obs ) (where now R α is a mapping from Y to X) will be a random variable, we have to use stochastic concepts of convergence, e.g. convergence in expectation. Other convergence concepts, in particular convergence in probability are also used frequently.
Definition 5.2 In the setting above, a family of operators R α : Y → X parameterized by a parameter α > 0 together with some parameter choice rule α : Again one may further ask not only for convergence, but even rates of convergence as t → ∞ on certain subsets K ⊂ X.

Regularization Methods
In this section we will discuss generalized Tikhonov regularization, which is given by finding the minimum of where S g obs is the data fidelity functional, which measures some kind of distance of F( f ) and the data g obs and causes the minimizerf α to still explain the data well, whereas R is the penalty functional which penalizes certain properties of the minimizer. This approach is called variational regularization, as our regularized solution is found by minimization of a functional (usually an integral functional).

Variational Regularization
We start with a probabilistic motivation of generalized Tikhonov regularization. From here until the end of this section we will consider the finite-dimensional setting X = R n , Y = R m and use boldface symbols to denote finite-dimensional vectors/mappings. We start from (5.1), where f ∈ R n , g ∈ R m and F : R n → R m is some injective function. Given the data g we want to find the solution f of F(f) = g, but recall that we cannot just apply the inverse F −1 as discussed in Sect. 5.1. Instead we might estimate f by maximizing the likelihood function L(f) = P(g|f), i.e. the probability that for a certain preimage f the data g will occur. If we assume that our data is normally distributed with covariance matrix σ 2 I , then we can rearrange the problem by using the monotonicity of the logarithm, as well as the fact, that neither additive nor multiplicative constants change the extremal point: This demonstrates the well-known fact that the maximum likelihood approach for Gaussian noise yields the least squares methods as first used by young Gauss to predict the path of the asteroid Ceres in 1801. However, asf ML = F −1 (g) for g in the range of F, this approach has no regularizing effect. In fact it is more reasonable to maximize P(f|g) instead of P(g|f), as our goal should be to find the solution f which is most likely to have caused the observation g, instead of just finding any f which causes the observation g with maximal probability. This leads to the Bayesian perspective on inverse problem with the characteristic feature that prior to the measurements a probability distribution (the so-called prior distribution) is assigned to the solution space X modeling our prior knowledge on f. By Bayes' theorem we have Estimating f by maximizing the posterior P(f|g) is called maximum a posteriori probability (MAP) estimate. To use this approach we have to model the prior P(f).
If we assume that f is normally distributed with mean f ∈ R n and covariance matrix τ 2 I , then we find is the standard (quadratic) Tikhonov functional, and therefore MAP and Tikhonov regularization coincide in this setting.
In photonic imaging the data is often given by photon counts, and these are typically Poisson distributed in the absence of read-out error (compare Chap. 4).
Recall that a random variable Z ∈ N 0 is called Poisson distributed with mean λ > 0, short Z ∼ Pois(λ), if P(Z = g) = e −λ λ g /(g!) for all g ∈ N 0 . Hence the negative log-likelihood function is given by If f has a probability distribution P( f ) = c exp(−R(f)/τ 2 ) this leads to generalized Tikhonov regularization with fidelity term S g = g − · 2 2 for normally distributed data, S g = KL(g, ·) for Poisson distributed data and the penalty term αR with regularization parameter α = τ −2 for Poisson data and α = σ 2 /τ 2 for Gaussian white noise.
Note that in the Bayesian setting above the regularization parameter α is uniquely determined by the prior distribution and the likelihood functional. However, often only qualitative prior knowledge on the solution is available, but the parameter τ is unknown. Then α has to be determined by a-posteriori parameter choice rules analogous to the discrepancy principle (5.6) for deterministic errors.
Let us discuss a few popular choices of the penalty functional R, which allows the incorporation of prior information on the solution or is simply the negative logarithm of the density of the prior in the above Bayesian setting. For a-priori known sparsity of the solution one should choose the sparsity enforcing penalty R(f) = f 1 . If f is an image with sharp edges, then the total variation seminorm is a good choice of R. However, we point out that for the total variation seminorm in a Bayesian setting there exists no straightforward useful infinite dimensional limit (see [16]). Bayesian prior modelling in infinite dimensional settings is often considerably more involved.
If f is a probability density, a frequent choice of the penalty functional is R(f) = KL(f, f 0 ), which naturally enforces nonnegativity of the solution because of the logarithm. Alternatively, more general inequality constraints N (f) ≤ 0 for some function N can be be incorporated into the penalty function by replacing R by

Implementation
In this paragraph we will discuss several possibilities to compute the minimizerf α of (5.9) for a linear forward operator denoted by F = T . In the case of quadratic Tikhonov regularization it follows from the first order optimality conditions that the Tikhonov functional J α has the unique minimizer for all α > 0. So in order to compute the regularized solution we have to solve the linear system of equations Solving this directly by for example Gauss-Jordan elimination requires O n 3 operations and we have to store the full matrix A. In imaging applications n is typically the number of pixels or voxels, which can be so large that storing A is impossible. Therefore, we have to resort to iterative methods which access A only via matrixvector products. Such matrix-vector products can often be implemented efficiently without setting up the matrix A, e.g. by the fast Fourier transform (FFT) or by solving a partial differential equation. As A is positive definite, the most common method for solving Af = b is the conjugate gradient (CG) method: If A = T * T + αI as in Tikhonov regularization, the stopping criterion s = 0 may be replaced by for some tolerance parameter TOL > 0. It can be shown that s l = b − Af l for all l. As A −1 ≤ 1/α, this guarantees the error bound f − f L ≤ TOL to the exact minimumf = A −1 b of the Tikhonov functional.
In the case of more general data fidelity terms S g and penalty terms R one can use a primal-dual algorithm suggested by Chambolle and Pock [3]. To formulate this algorithm, recall that for a functional F : R n → R ∪ {∞} the conjugate functional is given by If F is convex and continuous, then F * * = F. For more information on this and other basic notions of convex analysis used in the following we refer, e.g. to [17]. The algorithm is based on the saddle point formulation Note that an analytic computation of the maximum leads to the original problem (5.9) whereas a computation of the minimum leads to the dual problem The algorithm requires the computation of so called proximity operators (see Chap. 6). For a functional G : R n → R and a scalar λ > 0 the proximity operator prox G,λ : R n → R n is defined by For many popular choices of S g and R the proximity operator can either be calculated directly via a closed expression or there are efficient algorithms for their computation.
To give a simple example, for G(x) = 1 2 x 2 2 one can calculate the proximity operator just by (5.10) to be (1 + λ) −1 z.
One also needs to evaluate the proximity operator of the Fenchel convex conjugate G * , which can be done by Moreau's identity For constant parameters τ k X , τ k Y > 0 and θ k = 1 it has been shown [3] that f k converges to a solution of (5.9) and p k converges to a solution of the corresponding dual problem. Under certain assumptions on S g and R special choices for the parameters τ k X , τ k Y > 0 and θ k will speed up the convergence.
To compute the minimizer of (5.9) with additional inequality constraints one can apply semismooth Newton methods for which we refer to [19].

Iterative Regularization
Whereas in the previous subsection we have discussed iterative methods to compute the minimum of a generalized Tikhonov functional, in the following we will discuss iterative methods for the solution of F(f) = g without prior regularization. Typically the choice of the stopping index plays the role of the choice of the regularization parameter α. A motivation for the use of iterative regularization method is the fact that the Tikhonov functional is not convex in general for nonlinear operators. Therefore, it cannot be guaranteed that an approximation to the global minimum of the Tikhonov functional can be computed. For further information on iterative regularization methods we refer to the monograph [15].

Landweber Iteration
Landweber iteration can be derived as a method of steepest descent for the cost functional J (f) = 1 2 F(f) − g 2 2 . As the direction of steepest descent is given by the negative gradient −J (f) = −F (f) * (F(f) − g), this leads to the iteration formula with a step size parameter μ. This parameter should be chosen such that μ F (f) * F (f) < 1 for all f. Since a too small choice of μ slows down convergence considerably, it is advisable to compute the operator norm of F (f) * F (f) by a few iterations of the power method. Under certain conditions on the operator F it has been shown in [7] in a Hilbert space setting that Landweber iteration with the discrepancy principle as stopping rule is a regularization method in a sense analogous to Definition 5.1, i.e. the worst case error tends to 0 with the noise level for a sufficiently good initial guess.

Regularized Newton Methods
Although Landweber iteration often makes good progress in the first few iterations, asymptotic convergence is very slow. Faster convergence may be expected from Newton-type methods, which solve a linear system of equations or some minimization problem using the first order Taylor approximation (5.14) around a current iterate f k . Plugging this approximation into the quadratic Tikhonov functional and using the last iterate as initial guess leads to the Levenberg-Marquardtalgorithm For a convergence analysis we refer to [6]. The minimization problems can be solved efficiently by Algorithm 5.1 without the need to set up the full Jacobi matrices F (f k ) in each step. Newton-type methods converge considerably faster than Landweber iteration. In fact the number of Landweber steps which is necessary to achieve an accuracy comparable to k Newton steps increases exponentially with k (cf. [4]). On the other hand, each iteration step is more expensive. Which of the two methods is more efficient may depend on the size of the noise level. Newton-type methods are typically favorable for small noise levels. Plugging the first order Taylor approximation (5.14) into the generalized Tikhonov functional (5.8) leads to the iteration formula If S g obs (g) = 1 2 g − g obs 2 and R(f) = 1 2 f − f 0 , this leads the commonly used iteratively regularized Gauss-Newton method where compared to the Levenberg-Marquardt method the penalty term is replaced by α k 2 f − f 0 2 2 . Note that if S g obs and R are convex, a convex optimization problems has to be solved in each Newton step, which can be done, e.g., by Algorithm 5.2. For a convergence analysis of (5.15) including the case of Poisson data we refer to [12].

General Error Bounds for Variational Regularization
Under the deterministic noise model (5.4) we are now looking for error bounds of the form for the Tikhonov estimatorf α . In a second step the right hand side may be minimized over α. As discussed in Sect. 5.1.4 such estimates can only be obtained under additional conditions on f † , which are called source conditions. There are several forms of such conditions. Nowadays, starting with [8] such conditions are often formulated in the form of variational inequalities. For the sake of simplicity we confine ourselves to the case of Hilbert spaces with quadratic functionals R and S g obs here, although the concepts can be generalized with little additional effort to general convex R and S g obs . For a concave, monotonically increasing and continuous function ψ Such conditions are not easy to interpret at first sight, and we will come back to this in Sect. 5.3.2. However, they have been shown to be necessary for certain convergence rates of Tikhonov regularization and other regularization methods (see [11]), and sufficiency can be shown quite easily: (5.17), then the Tikhonov estimatorf α in (5.8) satisfies the error estimate where ψ * (t) := sup s≥0 st − ψ(t) denotes the conjugate function (see (5.11)).
Proof By definition off α and our noise model we have so together with our assumption (5.17) we find By the parallelogram law we have for all x, y, z ∈ X that Apply this with x = F(f α ), y = F( f † ) and z = g obs to find so that finally we have The two most commonly used types of functions ψ in the literature, are referred to as Hölder and logarithmic source function, respectively. For these functions we obtain Note that the two terms on the right hand side of (5.18) correspond to the error splitting (5.5). The following theorem gives an optimal choice of α balancing these two terms:

Theorem 5.4 If under the assumptions of Theorem 5.3 ψ is differentiable, the infimum of the right hand side of (5.18) is attained atᾱ(δ) if and only if
Proof Note from the definition of (−ψ) * that (−ψ) * (t * ) ≥ tt * + ψ(t) for all t ≥ 0, t * ∈ R. Further equality holds true if and only if t * = −ψ (t), as for this choice of t the concave functiont →tt * + ψ(t) attains its unique maximum at t and thus in particular tt * + ψ(t) ≥ (−ψ) * (t * ). Therefore we have with t = 4δ 2 and t * = − 1 4α that It can be shown that the same type of error bound can be obtained by a version of the discrepancy principle (5.6) which does not require knowledge of the function ψ describing abstract smoothness of the unknown solution f † [2]. This is an advantage in practice, because such knowledge is often unrealistic.

Connection to Stability Estimates
Variational source conditions (5.17) are closely related to so called stability estimates.
In fact if (5.17) holds true for all f † ∈ K ⊂ X, then all f 1 , f 2 ∈ K satisfy the stability estimate with the same function ψ, since one of the terms ± f 1 2 − f 2 2 will be nonpositive. There exists a considerable literature on such stability estimates (see e.g. [1,14]). However it is unclear if stability estimates also imply variational source conditions as two difficulties have to be overcome. Firstly the term f 2 − f † 2 might be negative and secondly one would have to extend the estimate from the set K to the whole space X.

General Strategy for the Verification of Variational Source Conditions
In general the rate at which the error of reconstruction methods converges to 0 as the noise level tends to 0 in inverse problems depends on two factors: The smoothness of the solution f † and the degree of ill-posedness of F. We will describe both in terms of a family of finite dimensional subspaces V n ⊂ X or the corresponding orthogonal projections P n : X → V n . The smoothness of f † will be measured by how fast the best approximations P n f † in V n converge to f † : Inequalities of this type are called Bernstein inequalities, and they are well studied for many types of subspaces V n such as spaces of polynomials, trigonometric polyno-mials, splines, or finite elements. We will illustrate this for the case of trigonometric polynomials below.
Concerning the degree of ill-posedness, recall that any linear mapping on a finite dimensional space is continuous. Therefore, a linear, injective operator T restricted to a finite dimensional space V n has a continuous inverse (T | V n ) −1 defined on T (V n ). However, the norm of these operators will grow with n, and the rate of growth may be used to measure the degree of ill-posedness. In the nonlinear case we may look at Lipschitz constants σ n such that P n f † − P n f ≤ σ n F(P n f † ) − F(P n f ) . However, to obtain optimal results it turns out that we need estimates of inner products of P n f † − P n f with f † . Moreover, on the right hand side we have to deal with The growth rate of σ n describes what we will call local degree of ill-posedness of F at f † .
Theorem 5.5 Let X and Y be Hilbert spaces and suppose that there exists a sequence of projection operators P n : X → X and sequences (κ n ) n∈N , (σ n ) n∈N of positive numbers such that (5.19) holds true for all n ∈ N. Then f † fulfills a variational source condition (5.17) with the concave, continuous, increasing function 20) which satisfies ψ(0) = 0.
Proof By straightforward computations we see that the variational source condition (5.17) has the equivalent form Using (5.19a), (5.19c), and the Cauchy-Schwarz inequality we get for each n ∈ N that Taking the infimum over the right hand side with respect to n ∈ N yields (5.20) with As ψ is defined by an infimum over concave and increasing functions, it is also increasing and concave. Moreover, (5.19b) implies ψ(0) = 0.

Example: Numerical Differentiation
Recall that the trigonometric monomials e n ( Therefore, the norm L 2 for k ∈ N 0 , but it also allows to measure non-integer degrees of smoothness of f .
We choose P n as the orthogonal projection Suppose that the kth distributional derivative belongs to L 2 ([0, 1]). Then H s which shows that (5.19a) and (5.19b) are satisfied with κ n = (2πn) −k . Moreover, we have that for s ∈ (0, 1], so (5.19c) is satisfied with σ n := (2πm) 1 It can be shown that this rate is optimal in the sense that there exists no reconstruction method R for which

Example: Fluorescence Microscopy
Similarly one can proceed for the example of fluorescence microscopy. As one has to work here with L 2 spaces rather than L 2 spaces, the Sobolev norm is defined by Assuming that the convolution kernel is a-times smoothing (a > 0) in the sense that F( f ) H a ∼ f H 0 , which is equivalent to the existence of two constant 0 < c < C such that the Fourier transform of the convolution kernel k fulfills one can show that f † H s < ∞ for s ∈ (0, a] implies a variational source condition with ψ(t) ∼ t s s+a and an error bound Again this estimate is optimal in the sense explained above.

Extensions
Variational source conditions with a given Hölder source function actually hold true on a slightly larger set. In the typical situation where the marker density of the investigated specimen is constant (or smooth) up to jumps, then it fulfills the same variational source condition with s = 1/2, although f † ∈ H s if and only if s < 1/2. The sets on which a variational souce condition is satisfied can be characterized in terms of Besov spaces B s 2,∞ , and bounded subsets of such spaces are also the largest sets on which a Hölder-type error bound like (5.21) and (5.22) are satisfied with uniform constants (see [11]).
In the case where the convolution kernel is infinitely smoothing, e.g. if the kernel is a Gaussian, then we cannot expect to get a variational source condition with a Hölder source function under Sobolev smoothness assumptions. Instead one obtains logarithmic source functions ψ log p as introduced above, which will again be optimal and lead to very slowly decaying error estimates as δ 0 [11]. Note that the rates (5.21) and (5.22) are restricted to smoothness indices s ∈ (0, 1] and s ∈ (0, a], respectively. This restriction to low order rates is a well-known shortcoming of variational source conditions. Higher order rates can be obtained by imposing a variational source condition on the dual problem (5.12), which can again be again by verified by Theorem 5.5 (see [5,18]).
With some small modifications the strategy in Theorem 5.5 can be extended to Banach space settings [9,21] and nonlinear forward operators F, in particular those arising in inverse scattering problems [10,20].

Error Bounds for Poisson Data
We already briefly discussed discrete versions of inverse problems with Poisson in Sect. 5.2.1. Such problems arise in many photonic imaging modalities such as fluorescence microscopy, coherent x-ray imaging, positron emission tomography, but also electron microscopy. In the following we briefly discuss a continuous setting for such problems.
We consider a forward operator F : X → Y that maps an unknown sample f † ∈ X a photon density g † ∈ Y = L 1 (M) generated by f † on some measurement manifold M ⊂ R d . The given data are modeled by a Poisson process with density tg † . Here {x 1 , . . . , x N } ⊂ M denote the positions of the detected photons and t > 0 can be interpreted as exposure time. Note that t ∼ E(N ), i.e. the exposure time is proportional to the number of detected photons. Now to discuss error bounds we first need some notion of noise level. But what is the "noise" in our setting? Our data G t do not belong to Y = L 1 (M) and the "noise" is not additive. However, it follows from the properties of Poisson processes that and the variance of 1 t G t , h is proportional to 1 t . This suggests that 1 √ t plays the role of the noise level. More precisely, it is possible to derive concentration inequalities where the distance function d is defined by a negative Besov norm (see [13,22] for similar results). As a reconstruction method we consider generalized Tikhonov regularization as in (5.8) with S G t given by a negative (quasi-)log-likelihood corresponding to the Poisson data. As discussed in Sect. 5.2 this amounts to taking the Kullback-Leibler distance as data fidelity term in the finite dimensional case, and in particular in the implementations of this method. (Sometimes a small shift is introduced in the Kullback-Leibler divergence to "'regularize" this term.) By assuming the VSC (5.17) with ψ( F(f † ) − F(f) 2 ) replaced by ψ KL(F(f † ), F(f))) and an optimal parameter choiceᾱ one can than show the following error estimate in expectation Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made. The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.