Inverse problems with inexact forward operator: iterative regularization and application in dynamic imaging

The classic regularization theory for solving inverse problems is built on the assumption that the forward operator perfectly represents the underlying physical model of the data acquisition. However, in many applications, for instance in microscopy or magnetic particle imaging, this is not the case. Another important example represent dynamic inverse problems, where changes of the searched-for quantity during data collection can be interpreted as model uncertainties. In this article, we propose a regularization strategy for linear inverse problems with inexact forward operator based on sequential subspace optimization methods (SESOP). In order to account for local modelling errors, we suggest to combine SESOP with the Kaczmarz’ method. We study convergence and regularization properties of the proposed method and discuss several practical realizations. Relevance and performance of our approach are evaluated at simulated data from dynamic computerized tomography with various dynamic scenarios.


Introduction
An inverse problem is concerned with extracting a searched-for quantity f ∈ X from measured data g ∈ Y where g = A f with forward operator A : X → Y modelling the data acquisition process. A characteristic property of an inverse problem is its ill-posedness: in general, a small noise in the measurement can cause a large error in the solution, and therefore, suitable regularization schemes are required.
Due to simplifications in the underlying physical model of the data acquisition process or due to calibration errors, the forward operator is usually not exactly known either. A prominent example arises for instance in the context of microscopic imaging, where the forward operator depends on the point spread function of the microscope which has to be determined experimentally. In magnetic particle imaging, the model-based determination of the system function is an unsolved problem which necessitates an experimental calibration step. Our main motivation to study inverse problems with inexact forward operator is to solve dynamic inverse problems, where changes of the searched-for quantity during data collection can be interpreted as model uncertainties.

Motivating example: dynamic inverse problems
In many practical applications, collecting the measurement g takes a considerable amount of time. In computerized tomography (CT), for instance, it takes time to rotate the radiation source around the investigated object. Within the standard formulation A f = g, this can be expressed by specifying the data variable, i.e., the equation more precisely reads with time interval [0, T] ⊂ R covering the time period required for collecting the data. Ω Y could be considered, e.g., as bounded subset of R m . To simplify the modelling step, A is typically derived under the assumption that the searched-for quantity f is stationary, i.e., independent of time. However, in many applications, the investigated quantity changes during the collection of the data, for instance in medical imaging due to patient or organ motion. In non-destructive testing, the dynamic scenario arises for example when visualizing driven liquid fronts [1] or performing elasticity experiments during the scan to determine material parameters [17]. The respective dynamic data then are inconsistent with the standard model A. Thus, simply applying a standard inversion procedure for A to such a dynamic data set causes motion artifacts in the reconstruction result which degrade the reconstruction quality. As a consequence, specific regularization methods taking into account the dynamic behaviour of the specimen have to be developed for such dynamic inverse problems.
Typically, the individual states of the specimen show a strong temporal correlation. They can be correlated to one reference configuration f via a deformation model Γ : [0, T] × R n → R n . In this case, the exact forward model for the dynamic case is given by the operator where A stat refers to the forward model of the static scenario. Thus, in the dynamic case, we have to solve the equation A Γ f (t, y) = g(t, y) for all (t, y) ∈ [0, T] × Ω Y , see, e.g., [12,18]. In practice, the exact deformation Γ, and therefore the exact dynamic forward operator A Γ , will in general be unknown. This either requires a delicate calibration step (socalled motion estimation) prior to the reconstruction [11,19,22] or the simultaneous recovery of Γ and f which is typically computationally demanding [2,3]. For specific applications and special deformations that preserve the underlying data acquisition geometry, exact analytic reconstruction methods have been derived, especially in CT, where this type of motion includes affine deformations, [6,7,10]. In this case, techniques for rebinning the measured data to make them feasible for standard reconstruction methods are proposed as well, [4,34].
To avoid the challenging task of an explicit motion estimation, several approaches incorporating less strong prior information on the solution have been studied as well. They typically invoke some spatio-temporal regularization term, such as, e.g., temporal smoothness [28,29]. In [5], a computational method in a Bayesian framework has been proposed, where spatio-temporal covariance matrices are considered. All-at-once and reduced iterative methods to solve time-dependent parameter identification problems are described and analyzed in [16,25]. Approaches for solving non-linear dynamic tomography problems have been studied in the context of electrical impedance tomography, for instance using an oppositional biogeography-based optimization technique to reconstruct organ boundaries [26].
In this article, we propose a different strategy which requires neither explicitly known deformation fields nor explicit spatio-temporal priors on the solution: instead of the exact but unknown operator A Γ , utilize a simplified model, such as the inexact operator A stat , and apply a regularization strategy that accounts for the respective model inexactness.

Inverse problems with inexact forward operator
The aforementioned examples illustrate that uncertainties in the forward model are highly relevant in practice. This inexactness needs to be taken into account within the reconstruction step in order to stably determine f from noisy data g δ with while the reconstruction procedure only uses an available approximate model A η with We refer to the model error (2) in the following as the global model inexactness.
In the literature, a few methods have been proposed to deal with an imperfection of the forward model. A well-known technique to compensate modelling errors in the reconstruction process is total least squares regularization [8]. An optimization approach based on partiallyordered spaces was proposed and evaluated for the deblurring problem in [21]. In [20], a reconstruction method for magnetic particle imaging was presented which computes simultaneously the searched-for quantity and deviations of the simplified discrete forward model. Uncertainties in the forward operator are often considered and dealt with in the context of statistical inverse problems, see, e.g., [33].
In this article, we propose a regularization strategy for linear inverse problems with inexact forward operator based on sequential subspace optimization methods. Sequential subspace optimization (or short SESOP) belongs to the class of iterative reconstruction techniques [9,15,23,31,32,35,36]. The underlying idea is to sequentially project from some initial value onto suitable subspaces that contain the searched-for solution. This method has been successfully applied, for example to solve an integral equation of the first kind [30] or in parameter identification [32,35], particularly in terahertz tomography [37]. The goal in our case is to choose these subspaces in accordance to the model inexactness such that the sequence of iterates generated by the imperfect model A η converges to a solution of A f = g.
The above framework can then be used to solve dynamic inverse problems, where the dynamic operator A Γ reprises the role of the unknown exact model A, i.e., we set A Γ A. As simplified operator, we could choose for instance the static forward operator, i.e., A stat A η , provided we have at least an estimate η on the overall model error.

Local inexactness and semi-discrete framework
In most cases, a global error such as (2) will quantify the model deviation too pessimistically. As an example, consider the case of periodic motion, such as respiratory or cardiac motion. Due to the periodicity, it might hold A f (t, ·) = A η f (t, ·) for certain time instances, i.e., locally, the model error is small. In addition, local deformations, e.g. cardiac motion, affect only a small region of the body and therefore, it can hold A f (t, y) = A η f (t, y) for some y ∈ Ω Y as well. Thus, we want to consider a framework that takes local model errors depending on the data variable into account.
For this purpose, we focus on a semi-discrete framework with a finite sequence of linear operators A k,l in this article. On the one hand, this framework accounts for the discrete nature of measured data in practice, see also, e.g., [27]. On the other hand, we can then consider a sequence of imperfect forward models A η k,l with This allows the level of inexactness, characterized by η k,l , to vary for each data point. Thus, such local error estimates will typically be less pessimistic than a global estimate (2) and further allow to incorporate prior knowledge on the motion (e.g. on periodicity). The indices k and l may refer, for example, to time instances and detector points in dynamic imaging.

Outline of the article
The article is organized as follows. Section 2 presents the basic principle of SESOP methods in the continuous setting with exact model operator and noise-free data. We then provide the mathematical framework of our semi-discrete setting. In particular, we suggest to combine SESOP with Kaczmarz' method and discuss the resulting advantages regarding practical applications. Section 2 concludes with the formulation of our regularized version of the SESOP-Kaczmarz method for inexact forward operators. Its convergence and regularization properties are then studied in section 3. In Section 4, we discuss possible variations of our proposed algorithms. Section 5 is devoted to the application in dynamic CT. In particular, our proposed algorithms are evaluated at different simulated CT data for various dynamic scenarios.

A SESOP-Kaczmarz method for linear inverse problems with inexact forward operator
This section is dedicated to introduce the basic tools that are needed to formulate the SESOP algorithms and their combinations with Kaczmarz' method. Detailed descriptions of SESOP techniques with an exact forward operator can be found in the literature, e.g., [30-32, 35, 36].

Basic principles of sequential subspace optimization
In the following, we provide a brief introduction to SESOP methods in the continuous setting with exact model operator and noise-free data. Let X, · X and Y, · Y denote Hilbert spaces equipped with their respective norms. By ·, · X we refer to the inner product on X. If there is no ambiguity, we will omit the indices on the norms and inner products to simplify the notation. Further, let A : X → Y be a bounded linear operator with for some Λ A > 0. By A * we denote the adjoint of A. The space is called the range of A and the space is called the null space of A. For given g ∈ Y, the set denotes the solution set of the operator equation A f = g.
The underlying idea of SESOP is to create a sequence of iterates f n ∈ X that converges to a solution f ∈ M A,g by sequentially projecting a starting value f 0 ∈ X onto suitable convex subsets of X. The following definition introduces specific types of subsets of X which play an important role in this context. An illustration of the concept of stripes is given in figure 1. For arbitrary w ∈ Y, the solution set M A,g is contained in the hyperplane Indeed, we see that for each z ∈ M A,g we have Next, we introduce the concept of metric projections, which play a crucial role for SESOP.

Definition 2.2.
Let C = ∅ be a convex subset of X. Then the metric projection P C (x) of x ∈ X onto C is defined by In Hilbert spaces, the metric projection P H(u,α) (x) of x ∈ X onto a hyperplane H(u, α) coincides with the orthogonal projection and can be expressed as The metric projection onto a halfspace or a stripe corresponds to the metric projection onto the respective bounding hyperplane. We now can formulate the nth iteration step (n > 0) of the SESOP method for solving A f = g, see, e.g., [32]: given a chosen finite index set I n and search directions A * w n,i with chosen w n,i ∈ Y, i ∈ I n , compute the iterate f n+1 as the metric projection of the previous iterate f n onto the intersection of hyperplanes In [31] it has been shown that the metric projection P H ( f) of f ∈ X onto the (non-empty) intersection wheret := (t i ) i∈I ∈ R |I| minimizes the function If {u i : i ∈ I} is linearly independent, the function h is strictly convex and thus,t is unique, see also [31]. Thus, the iterates (6) of the SESOP method are given by Under suitable conditions on the index set I n , the search directions A * w n,i and the parameters t n,i , the sequence { f n } n∈N converges to f ∈ M A,g , see [32,36]. We conclude this subsection by outlining how SESOP can be adapted to a regularization method in case only noisy data g δ with noise level δ (i.e. g − g δ δ) and an inexact version A η of the forward operator A with inexactness η (i.e. A − A η η) are given. The basic idea consists in replacing the hyperplanes within the projection step by stripes whose width is chosen in accordance to the noise and uncertainty level.
This idea is based on the property that since for all z ∈ M A,g with z ρ we have This approach has been studied so far for the case of noisy data and exact forward operator (i.e. η = 0). In particular, combining the respective projection algorithm with the discrepancy principle yields a finite stopping index n * and the respective iterate f n * constitutes a regularized solution of A f = g, see again [32,36]. One goal of this article is to extend the respective convergence and regularization results to the setting with inexact forward operator. In order to allow incorporating local estimates for the modelling errors, we focus in the following on the more general semi-discrete framework motivated in section 1.

A SESOP-Kaczmarz method for semi-discrete inverse problems
for some Λ A > 0. Their Hilbert space adjoints A * k,l as well as their range and null space are defined in analogy to the ones in the continuous setting. The solution set of the semi-discrete inverse problem is given by for some solution f + of (12). The results presented in this article can be directly transferred to this more general setting. However, to keep the notation as simple as possible, we consider in the following the setting where L is independent of k.
In order to solve the semi-discrete inverse problem (12), we propose to combine the SESOP method with Kaczmarz' method [24].

Remark 2.4.
A combination of Kaczmarz' method with SESOP allows to incorporate local properties of the inverse problem: on the one hand, SESOP is designed in such a way that the local properties of the inverse problem, i.e., particularly local modelling errors, can be easily included in the regularization. On the other hand, Kaczmarz' method allows to combine the respective subproblems.
We also remark that, in view of a comparison to the Landweber method, the additional computational effort due to solving the required optimization problem in the SESOP iteration is relatively low compared to the reduction in computation time due to this regulation of the step width and using multiple search directions, see also, e.g., [30,35,37] for a comparison.
To illustrate the basic idea, we first consider a system of linear operator equations with only one index set. The standard Kaczmarz' method reads in this case where [n] := n modK and f 0 ∈ X is an initial value. Here we use the notation from [13,14], where the Landweber-Kaczmarz method for nonlinear inverse problems is investigated.
The projection step can now be replaced by a metric projection onto the intersection H n := i∈I n H n,i of the hyperplanes H n,i : . The corresponding SESOP-Kaczmarz method then reads Here, t n := (t n,i ) i∈I n is calculated as a minimizer of the respective function h n from (8) with u n,i := A * [i] w n,i and α n,i := w n,i , g [i] . Each parameter t n,i can be interpreted as a regulation of the step width in the direction of the search direction u n,i .
Regarding our semi-discrete setting, the two discretization indices k and l allow for multiple versions of the basic SESOP-Kaczmarz iteration (14). For instance, the Kaczmarz iterations can be executed with respect to both k and l. Depending on the nature of the underlying inverse problem, it might be convenient to use search directions that are averaged with respect to one of the indices, for instance resulting from summing up all search directions with respect to the second discretization index.
Respectively, the roles of k and l could be interchanged. A detailed overview of these different versions can be found in section 4. Please note that the continuous case, i.e., for the operator equation A f = g, is included in the above framework since it corresponds to a setting with Throughout sections 2 and 3, we focus on the SESOP-Kaczmarz algorithm with searchdirections (15) for the semi-discrete case, since all the above mentioned versions can be treated in this framework, see the discussion in section 4.
In order to be able to provide a solution of A k,l f = g k,l , the hyperplanes associated to the search directions (15) have to include the solution set M sd A,g (13). According to the following lemma, this is indeed the case.

Lemma 2.5. For each iteration index n, the hyperplanes
with l ∈ L as well as the hyperplanes with w l n,i ∈ Y [i],l for all i ∈ I n , n ∈ N, and l ∈ L contain the solution set M sd A,g . Proof. The proof is analogous to (10).

Regularizing SESOP-Kaczmarz iteration for inexact forward operators
In this section, we propose a regularizing SESOP-Kaczmarz iteration that is suited for linear semi-discrete inverse problems (12) in Hilbert spaces with inexact forward operator and noisy data. Let A η k,l ∈ L(X, Y k,l ) denote an inexact version of the forward operator A k,l with inexactness Further, we consider noisy data g δ k,l ∈ Y k,l with noise level Remark 2.6. Please note that we label quantities that are influenced by the inexactness η = (η k,l ) k∈K,l∈L or the noise level δ = (δ k,l ) k∈K,l∈L by adding the index η or δ. For example, we assume that the forward operator A k,l is inexact with inexactness η k,l , yet we will only write In analogy to (13), we denote the set of solutions whose norm is bounded by ρ > 0 by For each element z ∈ M sd,ρ A,g , we obtain the estimate We propose to regularize the SESOP-Kaczmarz algorithm with search directions (15), see section 2.2, by replacing the hyperplanes H n,i by stripes with width in accordance to δ, respectively η, and by combining the algorithm with a discrete version of the discrepancy principle. The details are formulated in the following algorithm which is referred to as RESESOP-Kaczmarz. (20) by projecting the current iterate f η,δ n onto the intersection of the stripes with the parameters Stop iterating, as soon as f η,δ n = f η,δ n−K and n = 0 mod K. Remark 2.8. Algorithm 2.7 includes an adaption of the discrepancy principle for our semidiscrete setting, which considers the discrepancies of all subproblems separately to increase the accuracy of the reconstruction. Definition 2.9. The iterates f η,δ n , n ∈ N, in algorithm 2.7 are called auxiliary iterates. Since algorithm 2.7 is stopped when n = r · K for some integer r ∈ N (i.e., after a full sweep through all subproblems), we definef η,δ r := f η,δ rK and call {f η,δ r } r∈N the sequence of full iterates. Remark 2.10. With δ = 0 and η = 0, our RESESOP-Kaczmarz algorithm coincides with SESOP-Kaczmarz proposed in section 2.2 for the case of exact data and exact forward operator.

Convergence and regularization
In this section, we study convergence and regularization properties of our RESESOP-Kaczmarz algorithm 2.7. This requires in particular imposing suitable conditions on the index sets I η,δ n and the search directions.
Regarding the index set I η,δ n , we postulate that where N ∈ N\{0} is fixed.
Regarding the search directions, for each n ∈ N, we define with w η,δ,l n := w η,δ,l n,n := In other words, we define the search direction u η,δ n,n as the gradient of the functional [n],l 2 evaluated in the current iterate f η,δ n . The section is now organized as follows. We first prove in section 3.1, that the sequence of iterates for the exact forward operator and exact data converges to a solution f + ∈ M sd A,g . We then turn to the case of noisy data and inexact forward operator and show that the respective iteration scheme provides a regularized solution to the semi-discrete inverse problem A k,l f = g k,l , (k, l) ∈ K × L.

Convergence
In order to prove the convergence of algorithm 2.7 for an unperturbed forward operator and exact data, we establish in a first step the validity of a descent property.
holds for all z ∈ M sd A,g .

Proof.
For n ∈ N we setH n := H(u n , α n ) and Using f n+1 = P H n ( f n ) with H n := i∈I n H(u n,i , α n,i ) ⊂H n , we obtain Thanks to (11), we estimate The assertion follows with C := 1 (Λ A L) 2 . As a consequence of lemma 3.1, we obtain the descent property for the full iteratesf r := f rK , r ∈ N, i.e., the sequences { z − f n } n∈N and z −f r r∈N are monotonically decreasing. In addition, since z ∈ M sd A,g is fixed, we deduce from (24) that for a constant ρ > z − f 0 > 0. This yields the boundedness of the sequences { f n } n∈N and {f r } r∈N . The validity of these descent properties are the starting point for the following theorem.

Theorem 3.2.
The sequences { f n } n∈N and f r r∈N , generated by algorithm 2.7 with (22) and (23), each possess a subsequence that weakly converges to an element of M sd A,g . Proof. We consider the sequence { f n } n∈N of auxiliary iterates. From (24) we obtain for an arbitrary but fixed integer N > 0. We now let N → ∞. Since the right-hand side of the above estimate is independent of N, we deduce that the series ∞ n=0 w n 2 is absolutely convergent, and thus {w n } n∈N is a null sequence. Since { f n } n∈N is bounded according to (26), there exists a weakly converging subsequence of { f n } n∈N with weak limitf . The same holds for each subsequence, such that we have f n f . In particular, we havef r f . According to where we have used the continuity of A k,l and the weak lower semicontinuity of the norm, we havef ∈ M sd A,g . from an initial value f 0 ∈ X. If the optimization parameters |t n,i | t with i ∈ I n , n ∈ N are bounded by some t > 0, then { f n } n∈N and f r r∈N converge strongly to a solution f + , and The proof is given in appendix A.

Regularization
In the previous section we have shown that the SESOP-Kaczmarz iterates converge to a solution of the semi-discrete inverse problem (12) with exact data. The solution is the metric projection of the initial value onto the solution set M sd A,g . Based on these findings, we now want to proceed to noisy data, i.e., we show that the RESESOP-Kaczmarz iteration, i.e., algorithm 2.7, yields a regularized solution of a semidiscrete system (12) if only noisy data g δ k,l and inexact forward operators A η k,l , k ∈ K, l ∈ L, are given. We assume that the noise level can be estimated as in (19) and the inexactness by (18).
We first remark that if the set D η,δ j as defined in algorithm 2.7 is empty for all j = n, . . .  We proceed as in the case of the SESOP-Kaczmarz algorithm with (η, δ) = (0, 0) and show that the RESESOP-Kaczmarz algorithm yields a regularized solution f η,δ n * , i.e., we consider the iterates f η,δ n , n ∈ N. Due to definition 3.4, the properties of f η,δ r r∈N will be inherited from f η,δ n n∈N . Again we postulate that the conditions (22) and (23) are fulfilled. We can express the search direction u η,δ n := u η,δ n,n by The parameters of the stripes can be reformulated as These conditions assure that the descent property holds: and z − f η,δ which proves assertion (29). The descent property (30) follows in the same way as in the the unperturbed case, i.e., as in the proof of lemma 3.1, by using the representations of u η,δ n , α η,δ n and ξ η,δ n .
Proof. Let us assume that there is no finite stopping index r * , i.e., D η,δ,K rK = ∅ for all r ∈ N. Then, for each j ∈ D η,δ,K rK , there is at least one l ∈ L such that l ∈ D η,δ j and We set We plug relation (32) into the descent property from remark 3.6 and obtain In particular we have used our definition ofw η,δ,l j for the last estimate. Together with In a similar way as in the proof of theorem 3.2 we find that This is a contradiction to (32).
In the following we want to show that the RESESOP-Kaczmarz iteration is a regularization method in case the conditions (22) and for all i ∈ I η,δ n are fulfilled. We begin by showing that the iterates f η,δ n , n ∈ N, depend continuously on the noise level δ and the inexactness η.  (22), (33). By { f n } n∈N we denote the sequence generated by the SESOP-Kaczmarz algorithm 2.7 for (η, δ) = (0, 0) with (22), (27). If the set of search directions Proof. See appendix B.
We now have assembled all necessary statements to prove the main result for the RESESOP-Kaczmarz algorithm 2.7: if suitable conditions are fulfilled, the RESESOP-Kaczmarz algorithm yields a regularization method for semi-discrete linear inverse problems (12). Theorem 3.9. Let the conditions from lemma 3.8 be satisfied. We additionally postulate that the optimization parameters t η,δ n,i fulfill |t η,δ n,i | < t for all i ∈ I η,δ n , n ∈ N. Then the RESESOP-Kaczmarz algorithm 2.7 together with the discrepancy principle from definition 3.4 as a stopping criterion yields a regularized solution f η,δ r * (η,δ) ∈ X of the system (12) and The proof is given in appendix C.

Practical realizations of the RESESOP-Kaczmarz algorithm
The semi-discrete setting allows a range of modifications of our RESESOP-Kaczmarz algorithm 2.7. We recall that algorithm 2.7 first executes an averaging of the search directions u l [n],i : l ∈ L before executing a Kaczmarz loop with respect to the index k.
Evidently, the roles of k and l could be interchanged, resulting in an algorithmic version that performs the Kaczmarz projection with respect to the index l and averages regarding the index k. Alternatively, the Kaczmarz projection could be executed with respect to both k and l, Kaczmarz iteration regarding k and averaging regarding l K L Algorithm V4 Kaczmarz iteration regarding l and averaging regarding k L K omitting any averaging step. In the same spirit, it is possible to average with respect to both k and l.
The respective algorithms can be inferred from our general formulation in algorithm 2.7. In order to specify this, we introduce a separate notation for the index sets used in algorithm 2.7 for the Kaczmarz projection, K Alg , and for the averaging, L Alg , respectively. In case of algorithm 2.7, these sets coincide with the respective index sets K and L given by the semi-discrete inverse problem. With this additional notation, we can infer the modifications described in table 1.
The performance of these four versions of the RESESOP-Kaczmarz regularization is evaluated in section 5.1 at numerical examples from dynamic CT.
Furthermore, our framework includes the setting of an inverse problem which is discretized with respect to one data variable only. This is of particular interest in the context of dynamic inverse problems if we want to take into account model errors η k only depending on time. In this case, simply set L := {0} in our setting for semi-discrete inverse problems. In particular, in our RESESOP-Kaczmarz framework, algorithms V1 and V3 coincide (we refer to them as RESESOP-Kaczmarz iteration for a single discretization parameter), as well as algorithms V2 and V4 (referred to as averaged RESESOP iteration for a single discretization parameter in this case). Remark 4.1. Algorithm V2 with K = {0} and η = 0 has been used in [37] to solve a parameter identification problem in terahertz tomography, where the set L corresponds to the different detectors.
As remarked earlier, also the continuous setting with global model inexactness A − A η η ∈ [0, ∞) is covered by our framework with K = L = {0}.
In analogy to [30,32,35,36] we can derive a fast algorithm with two search directions per iteration, in which the new iterate is calculated by at most two metric projections onto (intersections of ) bounding hyperplanes. In a Hilbert space setting, these can be formulated explicitly without an evaluation of the underlying minimization problem. This is outlined in more detail in appendix D. For our numerical evaluation in section 5.1, we implemented algorithm 2.7 and the modifications V1-V4 with two search directions as given above.

Application in dynamic computerized tomography (CT)
In the following, we evaluate our regularization scheme for inverse problems with inexact forward operator at the example of dynamic CT.
In CT, a radiation source rotates around the studied object while emitting x-rays through the specimen to a detector panel. Depending on the density in the interior of the object, the photons of the x-ray beam get absorbed. According to the Beer-Lambert law, the observed intensity loss in 2D is modelled by the Radon transform which integrates the attenuation coefficient f of the studied tissue along straight lines (namely the travelling paths of the x-rays). In a semi-discrete setting, this model is given by where θ k ∈ S 1 , k ∈ K, characterizes the different source positions and s l ∈ [−1, 1], l ∈ L, denotes the affected detector points. The goal is to recover f from the measured CT data see [24] for a detailed overview.
In modern CT scanners, the time consuming process of data acquisition is the rotation of the x-ray source around the specimen. Hence, each source position-in our model each unit vector θ k -can be interpreted as a time instance. The standard model above is based on the assumption that the investigated object is stationary during the data collection. Thus, if the specimen shows a dynamic behavior during the collection of the data, it is no longer valid.
Correlating the individual states of the specimen to one reference configuration f via deformation fields Γ k : R 2 → R 2 , k ∈ K, the exact dynamic forward operator is given by While R k,l integrates along the straight lines C k,l := {x ∈ R 2 : s l = x T θ k }, the dynamic model R Γk,l integrates along curved lines If the deformation fields Γ k , k ∈ K are not explicitly known, neither is R Γk,l . Therefore, we propose to solve the dynamic inverse problem R Γk,l f = g k,l by applying our developed regularization strategy which uses an inexact version of R Γk,l , for instance the known static operator R k,l , along with estimates of the local model errors η k,l , i.e., R k,l − R Γk,l η k,l .

Numerical results
We evaluate our algorithms for two different numerical phantoms, one undergoing a periodic and affine deformation, the other performing a non-periodic and non-affine motion. The corresponding dynamic Radon data are computed for K = 300 source positions with L = 451 detector points using the composite trapezoidal rule and bilinear interpolation on a 400 × 400 grid. For the iterative reconstruction scheme, we choose the initial value f 0 = 0, and regarding the discrepancy principle, we use constantly τ = 1, 00 001 for all arising subproblems.  Periodic and affine deformation modelling respiratory and cardiac motion. As first example, we consider a chest phantom performing 30 breathing cycles and 100 heart beats during the data collection. The dynamic behavior of the phantom during one half of a breathing cycle is illustrated in figure 2. The local motion of the heart-without the effect of the breathing-is further highlighted in figure 3.
In terms of explicit formulas, this dynamic behavior is described as follows. We choose the initial state, shown in figure 2(a)), as reference configuration. Each ellipse, labelled consecutively by i = 1, . . . , 7 as in table 2, performs an affine deformation described by Γ i k : With the auxiliary variables a k := k mod 10,   Our goal is to recover an image of the initial state from the dynamic CT data. Since the dynamic forward operators R Γk,l are in general unknown, we only want to use the known static operator R k,l within the reconstruction step. Figure 4 shows the result of the classic SESOP algorithm that does not take the model inexactness into account. The naive use of the static forward operator yields a blurred and  deformed reconstruction, motivating that the model inexactness needs to be incorporated into the reconstruction step.
Therefore, we next apply our derived regularization technique for inverse problems with inexact forward operator. In a first step, we want to use a good approximation for the local model errors η k,l which we determine numerically as the discrepancy of computed static and dynamic data.
The reconstruction results of the individual algorithms presented in section 4 are shown in figure 5. Each of these algorithms provides an image of the interior structure at the reference time without motion artefacts. Thus, the suggested approach to treat the unknown dynamics as uncertainty in the forward model is indeed capable of compensating for the motion.
Drawing a comparison among the results of the different algorithms shows that considering local modelling errors η k,l depending on source as well as detector position provides a higher resolution and sharper images, see figures 5(c)-(f )) in comparison to figures 5(a) and (b)). This is to be expected since this approach incorporates the richest information on the dynamics. In particular, the local motion of the heart is more efficiently compensated for. We further observe that the algorithms involving an averaging part (figures 5(b) and (d), respectively) provide slightly smoother results than the true Kaczmarz algorithms (see figures 5(a) and (c)). When incorporating local model errors and combining Kaczmarz with an averaging part regarding the source position, i.e., regarding different time instances, can result in small artefacts in the region of local deformations (in our case in the upper part of the heart), see figures 5(d) and (f )).
In order to gain further insights into the performance of the four algorithms V1-V4 proposed for incorporating source and detector dependent model uncertainties, figure 6 displays the error of the generated iterates and the ground truth. With an increasing number of iterations, we clearly observe the monotonicity proved in proposition 3.5 and the consecutive remark. However, the slope and the magnitude slightly differ for each algorithm. This is due to the different number of subproblems that are considered for the averaging part. More precisely, the number of subproblems is    The more subproblems are considered, the faster the error decreases in the beginning (since more information are included in the first iteration steps). However, the averaging in each subproblem has an additional smoothing effect resulting in larger errors for higher iteration indices.
In order to test the stability of our method with respect to perturbations in the measurement, we add noise samples uniformly distributed in [−0.02, 0.02] to our dynamic Radon data. Again, we obtain good results with our algorithms as shown in figure 7. This illustrates that our approach can in fact handle both noise in the data and modelling errors. In particular, the additional smoothing property of the algorithms involving an inner summation part yields a better noise suppression than their respective counterparts.
So far, we used uncertainty levels computed numerically as the discrepancy of dynamic and static data. In practice, however, this is not feasible since a static data set will not be available. Determining good estimates for these error levels will be part of future research (see also the discussion in section 6). Nevertheless, we want to test the sensitivity of our method regarding overestimated uncertainty levels η k,l . For this purpose, we add a sample of noise uniformly distributed in [0.02, 0.06] to the computed values, which corresponds to an overestimation of the model errors by up to 5%. The respective reconstruction results shown in figure 8 illustrate that our method still reconstructs a good approximation to the actual reference state. Comparing the results for the different algorithms reveals that the algorithms of type V1 (V3) are more stable with respect to overestimations of the uncertainty.
Non-periodic and non-affine deformation. The dynamic behaviour considered so far corresponds to an affine motion of periodic nature. Due to the high periodicity, the phantom returns to its reference state multiple times during the data collection, and, as a consequence, R k,l − R Γk,l is small for these time instances k. Now, we want to study how our method performs in case of a non-periodic motion. Our test phantom and its (non-affine) dynamic behavior are illustrated in figure 9. The deformation is described by the motion functions We now aim to recover the state of the phantom at time instance k = 149 from the dynamic data. This is achieved by adapting the model errors η k,l such that η 149,l = 0 for all l. Our choice serves two purposes: on the one hand, we want to illustrate that we are not restricted to recovering the initial state. In addition, in case of a monotone motion, considering the initial state will typically be disadvantageous with respect to the magnitude of the model errors, which will increase monotonically up to the maximum value at the end of the scanning. This is illustrated in figure 10 (left) which shows the computed model errors with respect to the initial state as reference configuration. However, if we choose a state during the scanning, then the model errors will first decrease and then increase, see figure 10 (right). In particular, the maximum value is lower than for the initial state which corresponds to overall smaller stripe widths.
In order to recover the selected reference state from the dynamic data, we apply our method with the simplified static operator accounting for the respective modelling errors. Due to the strong non-linearity of the motion which affects the lower left part of the phantom much more than the upper right part, we incorporate modelling errors η k,l depending on both source position and detector point. The respective reconstruction results for our respective algorithms are shown in figure 11. A comparison with the reference state in figure 9(b)) reveals that our approach can also compensate for a non-periodic and highly non-linear dynamic behavior, whereas the naive approach ignoring the model errors results in a distorted visualization of the interior structure, see figure 12.

Conclusion
In this article, we presented a regularization strategy for linear inverse problems with inexact forward operator based on SESOP. The suggested combination with Kaczmarz' method in particular allows to account for local modelling errors. This plays an essential role when dealing with dynamic inverse problems since characteristics of the motion such as periodicity or locality can be incorporated in this way. Our numerical examples from dynamic CT illustrate that our iterative reconstruction method compensates for the object motion during the data acquisition without requiring the explicit deformation fields.
Studying suitable prior information on the motion in more detail will be subject to future research together with the application on real data. Instead of performing our approach with the static operator, it is further possible to use a dynamic forward operator with simplified motion model and to account for modelling and estimation errors regarding the deformation fields. In particular, we expect further insights regarding the relevant information on the model inexactness by taking into account specifics of the individual application. due to (a). Furthermore, using (b), we obtain f n m − f n j , f n j = f n m − f n j , f + + f n m − f n j , f n j − f + → f n m − f n j , f n j − f + for m → ∞.
We apply the iteration from algorithm 2.7, i.e., Together with ϑ j ∈ I ϑ j and w l ϑ j w ϑ j (by definition, see (25)) the above estimate yields whereC := CLNt. Due to the absolute convergence of n∈N W n 2 , the right-hand side of the above inequality converges to 0 as m → ∞. Altogether, this shows that f n j j∈N is indeed a Cauchy sequence, converging to somef ∈ X. The continuity of the norm and the forward operators A k,l , k ∈ K, l ∈ L, yield i.e.,f ∈ M sd A,g . As before, we conclude f n →f for n → ∞ andf r →f for r → ∞. It remains to show thatf = f + := P M sd We finally use proposition 3.7 from [31] for Hilbert spaces to obtaiñ i.e., the resultf of the SESOP-Kaczmarz iteration corresponds to the metric projection of the initial value f 0 onto the solution set M sd A,g .