A new class of accelerated regularization methods, with application to bioluminescence tomography

TL;DR

This paper introduces a new class of accelerated iterative regularization methods based on second order evolution equations, demonstrating their effectiveness in solving ill-posed problems and applying them to bioluminescence tomography.

Contribution

The paper develops a novel accelerated regularization framework using second order evolution equations and applies it to bioluminescence tomography, showing improved convergence and accuracy.

Findings

Proven regularization properties and convergence rates.

Numerical examples demonstrate enhanced accuracy and acceleration.

Comparison shows superiority over existing methods.

Abstract

In this paper we propose a new class of iterative regularization methods for solving ill-posed linear operator equations. The prototype of these iterative regularization methods is in the form of second order evolution equation with a linear vanishing damping term, which can be viewed not only as an extension of the asymptotical regularization, but also as a continuous analog of the Nesterov's acceleration scheme. New iterative regularization methods are derived from this continuous model in combination with damped symplectic numerical schemes. The regularization property as well as convergence rates and acceleration effects under the H\"older-type source conditions of both continuous and discretized methods are proven. The second part of this paper is concerned with the application of the newly developed accelerated iterative regularization methods to the diffusion-based bioluminescence tomography, which is modeled as an inverse source problem in elliptic partial differential equations with both Dirichlet and Neumann boundary data. A relaxed mathematical formulation is proposed so that the discrepancy principle can be applied to the iterative scheme without the usage of Sobolev embedding constants. Several numerical examples, as well as a comparison with the state-of-the-art methods, are given to show the accuracy and the acceleration effect of the new methods.

Tables

Table 1. Table 1: The iterative number k ∗ superscript 𝑘 k^{*} and the corresponding relative error L ​ 2 ​ E ​ r ​ r k ∗ 𝐿 2 𝐸 𝑟 subscript 𝑟 superscript 𝑘 L2Err_{k^{*}} vs τ 𝜏 \tau .

$\tau$	Example 1		Example 2
	L2Err${}_{k^{\ast }}$	$k^{\ast }$	L2Err${}_{k^{\ast }}$	$k^{\ast }$
$2^{-1}$	1.8066e-3	$N_{\max }$	1.6067e-2	$N_{\max }$
$1$	4.8745e-3	246	2.2110e-2	4384
$2$	5.1612e-3	245	2.4860e-2	3318
$2^2$	5.5573e-3	244	3.1812e-2	2226
$2^3$	6.6180e-3	242	3.9872e-2	1255
$2^4$	8.6874e-3	239	4.8587e-2	966
$2^5$	1.4807e-2	232	6.8973e-2	560
$2^6$	2.6672e-2	221	8.7679e-2	362
$2^7$	5.2232e-2	202	1.1497e-1	129

Table 2. Table 2: The iterative number k ∗ superscript 𝑘 k^{*} and the corresponding relative error L2Err k ∗ superscript 𝑘 {}_{k^{*}} vs Δ ​ t Δ 𝑡 \Delta t .

$\Delta t$	Example 1		Example 2
	L2Err${}_{k^{\ast }}$	$k^{\ast }$	L2Err${}_{k^{\ast }}$	$k^{\ast }$
$2^{-10}$	4.6816e-3	14964	8.4071e-2	$N_{\max }$
$2^{-9}$	4.6835e-3	7483	5.7224e-2	$N_{\max }$
$2^{-8}$	4.6852e-3	3743	4.1839e-2	37270
$2^{-7}$	4.6891e-3	1873	4.1849e-2	18629
$2^{-6}$	4.7268e-3	937	4.1866e-2	9309
$2^{-5}$	4.7560e-3	470	4.1894e-2	4650
$2^{-4}$	4.9548e-3	236	4.1945e-2	2321
$2^{-3}$	Divergence	-	4.1988e-2	1159
$2^{-2}$	Divergence	-	Divergence	-

Table 3. Table 3: The iterative number k ∗ superscript 𝑘 k^{*} and the corresponding relative error L2Err k ∗ superscript 𝑘 {}_{k^{*}} vs s 𝑠 s .

$s$	Example 1		Example 2
	L2Err${}_{k^{\ast }}$	$k^{\ast }$	L2Err${}_{k^{\ast }}$	$k^{\ast }$
-0.499	6.1991e-3	2506	2.5643e-2	1414
-0.4	9.5068e-3	651	2.3717e-2	4276
-0.3	1.0187e-2	229	3.1115e-2	3490
-0.2	6.3593e-3	666	2.6702e-2	1506
-0.1	7.6404e-3	244	3.8602e-2	703
0	5.2214e-3	681	3.8309e-2	708
$2^{-5}$	5.4304e-3	540	3.8893e-2	709
$2^{-4}$	6.4303e-3	256	3.9249e-2	710
$2^{-3}$	5.1272e-3	547	4.0054e-2	713
$2^{-2}$	5.0820e-3	413	3.9117e-2	798
$2^{-1}$	5.2929e-3	287	4.1006e-2	808
$1$	4.8956e-3	321	4.2226e-2	912
$2$	4.9548e-3	236	4.1988e-2	1159
$2^2$	4.7348e-3	349	4.1857e-2	1548
$2^3$	4.6814e-3	547	4.1702e-2	2137
$2^4$	4.6760e-3	884	4.1598e-2	2991
$2^5$	4.6744e-3	1316	4.1532e-2	4210

Table 4. Table 4: Comparison of different methods. The CPU time is measured in seconds.

${\delta }^{\prime }$	$0.5\%$		$1\%$		$5\%$
Example 1
Methods	L2Err${}_{k^{\ast }}$	$k^{\ast }$ (CPU)	L2Err${}_{k^{\ast }}$	$k^{\ast }$ (CPU)	L2Err${}_{k^{\ast }}$	$k^{\ast }$ (CPU)
Landweber	4.7418e-3	5079(149.61)	5.5450e-3	3735(109.83)	5.7679e-3	2787(88.44)
$\nu =0.25$	2.0463e-1	18594(541.84)	4.2972e-2	2703(80.31)	4.3583e-2	615(19.78)
$\nu =0.5$	1.1026e-2	1706(52.78)	8.0824e-3	314(9.28)	7.6839e-3	186(5.97)
$\nu =1$	4.6974e-3	218(6.78)	4.7379e-3	89(2.81)	5.0483e-3	86(2.72)
$\nu =2$	4.7277e-3	138(4.17)	5.2199e-3	135(4.86)	5.2410e-3	128(4.02)
Nesterov	4.7724e-3	787(23.08)	5.0624e-3	283(9.02)	1.0487e-3	152(4.92)
NSS	4.9850e-3	236(7.53)	6.1620e-3	232(6.97)	6.5021e-3	223(6.75)
ARM	4.9968e-3	236(7.78)	6.1991e-3	232(6.91)	6.5776e-3	223(6.78)
Example 2
Methods	L2Err${}_{k^{\ast }}$	$k^{\ast }$ (CPU)	L2Err${}_{k^{\ast }}$	$k^{\ast }$ (CPU)	L2Err${}_{k^{\ast }}$	$k^{\ast }$ (CPU)
Landweber	7.7560e-2	7534(236.28)	7.8878e-2	6857(217.75)	1.1761e-1	723(22.84)
$\nu =0.25$	3.7686e-2	955(31.28)	3.7667e-2	955(30.67)	7.5047e-2	132(4.61)
$\nu =0.5$	6.3875e-2	214(7.31)	6.7694e-2	179(5.95)	1.1228e-1	67(2.39)
$\nu =1$	7.5760e-2	190(6.25)	7.5604e-2	188(6.23)	1.1739e-1	47(1.73)
$\nu =2$	7.8135e-2	246(8.11)	7.9302e-2	236(7.95)	1.1921e-1	67(2.36)
Nesterov	7.2335e-2	365(11.55)	7.1938e-2	363(11.94)	1.1715e-1	8(2.75)
NSS	7.8257e-1	450(16.94)	7.9035e-2	436(13.97)	1.2095e-1	121(4.23)
ARM	7.8275e-2	450(14.22)	7.8956e-2	437(14.52)	1.212e-1	121(4.06)

Equations

$$Kf=y,$$

$$f_{k+1}^{\delta }=f_k^{\delta }+\Delta t{K}^{\ast }(y^{\delta }-K{f}_k^{\delta }),\Delta t\in (0,2/\Vert K^{\ast }K\Vert )(k=0,1,\mathrm{2...})$$

$${\dot{f}}^{\delta }(t)+K^{\ast }K{f}^{\delta }(t)=K^{\ast }{y}^{\delta },f(0)=f_0,$$

$${\ddot{f}}^{\delta }(t)+\frac{1+2s}{t}{\dot{f}}^{\delta }(t)+K^{\ast }K{f}^{\delta }(t)=K^{\ast }{y}^{\delta },f(0)=f_0,\dot{f}(0)=0,$$

$$\displaystyle\left\{\begin{array}[]{ll}-{\rm div}(D\nabla u)+\mu_{a}u=f\chi_{\Omega_{0}}\quad{\rm in\ }\Omega,\\ D\partial_{\nu}u=g_{2}\quad{\rm on}\ \Gamma\end{array}\right.$$

$$u=g_1\mathrm{on}\Gamma .$$

$$\displaystyle\left\{\begin{array}[]{ll}\langle\ddot{f}^{\delta}(t),u_{j}\rangle+\frac{1+2s}{t}\langle\dot{f}^{\delta}(t),u_{j}\rangle+\lambda^{2}_{j}\langle f^{\delta}(t),u_{j}\rangle=\lambda_{j}\langle y^{\delta},v_{j}\rangle,~{}j=1,2,\cdots\\ f^{\delta}(0)=f_{0},~{}\dot{f}^{\delta}(0)=0,\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}\ddot{\xi}_{j}(t)+\frac{1+2s}{t}\dot{\xi}_{j}(t)+\lambda^{2}_{j}\xi_{j}(t)=\lambda_{j}\langle y^{\delta},v_{j}\rangle,\\ \xi_{j}(0)=\langle f_{0},u_{j}\rangle,~{}\dot{\xi}_{j}(0)=0,\end{array}\right.\quad j=1,2,\cdots.$$

$${\xi }_j(t)=2^s\Gamma (s+1)\frac{J_s({\lambda }_{j}t)}{({\lambda }_{j}t{)}^s}⟨f_0,u_j⟩+\left(1-2^s\Gamma (s+1)\frac{J_s({\lambda }_{j}t)}{({\lambda }_{j}t{)}^s}\right){\lambda }_j^{-1}⟨y^{\delta },v_j⟩,$$

$$\begin{array}[]{ll}f^{\delta}(t)=\sum\limits_{j}2^{s}\Gamma(s+1)\frac{J_{s}(\lambda_{j}t)}{(\lambda_{j}t)^{s}}\langle f_{0},u_{j}\rangle u_{j}+\sum\limits_{j}\left(1-2^{s}\Gamma(s+1)\frac{J_{s}(\lambda_{j}t)}{(\lambda_{j}t)^{s}}\right)\lambda^{-1}_{j}\langle y^{\delta},v_{j}\rangle u_{j},\\ =:(1-K^{*}Kg(t,K^{*}K))f_{0}+g(t,K^{*}K)K^{*}y^{\delta},\end{array}$$

$$g(t,\lambda )=\frac{1-2^s\Gamma (s+1)\frac{J_s(\sqrt{\lambda }t)}{(\sqrt{\lambda }t{)}^s}}{\lambda }.$$

$$\underset{\delta \to 0}{\lim }{T}_{\ast }=\infty \text{ and }\underset{\delta \to 0}{\lim }\delta \cdot T_{\ast }=0,$$

$$r(t,\lambda )=1-\lambda g(t,\lambda )=2^s\Gamma (s+1)\frac{J_s(\sqrt{\lambda }t)}{(\sqrt{\lambda }t{)}^s}.$$

$$\displaystyle\left\{\begin{array}[]{ll}\ddot{r}(t,\lambda)+\frac{1+2s}{t}\dot{r}(t,\lambda)+\lambda r(t,\lambda)=0,\\ r(0,\lambda)=1,\quad\dot{r}(0,\lambda)=0.\end{array}\right.$$

$$\Vert f^{\delta }(t)-f^{†}\Vert \le \Vert f^{\delta }(t)-f(t)\Vert +\Vert f(t)-f^{†}\Vert$$

$$\le \delta \underset{\lambda >0}{\sup }\sqrt{\lambda }g(t,\lambda )+\Vert r(t,K^{\ast }K)(f_0-f^{†})\Vert$$

$$\dot{E}(t)=2\dot{r}(t,\lambda )[\ddot{r}(t,\lambda )+\lambda r(t,\lambda )]=-\frac{2(1+2s)}{t}{\dot{r}}^2(t,\lambda )\le 0,$$

$$|r(t,\lambda )|\le 1\text{ for all }\lambda >0,t\ge 0.$$

$$J_s(\sqrt{\lambda }t)=\sqrt{\frac{2}{\pi }}{\left(\sqrt{\lambda }t\right)}^{-\frac{1}{2}}\cos \left(\sqrt{\lambda }t-\frac{\pi (2s+1)}{4}\right)+\mathcal{O}\left(\frac{1}{\sqrt{\lambda }t}\right)\text{ as }t\to \infty ,$$

$$\Vert r(t,K^{\ast }K)(f_0-f^{†})\Vert \to 0\text{ as }t\to \infty .$$

$${\zeta }^{\prime }(\tau )=-2\tau \sum _{k=1}^{\infty }\left\{{ȷ}_{s,k}^{-2}\prod _{i\ne k}^{\infty }\left(1-\frac{{\tau }^2}{{ȷ}_{s,i}^2}\right)\right\}<0,$$

$$\zeta (\tau )\ge 1-\tau \ge 1-\frac{\tau }{2{\tau }_1},\text{ for }\tau \in (0,{\tau }_0].$$

$$\zeta (\tau )=r(t,\lambda )\ge -1\ge 1-\frac{2}{{\tau }_0}\tau \ge 1-\frac{\tau }{2{\tau }_1}\text{ for }\tau \ge {\tau }_0.$$

$$r(t,\lambda )\ge 1-\frac{\sqrt{\lambda }t}{2{\tau }_1}\text{ for all }\lambda >0,t\ge 0.$$

$$\sqrt{\lambda }g(t,\lambda )=\frac{1-r(t,\lambda )}{\sqrt{\lambda }}\le \frac{t}{2{\tau }_1},$$

$$f_0-f^{†}={\left(K^{\ast }K\right)}^{\mu }{v}_0\text{ with }\Vert v_0\Vert \le \rho .$$

$$\Vert f^{\delta }(T_{\ast })-f^{†}\Vert =\mathcal{O}\left({\delta }^{\frac{2\mu }{2\mu +1}}\right)\text{ as }\delta \to 0.$$

$$\Vert f^{\delta }(T_{\ast })-f^{†}\Vert =\mathcal{O}\left({\delta }^{\frac{2s+1}{2s+3}}\right)\text{ as }\delta \to 0.$$

$$\displaystyle\begin{array}[]{ll}\|r(t,K^{*}K)(f_{0}-f^{\dagger})\|=\|r(t,K^{*}K)\left(K^{*}K\right)^{\mu}v_{0}\|\leq\rho\sup_{\lambda}r(t,\lambda)\lambda^{\mu}\\ \quad\leq C_{0}2^{s}\Gamma(s+1)\rho\sup_{\lambda\in(0,\|K\|^{2}]}\min\left\{\lambda^{-\frac{1+2s}{4}}t^{-\frac{1+2s}{2}},1\right\}\lambda^{\mu}\\ \quad\leq C_{1}\cdot\left\{\begin{array}[]{l}t^{-2\mu},\textrm{~{}if~{}}\mu\in(0,\frac{1+2s}{4}]\\ t^{-\frac{1+2s}{2}},\textrm{~{}if~{}}\mu>\frac{1+2s}{4},\end{array}\right.\end{array}$$

$$\displaystyle\begin{array}[]{ll}\|f^{\delta}(T_{*})-f^{\dagger}\|\leq\left\{\begin{array}[]{l}\delta\frac{T_{*}}{2\tau_{1}}+C_{1}(T_{*})^{-2\mu}=(\frac{1}{2\tau_{1}}+C_{1})\delta^{\frac{2\mu}{2\mu+1}},\textrm{~{}if~{}}\mu\in(0,\frac{1+2s}{4}]\\ \delta\frac{T_{*}}{2\tau_{1}}+C_{1}(T_{*})^{-\frac{1+2s}{2}}=(\frac{1}{2\tau_{1}}+C_{1})\delta^{\frac{2s+1}{2s+3}},\textrm{~{}if~{}}\mu>\frac{1+2s}{4},\end{array}\right.\end{array}$$

$$(\lambda t^2{)}^{\mu }|r(t,\lambda )|\ge \gamma$$

Full text

The second part of this paper is concerned with the application of the newly developed accelerated iterative regularization methods to the diffusion-based bioluminescence tomography, which is modeled as an inverse source problem in elliptic partial differential equations with both Dirichlet and Neumann boundary data. A relaxed mathematical formulation is proposed so that the discrepancy principle can be applied to the iterative scheme without the usage of Sobolev embedding constants. Several numerical examples, as well as a comparison with the state-of-the-art methods, are given to show the accuracy and the acceleration effect of the new methods.

1 Introduction

In the first part of this paper we consider linear operator equations

[TABLE]

where $K$ is a compact linear operator acting between two infinite dimensional Hilbert spaces $Q$ and $Y$ such that the range $\mathcal{R}(K)$ of $K$ is an infinite dimensional subspace of $Y$. Then the range $\mathcal{R}(K)$ is a non-closed subset of $Y$. For simplicity, we denote by $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ in the sequel the inner products and norms for both Hilbert spaces $Q$ and $Y$. The non-closedness of the forward operator $K$ is typical for operator equations (1) which are models for linear inverse problems. More precisely, due to the compactness of $K$, the operator equation (1) is ill-posed of type II in the sense of Nashed (cf. [18]). As a consequence of this ill-posedness, a regularization method must be employed in order to obtain reasonable and stable approximate solutions to (1) if the measurement data contains noise. In this context, we consider iterative regularization methods and assume to know noisy data $y^{\delta}\in Y$ instead of the exact right-hand side $y\in\mathcal{R}(K)$, obeying the deterministic noise model $\|y^{\delta}-y\|\leq\delta$ with a priori known noise level $\delta>0$. The focus of our paper will be on studying a specific application in bioluminescence tomography, and we refer to Section 4 for details.

The dominant iterative regularization method for solving (1) should be the Landweber method, given by

[TABLE]

with some starting element $f_{0}\in Q$, where $K^{*}$ denotes the adjoint operator of $K$. The continuous analog to (2) as $\Delta t$ tends to zero is known as asymptotic regularization or Showalter’s method (see, e.g., [28, 29]). It is in the form of a first order evolution equation

[TABLE]

where an artificial scalar time $t$ is introduced. There must be chosen an appropriate finite stopping time $T_{*}=T_{*}(\delta)$ (a priori choice) or $T_{*}=T_{*}(\delta,y^{\delta})$ (a posteriori choice) in order to ensure the regularizing property $f^{\delta}(T_{*})\to f^{\dagger}$ as $\delta\to 0$. Here and later on, $f^{\dagger}$ represents the unique minimum-norm solution of (1). Moreover, it has been shown that by using Runge-Kutta integrators, all of the properties of asymptotic regularization (3) carry over to its numerical realization [25]. Hence, the continuous model (3) is of particular importance for studying the intrinsic properties of a broad class of general regularization methods for inverse problems, and can be used for the development of new iterative regularization algorithms by combining some appropriate numerical schemes. Inspired by this, the authors in [32] studied the second order asymptotical regularization with the fixed damping parameter.

However, a fatal defect for large-scale problems is the slow performance of the Landweber iteration (too many iterations required for optimal stopping) as well as of the (conventional and second order with a fixed damping parameter) asymptotical regularization methods, i.e. overly excessive stopping times $T_{*}$ are required for obtaining optimal convergence rates. Therefore, in practice, accelerating strategies are usually used. In so doing, the most commonly known methods are the $\nu$-method [7, § 6.3] and the Nesterov acceleration scheme [21]. Recently, the authors in [33] introduced the fractional order asymptotical regularization, and proved that the fractional order plays the role of acceleration. In this paper, we are interested in the following second order evolution equation with a linear vanishing damping term

[TABLE]

where $s>-1/2$ is a fixed number. One motivation to study (4) is that it can be viewed as an infinite dimensional extension of the Nesterov’s scheme in the sense that for all fixed $T>0$ ([27]): $\lim\limits_{\omega\to 0}\max\limits_{0\leq k\leq T/\sqrt{\omega}}\|f^{\delta}_{k}-f^{\delta}(k\sqrt{\omega})\|_{Q}=0$, where $f^{\delta}(\cdot)$ is the dynamical solution of (4) with $s\geq 1$, and $\{f^{\delta}_{k}\}_{k}$ is the sequence, generated by the Nesterov’s scheme with parameters $(\alpha,\omega)$, see formula (110) for details.

It should be noted that the second order dynamic (4) has recently been investigated in [4], where they have proven that the flow (1) with the vanishing initial data yields an optimal regularization method for the linear operator equation (1). In this paper, cf. Section 2, we focus on the acceleration effect in the sense of regularization theory of (1) with an arbitrary initial guess $f_{0}$. The main result of this paper regarding the discretized version of (4) is presented in Section 3, where we demonstrate that by using damped symplectic integrators, the regularization property and acceleration effect under optimal convergence rates of (4) carry over to its numerical realization. In Section 4, the developed accelerated iterative regularization methods, equipped with a posteriori stopping rule, are applied to a diffusion-based bioluminescence tomography, which can be formulated as:

Problem 1**.**

Given $g_{1}$ and $g_{2}$, find a bioluminescent source $f$ such that the solution $u$ of the boundary-value problem ($\Omega_{0}$ and $\Omega$ are bounded open domains in $\mathbb{R}^{n}$ ($n=2,3$))

[TABLE]

satisfies

[TABLE]

Some numerical examples, as well as a comparison with three well-known existing iterative regularization methods, are presented in Section 5. Concluding remarks are given in Section 6. Some proof details as well as some details on finite element discretization are postponed to the appendices.

2 Analysis of the continuous regularization method

2.1 Convergence analysis

We start with the convergence analysis of the continuous method (4) in the sense of regularization theory. Let $\{\lambda_{j};u_{j},v_{j}\}_{j=1}^{\infty}$ be the well-defined singular system for the compact linear operator $K$, i.e. we have $Ku_{j}=\lambda_{j}v_{j}$ and $K^{*}v_{j}=\lambda_{j}u_{j}$ with ordered singular values $\|K\|=\lambda_{1}\geq\lambda_{2}\geq\cdot\cdot\cdot\geq\lambda_{j}\geq\lambda_{j+1}\geq\cdot\cdot\cdot\to 0$ as $j\to\infty$. Since the eigenelements $\{u_{j}\}_{j=1}^{\infty}$ forms a orthogonal basis in $N(K)^{\bot}\subset Q$, any dynamical element $\hat{f}^{\delta}(t)\in Q$ has a decomposition $\hat{f}^{\delta}(t)=\sum_{j}\xi_{j}(t)u_{j}+\hat{f}_{0}(t)$, where $\hat{f}_{0}(t)\in N(K)$. As we are interested in an appropriate stable approximation of $f^{\dagger}$ – the minimum norm least square solution of (1), the designed regularized solution can be chosen in the form $f^{\delta}(t)=\sum_{j}\xi_{j}(t)u_{j}$ by noting that $\|f^{\delta}(t)\|\leq\|\hat{f}^{\delta}(t)\|$ for all $\hat{f}_{0}(t)\in N(K)$. Since $f^{\delta}(t)$ satisfies the dynamic (4), we have

[TABLE]

which implies the evolution equations for coefficient $\xi_{j}(t)$:

[TABLE]

Proposition 1**.**

Let $s>-1/2$ be a fixed number. Then, the differential equation (14) has a unique solution

[TABLE]

where $\Gamma(\cdot)$ and $J_{s}(\cdot)$ denote the Gamma function and the Bessel function of first kind of order $s$ respectively.

The proof of the above proposition can be found in Appendix A. By Proposition 1 and the decomposition $f^{\delta}(t)=\sum_{j}\xi_{j}(t)u_{j}$ we obtain the explicit formula for the solution of (11) as

[TABLE]

where (we identify $\lambda^{2}_{j}$ as $\lambda$)

[TABLE]

Theorem 1**.**

Let $f^{\delta}(t)$ be the dynamic solution of (4) with $s>-1/2$. Then, if the terminating time $T_{*}=T_{*}(\delta,y^{\delta})$ is chosen so that

[TABLE]

the approximate solution $f^{\delta}(T_{*})$ converges (strongly) to $f^{\dagger}$ as $\delta\to 0$.

Proof.

Let $f(t)$ be the solution of (4) with noise-free data, i.e., $f(t)=(1-K^{*}Kg(t,K^{*}K))f_{0}+g(t,K^{*}K)K^{*}y$. Furthermore, define the bias function by

[TABLE]

Obviously, $r(t,\lambda)$ is the unique solution to

[TABLE]

Then, with the help of the intermediate quantity $f(t)$ and bias function $r(t,\lambda)$, we obtain the well-known error estimates

[TABLE]

by noting that $y=Kf^{\dagger}$. Hence, to prove the convergence of the full regularization error, we have to show the convergence of both two terms in the right-hand side of (2.1).

Let’s first consider the estimate for $\|r(t,K^{*}K)(f_{0}-f^{\dagger})\|$. To this end, define the Lyapunov function of (21) by $E(t):=\dot{r}^{2}(t,\lambda)+\lambda r^{2}(t,\lambda)$. Since

[TABLE]

$E(t)$ is a non-increasing function, and consequently, we have $\lambda r^{2}(t,\lambda)\leq E(t)\leq E(0)=\lambda$, which implies

[TABLE]

On the other hand, by asymptotic [1, (9.2.1)]

[TABLE]

for any fixed $\lambda>0$, we obtain together with (23) and condition $s>-1/2$ that

[TABLE]

Now, consider the quality $\sqrt{\lambda}g(t,\lambda)$. To this end, we introduce the function $\zeta(\tau)=\zeta(\sqrt{\lambda}t)=r(t,\lambda)$. By the representation of the Bessel functions, see e.g. [1, (9.5.10)], we have $\zeta(\tau)=\prod^{\infty}_{k=1}\left(1-\frac{\tau^{2}}{\jmath^{2}_{s,k}}\right)$, where $\jmath_{s,k}$ denotes the $k$-th positive zero of Bessel functions $J_{s}$ (sorted in increasing order). Hence, for all $\tau\in(0,\jmath_{s,1})$,

[TABLE]

which means that $\zeta$ is monotonically decreasing on $(0,\jmath_{s,1})$. According to the the initial condition (21), we have $\zeta(0)=1$ and $\zeta^{\prime}(0)=0$, which implies together with the asymptotic (115) that a number $\tau_{0}\in(0,\jmath_{s,1})$ exists such that for all $\tau\in[0,\tau_{0}]$: $\zeta(\tau)\geq 1-\tau$. Setting $\tau_{1}:=\min\left\{1/2,\tau_{0}/4\right\}$, we obtain

[TABLE]

On the other hand, by (23), we deduce that

[TABLE]

Combine (26), (27) and the relations $\zeta(\tau)=\zeta(\sqrt{\lambda}t)=r(t,\lambda)$ to obtain

[TABLE]

Consequently, we have

[TABLE]

which implies that $\delta\sup_{\lambda>0}\sqrt{\lambda}g(t,\lambda)\leq\delta\frac{t}{2\tau_{1}}\to 0$ under the choice of terminating time in (17).

∎

2.2 Convergence rate and acceleration

The purpose of this subsection is to show that (4) with an appropriate terminating time yields an accelerated optimal regularization method. It is well-known that in order to prove the convergence rate for the approximate solution $f^{\delta}(t)$, additional smoothness assumptions on $f^{\dagger}$ in correspondence with the forward operator $K$ have to be fulfilled. For simplicity, in this paper, we only consider the Hödler type source conditions, i.e.

Assumption 1**.**

There exists an element $v_{0}$ and a number $\rho\geq 0$ such that

[TABLE]

By using the technique of the comparison of two qualifications, cf. [17, Def. 2] and [17, Prop. 3, Remark 5 and Lemma 2], the results in this subsection can be easily extended to general range-type source conditions such as the logarithmic source conditions.

Theorem 2**.**

*(A priori choice of the terminating time)

Let $f^{\delta}(t)$ be the solution of (4) with $s>-1/2$. Then, under the Assumption 1,*

• (i)

if $\mu\in(0,\frac{1+2s}{4}]$ and $T_{*}=\delta^{-\frac{1}{2\mu+1}}$, we have the order optimal convergence rate

[TABLE]

• (ii)

if $\mu>\frac{1+2s}{4}$ and $T_{*}=\delta^{-\frac{2}{2s+3}}$, we have the reduced convergence rate

[TABLE]

Proof.

According to (24), there exists a pair of numbers $(C_{0},T_{0})$ such that for all $t\geq T_{0}$: $J_{r}(\sqrt{\lambda}t)\leq C_{0}\lambda^{-1/4}t^{-1/2}$, which implies together with (18), (23), and Assumption 1 that

[TABLE]

where $C_{1}=C_{0}2^{s}\Gamma(s+1)\rho\max\left\{\|K\|^{2\mu-\frac{1+2s}{2}},1\right\}$. We complete the proof by the following inequalities

[TABLE]

due to (16), (28), (37) and the choice of $T_{*}$. ∎

Remark 1**.**

Denote by $\theta=\frac{4\mu}{1+2s}$. Then, assertion (ii) of Theorem 2 can be reformulated as follows: in the case $\mu>\frac{1+2s}{4}$ (or equivalently, $\theta>1$), if the terminating time is chosen by $T_{*}=\delta^{-\frac{\theta}{2\mu+\theta}}$, we have the convergence rate $\|f^{\delta}(T_{*})-f^{\dagger}\|=\mathcal{O}\left(\delta^{\frac{2\mu}{2\mu+\theta}}\right)$ as $\delta\to 0$. This means that the small choice of model parameter $s$ will reduce the accuracy of the estimated approximate solution of our regularization method (4). However, if the value of $s$ is not small, i.e. $\mu<s+1$, (4) still offers an accelerated regularization method by noting that $T_{*}=\delta^{-\frac{\theta}{2\mu+\theta}}<T_{Landweber}=\delta^{-\frac{2}{2\mu+1}}$.

The following proposition indicates that $\mu=\frac{1+2s}{4}$ cannot be a qualification in the sense of regularization theory [7].

Proposition 2**.**

Let $\mu\in(0,(1+s)/2)$. Then $\mu$ cannot be a qualification of method (4) under source conditions (29).

Proof.

We prove it by contradiction. According to [7, § 4.2], it is necessary to show that for any positive constants $c,\gamma$ there exists numbers $T_{1}$ and $\lambda\in[ct^{-2},\|K\|^{2}]$ such that the inequality

[TABLE]

cannot be hold for all $t\geq T_{1}$. Indeed, according to the asymptotic (24), there exist two constants $C$ and $T_{0}$ such that for all $t\geq T_{0}$,

[TABLE]

On the other hand, for

[TABLE]

there must exist an integer $k_{1}>\pi^{-1}(C2^{s}\Gamma(s+1)\gamma^{-1})^{2/(1+s-2\mu)}-1/2$ such that $k_{1}\pi+\pi/2\in[c,\|K\|^{2}t^{2}]$ for any $t\geq T_{1}$, and for $\lambda_{1}:=(k_{1}\pi+\pi/2)t^{-2}\in[ct^{-2},\|K\|^{2}]$ (according to (43)),

[TABLE]

which contradicts the inequlaity (42). ∎

Now, let us turn to the a posteriori choice of the terminating time $T_{*}=T_{*}(\delta,y^{\delta})$. We consider Morozov’s discrepancy principle as the most prominent version which exploits zeros of the discrepancy function

[TABLE]

where $\tau>1$ is a fixed parameter.

Lemma 1**.**

If $\|Kf_{0}-y^{\delta}\|>\tau\delta$, then $\chi(T)$ has at least one solution.

The proof of the above lemma can be found in Appendix B. If the function $\chi(T)$ has more than one root, we recommend selecting $T_{*}$ from the rule

[TABLE]

In other words, $T_{*}$ is the first time point for which the size of the residual $\|Kf^{\delta}(t)-y^{\delta}\|$ has about the order of the data error. By Lemma 1 such $T_{*}$ always exists. Furthermore, by the proof of Lemma 1 in Appendix B, it is easy to show that $\chi(T)$ is bounded by a monotonically decreasing function $\mathcal{E}(t)-\tau\delta\searrow-\tau\delta$ as $t\to\infty$. Hence, roughly speaking, the trend of $\chi(T)$ is to be a decreasing function, where oscillations may occur.

Theorem 3**.**

(A posteriori choice of the terminating time) Let the terminating time of $f^{\delta}(T_{*})$ of (4) be chosen as a root of the discrepancy function (44). Then, under Assumption 1,

• (i)

if $\mu\in(0,\frac{2s-1}{4}]$, we have the order optimal error estimates

[TABLE]

• (ii)

If $\mu>\frac{2s-1}{4}$, we have the reduced error estimates

[TABLE]

The proof of Theorem 3 can be done by a standard argument in regularization theory, and we omit it here. By Theorems 2 and 3, for the proposed continuous regularization method (4) the optimal convergence rates can be obtained with approximately the square root of iterations than would be needed for the conventional asymptotical regularization method [7, § 6.2], which means that our method is an accelerated regularization method. However, similar to the existing accelerated order-optimal regularization methods (e.g. $\nu$-method [7, § 6.3], Nesterov’s method [21], and fractional asymptotical regularization [33]), the proposed method (4) also shows a saturation phenomenon; i.e., the optimal convergence rate $\|f^{\delta}(T_{*})-f^{\dagger}\|=\mathcal{O}\left(\delta^{\frac{2\mu}{2\mu+1}}\right)$ and the asymptotic $T_{*}=\mathcal{O}\left(\delta^{-\frac{1}{2\mu+1}}\right)$ holds only for $\mu\in(0,\frac{2s+1}{4}]$ or $\mu\in(0,\frac{2s-1}{4}]$ according to the choice of the terminating time. Therefore, the choice of $s$ is crucial for applying the regularization method (4). The a priori knowledge of the degree of smoothness of unknown exact solution provides a low bound for the model parameter $s$. It should be noted that, though theoretically, the large value of $s$ can extend the region of optimal convergence rate and increase the convergence rate in the case of high smooth exact solution, in practice, too large value of $s$ may decrease the accuracy of the obtained approximate solution since the constant in the asymptotic $\mathcal{O}\left(\delta^{\frac{2\mu}{2\mu+1}}\right)$, e.g. $C_{1}$ in (37), blows up as $s\to\infty$. A detailed discussion will be given in Section 5.1.

3 A new class of accelerated iterative regularization methods

The evolution equation (4) with an appropriate numerical discretization scheme for the artificial time variable yields a concrete iterative method. This has motivated us to develop some novel iterative regularization methods based on the continuous method (4). The goal of this section is to realize this idea.

Just as with the Runge-Kutta integrators [25] or the exponential integrators [Hochbruck-1998] for numerically solving first order equations, the damped symplectic integrators are extremely attractive for solving (4), since the schemes are closely related to the canonical transformations [12], and the trajectories of the discretized second flow are usually more stable for its long-term performance. To this end, let us start with the simplest symplectic scheme – the symplectic Euler method, i.e.

[TABLE]

By elementary calculations, scheme (47) can express in the form of following three-term semi-iterative method

[TABLE]

with parameters $a_{k}=\frac{\Delta t_{k+1}}{\Delta t_{k}}\left(1-\Delta t_{k}\frac{1+2s}{t_{k}}\right)$ and $\omega_{k}=\Delta t_{k}\Delta t_{k+1}$. The adjoint scheme of (47), namely

[TABLE]

also shares the same recurrence form (48), but with parameters $a_{k}=\frac{\Delta t_{k}}{\Delta t_{k-1}}\left(1-\Delta t_{k-1}\frac{1+2s}{t_{k}}\right)$ and $\omega_{k}=\Delta t^{2}_{k}$. Obviously, when we fix the step size $\Delta t_{k}\equiv\Delta t$, iterations (47) and (49) coincide with each other . For the high order symplectic methods, we can consider the Störmer-Verlet scheme (with a constant step size $\Delta t$), which takes the form

[TABLE]

Surprisedly, the scheme (50) can also be rewritten in the form of (48), with parameters

[TABLE]

In the work, we consider the following modified Störmer-Verlet scheme

[TABLE]

where the third step in (51) is inspired by the Nesterov’s method. Here, parameters $a_{k}$ will be defined later in (53). The goal in this section is to prove that the scheme (51) with a posteriori iteration stopping rule yields an accelerated iterative regularization method.

Remark 2**.**

(a) At the beginning of iteration in (51), one can set $t_{0}:=t_{1}$ to avoid the singularity. (b) By adding $v^{k+1}=f^{k+1}+2\Delta ta_{k+1}q^{k+\frac{1}{2}}$ in the symplectic Euler method (50), one can obtain a new iteration scheme with a lower computational cost in the inner iterations. The convergence analysis is exactly the same as for the scheme (51). The behaviour of this method is similar to the Nesterov’s method, and thus, we omit it here.

Unlike the original the symplectic Euler method and Störmer-Verlet, the scheme (51) expresses the following recurrence form

[TABLE]

with parameters

[TABLE]

Without loss of generality, for $s>1/2$ define that

[TABLE]

Consequently, $\omega_{k}\geq\Delta t^{2}/2$ for all $k\in\mathbb{N}$.

As the Landweber iterates, according to (52), the iterates $f^{k}$ of (51) obviously belong to the Krylov subspace $\textrm{Span}\left\{K^{*}y^{\delta},\cdots,(K^{*}K)^{k-1}K^{*}y^{\delta}\right\}$. Therefore, the solution $f^{k}$ of (51) can be written as $f^{k}=g_{k}(K^{*}K)K^{*}y^{\delta}$, where $g_{k}$ is a polynomial of degree $k-1$, and the residual polynomials

[TABLE]

exhibit the following property.

Proposition 3**.**

Assume that $s>1/2$ and $\Delta t\in(0,\sqrt{2}/\|K\|)$. Then, the residual polynomials of scheme (51) satisfy the following inequality

[TABLE]

where $c_{1}=2^{3\mu}s^{2\mu}\Delta t^{-2\mu}$ and

[TABLE]

where $c_{2}=2^{3\mu}\Delta t^{-2\mu}(\mu-1/2)^{\mu-1/2}s$.

Proof.

This proof uses the technique in [21]. By using (52) and elementary calculations, the residual polynomials of (51) satisfy the recurrence relation

[TABLE]

which can be rewritten as

[TABLE]

where $a_{k}=\frac{\theta_{k}}{\theta_{k-1}}(1-\theta_{k-1})$, and

[TABLE]

Now, let us show that the sequence $(\theta_{k})_{k}$, defined in (59), satisfies the following inequality

[TABLE]

By definitions of $\omega_{k}$ and $\theta_{k}$ in (53), (54) and (59), we have under the assumption $s>1/2$ that

[TABLE]

which yields the inequality (60). The inequalities in the bracket of (61) are used when $k\geq k_{s}$.

Now, denote by $\hat{r}_{k}:=r_{k-1}+\frac{1}{\theta_{k-1}}(r_{k}-r_{k-1})$. Then, we derive together with (58) and (60) that

[TABLE]

by noting that $\theta_{0}=\hat{r}_{1}=1$ (as $r_{-1}=r_{0}\equiv 1$). Inequality (62) immediately yields

[TABLE]

as well as $|\hat{r}_{k+1}|\leq 1$. The latter inequality together with the recurrence (58) and initial data $r_{-1}=r_{0}\equiv 1$ implies

[TABLE]

If $\mu\in(0,1/2]$, (63), (64) and the definition of $\theta_{k}$ in (59) immediately gives

[TABLE]

If $\mu>1/2$, we obtain together with (63) that (let $k\geq 2$)

[TABLE]

∎

By inequality (63), it is not difficult to show that the following limit

[TABLE]

holds for all fixed $\lambda\in(0,\|K\|^{2}]$ and $\mu\geq 0$. Then, Based on the relation (65), Proposition 3 and standard argument for linear regularization theory, see e.g. [21, Theorems 3.1 and 4.1] or [7], we have the following convergence rate results.

Theorem 4**.**

Suppose that $s>1/2$ and $\Delta t\in(0,\sqrt{2}/\|K\|)$. Let $f^{k}$ be the approximate solution, generated by the scheme (51). Then, under Assumption 1,

• •

if $\mu\in(0,1/2]$ and $k^{*}=\mathcal{O}(\delta^{-\frac{1}{2\mu+1}})$, we have the convergence rate

[TABLE]

• •

If $\mu>1/2$ and $k^{*}=\mathcal{O}(\delta^{-\frac{1}{2\mu+3}})$, we have the convergence rate

[TABLE]

• •

For general positive $\mu$, if the iteration of (51) is terminated according to the discrepancy principle (with a fixed positive parameter $\tau$), i.e.

[TABLE]

then, it holds that

[TABLE]

We end this section by offering a few remarks.

Remark 3**.**

Unlike in the continuous situation, the a priori stopping rules in the first two assertions are not optimal. Consequently, the convergence rates in (66) and (67) are not optimal for our discretized regularization method (51). Indeed, similar to [21, Theorems 3.1], one can show that under the following a priori stopping rule (Obviously, it is not realizable, since $f^{\dagger}$ is not known):

[TABLE]

if $\mu\in(0,1/2]$, we have

[TABLE]

If $\mu>1/2$, we have

[TABLE]

Remark 4**.**

It is not difficult to show that Theorem 4 also holds for the scheme (51) with $\Delta t$ replaced by $\Delta t_{k}$ such that $(\Delta t_{k})_{k}\subset[\Delta t_{min},\sqrt{2}/\|K\|)$, where $\Delta t_{min}>0$ is a constant.

Remark 5**.**

It should be noted that the combination of the Nesterov’s method and the non-symplectic schemes for our second order asymptotical regularization (4) may also provide an accelerated iterative regularization method. For example, consider the following scheme (it is not a symplectic method as it belongs to explicit numerical scheme)

[TABLE]

where

[TABLE]

The above scheme can also be written in the form of (52), but with parameters

[TABLE]

It is not difficult to show that Proposition 3, and hence Theorem 4, also holds for the scheme (71). Consequently, iteration (71) also offers an accelerated iterative regularization method.

4 Application to the diffusion-based bioluminescence tomography (BLT)

4.1 Background of BLT and a reduced mathematical model

In the modern world, biomedical imaging has become extremely important not only for patient care but also for the study of biological structure and function, and for addressing fundamental questions in biomedicine. In molecular imaging, small animal organs and tissues are often labeled with reporter probes that generate detectable signals that can be tracked outside a living body. This technology has been widely used in clinical medicine for investigating tumorigenesis, cancer metastasis, cardiac diseases, etc. In comparison with traditional biomedical imaging approaches such as X-ray computed tomography, positron emission tomography and ultrasound and magnetic resonance imaging, optical molecular imaging has attracted considerable attention for its cost-effectiveness and performance as it directly reveals molecular and cellular activities sensitively [6]. Among various optical molecular imaging techniques, fluoresence molecular imaging [22] and bioluminescence imaging (BLI) [24] are among the most widely used in practice. In contrast with fluorescence imaging, there is no inherent tissue autofluorescence generated by external illumination in bioluminescence imaging, which makes it extremely sensitive. However, BLI is primarily qualitative and cannot provide information about the distribution of an in vivo bioluminescent source. For the problem of reconstructing an internal bioluminescent source from the measured bioluminescent signal on the external surface of a small animal, a quantitative prototype, termed as bioluminescence tomography (BLT), is introduced [13].

Bioluminescent photon propagation in biological tissue is governed by the radiative transfer equation (RTE) which has been utilized as the forward model for bioluminescence tomography [19]. However, the RTE is highly dimensional and presents a serious challenge for its accurate numerical simulations given the current level of development in computer software and hardware. Because the mean-free path of the photon is between 500 nm and 1000 nm in biological tissues, which is very small compared to the size of a typical object in this context, the predominant phenomenon in BLT is scattering, which provides a diffusion approximation of the RTE by the following reduced mathematical model [13]

[TABLE]

where $u$ denotes the (direction-averaged) photon density, $D=[3(\mu_{a}+\mu^{\prime}_{s})]^{-1}$ with $\mu_{a}$ and $\mu^{\prime}_{s}$ being the absorption and reduced scattering coefficients. The boundary $\Gamma$ of the domain $\Omega\subset\mathbb{R}^{n}$ ($n=2,3$) is assumed to be Lipschitz continuous. $\partial_{\nu}$ is the outward normal differentiation operator. $\Omega_{0}\subset\Omega$ is known as a permissible region of the source function, and $\chi$ is the indicator function such that $\chi_{\Omega_{0}}(x)=1$ for $x\in\Omega_{0}$, while $\chi_{\Omega_{0}}(x)=0$, when $x\not\in\Omega_{0}$. $g^{-}$ is an incoming flux on G and it vanishes when the imaging is implemented in a dark environment. $A=\frac{1+R(x)}{1-R(x)}$ with $R(x)\approx-1.4399\,\gamma(x)^{-2}+0.7099\,\gamma(x)^{-1}+0.6681+0.0636\,\gamma(x)$ and $\gamma(x)$ being the refractive index of the medium at $x\in\Gamma$. In the case when $\Omega$ is a unit circle centered at the origin, $\mu_{a}=0.04,\mu^{\prime}_{s}=1.5$ , and $A=3.2$ with refractive index $\gamma=1.3924$. In BLT, the measurement is the outgoing flux density on the boundary:

[TABLE]

If we denote by $g_{1}:=g^{-}+2A\,g$ and $g_{2}:=-g$, then the BLT problem (75)-(76) can be formulated as Problem 1, i.e., the problem (7)-(8). This inverse source problem has been intensively studied in [5, 9, 10, 11, 13, 14, 26, 30] and referenced therein. The essential methodology in these studies is to solve the inverse problem by a two-step strategy. The first step is to adopt Tikihonov variational regularization with a priori regularization parameter choice rule to overcome the ill-poseness of original inverse problem, and then solve the regularized PDE-controlled optimization problem by a numerical algorithm (usually we adopt an iterative method). The defects of these existing methods are: (a) Tikihonov regularization exhibits a “strong” saturation phenomenon, i.e., the optimal convergence rate is limited by $\mathcal{O}(\delta^{2/3})$ with respect to Hölder-type source condition and noise level $\delta$ of data. (b) The a priori stopping rule of regularization parameter is not realistic in practice as it requires some knowledge of the unknown exact solution. (c) Especially for large-scale inverse problems, variational regularization methods are time consuming. In order to overcome these shortcomings, we shall apply the developed accelerated iterative regularization method (51) for the fast solution of Problem 1. It should also be noted that, recently, by assuming the sourcewise representation of source function $f^{\dagger}$, the authors in [31] combined the coupled complex boundary method and the expanding compacts method to propose a new regularization method that can calculate a posteriori error estimate efficiently. However, no convergence rate can be derived for such a method.

4.2 Analysis of a mathematical formulation

The aim of this subsection is to reformulate the inverse source problem (7)-(8) as an abstract operator equation (in a relaxed weak form) so that we can adopt the developed accelerated iterative regularization method with the a posteriori stopping rule in the previous section. We start with the basic assumptions on the system parameter.

Assumption 2**.**

$D\in L^{\infty}(\Omega)$* and $D(x)\geq D_{0}$ for almost every $x\in\Omega$; $\mu_{a}\in L^{2}(\Omega)$ and $\mu_{a}(x)\geq\mu_{0}$ for almost every $x\in\Omega$. Here, $D_{0}$ and $\mu_{0}$ are two positive constants. Moreover, $\Gamma\subset\mathbb{R}^{n-1}$ is a open bounded set.*

Denote by

[TABLE]

where

[TABLE]

is the weighted $H^{1}(\Omega)$ norm. Moreover, we introduce the norm of the trace space $V^{1/2}(\Gamma)$ by

[TABLE]

where $\gamma_{0}:V\to V^{1/2}(\Gamma)$ denotes the standard trace operator. The space $V^{-1/2}(\Gamma)$ is defined as dual of $V^{1/2}(\Gamma)$, with the norm given

[TABLE]

It is not difficult to show that all of $V$, $V_{0}$, $V^{1/2}(\Gamma)$ and $V^{-1/2}(\Gamma)$ are Banach spaces, equipped with the corresponding norms (78), (79) and (80), respectively. We remark that if $D=\mu_{a}\equiv 1$, $V$, $V^{1/2}(\Gamma)$ and $V^{-1/2}(\Gamma)$ are reduced to the standard Sobolev spaces $H^{1}(\Omega)$, $H^{1/2}(\Gamma)$ and $H^{-1/2}(\Gamma)$, respectively. For simplicity, denote $Q_{0}=L^{2}(\Omega_{0})$, $Q=L^{2}(\Omega)$, and $Q_{\Gamma}=L^{2}(\Gamma)$. Set $V_{g_{1}}:=\{v\in V:v=g_{1}\,{\rm on}\ \Gamma\}$. Define

[TABLE]

Then $a(\cdot,\cdot)$ is symmetric, continuous and coercive on $V$. Therefore, by the Lax-Milgram Lemma ([8]), for any $f\in Q_{0}$, the problems

[TABLE]

and

[TABLE]

each have a unique solution. Moreover, a constant $c>0$ exists such that

[TABLE]

If we define

[TABLE]

we obtain that $u_{D}(f,g_{1})=u_{D}(f)+\widetilde{u}_{D}(g_{1})$ and $u_{N}(f,g_{2})=u_{N}(f)+\widetilde{u}_{N}(g_{2})$.

Define two operators $K_{D}$ and $K_{N}$ from $Q_{0}$ to $V$ by

[TABLE]

Furthermore, define

[TABLE]

It is easy to verify that for any $f\in Q_{0}$,

[TABLE]

Therefore, $K\,f=y$ means that $u_{D}(f,g_{1})=u_{N}(f,g_{2})$. In other words, the original BLT problem is equivalent to the following problem (in the sense of weak form): find $f\in Q_{0}$ such that

[TABLE]

Proposition 4**.**

The operator $K:Q_{0}\to V$ is compact.

The proof of Proposition 4 can be found in Appendix C. Now, let us consider the case with inexact measurement. Suppose that instead of exact boundary data $\{g_{1},g_{2}\}$, we are given noisy data $\{g^{\delta}_{1},g^{\delta}_{2}\}$ satisfying the following assumption.

Assumption 3**.**

Let $g_{1},g^{\delta}_{1}\in V^{1/2}(\Gamma)$ and $g_{2},g^{\delta}_{2}\in V^{-1/2}(\Gamma)$ such that

[TABLE]

where the noise level $\delta>0$ is known.

Remark 6**.**

In Assumption 3 the data space is assumed to be the trace spaces $V^{1/2}(\Gamma)\times V^{-1/2}(\Gamma)$, which is designed for noise-free boundary data of BLT problem. Consequently, the noise in Assumption 3 is not entirely random, and it also belongs to the considered data space $V^{1/2}(\Gamma)\times V^{-1/2}(\Gamma)$. However, in practice, the original measured data may hardly approximate true Dirichlet data in $V^{1/2}(\Gamma)$ as the noise usually exhibits a weaker regularity, e.g. $g^{\delta}_{1}$ usually only belongs to $L^{2}(\Gamma)$. In this case, one can use a smoothing technique to obtain a valid mollification $g^{\delta,\epsilon}_{1}$ of noisy Dirichlet data with small enough $\epsilon=\epsilon(\delta)>0$ such that

[TABLE]

In this paper, we use the mollification $g^{\delta,\epsilon}_{1}$, defined through the convolution smoother, i.e.

[TABLE]

where the mollifier $\eta_{\epsilon}$ is defined by

[TABLE]

with the positive constant $C$ chosen such that $\int_{\mathbb{R}^{n-1}}\eta(x)dx=1$. According to [2, Theorem 2.29], if ${\rm supp}(g^{\delta}_{1})\Subset\Gamma$ and ${\rm dist}({\rm supp}(g^{\delta}_{1}),\partial\Gamma)>\epsilon$, we have $g^{\delta,\epsilon}_{1}\in C^{\infty}_{0}(\Gamma)\subset V^{1/2}(\Gamma)$. It should be noted that when $\partial\Gamma=\emptyset$, e.g. $\Gamma$ is a circle or a sphere, ${\rm dist}({\rm supp}(g^{\delta}_{1}),\partial\Gamma)=+\infty$. Moreover, in the case $g^{\delta}_{1}\in L^{2}(\Gamma)$, $\|g^{\delta,\epsilon}_{1}-g^{\delta}_{1}\|_{L^{2}(\Gamma)}\to 0$ as $\epsilon\to 0$.

Proposition 5**.**

Under Assumption 3, it holds $\|y^{\delta}-y\|_{V}\leq\delta$, where $y^{\delta}=\widetilde{u}_{N}(g^{\delta}_{2})-\widetilde{u}_{D}(g^{\delta}_{1})$.

Proof.

Define $v_{D}:=\widetilde{u}_{D}(g^{\delta}_{1})-\widetilde{u}_{D}(g_{1})$ and $v_{N}:=\widetilde{u}_{N}(g^{\delta}_{2})-\widetilde{u}_{N}(g_{2})$. Then, $v_{D}$ and $v_{N}$ satisfy the following BVPs

[TABLE]

and

[TABLE]

Now, let us show that

[TABLE]

By the definition (79), we have

[TABLE]

On the other hand, according to equation (91), we have together with (78) that for any $v\in V$: $\langle v_{D},v\rangle_{V}=\int_{\Gamma}D\frac{\partial v_{D}}{\partial\mathbf{n}}vds$, which implies that

[TABLE]

By the trace theorem, we have $\left\|D\frac{\partial v_{D}}{\partial\mathbf{n}}\right\|_{V^{-1/2}(\Gamma)}\leq\|v_{D}\|_{V}$. Hence, we derive

[TABLE]

Combine (94) and (95) to obtain the first identity in (93).

Now, consider the second identity in (93). According to (92), for all $v\in V$, we have together with (78) that

[TABLE]

Set $v=v_{N}$ to get $\|v_{N}\|^{2}_{V}=\left|\int_{\Gamma}(g^{\delta}_{2}-g_{2})\gamma_{0}v_{N}ds\right|$, which gives

[TABLE]

The above inequality together with the trace inequality, i.e. $\|\gamma_{0}v_{N}\|_{V^{1/2}(\Gamma)}\leq\|v_{N}\|_{V}$, gives

[TABLE]

On the other hand, by using (96) we deduce that

[TABLE]

which implies the second identity of (93) by noting (97).

Finally, by using the definition of $y^{\delta}$ and identities (93), we complete the proof by following inequalities

[TABLE]

∎

Next we discuss the form of $K^{*}\,(K\,f-y^{\delta})$, which is used in our main algorithms (51) and (71), in the context of the BLT problem. To this end, denote by $K^{*}_{D}$ and $K^{*}_{N}$ the adjoint operators of $K_{D}$ and $K_{N}$ such that for any $v\in V$ and $f\in Q_{0}$:

[TABLE]

Then $K^{*}:V\rightarrow Q_{0}$ is such that $K^{*}=K^{*}_{D}-K^{*}_{N}$.

For any $f\in Q_{0}$, denote by $u_{DN}(f)=K\,f-y^{\delta}=u_{D}(f,g^{\delta}_{1})-u_{N}(f,g^{\delta}_{2})$, and define $w_{D}=w_{D}(u_{DN}(f))\in V_{0}$ and $w_{N}=w_{N}(u_{DN}(f))\in V$ the solutions of the adjoint variational problems

[TABLE]

and

[TABLE]

respectively. Then $K^{*}_{D}(K\,f-y^{\delta})=w_{D}(u_{DN}(f))|_{\Omega_{0}}$ and $K^{*}_{N}(K\,f-y^{\delta})=w_{N}(u_{DN}(f))|_{\Omega_{0}}$. Thus, we have

[TABLE]

Similarly, we can give a form of the source condition (29) for our BLT problem. In fact, (29) with $\mu=1$ reads that there exists an element $v_{*}\in Q_{0}$ such that

[TABLE]

where $w_{D}^{*}$ and $w_{N}^{*}$ are the solutions of (98) and (99), both with $u_{DN}$ being replaced by $u_{D}(v_{0})-u_{N}(v_{0})$, and $u_{D}(v_{0})=u_{D}(v_{0},0),u_{N}(v_{0})=u_{N}(v_{0},0)$.

5 Numerical experiments

In this section, we devote ourselves to presenting some numerical examples for demonstrating the effectiveness of the proposed accelerated iterative regularization method (51). We take the diffusion-based bioluminescence tomography considered in Section 4 as example. To that end, with the problem domain $\Omega$, parameters $\mu_{a},\mu^{\prime}_{s},A$, Robin data $g^{-}$, and a prescribed true source function $f^{*}$, we solve the forward BVP (75) to get $u^{*}$. A finite element method of solving (75) is briefly discussed in Appendix D.

The outgoing flux density and the Cauchy data on the boundary are

[TABLE]

Uniformly distributed noises with the relative error level $\delta^{\prime}$ are added to $g$ to get $g^{\delta}$

[TABLE]

where rand$(x)$ returns a pseudo-random value drawn from a uniform distribution on $[0,1]$. The corresponding noisy Cauchy data are $g^{\delta}_{1}=g^{-}+2A\,g^{\delta}$ and $g^{\delta}_{2}=-g^{\delta}$. Then the noise level of the measurement data is calculated by $\delta=\|y^{h,\delta}-y^{h}\|_{V}$, with $y^{h}=\widetilde{u}^{h}_{N}(g_{2})-\widetilde{u}^{h}_{D}(g_{1})$ and $y^{h,\delta}=\widetilde{u}^{h}_{N}(g^{\delta}_{2})-\widetilde{u}^{h}_{D}(g^{\delta}_{1})$. Here and later on, the superscript h (or the subscript h) denotes the linear finite element approximation of an element on a consistent triangulation, i.e. $\widetilde{u}^{h}_{N}$ and $\widetilde{u}^{h}_{D}$ are defined on the same triangulation with maximum triangle diameter $h$, see Appendix D for more details. Without loss of generality, in this section, let $\mu_{a}=0.04,\mu^{\prime}_{s}=1.5,D=1/[3(\mu_{a}+\mu^{\prime}_{s})]$, $A=3.2$, and $g^{-}=0$, which means the imaging is implemented in a dark environment.

Then, with the noisy data $g^{\delta}_{1}$ and $g^{\delta}_{2}$, properly chosen parameters, e.g. $s$ and $\Delta t$, approximate sources $f^{k}$ are computed by the proposed accelerated iterative regularization method (51). For the BLT problem, (51) is reduced to

[TABLE]

where $w^{k}_{D}$ and $w^{k}_{N}$ are the solutions of (98) and (99) respectively, both with $u_{DN}(f)$ replaced by $u_{DN}(f^{k})$. $w^{k+1}_{D}$ and $w^{k+1}_{N}$ have similar definitions. In the following, for the conciseness of the statements, we only consider the case that using Morozov’s discrepancy principle (68) to control the iterative procedure, namely that the iteration stops when $\|y^{h,\delta}-Af^{h,k}\|_{V}=\|u^{h}_{D}(f^{h,k},g^{\delta}_{1})-u^{h}_{N}(f^{h,k},g^{\delta}_{2})\|_{V}\leq\tau\delta$, where $u^{h}_{D}(f^{h,k},g^{\delta}_{1})$ and $u^{h}_{N}(f^{h,k},g^{\delta}_{2})$ are the finite element solutions of (82) and (83), both with $f$ being replaced by $f^{h,k}$, and with $g_{1}$ and $g_{2}$ being replaced by $g^{\delta}_{1}$ and $g^{\delta}_{2}$, respectively. Moreover, the initial guess $f_{0}$ is chosen so that the condition of Lemma 1 is satisfied: $\|y^{h,\delta}-Af^{h}_{0}\|_{V}=\|u^{h}_{D}(f^{h}_{0},g^{\delta}_{1})-u^{h}_{N}(f^{h}_{0},g^{\delta}_{2})\|_{V}>\tau\delta$, where $u^{h}_{D}(f_{0},g^{\delta}_{1})$ and $u^{h}_{N}(f^{h}_{0},g^{\delta}_{2})$ have similar definitions as $u^{h}_{D}(f^{h,k},g^{\delta}_{1})$ and $u^{h}_{N}(f^{h,k},g^{\delta}_{2})$ above.

We use $N_{\max}:=50000$ as the maximal number of iterations where the iteration (102) stops in all of simulations. To assess the accuracy of the approximate solutions, we define the finite element approximate $L^{2}$-norm relative error for an approximate solution $f^{k}$: L2Err${}_{k}:=\|f^{k}_{h}-f^{*}_{h}\|_{0,\Omega_{0}}/\|f^{*}_{h}\|_{0,\Omega_{0}}$. Obviously, $\|f^{k}_{h}-f^{*}_{h}\|_{0,\Omega_{0}}/\|f^{*}_{h}\|_{0,\Omega_{0}}\to\|f^{k}-f^{*}\|_{0,\Omega_{0}}/\|f^{*}\|_{0,\Omega_{0}}$ as $h\to 0$. All experiments in Subsection 5.1–5.2 are implemented for the following two examples:

Example 1: $\Omega:=\{(x_{1},x_{2})\in\mathbb{R}^{2}|\,x_{1}^{2}+x^{2}_{2}<1\}$, $f^{*}(x_{1},x_{2})=(1+x_{1}+x_{2})\chi_{\Omega_{0}}$ with $\Omega_{0}:=\{(x_{1},x_{2})\in\mathbb{R}^{2}|\,-0.5<x_{1},x_{2}<0.5\}$. The measurements are computed on a mesh with mesh size $h=0.01386$, 144929 nodes and 288768 elements.

Example 2: $\Omega$ is the same as Example 1, $f^{*}(x_{1},x_{2})=(1+x_{1}+x_{2})\chi_{\Omega_{1}}+e^{1+x_{1}+x_{2}}\chi_{\Omega_{2}}$ with $\Omega_{1}:=\{(x_{1},x_{2})\in\mathbb{R}^{2}|\,(x_{1}+0.5)^{2}+x^{2}_{2}<0.01\}$ and $\Omega_{2}:=\{(x_{1},x_{2})\in\mathbb{R}^{2}|\,(x_{1}-0.5)^{2}+x^{2}_{2}<0.01\}$. The measurements are computed on a mesh with $h=0.01228$, 156225 nodes and 311296 elements.

For Example 1, all approximate sources are reconstructed over a mesh with mesh size $h=0.0744$, 2325 nodes and 4512 elements. For Example 2, all approximate sources are reconstructed over a mesh with mesh size $h=0.0678$, 2505 nodes and 4864 elements.

5.1 Influence of parameters

The purpose of this subsection is to explore the dependence of the solution accuracy and the convergence speed on $\tau>0$, time step size $\Delta t$, model parameter $s$, and thus to give a guide on the choices for them in practice. For focusing on the effect of these parameters on the iteration (102), we fix $\delta^{\prime}=0.1\%$ in this subsection. Moreover, in the remaining part of this section, we simply set $f_{0}=0,q_{0}=0$.

We first investigate the influence of parameter $\tau$ on the convergence rate. For this purpose, we additionally set $s=2$, and $\Delta t=0.06$ for Example 1 or $\Delta t=0.125$ for Example 2. The detailed iterative numbers $k^{*}$ and the corresponding L2-norm relative errors ‘L2Err${}_{k^{*}}$’ for different values of $\tau$ are shown in Table 1, which shows that the smaller $\tau$ is, the more the iterative number for stopping (102) is. It is no surprise because the parameter $\tau$ does not involve the computation of the approximate solutions itself. It is used in the stop criterion and affects only the iterative number at which the iteration (102) stops. Table 1 indicates using a $\tau<1$ makes the iterative number increase dramatically. In the remaining experiments, we choose $\tau=1.1$ for Example 1 and $\tau=10$ for Example 2.

Now we investigate the influence of time step size $\Delta t$ on the solution accuracy and the convergence rate. To this end, set $s=2$, and $\tau=1.1$ for Example 1, $\tau=10$ for Example 2. The iterative numbers $k^{*}$ and the corresponding L2-norm relative errors ‘L2Err${}_{k^{*}}$’ are given in Table 2, which shows that the bigger the time step size $\Delta t$ is, the faster the iteration is. The iterative number halves when the $\Delta t$ doubles. However, our experiments suggest that $\Delta t$ should not be too big, e.g. $\Delta t\leq 0.0625$ in Example 1 and $\Delta t\leq 0.125$ in Example 2. Otherwise, the iteration will blow up as it breaks the consistency of the numerical scheme. In the remaining experiments, we choose $\Delta t=0.0625$ for Example 1 and $\Delta t=0.125$ for Example 2.

Finally, we discuss the influence of the model parameter $s$ on the solution accuracy and the convergence rate. In the experiments, set $\tau=1.1,\Delta t=0.0625$ for Example 1 and $\tau=10,\Delta t=0.125$ for Example 2. The required number of iterations $k^{*}$ and the corresponding relative error L2Err${}_{k^{*}}$ for different values of $s$ are given in Table 3, which indicates that in general, the value of $s$ is neither small nor big. Specifically, on one hand, small values of $s$, e.g. $s\leq 1$, usually bring the oscillation in solution accuracy during iterations (cf. Figures 1-2); on the other hand, large values of $s$ require more iterative number but with a little improvement on the solution accuracy. It is suggested that a value of $s$ near 2 produces satisfactory results in both solution accuracy and the iterative number for both Examples 1 and 2, which coincides with the empirical results about the optimal parameter choice for the Nesterov’s method. Therefore, in the remaining experiments, we set $s=2$.

5.2 Comparison with other methods

In this subsection, we compare the behaviors regarding the solution accuracy and the convergence rate between scheme (102), the non-symplectic scheme (71), two well-known acceleration methods: the Nesterov’s method, the $\nu$-method, and the Landweber method (2). For the BLT: Problem 1, the non-symplectic (71) has the form

[TABLE]

where $b_{k+1}$ is given in (72).

In our simulations, for the BLT problem, the $\nu$-method is defined by

[TABLE]

with $\mu_{1}=0,\omega_{1}=(4\nu+2)/(4\nu+1)$ and

[TABLE]

where $\omega>0$ is the weight. Note that $\omega=1$ in the conversional $\nu$-method (cf. [7, § 6.3]) as it works for a normalized operator equation $Kf=y$ with $\|K\|\leq 1$. For our BLT problem, $\omega$ in (106) plays the role of normalization, and it can be set as $\omega=\omega_{norm}(:=1/\|K^{*}K\|)$, which can be calculated by

[TABLE]

The dependence of stopping iterative number $k^{*}$ and the corresponding relative error L2Err${}_{k^{*}}$ on $\nu$ for different values of noise level $\delta$ are plotted under $\log$-scale in Figure 3, from which we can see that for both examples, a moderate value of $\nu$ have the best efficiency in both solution accuracy and iterative number: large values of $\nu$ do not improve the solution accuracy too much while require much more iterative numbers. In both examples, the $\nu$ of the best efficiency is near 1. In Table 4, we report relative error L2Err${}_{k^{*}}$ as well as the corresponding iterations numbers of $\nu$-method with $\nu=1/4$, $\nu=1/2$, $\nu=1$ and $\nu=2$.

The Nesterov’s method is defined by ([21])

[TABLE]

with $\alpha\geq 3$ and $0<\omega\leq\omega_{norm}$. In all simulations, we choose $\alpha=3$. For Example 1, $\omega_{norm}\approx 0.005422264152263$, we set $\omega=0.005$; for Example 2, $\omega_{norm}\approx 0.021370788062004$, we set $\omega=0.02$.

The Landweber method (2) has the form

[TABLE]

with $0<\Delta t<2\omega_{norm}$. For Example 1, we set $\Delta t=2\times 0.005=0.01$; for Example 2, we set $\Delta t=2\times 0.02=0.04$.

As suggested by Subsection 5.1, in our methods (102) (termed as “ARM”) and (103) (termed as “NSS”), we set $s=2$, $\Delta t=0.0625$ for Example 1 and $\Delta t=0.125$ for Example 2. In all methods, the initial guess $f_{0}=0$ and the iterations stop when $\|y^{h,\delta}-Af^{h,k}\|_{V}=\|u^{h}_{D}(f^{h,k},g^{\delta}_{1})-u^{h}_{N}(f^{h,k},g^{\delta}_{2})\|_{V}\leq\tau\delta$ with $\tau=1.1>1$ for Example 1 and $\tau=10>1$ for Example 2. We note that when the noisy level is large, bigger $\tau$ is suggested so that iterations can stop before the solution accuracy gets worse.

The results of the simulations are presented in Table 4, from which we conclude that, with properly chosen parameters, all the above mentioned methods are stable and can produce satisfactory solutions. Moreover, on one hand, compared with the conventional Landweber method, all of the other methods produce better accuracy with considerably fewer iterations; on the other hand, for the given model problems, the proposed ARM and NSS are comparable to well-known accelerated regularization methods, i.e. the $\nu$-method and the Nesterov’s method, in both solution accuracy and convergence rate.

We note that inverse source problems with only one measurement on the boundary have infinite solutions. In the context of the BLT problem, one cannot distinguish between a strong source over a small region and a weak source over a large region. In this paper, we are interested in the minimum norm solution, which is always unique for our linear inverse problems. However, in practice, we don’t know which solution is the minimum norm solution. Therefore, in all the above experiments, we suppose that the phantom, which has been used to generate data, is just the minimum norm solution. In order to make this assumption acceptable, we are assumed to know exactly the positions of sources (i.e. the geometry of $\Omega_{0}$). With the help of appropriately selected $\Omega_{0}$ (e.g. the one we used in our examples), we found that in the noise-free case, for large $k$ (e.g. $k=1000$), $f_{k}$ is very close (in the accuracy L2Err${}_{k}<$ 1e-5) to the input source $f^{*}$. Since $f_{k}\to f^{\dagger}$ as $k\to\infty$, we can assume that $f^{\dagger}=f^{*}$.

Of course, one can also apply the proposed method to other linear inverse problems and compare the behavior of different methods. Another well-known linear inverse problem is the Cauchy problem of finding $(\phi,t)$ on unaccessible boundary $\Gamma_{u}$ from the Cauchy data $(\Phi,T)$ on accessible boundary $\Gamma_{a}$ such that the following relations hold:

[TABLE]

In contrast to the BLT problem, the Cauchy problem above admits solution uniqueness, provided a solution exists. We can expect good behavior of the proposed method for this Cauchy problem. However, for the conciseness of the paper, we omit these numerical results.

6 Conclusions

In this paper we have proposed a new class of accelerated regularization methods for solving ill-posed linear operator equations. A series of theoretical results including limiting behavior and convergence rates are proved. Moreover, as an application of the proposed method, in this paper, a model problem arising from bioluminescence tomography is discussed in detail. Since the proposed methods are comparable to the Nesterov’s acceleration method and the $\nu$-method about the convergence rate and the solution accuracy, they are promising approaches which merit further theoretical and numerical development as well as more extensive comparison to state-of-the-art methods. Similar to Nesterov’s acceleration method [20, 23] or it’s modified versions [15, 16], the introduced iterative regularization methods can also be used to solve to some non-linear ill-posed problems. However, for performing a rigorous theoretical analysis, the concept of acceleration in the sense of regularization theory should be extended so that it can be used for evaluating general non-linear regularization methods.

Acknowledgement

We express our gratitude to the anonymous referees whose valuable comments and suggestions allowed us to eliminate weak points of the manuscript and thus to improve the paper.

The work of R. Gong is supported by the Natural Science Foundation of China (No. 11401304, 11971230) and the Fundamental Research Funds for the Central Universities (No. NS2018047). The work of B. Hofmann is supported by the German Research Foundation (DFG-grant HO 1454/12-1), and the work of Y. Zhang is supported by the Alexander von Humboldt foundation through a postdoctoral researcher fellowship.

Appendix A: Proof of Proposition 1

For simplicity, we only consider the case of positive $s$. For the case $s\in(-1/2,0)$, we refer to the similar result, presented in [4, Lemma 6.1]. The general solution to (14) is

[TABLE]

where $Y_{s}$ denotes the Bessel functions of second kind. In order to determine the constants $C^{j}_{1,s}$ and $C^{j}_{2,s}$ from the initial conditions, we distinguish three different cases: (i) $s=0$, (ii) $s\in\mathbb{N}$, and (iii) $s\not\in\mathbb{N}$. We show that for all of three cases $C^{j}_{2,s}=0$ according to the boundedness of initial data. In case (i), by using the divergence behaviour $Y_{0}(\lambda_{j}t)=\mathcal{O}(\log(\lambda_{j}t))$ as $t\to 0$ [1, (9.1.12)], $C^{j}_{2,0}$ must be zero. For case (ii), the asymptotic, cf. [1, (9.1.11)], $Y_{s}(\lambda_{j}t)=\mathcal{O}((\lambda_{j}t)^{-s})$ as $t\to 0$ implies $\frac{Y_{s}(\lambda_{j}t)}{(\lambda_{j}t)^{s}}=\mathcal{O}((\lambda_{j}t)^{-2s})$. Therefore, $C^{j}_{2,s}=0$ for all $s\in\mathbb{N}$. Now, consider the last case. According to the asymptotic $J_{-s}(\lambda_{j}t)=\mathcal{O}((\lambda_{j}t)^{-s})$ as $t\to 0$, cf. [1, (9.1.10)], $C^{j}_{2,s}=0$ for all $s\not\in\mathbb{N}$.

By the above analysis, the general solution to (14) bounded initial data should be

[TABLE]

By the initial data $f_{0}=\sum_{j}\langle f_{0},u_{j}\rangle u_{j}$ and the limit $\lim_{t\to 0}\frac{J_{s}(\lambda_{j}t)}{(\lambda_{j}t)^{s}}=\frac{1}{2^{s}\Gamma(s+1)}$, we conclude that

[TABLE]

which gives the desired formula for $\xi_{j}(t)$.

Finally, check that $\dot{\xi}_{j}(0)=0$. It can be done by the following limit

[TABLE]

by noting that [1, (9.1.10)]

[TABLE]

Appendix B: Proof of Lemma 1

This proof uses the technique in [3]. Consider the Lyapunov function of (4) by $\mathcal{E}(t)=\frac{1}{2}\|\dot{f}^{\delta}(t)\|^{2}+\|Kf^{\delta}(t)-y^{\delta}\|^{2}$. It is easy to show that

[TABLE]

Hence, $\mathcal{E}(t)$ is non-increasing, and $\mathcal{E}(\infty):=\lim_{t\to\infty}\mathcal{E}(t)$ exists by noting that $\mathcal{E}(t)\geq 0$ for all $t$. Now, consider the function $e(t)=\frac{1}{2}\|f^{\delta}(t)-f^{*}\|^{2}$. It is not difficult to obtain

[TABLE]

Divide this expression by $t$ to obtain

[TABLE]

Integrating the above inequality from $1$ to $t$ and using integration by parts for $\ddot{e}(t)$, we obtain

[TABLE]

On the one hand, using the integration by parts and the positivity of functional $e(\cdot)$, we have

[TABLE]

On the other hand, relation (116) gives

[TABLE]

Combine (118)-(120) to get

[TABLE]

where $C(1)=\dot{e}(1)+2(s+1)e(1)+\frac{3\mathcal{E}(1)}{2(1+2s)}$ collects the constant terms. Therefore, for any $T\geq t>1$, we have

[TABLE]

by noting the non-increasing of Lyapunov function $\mathcal{E}(t)$. Rewrite (123) as $\mathcal{E}(T)\left(\ln(t)+\frac{3}{2(1+2s)}\right)\leq C(1)-\frac{\dot{e}(t)}{t}$, and then integrate it from $t=1$ to $t=T$ to derive

[TABLE]

Moreover, using the integration by parts and the positivity of functional $e(\cdot)$, we have

[TABLE]

By combining (124) and (125), we deduce that

[TABLE]

where $C_{1}=\frac{3}{2(1+2s)}-1$, $C_{2}=-C_{1}$ and $C_{3}=e(1)-C(1)$ are three constants.

Inequality (126) immediately yields $\mathcal{E}(\infty)\leq 0$. By the non-negativity of Lyapunov function $\mathcal{E}(\cdot)$, we conclude

[TABLE]

The continuity of $\chi(T)$ is obvious as our problem is linear. Hence, from (127) and the assumption of the lemma, we conclude that

[TABLE]

which implies the existence of the root of $\chi(T)$.

Appendix C: Proof of Proposition 4

Let $\{f^{n}\}_{n}\subset Q_{0}$ be bounded. Then there is a subsequence, denoted again by $\{f^{n}\}_{n}$, which converges weakly in $Q_{0}$ to some element $f^{*}\in Q_{0}$ because of the reflexivity of space $Q_{0}$. Let $u^{n}_{D}=u_{D}(f^{n})$, $u^{n}_{N}=u_{N}(f^{n})$, i.e., $u_{D}^{n}\in V_{0},u_{N}^{n}\in V$, and

[TABLE]

Then $\{u_{D}^{n}\}_{n}$ and $\{u_{N}^{n}\}_{n}$ are bounded in $V$ from the properties (84) and (85). Hence, we can extract two further subsequences, denoted again by $\{u_{D}^{n}\}_{n}$ and $\{u_{N}^{n}\}_{n}$, which converge weakly in $V$ and strongly in $Q$ to $u_{D}^{*}\in V_{0}$ and $u_{N}^{*}\in V$, respectively. Let $n\to\infty$ in (128) and (129) to get $u_{D}^{*}=u_{D}(f^{*})$ and $u_{N}^{*}=u_{N}(f^{*})$. Strong convergence of $\{u_{D}^{n}\}_{n}$ to $u_{D}^{*}$ in $V$ follows from

[TABLE]

as $n\to\infty$. Similarly, $u_{N}^{n}\to u^{*}_{N}$ as $n\to\infty$.

Denote $g^{n}=K\,f^{n}$. Then $\{g^{n}\}_{n}$ is bounded in $V$. Repeating the above argument, we conclude that there exists an element $s^{*}\in V$ such that

[TABLE]

Therefore, $\forall\ v\in V$,

[TABLE]

Thus, we have $g^{*}=Kf^{*}$. Consequently, strong convergence of $g^{n}$ to $g^{*}$ in $V$ follows from

[TABLE]

as $n\rightarrow\infty$, and the proof is completed.

Appendix D: Finite element discretization of boundary value problems

In this appendix, we discuss the numerical implementations of (82) and (83) by standard finite element method. We use linear finite element space for an approximation of the light source space $Q_{0}$. Specifically, let $\{\mathcal{T}_{0,H}\}_{H}$ be a regular family of triangulations over domains $\overline{\Omega}_{0}\subset\overline{\Omega}$ with meshsize $H>0$. For each triangulation $\mathcal{T}_{0,H}=\{K_{H}\}$, define finite element space $Q^{H}_{0}=\{q\in C(\overline{\Omega}_{0})\mid q|_{K_{H}}\in\mathcal{P}_{1}(K),\ \forall\,K_{H}\in\mathcal{T}_{0,H}\}$, where $\mathcal{P}_{k}$ represents the space of all polynomials of degree no greater than $k$. Let $\{\mathcal{T}_{h}\}_{h}$ be a regular family of triangulations over domains $\overline{\Omega}\subset\mathbb{R}^{d}$ with a mesh size $h>0$. For each triangulation $\mathcal{T}_{h}=\{K_{h}\}$, define finite element spaces $V^{h}$ and $V_{0}^{h}$ as follows.

[TABLE]

Moreover, we use the symbol $g^{\delta}_{1}+V_{0}^{h}$ for the set

[TABLE]

For each $f\in Q_{0}$, the finite element discretization of (82) and (83) read

[TABLE]

Similar to the continuous case, we use the symbols $u^{h}_{D}(f)$, $\widetilde{u}^{h}_{D}(g^{\delta}_{1})$, $u^{h}_{N}(f)$ and $\widetilde{u}^{h}_{N}(g^{\delta}_{2})$ for $u^{h}_{D}(f,0)$, $\widetilde{u}^{h}_{D}(0,g^{\delta}_{1})$, $u^{h}_{N}(f,0)$ and $u^{h}_{N}(0,g^{\delta}_{2})$, respectively.

Suppose that $\mathcal{T}_{0,H}$ and $\mathcal{T}_{h}$ are consistent, i.e., the triangulation $\mathcal{T}_{0,H}$ is a restriction of the triangulation $\mathcal{T}_{h}$ on $\overline{\Omega}_{0}$, and let $n_{0}$ and $n$ be the numbers of nodes of the triangulations $\mathcal{T}_{0,H}$ and $\mathcal{T}_{h}$. Denote $\varphi_{i}(x)\in V^{h}$, $1\leq i\leq n$, be the node basis functions of the finite element space $V^{h}$ associated with grid nodes $x_{i}\in\overline{\Omega}$. Let $x_{i_{j}}\in\overline{\Omega}_{0},1\leq j\leq n_{0}$ be the nodes of $\mathcal{T}_{0,H}$, and $\varphi_{i_{j}}(x)\in V^{h}$ the corresponding basis functions. Then, the approximate source function $f^{H}$ of $f$ can be expressed by $f^{H}=\sum^{n}_{j=1}f_{j}^{H}\varphi_{i_{j}}$ with $f_{j}^{H}=f(x_{i_{j}})$. For the problems (130) and (131), the solutions $u_{D}^{h}\in g^{\delta}_{1}+V_{0}^{h}$ and $u_{N}^{h}\in V^{h}$ can be expanded by $u_{D}^{h}=\sum^{n}_{i=1}u_{D,i}\varphi_{i}$ and $u_{N}^{h}=\sum^{n}_{i=1}u_{N,i}\varphi_{i}$, respectively, where $u_{D,i}=u_{D}^{h}(x_{i})$ and $u_{N,i}=u_{N}^{h}(x_{i})$.

Denote $I=\{1,2,\cdot\cdot\cdot,n\}$, $I_{0}=\{1,2,\cdot\cdot\cdot,n_{0}\}$, $I_{b}=\{i\in I|x_{i}\in\Gamma\}$, and define

[TABLE]

In the following, we use the same symbol for a finite element function and its vector representation associated with the given finite element basis functions. Then, the finite element solutions $u_{D}^{k}$ and $u_{N}^{k}$ of the forward problems (130) and (131) corresponding to the source $f^{k}$, can be calculated by

[TABLE]

Similarly, for the discretization of the quality $K^{*}\,(K\,f-y^{\delta})$, define

[TABLE]

Then the finite element approximation of $K^{*}\,(K\,f^{k}-y^{\delta})=w^{k}_{D}-w^{k}_{N}$ can be calculated through

[TABLE]

Note that (132) and (133) are the finite element discretization of the adjoint problems (98) and (99) with $u_{DN}$ being replaced by $u^{k}_{D}-u^{k}_{N}$.