New Pair of Primal Dual Algorithms for Bregman Iterated Variational Regularization

TL;DR

This paper introduces two new primal-dual algorithms for variational regularization using Bregman distances, providing convergence rates and stability analysis, and demonstrates their effectiveness on tomographic image reconstruction problems.

Contribution

The work develops novel primal-dual iterative algorithms for Bregman-penalized variational regularization, with explicit parameter choices and convergence analysis under variational source conditions.

Findings

Algorithms effectively reconstruct 3D atmospheric tomography images.

Methods demonstrate stability and convergence in both continuous and iterative senses.

Successful application to 2D and 3D tomographic problems confirms practical utility.

Abstract

Primal-dual splitting involving proximity operators in order to be able to find some approximation to the minimizer for a general form of Tikhonov type functional is in the focus of this work. This approximation is produced by a pair of iterative variational regularization procedures. Under the assumption of some variational source condition (VSC), total error estimation both in the iterative sense and in the continuous sense has been analysed separately. Rates of convergence will be obtained in terms some concave and positive definite index function. Of the choice of the penalty term, we are interested in Bregman distance penalization associated with the non-smooth total variation (TV) functional. Furthermore, following up the lower and bounds defined for the regularization parameter, some deterministic choice of the regularization parameter is given explicitly. It is in the emphasis of this work that the regularization parameter obeys Morozov`s discrepancy principle (MDP) in order for the stability analysis of regularized solution. In the computerized environment, the algorithms are verified as iterative regularization methods by applying it to an atmospheric tomography problem named as GPS-Tomography. Apart from this 3-D tomographic inverse problem, we also apply the algorithms to some 2-D conventional tomographic image reconstruction problems in order to be able test algorithms` capability of capturing the details and observe that algorithms behave as iterative regularization procedures.

Equations

$$\delta \xi +Tu=v^{\delta },$$

$$h(u)=1_{\Omega}(u):=\left\{\begin{array}[]{rcl}0&\mbox{, for}&u\in\Omega\subset{\mathcal{X}}\\ \infty&\mbox{, for}&u\notin\Omega\subset{\mathcal{X}}.\end{array}\right.$$

$$||T||:=\underset{1\le k\le M}{\max }\left\{\sqrt{{\sigma }_k}\right\}.$$

$$||\lambda u_1+(1-\lambda )u_2|{|}^2=\lambda ||u_1|{|}^2+(1-\lambda )||u_2|{|}^2-\lambda (1-\lambda )||u_1-u_2|{|}^2.$$

$$||T^{T}T(u_1-u_2)|{|}^2$$

$$-⟨u_1-u_2,T^{T}T(u_1-u_2)⟩\le -\frac{1}{||T|{|}^2}||T^{T}T(u_1-u_2)|{|}^2$$

$$\partial f(x^{\prime }):=\left\{\eta \in {\mathcal{X}}^{\ast }:f(x)-f(x^{\prime })\ge ⟨\eta ,x-x^{\prime }⟩\text{mbox}forallx\in \mathcal{X}\right\}.$$

$$D_{\mathcal{J}}:$$

$$D_{\mathcal{J}}(u,u^{\ast })\ne D_{\mathcal{J}}(u^{\ast },u),$$

$$D_{\mathcal{J}}(u,u^{\ast })\ge 0.$$

$$F_{\alpha }:$$

$$TV(u,\Omega )=\mathcal{J}(u):={\int }_{\Omega }|\nabla u(x)|dx\approxeq \sum _i|{\nabla }_{i}u|,$$

$$J(u)=g(Du)\text{mbox}where,g(\cdot )=||\cdot |{|}_1\text{mbox}withD(\cdot )=\nabla (\cdot ).$$

$$\partial J(u)=D^{\ast }\partial g(Du).$$

$$u_i=R_I(v^{\delta },w_{i-1},\Gamma ),$$

$$||T{u}_{i^{\ast }}-v^{\delta }||\le \tau \delta <||T{u}_i-v^{\delta }||$$

$$||T{u}_{i^{\ast }}^{\delta }-T{u}^{†}||\le (\tau +1)\delta ,$$

$$\tau \delta \le ||T{u}_i-v^{\delta }||$$

$${\mathrm{prox}}_{\mathcal{J}}(\tilde{u}):={\arg \min }_{u\in {\mathbb{R}}^N}\mathcal{J}(u)+\frac{1}{2}||u-\tilde{u}|{|}^2.$$

$$E:$$

$$\sigma ||u-u^{†}||\le \mathcal{J}(u)-\mathcal{J}(u^{†})+\Psi \left(||Tu-T{u}^{†}||\right)\text{mbox},forallu\in \mathcal{X}.$$

$$0\in \partial F_{{\alpha }_{\ast }}(u_{{\alpha }_{\ast }}^{\delta },v^{\delta }),$$

$$0=T^T(T{u}_{{\alpha }_{\ast }}^{\delta }-v^{\delta })+{\alpha }_{\ast }\partial \mathcal{J}(u_{{\alpha }_{\ast }}^{\delta })-{\alpha }_{\ast }\partial \mathcal{J}(u_0)+\hat{z}\text{mbox},where\hat{z}\in \partial h(u_{{\alpha }_{\ast }}^{\delta }).$$

$$0=T^T(T{u}_{{\alpha }_{\ast }}^{\delta }-v^{\delta })+{\alpha }_{\ast }{D}^T{\hat{w}}_{{\alpha }_{\ast }}^{\delta }-{\alpha }_{\ast }{D}^T{\hat{w}}_0+\hat{z},$$

$$\displaystyle\left\{\begin{array}[]{rcl}u_{\alpha}^{\delta}&=&\mathrm{prox}_{\mu h}\left[u_{\alpha}^{\delta}-\mu\left(T^{T}(Tu_{\alpha}^{\delta}-v^{\delta})+\alpha_{\ast}D^{T}(\hat{w}_{\alpha}^{\delta}-\hat{w}_{0})\right)\right]\\ \hat{w}_{\alpha}^{\delta}&=&\mathrm{prox}_{\nu g^{\ast}}\left(\hat{w}_{\alpha}^{\delta}+\nu Du_{\alpha}^{\delta}\right),\end{array}\right.$$

$$||u_i-u^{†}||\le ||u_i-u_{\alpha }^{\delta }||+||u_{\alpha }^{\delta }-u^{†}||$$

$$||x^{+}-y|{|}^2\le ||x^{-}-y|{|}^2-||x^{+}-x^{-}|{|}^2+2⟨x^{+}-y,\Delta ⟩+2g(y)-2g(x^{+}).$$

$$||u_{i+1}-u_{\alpha }^{\delta }|{|}^2\le ||u_i-u_{\alpha }^{\delta }|{|}^2-||u_{i+1}-u_i|{|}^2$$

$$||u_{\alpha }^{\delta }-u_{i+1}|{|}^2\le ||u_{\alpha }^{\delta }-u_{i+1}|{|}^2-||u_{\alpha }^{\delta }-u_{\alpha }^{\delta }||$$

$$||u_{i+1}-u_{\alpha }^{\delta }|{|}^2\le ||u_i-u_{\alpha }^{\delta }|{|}^2-||u_{i+1}-u_i|{|}^2$$

Full text

Under the assumption of some variational source condition (VSC), total error estimation both in the iterative sense and in the continuous sense has been analysed separately. Rates of convergence will be obtained in terms some concave and positive definite index function. Of the choice of the penalty term, we are interested in Bregman distance penalization associated with the non-smooth total variation (TV) functional. Furthermore, following up the lower and bounds defined for the regularization parameter, some deterministic choice of the regularization parameter is given explicitly. It is in the emphasis of this work that the regularization parameter obeys Morozov‘s discrepancy principle (MDP) in order for the stability analysis of regularized solution.

In the computerized environment, the algorithms are verified as iterative regularization methods by applying it to an atmospheric tomography problem named as GPS-Tomography. Apart from this 3-D tomographic inverse problem, we also apply the algorithms to some 2-D conventional tomographic image reconstruction problems in order to be able test algorithms‘ capability of capturing the details and observe that algorithms behave as iterative regularization procedures.

Keywords. iterative regularization, primal dual algorithm, Bregman iteration, total variation

\MakePerPage

footnote

1 Introduction

Following up a recent work in [3], we introduce the primal dual algorithms in the simpler forms without nested loops. Both algorithms are to be used for the purpose of finding iteratively regularized approximations for the inverse ill-posed problems. Stability analyses are reduced to be in the emphasis of the stability of the iteratively regularized approximations of the regularized minimizer for some general form of Tikhonov functional.

In general terms, regularization theory deals with approximation of some ill-posed inverse problem by a family of parametrized well-posed problems. Traditional quadratic-Tikhonov regularization [76, 77] has been well established and analyzed [34]. This work focuses on the analysis of some iterative regularization procedures.

Proximal mapping algorithms and Bregman iterated regularization procedures have been commonly known. Here, Bregman distance plays the role of penalization in the corresponding minimization problems. Having a look at the literature, firstly in the work [64], Bregman iteration has been proposed for providing solution to the basis pursuit problem. In a recent study by Sprung and Hohage et al., 2017, [74], authors have investigated convergence rates results for the such objective functionals with the Bregman distance as penalty term. Having optimization algorithms in the field of inverse problems as iterative regularization method has also become popular. Authors in [37] have proposed some primal-dual algorithm, wherein the convergence has been studied for the given noiseless measurement data. We consider linear, inverse ill-posed problems in the general form. Stability of the algorithms will be developed in the context of convex variational regularization and be verified in the Hadamard sense. Main results of our work are derived in the case of noisy measurement and are in the best interest of variational regularization theory.

2 Notations and Mathematical Setting

Over the finite dimensional Hilbert spaces ${\mathcal{X}}=\mathbb{R}^{N}$ and ${\mathcal{Y}}=\mathbb{R}^{M}$, let us be given some linear, injective, forward operator $T:{\mathcal{X}}\rightarrow{\mathcal{Y}}.$ In this work, we concentrate on the numerical solution for the linear inverse ill-posed problem of the form

[TABLE]

with some iterative regularization procedure. Here, the given noisy data $v^{\delta}\in{\mathcal{Y}}$ and the noise model is denoted by $\xi$ with the noise magnitude $\delta.$

Non-negativity constraint on our targeted data $u$ is imposed. Then, the constraint domain is treated as the indicator function $h:{\mathcal{X}}\longrightarrow\{0,1\}$ defined by

[TABLE]

Throughout the work, unless otherwise stated, the notation $||\cdot||$ without any subscript will be used as to denote the usual Euclidean norm. Let $\sigma(T^{T}T)$ be the spectrum of $T^{T}T$ is the set of those $\sigma_{k}\in\mathbb{R}^{N}.$ Then, for the finite dimensional forward operator $T:\mathbb{R}^{N}\rightarrow\mathbb{R}^{M}$ where $M<N,$ we define

[TABLE]

Below, we give two norm estimations that will be in use of our mathematical development. For some $u_{1},u_{2}\in{\mathcal{X}}$ and $\lambda\in\mathbb{R},$ the following equality holds, [75, Eq. (2.1)],

[TABLE]

Also, nonexpansiveness of the misfit term provides,

[TABLE]

which implies

[TABLE]

For some function $f:{\mathcal{X}}\rightarrow{\mathcal{Y}}$ and some point $x$ in the domain of $f,$ the subdifferential of $f$ at $x^{\prime},$ denotes $\partial f(x^{\prime})$ is defined by

[TABLE]

Definition 2.1**.**

[Generalized Bregman Distances] Let ${\mathcal{J}}:{\mathcal{X}}\rightarrow\mathbb{R}_{+}\cup\{\infty\}$ be a convex functional with the subgradient $q^{\ast}\in\partial{\mathcal{J}}(u^{\ast}).$ Then, for $u,u^{\ast}\in{\mathcal{X}},$ Bregman distance associated with the functional $J$ is defined by

[TABLE]

It is well known that the Bregman distance does not satisfy symmetry,

[TABLE]

and for the defined convex functional ${\mathcal{J}}$

[TABLE]

With all these tools stated, we consider the following objective functional,

[TABLE]

with some initial estimation $u_{0}\in{\mathcal{X}}.$ In particular, we associate the Bregman distance penalty term with the total variation (TV) functional defined by

[TABLE]

which is, in 3D case, $i=(i_{x},i_{y},i_{z}).$ For the sake of following the further calculations easily in the future developments of this work, we introduce TV in the composite form

[TABLE]

Thus,

[TABLE]

2.1 Overview on the iterative regularization and the choice of the regularization parameter

By an iterative procedure involving some iteration operator $R_{I},$ [14, Ch. 6], we aim to construct some approximation to the given inverse ill-posed problem (2.1)

[TABLE]

where $w^{i-1}$ is the collection of dual variables used during $i-1$ iteration steps, and $\Gamma$ is the auxiliary parameters such as step-size, relaxation parameter, regularization parameter.

In the iterative regularization procedures, discrepancy principles act as the stopping rules for the corresponding algorithms, [14, Section 6].

Definition 2.2**.**

[Morozov‘s Discrepancy Principle (MDP), [14, Def. 6.1]]

Given deterministic noise model $||v^{\dagger}-v^{\delta}||\leq\delta,$ if we choose $\tau>1$ and $i^{\ast}=i^{\ast}(\delta,v^{\delta})$ such that

[TABLE]

is satisfied for $u_{i^{\ast}}=R_{I}(v^{\delta},w_{i^{\ast}-1},\vec{\alpha})$ and $u_{i}=R_{I}(v^{\delta},w_{i-1},\vec{\alpha})$ for all $i<i^{\ast},$ then $u_{i^{\ast}}$ is said to satisfy Morozov‘s discrepancy principle.

Following up MDP, some immediate consequences can be given below,

[TABLE]

likewise,

[TABLE]

Our primal-dual splitting algorithms involve proximal mapping that is defined below.

Definition 2.3**.**

**[Proximal mapping]

Let ${\mathcal{J}}:\mathbb{R}^{N}\rightarrow\mathbb{R}\cup\{\infty\}$ be a proper, convex, lower-semicontinious function. Then $\mathrm{prox}_{{\mathcal{J}}}$ is defined as the unique minimizer**

[TABLE]

Measuring the deviation of the regularized solution $u_{\alpha}^{\delta}$ from the minimum norm solution $u^{\dagger}$ by the a priori and a posteriori strategies for the choice of the regularization parameter in Banach spaces with the VSC has been widely studied,

The objective of convergence and convergence rates results in the regularization theory is to be able find some stable bound for the total error estimation function defined by

[TABLE]

where the coefficient $\Lambda\in\mathbb{R}_{+}$ depends on the functional properties of the data on the pre-image space. Such estimation requires the knowledge of the smoothness of the minimum norm solution $u^{\dagger}$ which is some source condition in the form of varitional inequality, [15]****, [35], [38], [46, Eq. (1.4)], [45], [52, Section 4] [71, Theorem 2.60 - (g), Subsection 3.2.4]. Convergence and convergence rates results, or in other words the total error estimation, are derived in terms of a concave, monotonically increasing, positive definite index function that is a part of the VSC expression. Following up the work [3], the form of VSC that will be used in our analysis is given below. Also in the aforementioned works that have been dedicated to variational regularization, derivation of noise dependent error estimation following from VSC has been well explained. In our work, any data that is in the constraint domain is assumed to satisfy VSC.

Assumption 2.4**.**

[Variational Source Condition] Let $T:{\mathcal{X}}\rightarrow{\mathcal{Y}}$ be linear, injective forward operator and $v^{\dagger}\in\mathrm{range}(T).$ There exists some constant $\sigma\in(0,1]$ and a concave, monotonically increasing index function $\Psi$ with $\Psi(0)=0$ and $\Psi:[0,\infty)\rightarrow[0,\infty)$ such that for $q^{\dagger}\in\partial{\mathcal{J}}(u^{\dagger})$ the minimum norm solution $u^{\dagger}\in\mathrm{BV}(\Omega)$ satisfies

[TABLE]

3 Subdifferential Characterization

Both algorithms will evolve from the subdifferential characterization of the regularized minimizer of the objective functional (2.7). Then, by its first order optimality condition,

[TABLE]

which implies,

[TABLE]

Furthermore, recall the settings in (2.9) and (2.10) to represent (3.1) in the following form

[TABLE]

where $\hat{w}_{\alpha}^{\delta}\in\partial g(Du_{\alpha}^{\delta})$ and likewise $\hat{w}_{0}\in\partial g(Du_{0}).$

Theorem 3.1**.**

[3, Theorem 4.1]** [Subgradient characterization of the regularized solution]

For any positive valued $\alpha,\nu,\mu,$ the regularized minimizer $u_{\alpha}^{\delta}$ of the objective functional (2.7) is characterized by

[TABLE]

with $\hat{w}_{0}=\partial||Du_{0}||_{1}.$

4 The Primal-Dual Algorithms

We introduce a pair of algorithms involving proximal mappings. Both algorithms aim to provide approximation for the regularized minimizer of the objective functional (2.7). Algorithm 1 can be interpreted as a direct discrete form of the subdifferential characterization given above. However, Algorithm 2 is endowed with some projected convex extrapolation on a line segment due to the choice of the relaxation parameter $\lambda$.

Further than investigating whether the regularized iterations are better approximations, we are also interested in the convergences of those approximations towards the minimum norm solution.

5 Stability Analysis of the Algorithms

Iterative total error estimation can be decomposed in the following form

[TABLE]

The term on the far right has been very well analysed in the context of variational regularization. In what follows, we will focus on the term on the left hand side and the first term on the right hand side. Before stating that the $u_{i}$ is the approximation to the regularized minimizer $u_{\alpha}^{\delta}$ of the objective functional, necessary and sufficient conditions for the boundedness of $||u_{i}-u_{\alpha}^{\delta}||$ must be established. To this end, we shall give a pair of observations on the iterative approximations produced by the both algorithms. The assertions in the following formulations describe how the parameters in the algorithms must be chosen. Further assumption for the theoretical developments is the initial guess. In order to overcome the mathematical difficulties, it is always assumed that the initial guess of the objective functional and the initial guess for the algorithms are the same.

5.1 Iterative approximations of the regularized minimizer

The primary tool to study stability of the iteratively regularized approximations that are produced by the proximal gradient algorithms is formulated below.

Property 5.1**.**

[55*, Lemma 1]**

If $x^{+}=\mathrm{prox}_{g}(x^{-}+\Delta)$, then for any $y\in\mathbb{R}^{N},$*

[TABLE]

Before stability analysis, we give some observations on the primal variables produced by the both algorithms.

Proposition 5.2**.**

Let the step-length $\mu_{i}$ satisfy $\mu_{i}\leq\frac{2}{||T||^{2}}.$ Furthermore, let the iterative regularization parameter be chosen as $\alpha_{i}=i(\delta,v^{\delta})$ and $\nu_{i}\leq\frac{i(\delta,v^{\delta})^{2}}{\mu_{i}^{2}}.$ Then, iteratively regularized primal variable $u_{i+1}$ that is produced by Algorithm 1 is a better approximation of $u_{\alpha}^{\delta}$ than $u_{i}$ for each $i=0,1,2,\cdots$.

Proof.

Following first two error estimations between the final update $u_{i+1}$ of the Algorithm 1 and the regularized minimizer $u_{\alpha}^{\delta}$ of the objective functional (2.7) according to the Property 5.1 is given by

[TABLE]

Likewise, still by Property 5.1,

[TABLE]

After the necessary simplifications, these both estimations in total return,

[TABLE]

Some useful upper bound for the inner product on the right hand side of the 1st line is given below,

[TABLE]

where we have used (2.4). Also, total of the inner products on the 2nd and the 3rd lines of (5.4) can be given in a simpler form. Thus,

[TABLE]

Analagous estimations on the dual iterative variable $w_{i+1}$ and the dual variable $w_{\alpha}^{\delta}$ can be derived as follows,

[TABLE]

We rewrite the inner product on the right hand side

[TABLE]

Now, $\nu_{i}$ times (LABEL:primal_final_estI) plus $\mu_{i}\alpha_{i}$ times (5.7) with taking into account (5.8) will result in further term reduction,

[TABLE]

where we have also dropped the term $\mu_{i}\alpha_{i}||w_{i+1}-w_{\alpha}^{\delta}||^{2}$ from the left hand side since the boundedness of the error estimation for the primal variables are in the interest of this result. Also, the inner products, after using $\langle\cdot,D^{T}(\cdot)\rangle=\langle D(\cdot),\cdot\rangle$ are bounded in the following ways,

[TABLE]

In the light of these bounds, after multiplying both sides by $\frac{1}{\nu_{i}}$ of (5.9),

[TABLE]

Hence, asserted parameter choices will yield the result. ∎

As for the Algorithm 2, the similar observation is also formulated below.

Proposition 5.3**.**

Let the relaxation parameter of the step 8 of the Algorithm 2 be $\lambda\in(1,2)$ and the step-length be $\mu_{i}\leq\frac{2}{||T||^{2}}.$ Furthermore, let the iterative regularization parameter be chosen as $\alpha_{i}=i(\delta,v^{\delta})$ and $\nu_{i}\leq\frac{i(\delta,v^{\delta})^{2}}{\mu_{i}^{2}}.$ Then, the iteratively regularized primal variable $u_{i+1}$ is a better approximation of $u_{\alpha^{\ast}}^{\delta}$ than $u_{i}$ for each $i=0,1,2,\cdots$.

Proof.

We begin with applying the equality given in (2.3) to the step 8 of the algorithm,

[TABLE]

Now, by Property 5.1 the error estimation between the primal variable $\hat{u}_{i+1}$ and the regularized minimizer $u_{\alpha}^{\delta}$ is given,

[TABLE]

Still by Property 5.1, analagous error estimation is given below

[TABLE]

which is in other words, after multiplying both sides by $\frac{\mu_{i}}{\mu},$

[TABLE]

Recall that we consider the constant valued initial guess $u_{0}$ which makes the dual variable $w_{0}$ zero valued by its definition. Keeping this in mind and summing up (5.12) and (5.13) will return

[TABLE]

Here, the first inner product on the first line of the right hand side can be bounded similar to the proof above. So, one can quickly write down that

[TABLE]

Then, (LABEL:primal_var_total0) is reduced to

[TABLE]

This estimation together with (5.11), after putting the term $||\hat{u}_{i+1}-u_{\alpha}^{\delta}||^{2}$ on the other right hand side, in total

[TABLE]

We now apply Property 5.1 for the dual variables as in the proof of Proposition 5.2,

[TABLE]

We arrive at the following result after summing $\nu_{i}$ times (LABEL:primal_var_total2) with $\mu_{i}\alpha_{i}$ times (5.16),

[TABLE]

As in the previous proof, the inner products are bounded by

[TABLE]

We now multiply both sides by $\frac{1}{\nu_{i}}$ of (5.17) with taking into account the bounds above. So that we obtain,

[TABLE]

Hence, the result is a natural consequence of the choices of the asserted parameters. ∎

5.2 Convergence of the Iteratively Regularized Approximation Against the Minimum Norm Solution

This section is dedicated to analyse the convergence of the iterative regularized approximation towards the minimum norm solution.

For each algorithm, we will establish convergence results in the Hadamard sense on each iterative step $i=1,2,\cdots.$ Following up these results, cumulative error estimations for both algorithms will be formulated. In these estimations, exact choices of the step-length $\mu$ and the relaxation parameter $\lambda_{i}$ will be conveyed.

Theorem 5.4**.**

Let the iterative regularization parameter be $\alpha_{i}(\delta,v^{\delta})=\frac{1}{i(\delta,v^{\delta})}$ and $\mu_{i}\leq\frac{2}{||T||^{2}}.$ Then the convergence of $u_{i+1}$ produced by Algorithm 1 to the minimum norm solution $u^{\dagger}$ of the linear inverse problem $Tu^{\dagger}=v^{\dagger}$ which satisfies the VSC (2.16), with the given deterministic noise model $||v^{\dagger}-v^{\delta}||\leq\delta$ is satisfied as $\delta\rightarrow 0$ and $i\rightarrow\infty.$

Proof.

According to Property 5.1,

[TABLE]

By means of the deterministic noise model $||v^{\dagger}-v^{\delta}||\leq\delta,$ consequence of the MDP (2.14) and the condition on the step-length $\mu_{i}\leq\frac{2}{||T||^{2}},$ the first inner product on the right hand side is bounded as follows,

[TABLE]

Note that we have also used Assumption 2.4 for the term $u_{i}$ since by definition $u_{i}\in BV(\Omega).$ Thus, if we plug this estiamation in (5.19), then we obtain

[TABLE]

Once more, the inner product on the right hand side can suitably be bounded by,

[TABLE]

Hence, again by the condition on the step-length $\mu_{i}\leq\frac{2}{||T||^{2}}$ with $\alpha_{i}(\delta,v^{\delta})=\frac{1}{i(\delta,v^{\delta})},$ the following form of (5.20) yield the result,

[TABLE]

∎

Although in Theorem 5.4, the step-length has been given in dynamical sense, below we see that in order to obtain desired convergence we must guarantee that some terms are Cesáro summable only when $\mu_{i}$ is rather chosen fixed.

Theorem 5.5**.**

Let the initial guess $u_{0}$ be of $\mathrm{BV}(\Omega)$ and the dynamical regulariazation parameter that satisfies MDP be $\alpha_{i}(\delta,v^{\delta})=\frac{1}{i(\delta,v^{\delta})}.$ If the step-length is chosen as $\mu=\frac{1}{2i^{\ast}||T||^{2}},$ then after $i^{\ast}$ times iteration of Algorithm 1 and in the light of Assumption 2.4 convergence of the iteratively regularized approximation towards the minimum norm solution is satisfied, i.e., $||u_{i^{\ast}}-u^{\dagger}||\rightarrow 0$ as $i^{\ast}\rightarrow\infty$ whilst $\delta\rightarrow 0.$

Proof.

If we iterate the estimation (5.19) from $i=0$ to $i=i^{\ast}-1$ and sum-up over $i,$ after simplifications all the necessary terms, we then obtain the following cumulative error estimation,

[TABLE]

Let us begin with rewriting the first inner product on the right hand side,

[TABLE]

With the inclusion of MDP, VSC and the determenistic noise error, each piece of the right hand side will be analysed separately after applying the Cauchy-Schwartz inequality.

[TABLE]

[TABLE]

[TABLE]

We will discuss the convergence of the sequences on the first two lines as $i^{\ast}\rightarrow\infty$ in the sense of Cesáro summation. Boundedness of the term $||u_{i+1}-u^{\dagger}||$ has been discussed already above. Of the convergence conditions formulated in Theorem 5.4 is the choice of the step-length $\mu.$ In the light of this condition, if we fix the step-length $\mu=\frac{1}{2i^{\ast}||T||^{2}},$ then

[TABLE]

Thus the sequence on the right hand side is Cesáro summable as $i^{\ast}\rightarrow\infty.$ Following up this estimation, bounds for rest of the terms read

[TABLE]

and likewise,

[TABLE]

Also, the last inner product can be bounded again by using Cauch-Schwartz inequality,

[TABLE]

Now, if the negative term drops and the Assumption 2.4 is taken into account,

[TABLE]

∎

Theorem 5.6**.**

Let the iterative regularization parameter $\alpha_{i}(\delta,v^{\delta})=\frac{1}{i(\delta,v^{\delta})}$ be chosen according to (MDP), the relaxation parameter $\lambda\in(1,2)$ of the step 8 of the Algorithm 2 and $\mu_{i}\leq\frac{2}{||T||^{2}}.$ Then the convergence of $u_{i+1}$ produced by Algorithm 2 to the minimum norm solution $u^{\dagger}$ of the linear inverse problem $Tu^{\dagger}=v^{\dagger}$ which satisfies the VSC (2.16), with the given deterministic noise model $||v^{\dagger}-v^{\delta}||\leq\delta$ is satisfied as $\delta\rightarrow 0$ and $i\rightarrow\infty.$

Proof.

According to equality (2.3), the error estimation between the final update $u_{i+1}$ and the minimum norm solution $u^{\dagger}$ is

[TABLE]

Also by Property 5.1, some error estimation between the primal variable and the minimum norm solution

[TABLE]

Similar to the proof above, the inner products on the right hand side can easily be handled by replacing $u_{i+1}$ with $\hat{u}_{i+1}.$ On account of $\mu_{i}\leq\frac{2}{||T||^{2}}$ and by definition $\hat{u}_{i+1}$ and $u_{i}$ are from the constraint domain $\mathrm{BV}(\Omega)$, we then obtain,

[TABLE]

Also, likewise, as for the second inner product on the right hand side of (5.25),

[TABLE]

Plugging these both estimations into (5.25), summing up with (5.24) and since the variables $u_{i}$ and $\hat{u}_{i+1}$ are of $\mathrm{BV(\Omega)},$ the following estimation yields the result,

[TABLE]

∎

Now, we are also concerned about the exact choice of the relaxation parameter $\lambda\in(1,2)$ in Algorithm 2. As in the proof of Theorem 5.5, Cesáro summation of some terms are guaranteed depending on the choice of $\lambda$ that is formulated below.

Theorem 5.7**.**

Let the initial guess $u_{0}$ be of $\mathrm{BV}(\Omega)$ and the dynamical regulariazation parameter that satisfies MDP be $\alpha_{i}(\delta,v^{\delta})=\frac{1}{i(\delta,v^{\delta})}.$ If the relaxation parameter $\lambda=\frac{2}{i^{\ast}}$ and step-length $\mu=\frac{1}{2i^{\ast}||T||^{2}},$ then $||u_{i^{\ast}}-u^{\dagger}||\rightarrow 0$ as $i^{\ast}\rightarrow\infty$ whilst $\delta\rightarrow 0.$

Proof.

Let us iterate the equality (5.24) from $i=0$ to $i=i^{\ast}-1$ and sum up over $i,$

[TABLE]

Note that we have pulled the term2 $||u_{0}-u^{\dagger}||^{2}$ and $\lambda||u_{i^{\ast}}-u^{\dagger}||^{2}$ out of the sequences $(1-\lambda)\sum_{i=1}^{i=i^{\ast}-1}||u_{i}-u^{\dagger}||^{2}$ and $\lambda\sum_{i=0}^{i=i^{\ast}-2}||\hat{u}_{i+1}-u^{\dagger}||^{2}$ respectively, so that we can make use of these individuals in the next step. Then, we repeat the same for the estimation (5.25),

[TABLE]

In this last estimation, analagous sequences of inner products on the right hand side have already been bounded in the proof of Theorem 5.5. In order to avoid dublicating the calculations, one can replace $u_{i+1}$ by $\hat{u}_{i+1}.$ Then quick adaptation of those estimations above must reveal

[TABLE]

Also likewise,

[TABLE]

If we sum (LABEL:AlgII_cumulative1) and (5.27) with the consideration of both (5.28) and (5.29), after the necessary algebraic arrangment, we then arrive at,

[TABLE]

With the given choice of the of the relaxation parameter $\lambda$ all the terms on the first and the second lines are either bounded or negative. Furthermore, the initial guess $u_{0}$ and $\hat{u}_{i^{\ast}}$ are from $\mathrm{BV}(\Omega).$ Hence, all these aforementioned assumptions yield the assertion. ∎

6 Numerical Results; Behaviours of the Algorithms From Regularization Aspect

In this section, we will put the algorithms into tests in order to observe they behave as iterative regularization prcedures. To this end, the followings in our numerical tests will be analysed in the computerized environment;

Iterative error values on the image space, 2. 2.

Iterative error values on the pre-image space, 3. 3.

Iterative error values of both algorithms against each other, 4. 4.

Iterative error values with the different noise amount input on the measured data.

In what follows, the data on the pre-image space are painted drawings, see Acknowledgement. The measured data is the sinogram tomographic measurement applied on the image, [9].

Firstly, we will provide some simple benchmark on the efficiency of each algorithm against each other. The mathematical developlment has proved that Algorithm 2 produce

As can be seen in the Figures 2 and 3, it is verified that the less noise amount in the measured data provide less error estimation both on the image and the pre-image spaces.

7 An Atmospheric Tomography Problem: GPS Tomography

One important predictor in meteorology is the humidity of the atmosphere. This is estimated by fan-beam measurements between satellite transmitters and land-based receivers. The measurements are sparse and fluctuate randomly with receiver availability. The task is to reconstruct from these measurements the 3-dimensional, spatially varying index of refraction of the atmosphere, from which the relative humidity can be inferred.

GPS-tomography involves the reconstruction of some quantity (e.g. humidity), pointwise within a volume, from geodesic X-ray measurements transmitted by nonuniformly distributed transducers (satellites). These measurements are collected by nonuniformly distributed receivers on the ground (ground stations). As in the conventional tomography, the task here is the reconstruction of the density volume profile of a layer in the atmosphere from a set of line integrals that are of fan-beam projections.

Function reconstruction from its measured line integrals was firstly proposed and solved in [66]. Profound mathematical and numerical aspects of the computerized tomography have been studied in [59, 61]. Measurement from the Radon transform is obtained by integrating some integrable function over the hyperplanes in $\mathbb{R}^{N}.$ The ray transform, on the other hand, produces measurement by integrating the function over straight lines. It is known that in the two dimensional tomography, general Radon and ray transformations coincide, [61, p. 17].

In the discretized form of the problem, it is assumed that each station receives equal number of signals transmitted by the satellites. Also for the sake of simplicity, we ignore any deviations from the shortest path between transmitters and receivers due to atmospheric refractivity. The received signal is then modelled as a line integral along the shortest path between the satellites and the ground stations.

Peculiar to this problem, reconstructions by Kalman filtering and ART have been widely applied, [13, 57, 65, 84]. Different from these conventional numerical reconstruction methods, a quasi-Newton approach, which is limited memory BFGS (L-BFGS), has also been proposed to obtain the optimal regularized solution, [2]. The L-BFGS algorithm has been also applied for atmospheric imaging wherby the forward problem has been modelled as a phase retrieval problem, see [78].

7.1 Physical Problem

This is an inverse problem with incomplete data. It is well known that the incompleteness of data causes nonuniqueness issue in inverse problems, [61, p. 144]. Thorough implementation and inversion of geodesic X-ray transform has been studied in [58]. The model of this problem is also widely known. We refer readers to [2] how the compact support assumption on the targeted data has been established for satisfying unique solvibility of the problem.

7.2 Numerical Results; Response of the algorithms to GPS-Tomography

Above section has already been dedicated to understand that the algorithms are iterative regularization procedures, furthermore it has been demonstrated that Algorithm 2 provide better results. Thus in this section, we will rather present the numerical results of Algorithm 2 when it is applied to the GPS-Tomography problem. Specifically in this problem, recent novel reserch works [81, 83] indicate the importance of the number and the sources of the measured data. In this work, we are only able to emphasize that how number of the measurement data, e.g. the signals, affects the quality of the reconstructed data. In Figure 4, the numerical results have been produced when all the parameters are chosen as fixed coefficient. Figure 4 display the results when the parameters are dynamically chosen. Having a comparison two results against each other, dynamical choice of the parameters allow us to obtain better results with much less iteration steps during the procedure. Regarding the impact of the number of the measured data on the reconstruction, Figures 6 and 7 display less quality on the reconstructed data when the problem is defined as underdetermined.

In the Figure 8, we demonstrate the importance of the parameter choices. As in the numerical results, with the dynamically chosen parameters, one can observe less error estimation at the end of the procedures.

8 Discussion and Future Prospects

Having the initial guess as constant in the mathematical analysis may not reveal the purpose of having Bregman distance as penalizer in our objective functional (2.7). Therefore, it could be worthwhile to consider the following iteration scheme

[TABLE]

Then, as a result of some quick calculations, the subdifferential characterization of the iterative minimizer is given by

[TABLE]

An algorithm solving this system will definitely convey minimization impact when the iterative form of Bregman distance is introduced as penalty term.

Acknowledgement

The author is indepted to Ignace Loris for the valuable discussions throughout the work. Also, the artist Bengü Uğurlu Öğretmen is deeply recognized and appreciated for sharing the contemporary artwork.