Fast matrix exponential-based quasi-boundary value methods for inverse space-dependent source problems

1.   Introduction

Inverse source problems arise often in real-world applications, such as localizing unknown groundwater contaminant sources, geophysical prospecting, crack identification and pollutant detection. Let $T > 0$ and $\Omega\subset {\mathbb{R}}^d (d = 1, 2, 3)$ be an open and bounded domain with a piecewise smooth boundary $\partial\Omega$. We consider the inverse source problem (ISP) ^[6,10,37] of reconstructing the unknown space-dependent source term $f\in L^2(\Omega)$ from the final condition $g = u(\cdot, T) \in H_0^1(\Omega)$, according to a heat equation with homogeneous Dirichlet boundary conditions

where we assume zero initial condition $\phi\equiv 0$ for simplicity. The general case with nonzero initial condition can be treated similarly; see ^[20] for related discussion. In practice, the exact final condition $g$ is unknown and it is available as a noisy measurement $g_\delta\in L^2(\Omega)$ satisfying $\|g-g_\delta\|_{2}\le \delta$ for a noise level $\delta > 0$. This leads to an ill-posed inverse problem that requires regularization ^{[12,24,26,30]}.

There were many research works on the above inverse source problem with the source term $f$ being of various a priori form. For $f$ that depends on the state function $u$, the problem was investigated in ^[15,16]. For $f$ that is a function of both time and space variables but is additive or separable, we refer to ^{[32,39,40,53]}. For $f$ that depends on space or time variable only, many regularization methods have been developed, including Fourier method ^[11], quasi-reversibility method ^[10], quasi-boundary value method ^[49], simplified Tikhonov regularization method ^[48], the boundary element method ^[14], the method of fundamental solutions ^[1,46,47] and the finite element method ^[41]. Some iterative algorithms can be found in ^{[22,23,51,52]}. Regularization methods allowing efficient implementations are desirable in practical use. Besides the above mentioned ISPs for PDEs based on ordinary integer-order derivatives, there are also several recent works on solving ISPs in the framework of time-fractional PDEs, to name just a few ^{[3,8,17,21,35,42,43,44,50]}.

Being quite different from the standard Tikhonov regularization, the quasi-boundary value method ^[49] and its modified version ^[45] have been established as effective ways for regularizing such inverse source problems. They are shown to achieve an optimal order convergence rate under suitable regularity assumptions on the to-be-recovered source term. The main objective of this study is to develop optimal error estimates that cover in a unified way several source reconstruction methods in use, such as the simplified Tikhonov method, the quasi-boundary value method and its modified version. Furthermore, having as motivation the fact that the large-scale sparse linear systems resulting from the all-at-once space-time discretization of the obtained quasi-boundary value problem are expensive to solve, our second objective aims to present computationally more efficient matrix exponential based algorithms that eliminate the time variable and hence reduce the discretized system sizes and overall CPU times. In our recent work ^[20], we have designed efficient diagonalization-based parallel-in-time algorithms to solve such structured linear systems, which achieve a dramatic speedup in CPU times when compared with MATLAB built-in backslash solver.

The remaining of this paper is organized as follows. In the next Section 2, we give the explicit solutions in series form for the direct and inverse problems. The unified convergence analysis of the general QBVM regularization model is presented in Section 3. Two different implementations based on finite difference discretization are described in Section 4, where the proposed matrix exponential method was not studied in the literature of inverse problems. Several 1D and 2D numerical examples are reported in Section 5 to compare the discussed algorithms. Finally, Section 6 concludes the paper with some ideas for future work.

2.   Explicit solution in series form

In this section, we present the explicit solutions in series form corresponding to the direct and inverse problems, which lays the foundation for the subsequent convergence analysis.

2.1.   Direct forward problem

The forward problem consists in the following initial boundary-value problem:

with $f\in L^2(\Omega)$. We denote $< ., . >$ and $\|\cdot\|$ as the standard inner product and respective norm in $L^2(\Omega)$. Let $\{\lambda_k\}_{k\in N}$ be the eigenvalues of negative Laplacian operator $(-\Delta)$ with associated eigenfunctions $\{w_k\}_{k\in \mathbb{N}}$, $w_k\in H_0^1(\Omega)$, in the sense that

where $a:H_0^1(\Omega)\times H_0^1(\Omega)\rightarrow \mathbb{R}$ is given by $a(u, \eta) = \int_\Omega \nabla u\cdot \nabla \eta dx.$ It is well-known that $0 < \lambda_1\leq \lambda_2\leq \lambda_3\leq \ldots$, with $\lambda_k \rightarrow \infty$ as $k\rightarrow \infty$ and that the vector subspace generated by $\{w_k\}_{k\in \mathbb{N}}$ forms an orthonormal basis for $L^2(\Omega)$ [13,Theorem 1,pp. 335]. Classical Fourier method shows that the solution of the forward problem (2.1) in series form is given by

which converges in $L^2(0, T;H_0^1(\Omega))\cap C([0, T];L^2(\Omega))$ and satisfies Eq (2.1) in the generalized sense [31,pp. 128–131]. This explicit series solution Eq (2.3) will be used in convergence analysis.

2.2.   Inverse Source Problem

Given $g\in L^2(\Omega)$, we aim to recover the source term $f\in L^2(\Omega)$ using as input data final time measurements $u(., T)$ satisfying

To this end, we introduce the linear operator $B:L^2(\Omega)\rightarrow L^2(\Omega)$ defined by

Clearly, from Eq (2.3) for the exact source function $f$ we have $B(f) = u(., T)$. As $u\in C([0, T];L^2(\Omega))$, $B$ is well defined and it is a linear, compact, self-adjoint positive and one-to-one operator. Let $\{\sigma_k, w_k\}_{k = 1}^\infty$ be its singular system [26,Appendix A.6], where $\sigma_k = (1-e^{-\lambda_k T})\lambda^{-1}_k > 0$ is called singular value. In particular, we have $B(w_k(x)) = \sigma_k w_k(x)$ and $B^* = B$.

From here on we will focus on methods to recover $f$ from the final time data given in Eq (2.4). We start by noting from Eq (2.3) and the additional condition (2.4) that for $g\in L^2(\Omega)$ we have

where $g_k = < g, w_k >$. Thus, recovering the source $f\in L^2(\Omega)$ from additional data $g\in L^2(\Omega)$ amounts to solve the abstract operator equation

Obviously, if $g\in \mathcal{R}(B)$, Eq (2.7) is uniquely solvable as $q\in{\cal R}(B)$ and $B$ is a positive operator. Moreover, from Eq (2.6) we can easily obtain that the Fourier coefficient $f_k$ of the exact source function $f$ is given by

so that

i.e., a closed form to the source recovering problem can be expressed as

where $g_k = < g, w_k >$ and $\sigma_k = (1-e^{-\lambda_k T})\lambda^{-1}_k$. It is worth noting that as $B$ is a compact operator, its inverse $B^{-1}:\mathcal{R}(B)\rightarrow L^2(\Omega)$ is unbounded so that the linear equation (2.7) is ill-posed. This means that small errors in the input data $g$ can result in arbitrarily large perturbations in the computed solution via the series solution Eq (2.8). Usually, the available final time measurement is contaminated with noise. Then, even if the exact data $g$ belongs to $\mathcal{R}(B)$, we cannot expect the same for the noisy data $g_\delta$ and, consequently, the series solution Eq (2.8) with noisy data $g_\delta$ instead of $g$ will be divergent or inaccurate. In the following section, we will introduce effective regularization techniques to address the ill-posedness in practical computation, and their implementation will be explained in Section 4.

3.   Convergence analysis of QBVM by filtering regularization

In this section we will take advantage of the singular system of $B$ to build stable approximate solutions to the source recovery problem. As usual, let us assume that $g_\delta\in L^2(\Omega)$ such that

for some $\delta > 0$. Hereafter $\|\cdot\|$ denotes the standard $L^2(\Omega)$ norm. With assumed zero initial condition, define the recovered source term corresponding to the given noisy final data $g_\delta$:

This expansion not only illustrates the influence of the errors in $g$ but also suggests trying to construct stable approximations for the source function $f$ by damping or filtering out the factors $1/\sigma_k$. This can be achieved by introducing approximate solutions given by

parameterized by a positive regularization parameter $\alpha$, where ${\boldsymbol{\digamma}}_{\alpha}: (0, \|B\|]\rightarrow \mathbb{R}$ is a bounded function referred to as filter function or filter factor, such that

Here $\|B\|$ denotes the operator norm of $B$ with respect to $L^2(\Omega)$ norm, that is $\|B\| = \sup_{\|v\| = 1} \|B(v)\|.$ It is known from regularization theory ^[26] that if ${\boldsymbol{\digamma}}_\alpha$ satisfies i)–iii) in Eq (3.3) and $\alpha = \alpha(\delta)$ is chosen by a priori or a posteriori choice rules such that $\alpha\rightarrow 0$ as $\delta\rightarrow 0$, then $f_\alpha^\delta \rightarrow f$ as $\delta\rightarrow 0$.

For a general operator equation $Af = g_\delta$ with a compact operator $A$, the Tikhonov regularization method computes an approximate solution $f_\alpha^\delta$ by minimizing the regularized functional

where $\alpha > 0$ is the regularization parameter. A straightforward calculation shows that the unique minimizer of $J_\alpha$ is given by

where $A^*$ denotes the adjoint of $A$. Let $\{\sigma_k, \xi_k, \eta_k\}_{k = 1}^\infty$ be the singular system ^[26] of $A$ satisfying $A \xi_k = \sigma_k \eta_k$ and $A^*\eta_k = \sigma_k \xi_k$, one has

Alternatively, in the case of positive and self-adjoint operator $B$ with $B^* = B$, a simpler approximation solution $f_\alpha$ can be obtained by minimizing the regularized quadratic functional

which is equivalent to solving a simpler regularized equation

With a positive self-adjoint $B$, the singular system of $B$ reduces to $\{\sigma_k, w_k\}_{k = 1}^\infty$ with $w_k = \xi_k = \eta_k$, we have

This method of constructing regularized solution as Eq (3.7) is referred to as Lavrentiev regularization ^[28,34,38] or, simplified regularization. Notice that the Tikhonov regularization applied to Eq (2.7) with a self-adjoint operator $B$ results in a more complicated approximation $f_\alpha^\delta = (B^2 +\alpha I)^{-1}Bg_\delta$, which explains why it is also referred to as simplified regularization.

Both Eqs (3.5) and (3.7) are particular cases of filtering regularization methods of the form

with filter factors ${\boldsymbol{\digamma}}_\alpha(\sigma_k)$. For Tikhonov and Lavrentiev regularization the filter factors ${\boldsymbol{\digamma}}^{\rm Tik}_{\alpha}(\sigma)$ and ${\boldsymbol{\digamma}}^{\rm Lav}_{\alpha}(\sigma)$ are given by

respectively. It is seen that for singular values $\sigma_k$ much smaller than $\sqrt{\alpha}$ the filter factors ${\boldsymbol{\digamma}}_\alpha^{\rm Tik}(\sigma_k)$ are small, and thus the corresponding components in Eq (3.5) are damped or filtered. In this case, the filter factors are approximately proportional to $\sigma_k^2$ and so we can use $\alpha$ to control the filtering of potentially increasing ratios ${ < g_\delta, w_k > }/{\sigma_k}$. That is the amount of regularization or filtering depends on a judicious choice of the regularization parameter $\alpha.$

We consider the generalized QBVM approximations $f_{\alpha, \beta}^\delta$ to the source function $f$, within the framework of filtering methods, with $f_{\alpha, \beta}^\delta$ being defined as the solution of the regularized model

which reduces to the QBVM ^[49] if $\beta = 0$ and the MQBVM ^[45] if $\alpha = 0$. Following the discussion in ^[20], the solution of Eq (3.9) in series form can be expressed as

and that $f_{\alpha, \beta}^\delta \rightarrow f$ as $\alpha, \beta \rightarrow 0$, with the exact source function $f$ being written as

thus indicating that $\alpha$ and $\beta$ play the role of regularization parameters. We will study the convergence properties of $f_{\alpha, \beta}^\delta$ within the framework of filtering regularization methods. In particular, we will investigate the error estimates $\|f-f_{\alpha, \beta}\|$ for suitable choices of the regularization parameters and focus on efficient numerical implementation in the practical case of using discrete data. For clarity, we split our discussion into 3 different cases: a) $\beta = 0$; b) $\alpha = 0$; and c) $\beta\ne 0$, $\alpha\ne0$.

Case a) Since $\sigma_k = (1-e^{-\lambda_kT})/\lambda_k$, when $\beta = 0$ we obtain the regularized approximation

and the approximate solution is nothing but the solution obtained by the well-studied Lavrentiev regularization method. This interesting equivalence between QBVM and Lavrentiev regularization seems to be not known in literature. In particular, provided the regularization parameter is chosen either a priori or a posteriori by proper choice rules, it is known that $\|f-f_{\alpha, 0}\|\rightarrow 0$ as $\alpha \rightarrow 0$. Moreover, the best convergence rate of $f_{\alpha, 0}$ is $O(\delta^{1/2})$ ^[33] and this rate cannot improved for a compact $B$ with non-closed range.

Case b) With $\alpha = 0$ in Eq (3.10), we can get the approximate solution

with filter factors given by

Since

it is clear that conditions i)–iii) in Eq (3.3) hold with $c(\beta) = C_1/\sqrt{\beta}$ and $C_1 = 1/2\sqrt{1-e^{-\lambda_1T}}$.

Estimates on the error norm $\|e_\beta^\delta\| : = \|f-f_{0, \beta}^\delta \|$ can be obtained based on both the filtering properties of filter factors and a priori assumptions on the exact solution $f$. Recall that for a general operator equations $Kf = g_\delta$ with a compact $K$ having non-closed range ${\cal R}(K)$ and the Moore–Penrose inverse $K^\dagger$, under the assumption that the exact solution satisfies the so-called source condition

it is known that ^[9,12]

where the index $\mu_0$ denotes the qualification of the regularization method ^[9,12]. For our following analysis let us introduce $r_\beta(\sigma_k) : = 1-{\boldsymbol{\digamma}}_\beta(\sigma_k)$, the regularized solution for exact data,

and then decompose the error $e_\beta^\delta = f-f_{0, \beta}^\delta$ as the sum of regularization and noise errors,

where

and

In our context, the above source condition reads $f\in {\cal R}(B^{2\mu})$ and is equivalent to

with some $z\in L^2(\Omega)$ such that $B^{2\mu}z = f$. Then the regularization error norm satisfies

But since

and hence

to bound the regularization error norm we have to bound the function

It is easy to see that this function attains its maximum at $\sigma = \sqrt{\mu \beta/(1-\mu)}$ for $0 < \mu < 1$ and therefore we have

with $C = \mu^\mu(1-\mu)^{1-\mu}.$ For $\mu \ge 1$, $h_\mu(\sigma)$ is strictly increasing and attains its maximum at $\sigma = \sigma_1$. Thus for $\mu\ge 1,$ $h_\mu(\sigma)\le\sigma_1^{(2\mu-2)}\beta.$ Consequently, for proper $C_2$ the regularization error norm can be bounded as

Further, for the noise error norm, since

using Eq (3.14) we obtain

Now Eqs (3.23) and (3.24) together imply

and the a priori selection parameter rule $\beta = \left (\dfrac{\delta}{\|z\|}\right)^\frac{2}{2\mu+1}$ delivers the estimate

The error estimate we just derived is essentially the same estimate obtained by applying Tikhonov regularization to an operator equation that involves a compact operator. The reason behind this is that our estimate depends on Eq (3.14) and Eqs (3.21)–(3.24) which in turn depend heavily on the properties of the Tikhonov filter. Based on this observation, an error estimate with the regularization parameter $\beta$ chosen by Morozov discrepancy principle can also be deduced. This essentially leads to an a posteriori regularization parameter choice rule of choosing $\beta$ such that

for some given $\rho > 1$. Based on Eq (3.14) and Eqs (3.21)–(3.24) an estimate of the form Eq (3.25) can be derived which holds for $0 < \mu\le 1/2,$ see [12,Thm. 4.17] or [9,Appendix C] for an illustrative analysis in finite dimension. Error estimates obtained by other means can also be found in ^[45].

Case c) If $\alpha\ne 0$ and $\beta\ne 0$, the approximation can be described in terms of filter factors as

with

which describes a two-parameter regularization problem. This is a difficult problem that will not be fully addressed in this work. Instead, we prefer to only discuss the case where $\alpha = \alpha(\beta)$, e.g., $\alpha = 2\sqrt{\beta}$. For this, as before, we decompose

and then estimate regularization and noise errors separately. For both choices of $\alpha$ we see that

and therefore

From here on we will focus on the choice $\alpha = 2\sqrt{\beta}$. Following in the same lines as in Case b let us introduce $r_{2\sqrt{\beta}, \beta}(\sigma_k) = 1- {\boldsymbol{\digamma}}_{2\sqrt{\beta}, \beta}(\sigma_k)$ and then note that since in this case

and since

we need to analyze the function

In fact, elementary calculations show that the derivative of this function is

and that its critical points are

This shows that the only positive critical point occurs when $0 < \mu < 1/2$ and that it is

Since the quadratic function between brackets in Eq (3.30) opens downward and the derivative changes sign from positive on $[0, s^*]$ to negative for $s > s^*,$ we conclude that $g_\mu$ attains its maximum at $s = s^*.$ Further, similar to the analysis of $h_\mu$, $g_\mu$ increases for $\mu\ge 1/2$ and $g_\mu$ attains its largest values at $s = \sigma_1$. Hence, for proper $C_4$ we have $g_{\mu}(\sigma)\le C_4\beta^{\mu}, \; 0 < \mu\le 1/2, \; 0\le \sigma\le \sigma_1.$ Using this in Eq (3.29) it is readily seen that the regularization error norm can be bounded as

Then, Eqs (3.28) and (3.31) imply that

and an a priori choice rule of $\beta$ as the one in Case b leads to

To summarize, we have obtained the following convergence results.

Theorem 3.1. Let $f_{\alpha, \beta}^\delta$ be the general QBVM approximation as given by Eq (3.9) and $f$ be the exact source function. We have the following estimates:

(a) With $\beta = 0$, if $f\in {\cal R}(B^{2\mu})$ and choosing $\alpha = \left (\dfrac{\delta}{\|z\|}\right)^{\frac{1}{2\mu+1}}$, there holds [33,Corollary 3.3]

(b) With $\alpha = 0$, $f\in {\cal R}(B^{2\mu})$ and choosing $\beta = \left (\dfrac{\delta}{\|z\|}\right)^{\frac{2}{2\mu+1}}$, there holds

(c) With $\alpha = 2\sqrt{\beta} > 0$, if $f\in {\cal R}(B^{2\mu})$ and choosing $\beta = \left (\dfrac{\delta}{\|z\|}\right)^{\frac{2}{2\mu+1}}$, there holds

Clearly, Case b gives a faster convergence rate than Cases a and c whenever $\mu > 1/2$.

4.   Two different finite difference implementations

We will use a second-order center finite difference scheme in space and a first-order backward Euler scheme in time for full discretization of the continuous QBVM model. More specifically, let $h > 0$ be the uniform spatial step size and $\tau = T/n$ be the uniform time step size. Let $\Delta_h\in {\mathbb{R}}^{m\times m}$ denotes the discrete Laplacian matrix and $I_h\in {\mathbb{R}}^{m\times m}$ be an identity matrix. Here we will describe the numerical scheme with a general initial condition $\phi$ for the purpose of better practical use.

4.1.   An all-at-once scheme for QBVM

Let $\phi_h$ and $g_{\delta, h}$ denotes function values of $\phi$ and $g_{\delta}$ over all spatial grids in lexicographical order, respectively. Let $f_h$ and $u^j$ denotes the finite difference approximation of $f$ and $u(\cdot, j \tau)$ over all spatial grids with the initial condition given by $u^0 = u(\cdot, 0) = \phi_h$. The full discretization of Eq (3.9) reads

which can be reformulated into a large-scale nonsymmetric sparse linear system

where

For large $m$ and $n$, the above all-at-once sparse linear system in Eq (4.2) is very costly to solve by direct solvers, including MATLAB's build-in highly optimized backslash ('$\backslash$') sparse direct solver. In our recent work ^[20], a novel diagonalization-based parallel-in-time algorithm was proposed to speed up the direct inversion of $S$, where the special choice of $\alpha$ and $\beta$ are crucial to the stability of diagonalization. Notice that in the original inverse source problem we are only interested in recovering $f$ from the final measurement $u(\cdot, T)$ and the actual values of $u^j, j = 1, 2, \ldots, n$ are unnecessary to compute and store. This drawback of the all-at-once scheme for QBVM can be addressed by the following matrix exponential-based implementation with a much better computational efficiency, where the intermediate values $u^j, j = 1, 2, \ldots, n$ are not computed and stored anymore. It is also possible to incorporate other finite difference schemes in time.

4.2.   Matrix exponential-based implementation for QBVM

Different from the above all-at-once full discretization in Eq (4.1), the semi-discretization in space of Eq (3.9) together with the initial condition can be written as linear ODE system with constant coefficient matrix

where $A_h = -\Delta_h$ denotes the negative discretized Laplacian after enforcing Dirichlet boundary conditions, and since ${f_h}$ is independent of time $t$. With general initial condition $\phi$, the explicit solution of Eq (4.3) is

With the noisy final measurement ${{{\boldsymbol{v}}}}(T) = g_{\delta, h}$, it follows from Eq (4.4) that (compare to Eq (2.7))

which can be solved for $f_h$. In a similar manner, the above matrix exponential form of QBVM regularization model (3.9) leads to

which gives a Lavrentiev-type regularized system (avoided computing the dense matrix $A^{-1}$)

Here the suitable choice of regularization parameters $(\alpha, \beta)$ depends on the noise level $\delta$. Notice that the matrix exponential $e^{-A_h T}$ is in general very expensive (with $O(m^3)$ complexity) to compute exactly, but the matrix exponential-vector product $e^{-A_hT}f_h$ with a sparse matrix $A_h$ can be approximately computed very efficiently ^[2]. This suggests us to solve the regularized linear system (4.7) by iterative Krylov subspace methods (such as PCG) that only makes use of matrix-vector products. The above matrix exponential formulation is computationally attractive (especially for 2D/3D problems), since it eliminates the time variable as in the all-at-once scheme giving large-scale sparse linear systems for QBVM. The same idea of matrix exponential operator can also be generalized to nonlinear cases with exponential integrators ^{[18,19,27,29]}.

Alternatively, if the eigenpairs of $A_h = -\Delta_h$ are simple to construct or compute analytically, then one can avoid the direct computation of the matrix exponential $e^{-A_h T}$. More specifically, using an eigendecomposition of $A_h = -\Delta_h$, i.e., $A_h = P\Lambda P^T,$ with eigenvalues $\hat\lambda_k$, $0 < \hat\lambda_1 < \hat\lambda_2 < \cdots < \hat\lambda_m$, and associated orthonormal eigenvectors $p_k$ from the $k$-th column of $P$, the solution $f_h$ of Eq (4.5) reads

which matches with the truncated exact series solution given in Eq (2.8). Analogously, in the case of Eq (4.7), using eigenpairs of $A_h$, the regularized solution reads

which matches with the truncated exact series solution given in Eq (3.10). Nevertheless, these explicit solutions in Eq (4.9) are of more theoretical use, since the eigenpairs of $A_h$ are in general difficult to obtain, except for special cases with uniform grids and simple boundary conditions.

Compared with the all-at-once sparse linear system (4.2) for the above QBVM, the size of the system matrix $L_{\alpha, \beta}$ is much smaller and hence computationally cheaper to solve. But it involves the matrix exponential $e^{-AT}$ that requires extra costs for its computation. For 1D examples, we explicitly construct the matrix exponential $e^{-AT}$ by the MATLAB function expm, which is very efficient for matrix $A$ of a small size. Nevertheless, for 2D example it is much more efficient to use PCG as the iterative system solver, where the matrix exponential times a vector is approximated by the MATLAB function expmv ^[2] The MATLAB codes for implementing the above methods are available online at the GitHub link: https://github.com/junliu2050/MatExpQBVM.

5.   Numerical examples

In this section, we present some numerical examples to illustrate the computational efficiency of our proposed methods. All simulations are implemented in serial with MATLAB on a Dell Precision 5820 Workstation with Intel(R) Core(TM) i9-10900X CPU@3.70GHz CPU and 64GB RAM, where CPU times (in seconds) are estimated by the timing functions tic/toc. For solving the sparse linear systems from all-at-once scheme for QBVM, we use MATLAB's backslash sparse direct solver, which runs very fast for several thousands (but not millions) of unknowns. To avoid inverse crimes, given an exact source $f$ we solve the forward (direct) problem by the Crank-Nicolson (different from backward Euler) time-stepping scheme in time to compute $g_h$, and then generate the noisy final condition measurement by adding random noise according to

where $\epsilon > 0$ controls the noise level and ${\rm{rand}}(-1, 1)$ denotes random noise uniformly distributed within $[-1, 1]$. We then compute the estimated noise bound $\delta: = \|g_{\delta, h}-g_h\|_{2, h}$ in discrete $L_2$ norm.

Since $\|z\|$ and $\mu$ in Theorem 3.1 are usually unknown, we select more practical regularization parameters as follows. For the first Case (a) with $\beta = 0$, inspired by the QBVM ^[49], the choice of the regularization parameter $\alpha = \sqrt{\delta}$ seems to work well. For the second Case b with $\alpha = 0$, inspired by the MQBVM ^[45], the choice of the regularization parameter $\beta = \delta$ seems to work well. For the third Case c with $\alpha > 0$ and $\beta > 0$, based on the above discussion we can choose $\alpha = \sqrt{\delta}$ and $\beta = \alpha^2/4 = \delta/4$ to get a similar convergence rate as observed in Case a. Hence, the Case c is of less practical use since it gives slower convergence rate than the Case b with similar computational costs. Therefore, for brevity we choose to not compare the Case c in our examples. As studied in ^[43], it is more practical to use an a posteriori regularization parameter choice rule Eq (3.26) that does not require the knowledge of $\|z\|$ and $\mu$. However, it can be more expensive since it needs to iteratively solve the nonlinear equation (3.26) for determining the regularization parameter, which also necessitates further separate convergence analysis. In particular, the QBVM based on all-at-once scheme seems to be too costly to solve multiple times.

For solving $L_{\alpha, \beta} f_h = {{\boldsymbol{r}}}_h$ in matrix exponential-based implementation, we also use MATLAB's backslash sparse direct solver for 1D examples since $e^{-AT}$ can be explicitly computed very fast. For 2D examples, however, we will use preconditioned conjugate gradient (PCG) iterative solver (with a stopping tolerance $tol = 10^{-3}$) since it only requires to compute or approximate the matrix exponential times a vector $e^{-AT}v$ without constructing $e^{-AT}$, which can be performed very efficiently by using just matrix-vector product with the sparse matrix $A$; see ^[2] for more technical details. To the best of our knowledge, the application of existing fast matrix exponential based algorithms to our considered inverse source problems was not discussed in the literature; see ^[25] for a fast structured preconditioner when apply GMRES method to solve the system like Eq (4.2).

Example 1. Choose $\Omega = (0, \pi), T = 1$, $\phi(x) = 0$, and a smooth source function

Figure 1 compares the reconstructed source function by 4 different regularization methods: QBVM and MQBVM based on all-at-once scheme, MatExp-a (Case a) and MatExp-b (Case b) based on matrix exponential implementation. Clearly, MatExp-a and MatExp-b shows very similar approximation accuracy (measured as 'error' in $L^2$ norm corresponding to the smallest $\delta$) as QBVM and MQBVM, respectively, but costs much less CPU times (shown in titles) due to the elimination of time variable by using matrix exponential. For example, with $h = \pi/1024$, MQBVM costs about 20 seconds while MatExp-b takes only 0.05 second, giving about 40 times speedup. Notice that both MQBVM and MatExp-b provide about two times more accurate reconstruction in this case with a smooth source term, which matches with our convergence rates estimate in Theorem 3.1.

Example 2. Choose $\Omega = (0, \pi), T = 1$, $\phi(x) = 0$, and a non-smooth source function

Figure 2 compares the reconstructed source function by 4 different regularization methods: QBVM, MQBVM, MatExp-a (Case a), and MatExp-b (Case b). Again, both MatExp-a and MatExp-b are significantly faster than QBVM and MQBVM, while both MQBVM and MatExp-b deliver smoother recovery with a better accuracy due to the introduced Laplacian regularization term $\beta\Delta f$ that penalizes undesirable oscillation as observed in QBVM and MatExp-a.

Example 3. Choose $\Omega = (0, \pi), T = 1$, $\phi(x) = 0$, and a discontinuous source function

Figure 3 compares the reconstructed source function by 4 different regularization methods: QBVM, MQBVM, MatExp-a (Case a), and MatExp-b (Case b). Similarly, both MatExp-a and MatExp-b are significantly faster than QBVM and MQBVM, while both MQBVM and MatExp-b deliver smooth source terms that largely fit the discontinuous pattern of $f$. Nevertheless, in this case both MQBVM and MatExp-b achieve a comparable accuracy, which is expected from our convergence analysis since a discontinuous $f$ has a much lower regularity. To accurately capture the discontinuous jumps, it requires different regularization techniques which will be studied in our future work.

Example 4. Choose $\Omega = (0, \pi)^2, T = 1$, $\phi(x, y) = 0$, and a smooth source function

Figure 4 compares the reconstructed source function by 4 different regularization methods: QBVM, MQBVM, MatExp-a (Case a), and MatExp-b (Case b). Similarly, MatExp-a and MatExp-b shows a comparable approximation accuracy as QBVM and MQBVM, respectively, but costs much less CPU times (e.g., reduced by over 30 times from about 120 seconds by MQBVM to 4 seconds by MatExp-b). The speedup of CPU times will become more significant for a fine mesh.

6.   Conclusion

The classical quasi-boundary value method (QBVM) and its variants are widely used for regularizing inverse space-dependent source problems, which upon full space-time finite difference discretization lead to large-scale ill-conditioned nonsymmetric sparse linear systems that are costly to solve. In this paper we first investigate the convergence rates of the general QBVM regularization model and then propose to integrate matrix exponential algorithms to eliminate the time variable so that the corresponding discretized linear systems are of smaller sizes and hence cheaper to solve. Both 1D and 2D examples show our proposed matrix exponential based algorithms can achieve a comparable accuracy with significantly faster CPU times. Compared with our recent work ^[20] utilizing parallel-in-time algorithms, our proposed matrix exponential based methods also have the advantage of using less memory storage.

The QBVM regularization approaches can not accurately capture the corners or jumps in less regular source term, which requires advanced regularization techniques, such as the widely used nonlinear total variation-based regularization ^[7,36]. The above used choice of $\alpha = \sqrt{\delta}$ or $\beta = \delta$ depends on the noise level $\delta$, which may not be available in practice. Hence, it is interesting to generalize the improved maximum product criterion (IMPC) techniques ^[4,5] to estimate the effective regularization parameters $\alpha$ and $\beta$ without the exact knowledge of noise level.

Acknowledgments

The authors appreciate the editor and two anonymous referees for their revision suggestions that have significantly improved the quality of this paper.

Conflict of interest

The authors declare there is no conflict of interest.