# A data-driven model reduction method for parabolic inverse source problems
and its convergence analysis

Zhongjian Wang  Department of Statistics and CCAM, The University of Chicago,
Chicago, IL 60637, USA. (zhongjian@statistics.uchicago.edu).  Wenlong Zhang
Corresponding author. Department of Mathematics, Southern University of
Science and Technology (SUSTech), 1088 Xueyuan Boulevard, University Town of
Shenzhen, Xili, Nanshan, Shenzhen, Guangdong Province, P.R.China.
(zhangwl@sustech.edu.cn).  Zhiwen Zhang Corresponding author. Department of
Mathematics, The University of Hong Kong, Pokfulam Road, Hong Kong SAR,
P.R.China. (zhangzw@hku.hk).

###### Abstract

In this paper, we propose a data-driven model reduction method to solve
parabolic inverse source problems efficiently. Our method consists of offline
and online stages. In the off-line stage, we explore the low-dimensional
structures in the solution space of the parabolic partial differential
equations (PDEs) in the forward problem with a given class of source functions
and construct a small number of proper orthogonal decomposition (POD) basis
functions to achieve significant dimension reduction. Equipped with the POD
basis functions, we can solve the forward problem extremely fast in the online
stage. Thus, we develop a fast algorithm to solve the optimization problem in
the parabolic inverse source problems, which is referred to as the POD
algorithm in this paper. Under a weak regularity assumption on the solution of
the parabolic PDEs, we prove the convergence of the POD algorithm in solving
the forward parabolic PDEs. In addition, we obtain the error estimate of the
POD algorithm for parabolic inverse source problems. Finally, we present
numerical examples to demonstrate the accuracy and efficiency of the proposed
method. Our numerical results show that the POD algorithm provides
considerable computational savings over the finite element method.

AMS subject classification: 35R30, 65J20, 65M12, 65N21, 78M34.

Keywords: Parabolic inverse source problems; regularization method; data-
driven model reduction; proper orthogonal decomposition (POD); stochastic
error estimate; optimal regularization parameter.

##  1 Introduction

This paper presents a data-driven model reduction method to solve parabolic
inverse source problems and studies the convergence analysis of the proposed
method. Inverse problems are very important in physics, engineering, and
bioengineering. The inverse source problems, which seeks reconstruction of
source from final time observation, have attached much attention of the
researchers over the past decades, see an introduction and references in [21].
They have been widely studied in the literature, applied to many physical and
engineering source identification problems, e.g. migration of groundwater,
groundwater pollution detection, control of pollution source and environmental
protection [8, 12, 15, 11, 20] and references therein. The accurate recovery
pollutant source is crucial to environmental safeguard in cities with high
populations [9]. The estimation of the strength of acoustic sources from
measurements can be found in e.g. [20, 13, 27]. The inverse source problems
that arise from PDEs are ill-posed in the sense of Isakov and Hadamard [20,
17], since the eigenvalues of the elliptic operators decay exponentially fast,
especially the lack of stability with respect to the uncertainty in the
measurement data is the most difficult challenge for numerical inversions.
Namely, a small change in the data may lead to a significant difference in the
reconstructed source strength. Due to the important applications of inverse
source problems, numerical methods have been widely explored [8, 7, 26, 14,
15] and references can be found therein.

In this paper, we will work on a very practically physical scenario, assuming
that the observational measurement is collected point-wisely over a set of
distributed sensors located at $\\{x_{i}\\}^{n}_{i=1}$ over the physical
domain [11, 20, 25, 28]. At each sensor, independently additional noise or
random error will be considered due to the uncertainty of natural noise,
measurement errors, and the other uncertainty of the model itself. We will
take one of the realistic approaches for such inverse problems, i.e.
optimizing the mean-square error with a proper Tikhonov type regularizations
[7, 15, 35]. Classical methods such as regression methods [15], linear and
nonlinear programming methods [15], linear and nonlinear conjugate gradient
methods [1, 35], Newton-type methods, etc. can be applied during optimization.
For the inverse parabolic source problem, one usually uses iteration
optimization methods to find the true source term [7, 10, 22]. In each
iteration, one has to solve the forward parabolic equation one or two times.
However, as the sizes of discrete problems grow (e.g. finite element method
(FEM) or finite difference method (FDM)), the computation time will grow
rapidly, especially for such evolution problems. As a consequence, the
computation of the forward equation will cost the most time throughout the
whole procedure. This motivates us to develop efficient model reduction
methods to address this issue.

One of the successful model reduction ideas in solving evolution problems is
the proper orthogonal decomposition (POD) method [31, 5]. The POD method uses
the data from an accurate numerical simulation and extracts the most energetic
modes in the system by using the singular value decomposition. This approach
generates low-dimensional structures that can approximate the solutions of the
evolution problems with high accuracy. The POD method has been used to solve
many types of PDEs, including linear parabolic equations [34], Navier‐Stokes
equations [23], viscous G-equations [16], Hamilton–Jacobi–Bellman (HJB)
equations [24], and optimal control problem [2]. The interested reader is
referred to [29, 3, 18] for a comprehensive introduction of the model
reduction methods.

Since the POD method can significantly accelerate the computation of the
forward problem compared with the traditional methods, e.g. FEM and FDM, we
apply the POD method to solve parabolic inverse source problems with scattered
measured data at the final time. The key idea is to exploit (model-based) and
construct (data-based) the intrinsic approximate low-dimensional structures of
the underlying problem that consists of two components. First, we have a
training component that computes a set of data-driven basis functions (i.e.
the POD basis functions) to achieve significant dimension reduction in the
solution space. Following up, there is a fast-solving component that solves
the optimization problem using the constructed POD basis functions in each
iteration. Hence we achieve an effective data-driven model reduction method in
solving the parabolic inverse source problems and overcome the typical
computational bottleneck in the traditional methods, e.g. FEM and FDM.

Then, we study the convergence analysis of the proposed method. In [7], a
general theory of the stochastic convergence of numerical method has been
established for a certain type of inverse problem. Based on this framework, we
explore the convergence theory in the POD settings. Specifically, we will
prove a relatively weaker convergence of the POD method when the source term
only belongs to $L^{2}(\Omega)$, where $\Omega$ is the computational domain.
Unlike the traditional analysis of the POD method [23] or FEM convergence, we
do not assume the higher regularity for parabolic PDE solution $u$, i.e.
$u_{tt}$ to be bounded in $L^{2}(\Omega)$, which is quite strict in many
cases. Based on our analysis, we derive the stochastic convergence when
applying the POD method to the parabolic inverse source problem with uncertain
data. Our analysis in Theorem 4.3 shows that the optimal error of the Tikhonov
type least-square optimization problem depends on the noise level, the number
of sensors, and without any higher source regularization but the
$L^{2}(\Omega)$ norm of the source term. A self-consistent algorithm to get
the optimal smoothing parameter is also established using the POD method to
solve the inverse parabolic source problem, motivated by the recent study in
[6, 7]. Finally, we conduct numerical experiments to demonstrate the accuracy
and efficiency of the proposed method. Several kinds of source functions are
involved. Comparing to original optimization with FEM basis functions, the
proposed POD method significantly reduces the computational cost and maintains
the same accuracy as the FEM.

The rest of the paper is organized as follows. In Section 2, we introduce the
setting of the parabolic inverse source problems. In Section 3, we develop the
data-driven model reduction method for solving parabolic inverse source
problems. In Section 4, we provide the error estimate for the proposed method.
In Section 5, we present numerical results to demonstrate the accuracy of our
method. Finally, concluding remarks are made in Section 6.

##  2 Parabolic inverse source problems

To start with, we consider a parabolic equation,

$\left\\{\begin{aligned} u_{t}+\mathcal{L}u&=f(x)&\mbox{in
}\Omega\times(0,T),\\\ u(x,t)&=0&\mbox{on }\partial\Omega\times(0,T),\\\
u(x,0)&=g(x)&\mbox{in }\Omega\,,\end{aligned}\right.$ (2.1)

where $\Omega\subset\mathbb{R}^{d}$ $(d=1,2,3)$ is a bounded domain with
$C^{2}$ boundary or a convex domain satisfying the uniform cone condition,
$\mathcal{L}$ is a second-order elliptic operator given by
$\mathcal{L}u=-\nabla\cdot(a(x)\nabla u)+c(x)u$, and $g(x)$ is the initial
value. We assume the elliptic operator $\mathcal{L}$ is uniform elliptic, i.e.
there exist $a_{\min},a_{\max}>0$, such that $a_{\min}<a(x)<a_{\max}$ for all
$x\in\Omega$. Moreover, we assume $a(x)\in C^{1}(\bar{\Omega})$, $c(x)\in
C(\bar{\Omega})$ and $c(x)\geq 0$.

Let $u$ be the solution of the parabolic equation (2.1). We define the forward
operator $\mathcal{S}:$ $L^{2}(\Omega)\rightarrow H^{2}(\Omega)$ by
$\mathcal{S}f=u(\cdot,T)$. The forward problem is to compute the solution
$u(\cdot,t)$ for $t>0$ with known source term $f(x)$ and known initial
condition $g(x)$.

In the parabolic inverse source problem, $f(x)\in L^{2}(\Omega)$ is an unknown
source term that we want to reconstruct based on the final time measurement
$u(\cdot,T)$. To be specific, we will focus on a very practically physical
scenario, where we assume that the observational measurement is collected
point-wisely over a set of distributed sensors located at
$\\{x_{i}\\}^{n}_{i=1}$ over the physical domain $\Omega$ (see e.g.[11, 20,
25, 28]). We also assume the measurement data is always blurred by noise and
takes the form $m_{i}=\mathcal{S}f^{*}(x_{i})+e_{i}$, $i=1,\cdots,n$, where
$f^{*}\in L^{2}(\Omega)$ is the true source term and $\\{e_{i}\\}^{n}_{i=1}$
are independent and identically distributed (i.i.d.) random variables on a
proper probability space ($\mathfrak{X},\mathcal{F},\mathbb{P})$. From [7] and
the analysis therein, we know that $\|u\|_{C([0,T];H^{2}(\Omega))}\leq
C\|f\|_{L^{2}(\Omega)}$. According to the embedding theorem of Sobolev spaces,
we know that $H^{2}(\Omega)$ is continuously embedded into $C(\bar{\Omega})$
so that $\mathcal{S}f^{*}(x)$ is well defined point-wisely for all
$x\in\Omega$. Without loss of generality, we assume that the scattered
locations $\\{x_{i}\\}_{i=1}^{n}$ are uniformly distributed in $\Omega$, i.e.,
there exists a constant $B>0$ such that ${d_{\max}}/{d_{\min}}\leq B$, where
${d_{\max}}$ and ${d_{\min}}$ are defined by

$\displaystyle d_{\max}=\mathop{\rm sup}\limits_{x\in\Omega}\mathop{\rm
inf}\limits_{1\leq i\leq
n}|x-x_{i}|~{}~{}~{}\mbox{and}~{}~{}~{}d_{\min}=\mathop{\rm inf}\limits_{1\leq
i\neq j\leq n}|x_{i}-x_{j}|.$ (2.2)

We denote the inner product between the measurement data and any function
$v\in C(\bar{\Omega})$ by $(m,v)_{n}=\frac{1}{n}\sum^{n}_{i=1}m_{i}v(x_{i})$.
Moreover, we denote the inner product between two functions by
$(u,v)_{n}=\frac{1}{n}\sum^{n}_{i=1}u(x_{i})v(x_{i})$ for any $u,v\in
C(\bar{\Omega})$ and the empirical norm
$\|u\|_{n}=\big{(}\frac{1}{n}\sum_{i=1}^{n}u^{2}(x_{i})\big{)}^{1/2}$ for any
$u\in C(\bar{\Omega})$.

Equipped with these definitions, we can define the parabolic inverse source
problem as to recover the unknown source term $f^{*}$ from the noisy final
time measurement data $m_{i}=\mathcal{S}f^{*}(x_{i})+e_{i}$, $i=1,...,n$. We
will use the regularization method to solve this parabolic inverse source
problem. Specifically, we will look for an approximate solution of the true
source term $f^{*}$ by solving the following least-squares regularized
minimization problem:

$\displaystyle\mathop{\rm min}\limits_{f\in
X}\|Sf-m\|^{2}_{n}+\lambda_{n}\|f\|_{L^{2}(\Omega)}^{2}\,,$ (2.3)

where $\lambda_{n}$ is the regularization parameter.

We first recall an existing work [7], where the optimal stochastic convergence
of regularized solutions and finite element solutions to time-dependent
parabolic inverse source problems (2.1) have been studied. The following
result represents the stochastic convergences corresponding to the random
variables with bounded variance.

######  Proposition 2.1 (Theorem 3.4 in [7])

Suppose $\\{e_{i}\\}^{n}_{i=1}$ are independent random variables satisfying
$\mathbb{E}[e_{i}]=0$ and $\mathbb{E}[e^{2}_{i}]\leq\sigma^{2}$. Let $f_{n}\in
L^{2}(\Omega)$ be the solution of the least-squares regularized minimization
problem (2.3). Then there exist constants $\lambda_{0}>0$ and $C>0$ such that
the following estimates hold for any $0<\lambda_{n}\leq\lambda_{0}$:

$\displaystyle\mathbb{E}\big{[}\|\mathcal{S}f_{n}-\mathcal{S}f^{*}\|^{2}_{n}\big{]}\leq
C\lambda_{n}\|f^{*}\|^{2}_{L^{2}(\Omega)}+\frac{C\sigma^{2}}{n\lambda^{d/4}_{n}},$
(2.4)
$\displaystyle\mathbb{E}\big{[}\|f_{n}-f^{*}\|^{2}_{L^{2}(\Omega)}\big{]}\leq
C\|f^{*}\|^{2}_{L^{2}(\Omega)}+\frac{C\sigma^{2}}{n\lambda^{1+d/4}_{n}},$
(2.5)
$\displaystyle\mathbb{E}\big{[}\|f_{n}-f^{*}\|^{2}_{H^{-1}(\Omega)}\big{]}\leq
C\lambda_{n}^{1/2}\|f^{*}\|^{2}_{L^{2}(\Omega)}+\frac{C\sigma^{2}}{n\lambda^{1/2+d/4}_{n}}.$
(2.6)

A stronger stochastic convergence when $\\{e_{i}\\}^{n}_{i=1}$ are independent
Gaussian random with variance $\sigma^{2}$ can also be found in [7]. Our
convergence analysis below also applies to this case. From Proposition 2.1, we
have an optimal choice of the regularization parameter $\lambda_{n}$ in (2.3)
as follows:

$\displaystyle\lambda_{n}^{1/2+d/8}=O\big{(}\sigma
n^{-1/2}\|f^{*}\|_{L^{2}(\Omega)}^{-1}\big{)}.$ (2.7)

It is an a priori estimate with knowing the noise level $\sigma$ and the
knowledge of the true source term $f^{*}$, yet it is our goal to recover the
true source term $f^{*}$ with unknown $\sigma$. The noise level is difficult
to measure in many circumstances. In [7], the authors proposed a self-
consistent algorithm to compute the optimal $\lambda_{n}$ without any
knowledge of the noise level $\sigma$ and the true source term $f^{*}$. Apply
the POD method to this algorithm, we will also successfully determine the
optimal parameter $\lambda_{n}$ iteratively. We find the proposed POD method
maintains the same accuracy as the FEM. Details can be found at the end of
Section 5.1.

An effective way to approximate the optimal control problem (2.3) with a
proper regularization parameter $\lambda_{n}$ is the iteration method. For
each iteration, one has to solve the forward problem (2.1) and its adjoint
problem for at least one time. In [7], the FEM was used to approximate (2.3)
and the optimal second-order convergence with respect to space has been
proved. But as the DOF grows for a discrete method such as FEM and FDM, the
iteration method will cost too much computational time. Therefore, we need to
design numerical methods that allow us to efficiently and accurately solve
(2.3).

##  3 The data-driven model reduction method

In this section, we will use the POD method to accelerate the inverse problem
computation. We first construct the POD basis functions from the snapshot
solutions of the parabolic equation (2.1) with some known type of source
function. Then, we solve the optimization problem in the inverse problems with
the constructed POD basis.

###  3.1 Construction of the POD basis functions

Assume that $u\in H^{1}_{0}(\Omega)$ is the solution to the following weak
formulation of the parabolic equation (2.1)

$\displaystyle(\partial_{t}u,\psi)+a(u,\psi)=(f,\psi),\quad\forall\psi\in
H^{1}_{0}(\Omega),\quad t\in[0,T],$ (3.1)

where $\Omega\subset\mathbb{R}^{d}$ and $a(\cdot,\cdot)$ is a bilinear form on
$H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)$ that is defined according to the
elliptic operator $\mathcal{L}$. Given a set of solutions at different time
instances
$\big{\\{}u(\cdot,t_{0}),u(\cdot,t_{1}),\ldots,u(\cdot,t_{m})\big{\\}}$ where
$t_{k}=k\Delta t$ with $\Delta t=\frac{T}{m}$, we first get the solution
snapshopts $\big{\\{}y_{1},\ldots,y_{m+1},y_{m+2},\ldots,y_{2m+1}\big{\\}}$,
where $y_{k}=u(\cdot,t_{k-1})$, $k=1,\ldots,m+1$, and
$y_{k}=\overline{\partial}u(\cdot,t_{k-m-1})$, $k=m+2,\ldots,2m+1$ with
$\overline{\partial}u(\cdot,t_{k})=\frac{u(\cdot,t_{k})-u(\cdot,t_{k-1})}{\Delta
t}$, $k=1,\ldots,m$.

Then, the POD basis functions $\\{\psi_{k}\\}_{k=1}^{{N_{pod}}}$ are built
from the solution snapshopts by minimizing the following error:

$\displaystyle\frac{1}{2m+1}\Big{(}\sum_{j=0}^{m}\|u(t_{j})-\sum_{k=1}^{{N_{pod}}}(u(t_{j}),\psi_{k})_{L^{2}(\Omega)}\psi_{k}\|_{L^{2}(\Omega)}^{2}+\sum_{j=1}^{m}\|\overline{\partial}u(t_{j})-\sum_{k=1}^{{N_{pod}}}(\overline{\partial}u(t_{j}),\psi_{k})_{L^{2}(\Omega)}\psi_{k}\|_{L^{2}(\Omega)}^{2}\Big{)}$
(3.2)

subject to the constraints that
$\big{(}\psi_{k_{1}}(\cdot),\psi_{k_{2}}(\cdot)\big{)}_{L^{2}(\Omega)}=\delta_{k_{1}k_{2}}$,
$1\leq k_{1},k_{2}\leq l$, where $\delta_{k_{1}k_{2}}=1$ if $k_{1}=k_{2}$,
otherwise $\delta_{k_{1}k_{2}}=0$.

Using the method of snapshot proposed by Sirovich [32], we know that the
optimization problem (3.2) can be reduced to an eigenvalue problem:

$Kv=\mu v,$ (3.3)

where the correlation matrix $K$ is computed from the solution snapshots
$\\{y_{1},y_{2},\ldots,y_{2m+1}\\}$ with entries
$K_{ij}=(y_{i},y_{j})_{L^{2}(\Omega)}$, $i,j=1,\ldots,2m+1$, and is symmetric
and semi-positive definite. We sort the eigenvalues in a decreasing order as
$\mu_{1}\geq\mu_{2}\geq...\geq\mu_{2m+1}>0$ and the corresponding eigenvectors
are denoted by $v_{k}$, $k=1,...,2m+1$. It can be shown that if the POD basis
functions are constructed by

$\varphi_{k}(\cdot)=\frac{1}{\sqrt{\mu_{k}}}\sum_{j=1}^{2m+1}(v_{k})_{j}u(\cdot,t_{j}),\quad
1\leq k\leq 2m+1,$ (3.4)

where $(v_{k})_{j}$ is the $j$-th component of the eigenvector $v_{k}$, they
will minimize the error (3.2). This result as well as the error formula were
proved in [19].

######  Proposition 3.1 (Sec. 3.3.2, [19] or p. 502, [4])

Let $\mu_{1}\geq\mu_{2}\geq...\geq\mu_{2m+1}>0$ denote the positive
eigenvalues of $K$ in (3.3). Then $\\{\psi_{k}\\}_{k=1}^{N_{pod}}$ constructed
according to (3.4) is the set of POD basis functions, and we have the
following error formula:

$\frac{\sum_{i=1}^{2m+1}\left|\left|y_{i}-\sum_{k=1}^{{N_{pod}}}\big{(}y_{i},\psi_{k}(\cdot)\big{)}_{L^{2}(\Omega}\psi_{k}(\cdot)\right|\right|_{L^{2}(\Omega)}^{2}}{\sum_{i=1}^{2m+1}\left|\left|y_{i}\right|\right|_{L^{2}(\Omega)}^{2}}=\frac{\sum_{k={N_{pod}}+1}^{2m+1}\mu_{k}}{\sum_{k=1}^{2m+1}\mu_{k}},$
(3.5)

where the number ${N_{pod}}$ is determined according to the ratio
$\rho=\frac{\sum_{k={N_{pod}}+1}^{2m+1}\mu_{k}}{\sum_{k=1}^{2m+1}\mu_{k}}$.

In practice, we shall make use of the decay property of eigenvalues in
$\mu_{k}$ and choose the first $m$ dominant eigenvalues such that the ratio
$\rho$ is small enough to achieve an expected accuracy, for instance
$\rho=1\%$. One would prefer the eigenvalues decays as fast as possible so
that one can ensure high accuracy with few POD basis functions.

###  3.2 A Fast algorithm for solving parabolic inverse source problems

We now propose the fast algorithm for solving parabolic inverse source
problems based on some given discretization, including our proposed POD
method. We first define the functional $\mathcal{J}$ to be

$\displaystyle\mathcal{J}[f]=\|\mathcal{S}f-m\|^{2}_{n}+\lambda_{n}\|f\|_{L^{2}(\Omega)}^{2}.$
(3.6)

Then, the least-squares regularized minimization problem (2.3) becomes to
solve the following misfit functional:

$\displaystyle\mathop{\rm min}\limits_{f\in L^{2}(\Omega)}\mathcal{J}[f].$
(3.7)

###### Lemma 3.2

The misfit functional $\mathcal{J}[f]$ is Fr$\acute{e}$chet-differentiable.

Proof. From the definition of Fr$\acute{e}$chet-differentiable, we need to
compute

$\displaystyle d\mathcal{J}[f](v)$ $\displaystyle=\lim_{t\rightarrow
0}\frac{\mathcal{J}[f+tv]-\mathcal{J}[f]}{t}$
$\displaystyle=(\mathcal{S}f-m,\mathcal{S}v)_{n}+\lambda_{n}(f,v)$
$\displaystyle=\big{(}\mathcal{S}^{*}(\mathcal{S}f-m),v\big{)}+\lambda_{n}(f,v)$
$\displaystyle=\big{(}\mathcal{S}^{*}(\mathcal{S}f-m)+\lambda_{n}f,v\big{)},$
(3.8)

where $v\in L^{2}(\Omega)$. In (3.2), the second equality is easily to get
from the quadratic form of the misfit functional $\mathcal{J}[f]$ and
$\mathcal{S}^{*}$ is the adjoint operator of $\mathcal{S}$ in the third
equality. Thus, one can directly obtain that

$\displaystyle d\mathcal{J}[f]=\mathcal{S}^{*}(\mathcal{S}f-m)+\lambda_{n}f.$
(3.9)

The formula (3.9) in Lemma 3.2 allows us to apply the gradient descent method
to minimize the discrepancy functional $\mathcal{J}[f]$. Let $f_{0}$ be the
initial guess and $f_{k}$ denote the the solution of the least-squares
regularized minimization problem (2.3) at the $k$-th iteration step, we update
the iterative solution by

$\displaystyle f_{k+1}=f_{k}-\alpha d\mathcal{J}[f_{k}],\quad\forall
k\in\mathbb{N},$ (3.10)

where $\alpha$ is the step size.

Given a fully discrete scheme with solution space $V_{dis}$, we define an
approximate forward operator $\mathcal{S}_{dis}$: $L^{2}(\Omega)\rightarrow
V_{dis}$, such that, given any source function $f$, $\mathcal{S}_{dis}f$ gives
the numerical solution at final time. In such discrete setting, we turn to
solve the following misfit functional:

$\displaystyle\mathop{\rm min}\limits_{f\in\Psi}\mathcal{J}_{dis}[f],$ (3.11)

where the functional
$\mathcal{J}_{dis}[f]=\|\mathcal{S}_{dis}f-m\|^{2}_{n}+\lambda_{n}\|f\|_{L^{2}(\Omega)}^{2}$.
We can compute the Fr$\acute{e}$chet derivative of $\mathcal{J}_{dis}[f]$ and
obtain the following iterative scheme:

$\displaystyle f_{k+1}=f_{k}-\alpha d\mathcal{J}_{dis}[f_{k}],\quad\forall
k\in\mathbb{N},$ (3.12)

where $\alpha$ is the step size,
$d\mathcal{J}_{dis}[f]=\mathcal{S}_{dis}^{*}(\mathcal{S}_{dis}f-m)+\lambda_{n}f$,
and $f_{0}$ is some initial guess.

The iterative scheme (3.12) works with both FEM-based Galerkin method and POD-
based Galerkin method. However, we emphasize that the DOF of the POD basis
functions is much smaller than that of the FEM space. Thus, the POD-based
Galerkin method can provide significant computational savings in solving the
parabolic equation (2.1) than the FEM, which allows us to quickly compute the
iterative scheme (3.12). Therefore, we obtain a fast algorithm in solving
parabolic inverse source problem.

###  3.3 Complete algorithm

Due to the nature of the inverse source problem, the source term is unknown.
We assume it to be a sample from some random space. Such random space may not
have a closed form or a finite-dimensional parameterization, so we further
assume that we have $N_{f}$ realizations, $\\{f_{l}\\}_{l=1,\cdots,N_{f}}$ of
the source term. In other words, $\\{f_{l}\\}_{l=1,\cdots,N_{f}}$ are some
possible ground truth source function in the inverse problem. Hence, deviating
from the classic POD algorithm in solving parabolic equation introduced in
Section 3.1, we include snapshots from each proposed source function into the
minimization (3.2). Estimation in (3.5) still holds as it is only an algebraic
result.

Finally, we summarize the proposed data-driven model reduction method for
solving parabolic inverse source problems in Algorithm 1, where the notations
have been defined before.

1 Input: Observation data $m$. Proposed source function
$\\{f_{l}\\}_{l=1,\cdots,N_{f}}$, error thresholds for constructing POD basis
functions $\epsilon>0$ and for the optimization problem $tol$, computational
time $T$, time step of the FEM $\Delta t$, $N_{T}=\lceil T/\Delta t\rceil$,
mesh size of the FEM $h$, one step finite element solver for the forward
problem $\mathcal{O}_{FEM}$, where
$\mathcal{S}_{FEM}f=(\mathcal{O}_{FEM}(f))^{N_{T}}u_{0}$, and step size in the
gradient descent method $\alpha$.

2 for  _$l=1:N_{f}$_ do

3 Solve the forward problem and store the snapshots
$(\mathcal{O}_{FEM}(f_{l}))^{i}u_{0}$ for $i=1,\cdots,N_{T}$.

4

5 end for

6Concatenate all snapshots as $\\{S_{j}\\}_{j=1,\cdots,N_{S}=N_{T}\times
N_{f}}=\cup_{i=1,\cdots,N_{T},l=1,\cdots,N_{f}}\\{(\mathcal{O}_{FEM}(f_{l}))^{i}u_{0}\\}$;

7 compute covariance matrix $K=(K_{ij})$ where
$K_{ij}=\frac{1}{N_{S}}(S_{i},S_{j})_{L^{2}(\Omega)}$;

8 compute the SVD of $K=(K_{ij})$ and denote the eigenvalues in a decreasing
order as $\mu_{1}\geq\mu_{2}\geq...\geq\mu_{N_{S}}>0$ and the corresponding
eigenvectors are denoted by $v_{k}$;

9 find minimal $N_{pod}$ such that
$\frac{\sum_{k=N_{pod}+1}^{N_{S}}\mu_{k}}{\sum_{k=1}^{N_{S}}\mu_{k}}<\epsilon$;

10 construct POD basis by
$\varphi_{k}(\cdot)=\frac{1}{\sqrt{\mu_{k}}}\sum_{j=1}^{N_{S}}(v_{k})_{j}S_{j},\quad
1\leq k\leq N_{pod}$;

11 construct Galerkin type solver for forward problem on POD basis
$\\{\varphi_{k}\\}_{k=1,\cdots,N_{pod}}$ as $S_{POD}$;

12 Set $f$ as an initial guess $f=f_{0}$.

13 while
_$\mathcal{J}_{pod}[f_{k}]=\|\mathcal{S}_{pod}f_{k}-m\|^{2}_{n}+\lambda_{n}\|f_{k}\|_{L^{2}(\Omega)}^{2}
>tol$ _ do

14 update $f_{k+1}\leftarrow f_{k}-\alpha d\mathcal{J}_{pod}[f_{k}]$ where
$d\mathcal{J}_{pod}[f_{k}]=\mathcal{S}_{pod}^{*}(\mathcal{S}_{pod}f_{k}-m)+\lambda_{n}f_{k}$.

15 end while

Output: computed source term $f$.

Algorithm 1 A fast algorithm for solving parabolic inverse source problem

##  4 Convergence analysis

In this section, we will first present the general discrete approximation to
the optimization problem (2.3). Recall that $V_{dis}\subset X$ be the discrete
function space with dimension $N_{dis}$ and $\mathcal{S}_{dis}$ be the
discrete approximation of the operator $\mathcal{S}:L^{2}(\Omega)\to
H^{2}(\Omega)$. We make the following assumptions on the discrete function
space $V_{dis}$ and the discretized operator $\mathcal{S}_{dis}$.

###### Assumption 4.1

(1) There exists an error estimate term $e_{dis}$ such that the discrete
operator $S_{dis}$ satisfies

$\displaystyle\|\mathcal{S}f-\mathcal{S}_{dis}f\|^{2}_{n}\leq
Ce_{dis}\|f\|^{2}_{L^{2}(\Omega)}.$ (4.1)

(2) For any $f\in L^{2}(\Omega)$, there exists $v\in V_{dis}$ such that,

$\displaystyle\lambda_{n}\|f-v\|^{2}_{L^{2}(\Omega)}+\|\mathcal{S}_{dis}(f-v)\|^{2}_{n}\leq
C(\lambda_{n}+e_{dis})\|f\|^{2}_{L^{2}(\Omega)}.$ (4.2)

In the general discrete approximation of the problem (2.3), we will solve the
least-squares regularized minimization problem:

$\displaystyle\mathop{\rm min}\limits_{f_{dis}\in
V_{dis}}\|\mathcal{S}_{dis}f_{dis}-m\|^{2}_{n}+\lambda_{n}\|f_{dis}\|_{L^{2}(\Omega)}^{2}.$
(4.3)

The following proposition gives the convergence analysis for the discretized
optimization schemes in solving the parabolic inverse source problems, which
is proved in [7].

######  Proposition 4.2 (Theorem 3.10 in [7])

Let $\\{e_{i}\\}^{n}_{i=1}$ are independent random variables satisfying
$\mathbb{E}[e_{i}]=0$ and $\mathbb{E}[e^{2}_{i}]\leq\sigma^{2}$ for
$i=1,\cdots,n$ and $f_{dis}$ be the solution of (4.3). Then there exist
constants $\lambda_{0}>0$ and $C>0$, such that for any
$\lambda_{n}\leq\lambda_{0}$ the following estimates hold true:

$\displaystyle\mathbb{E}\big{[}\|\mathcal{S}f^{*}-\mathcal{S}_{dis}f_{dis}\|_{n}^{2}\big{]}\leq
C(\lambda_{n}+e_{dis})\|f^{*}\|^{2}_{L^{2}(\Omega)}+C\left(1+\frac{e_{dis}}{\lambda_{n}}+\frac{N_{dis}e_{dis}}{\lambda^{1-d/4}_{n}}\right)\frac{\sigma^{2}}{n\lambda_{n}^{d/4}},$
(4.4)
$\displaystyle\mathbb{E}\big{[}\|f^{*}-f_{dis}\|_{L^{2}(\Omega)}^{2}\big{]}\leq
C\frac{\lambda_{n}+e_{dis}}{\lambda_{n}}\|f^{*}\|^{2}_{L^{2}(\Omega)}+C\left(1+\frac{e_{dis}}{\lambda_{n}}+\frac{N_{dis}e_{dis}}{\lambda^{1-d/4}_{n}}\right)\frac{\sigma^{2}}{n\lambda_{n}^{1+d/4}},$
(4.5)

and

$\displaystyle\mathbb{E}\big{[}\|f^{*}-f_{dis}\|_{H^{-1}(\Omega)}^{2}\big{]}$
$\displaystyle\leq
C(\lambda^{1/2}_{n}+e^{1/2}_{dis})\frac{\lambda_{n}+e_{dis}}{\lambda_{n}}\|f^{*}\|^{2}_{L^{2}(\Omega)}$
$\displaystyle+C(\lambda^{1/2}_{n}+e^{1/2}_{dis})\left(1+\frac{e_{dis}}{\lambda_{n}}+\frac{N_{dis}e_{dis}}{\lambda^{1-d/4}_{n}}\right)\frac{\sigma^{2}}{n\lambda_{n}^{1+d/4}}.$
(4.6)

Now we study the error analysis of the proposed method in solving the
parabolic inverse source problem, where we used the backward Euler scheme to
discretize the temporal derivative. Let $u$ be the solution to (3.1) and
$\hat{u}$ is the solution computed by semi-discrete scheme based on FEM basis
functions. We also assume the POD basis functions are constructed by snapshots
of $\hat{u}$ and their temporal finite differences. The FEM space and POD
space are denoted by $V_{h}$ and $\Psi$, respectively. The fully discrete
scheme is constructed on $\Psi$ and the solution is denoted by $U_{k}$ for
$k=1\cdots m$ with $m=\frac{T}{\Delta t}$. To be precise, we seek $U_{k}$ such
that

$(\bar{\partial}U_{k},\psi)+a(U_{k},\psi)=(f,\psi),\quad\forall\psi\in\Psi.$
(4.7)

Now we denote the solution operator from source term $f$ to the final time
solution $U_{m}$ as $\mathcal{S}_{pod}$, such that $\mathcal{S}_{pod}f=U_{m}$.

###### Theorem 4.3

Let $f_{pod}$ be the solution of (4.3), then there exist constants
$\lambda_{0}>0$ and $C>0$, such that for any $\lambda_{n}\leq\lambda_{0}$ the
following estimates hold true:

$\displaystyle\mathbb{E}\big{[}\|\mathcal{S}f^{*}-\mathcal{S}_{pod}f_{pod}\|_{n}^{2}\big{]}\leq
C(\lambda_{n}+e_{pod})\|f^{*}\|^{2}_{L^{2}(\Omega)}+C\left(1+\frac{e_{pod}}{\lambda_{n}}+\frac{{N_{pod}}e_{pod}}{\lambda^{1-d/4}_{n}}\right)\frac{\sigma^{2}}{n\lambda_{n}^{d/4}},$
(4.8)
$\displaystyle\mathbb{E}\big{[}\|f^{*}-f_{pod}\|_{L^{2}(\Omega)}^{2}\big{]}\leq
C\frac{\lambda_{n}+e_{pod}}{\lambda_{n}}\|f^{*}\|^{2}_{L^{2}(\Omega)}+C\left(1+\frac{e_{pod}}{\lambda_{n}}+\frac{{N_{pod}}e_{pod}}{\lambda^{1-d/4}_{n}}\right)\frac{\sigma^{2}}{n\lambda_{n}^{1+d/4}},$
(4.9)

and

$\displaystyle\mathbb{E}\big{[}\|f^{*}-f_{pod}\|_{H^{-1}(\Omega)}^{2}\big{]}$
$\displaystyle\leq
C(\lambda^{1/2}_{n}+e^{1/2}_{pod})\frac{\lambda_{n}+e_{pod}}{\lambda_{n}}\|f^{*}\|^{2}_{L^{2}(\Omega)}$
$\displaystyle+C(\lambda^{1/2}_{n}+e^{1/2}_{pod})\left(1+\frac{e_{pod}}{\lambda_{n}}+\frac{{N_{pod}}e_{pod}}{\lambda^{1-d/4}_{n}}\right)\frac{\sigma^{2}}{n\lambda_{n}^{1+d/4}},$
(4.10)

where

$e_{pod}=\left(h^{2}+\Delta t|\ln(\Delta t)|+\sqrt{\frac{T}{\Delta
t}\rho}\right)^{2}.$ (4.11)

Proof. To simplify the notation in the proof, let $\|\cdot\|$ denote the norm
in $L^{2}(\Omega)$, unless otherwise specified. For $k=1,...,m$ , we consider
the decomposition as follows:

$\displaystyle
U_{k}-\hat{u}\left(t_{k}\right)=U_{k}-\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right)+\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right)-\hat{u}\left(t_{k}\right)=\vartheta_{k}+\varrho_{k},$
(4.12)

where $\vartheta_{k}=U_{k}-\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right)$,
$\varrho_{k}=\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right)-\hat{u}\left(t_{k}\right)$,
and $\mathcal{P}^{{N_{pod}}}$ denotes the Ritz-projection to $\Psi$ under the
bilinear form $a(\cdot,\cdot)$. The projection error $\varrho_{m}$ can be
bounded by (3.5). In addition, by applying classic result of parabolic
equation, i.e., $\|u(\cdot,t_{j})\|\leq C\|f\|$, $\forall j=1,\cdots,m$ and
$\|u_{t}(\cdot,t_{j})\|\leq C\|f\|$, $\forall j=1,\cdots,m$, we know that,

$\displaystyle\varrho_{m}\leq C\sqrt{(2m+1)\rho}\|f\|.$ (4.13)

Now we only need to estimate the term $\vartheta_{m}$. According to (4.7), we
have

$\displaystyle\begin{aligned}
\left(\bar{\partial}\vartheta_{k},\psi\right)+a\left(\vartheta_{k},\psi\right)=&\left(\bar{\partial}U_{k},\psi\right)+a\left(U_{k},\psi\right)-\left(\bar{\partial}\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right),\psi\right)-a\left(\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right),\psi\right)\\\
=&\left(f\left(t_{k}\right),\psi\right)-\left(\bar{\partial}\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right),\psi\right)-a\left(\hat{u}\left(t_{k}\right),\psi\right)\\\
=&\left(v_{k},\psi\right),\end{aligned}$ (4.14)

where

$\displaystyle
v_{k}=\hat{u}_{t}\left(t_{k}\right)-\bar{\partial}\mathcal{P}^{{N_{pod}}}\hat{u}\left(t_{k}\right)=\hat{u}_{t}\left(t_{k}\right)-\bar{\partial}\hat{u}\left(t_{k}\right)+\bar{\partial}\hat{u}\left(t_{k}\right)-\bar{\partial}P^{{N_{pod}}}\hat{u}\left(t_{k}\right).$

We define
$w_{k}=\hat{u}_{t}\left(t_{k}\right)-\bar{\partial}\hat{u}\left(t_{k}\right)$
and
$z_{k}=\bar{\partial}\hat{u}\left(t_{k}\right)-\mathcal{P}^{{N_{pod}}}\bar{\partial}\hat{u}\left(t_{k}\right)$.
Substituting $\psi=\vartheta_{k}$ into (4.14), we derive that

$\displaystyle\left\|\vartheta_{k}\right\|^{2}-\left(\vartheta_{k},\vartheta_{k-1}\right)+\Delta
ta\left(\vartheta_{k},\vartheta_{k}\right)\leq\Delta
t\left\|v_{k}\right\|\left\|\vartheta_{k}\right\|$ (4.15)

and therefore,

$\displaystyle\left\|\vartheta_{k}\right\|\leq\frac{1}{1+\frac{a_{min}}{\alpha}\Delta
t}\left(\left\|\vartheta_{k-1}\right\|+\Delta t\left\|v_{k}\right\|\right),$
(4.16)

where $a_{min}$ is the minimum of the coefficient $a(x)$ in the elliptic
operator and $\alpha$ is a constant such that

$\displaystyle\|\phi\|^{2}_{L_{2}}\leq\alpha\|\phi\|_{H^{1}}\quad\forall\phi\in
V_{h}.$ (4.17)

Taking $\gamma=\frac{a_{min}}{\alpha}$, we have

$\displaystyle\left\|\vartheta_{m}\right\|\leq\left(\frac{1}{1+\gamma\Delta
t}\right)^{m}\left\|\vartheta_{0}\right\|+\Delta
t\sum_{j=1}^{m}\left(\frac{1}{1+\gamma\Delta
t}\right)^{m-j+1}\left\|v_{j}\right\|.$ (4.18)

As $\|v_{j}\|\leq\|w_{k}\|+\|z_{k}\|$, so (4.18) is taken into two parts.

$\displaystyle\Delta t\sum_{j=1}^{m}\xi^{m-j+1}\left\|z_{j}\right\|\leq$
$\displaystyle\Delta
t\sqrt{(\sum_{j=1}^{m}\xi^{2(m-j+1)})(\sum_{j=1}^{m}\|z_{k}\|^{2})},$ (4.19)

where $\xi$ denotes $\frac{1}{1+\gamma\Delta t}$. Using the facts that
$\left(\frac{1}{\xi}\right)^{2m}=(1+\gamma\Delta t)^{2m}=\left(1+\frac{2\gamma
T}{2m}\right)^{2m}\leq e^{2\gamma T}$ and $\xi^{-2}-1=(1+\gamma\Delta
t)^{2}-1\geq 2\gamma\Delta t$, we can obtain that

$\displaystyle\Delta
t\sqrt{(\sum_{j=1}^{m}\xi^{2(m-j+1)})(\sum_{j=1}^{m}\|z_{k}\|^{2})}=$
$\displaystyle
T\sqrt{(\frac{1}{m}\sum_{j=1}^{m}\xi^{2(m-j+1)})(\frac{1}{m}\sum_{j=1}^{m}\|z_{j}\|^{2})}$
$\displaystyle\leq$ $\displaystyle
T\sqrt{(\frac{1}{m}\frac{1-\xi^{2m}}{\xi^{-2}-1})(\frac{1}{m}\sum_{j=1}^{m}\|z_{j}\|^{2})}$
$\displaystyle\leq$ $\displaystyle T\sqrt{\frac{1-e^{-2\gamma T}}{2\gamma
T}(\frac{1}{m}\sum_{j=1}^{m}\|z_{j}\|^{2})}.$ (4.20)

By Corollary 4 in [23], we have that
$\frac{1}{m}\sum_{j=1}^{m}\|z_{j}\|^{2}\leq C\rho\|f\|^{2}$, which implies
that

$\displaystyle\Delta t\sum_{j=1}^{m}\xi^{m-j+1}\left\|z_{j}\right\|\leq
C(T)\sqrt{\rho}\|f\|.$ (4.21)

Now we turn to estimate the term $\Delta
t\sum_{j=1}^{m}\xi^{m-j+1}\left\|w_{j}\right\|$. By definition,

$\displaystyle\begin{aligned}
w_{j}&=u_{t}\left(t_{j}\right)-\frac{u\left(t_{j}\right)-u\left(t_{j-1}\right)}{\Delta
t}\\\ &=\frac{1}{\Delta t}\left(\Delta
tu_{t}\left(t_{j}\right)-\left(u\left(t_{j}\right)-u\left(t_{j-1}\right)\right)\right)\\\
&=\frac{1}{\Delta
t}\int_{t_{j-1}}^{t_{j}}\left(s-t_{j-1}\right)u_{tt}(s)ds.\end{aligned}$
(4.22)

By the analytic semigroup theory [30] or Lemma 3.6 in [7], we have the
following regularity estimate,

$\displaystyle\|u_{tt}(s)\|\leq\frac{C}{s}\|f\|,$ (4.23)

which implies,

$\displaystyle\|w_{j}\|\leq\frac{1}{\Delta
t}\int_{t_{k-1}}^{t_{j}}\left(s-t_{k-1}\right)\frac{1}{s}\|f\|ds.$ (4.24)

When $j=1$,

$\displaystyle\|w_{1}\|\leq\frac{1}{\Delta t}\int_{0}^{\Delta
t}\|f\|ds=\|f\|,$ (4.25)

and when $j\geq 2$,

$\displaystyle\|w_{k}\|\leq\frac{1}{\Delta
t}\int_{t_{j-1}}^{t_{j}}\left(s-t_{j-1}\right)\frac{1}{t_{j-1}}\|f\|ds=\frac{1}{2}\|f\|\frac{1}{(j-1)}\leq\frac{\|f\|}{j}.$
(4.26)

Therefore, we obtain that

$\displaystyle\Delta t\sum_{j=1}^{m}\xi^{m-j+1}\left\|w_{j}\right\|\leq\Delta
t\sum_{j=1}^{m}\xi^{m-j+1}\frac{\|f\|}{j}\leq\Delta
t\sum_{j=1}^{m}\frac{\|f\|}{j}\leq C\Delta t|\ln(\Delta t)|\|f\|.$ (4.27)

Summing up (4.21), (4.27) and using the classic parabolic FEM theory (e.g.
Theorem 1.2 in [33]), we have that

$\displaystyle\|\hat{u}\left(t_{m}\right)-u\left(t_{m}\right)\|\leq
Ch^{2}\|f\|.$ (4.28)

Therefore, we can prove that

$\displaystyle\|U_{m}-u(T)\|\leq C\left(h^{2}+\Delta t|\ln(\Delta
t)|+\sqrt{\frac{T}{\Delta t}\rho}\right)\|f\|.$ (4.29)

Finally, we apply the Proposition 4.2 by taking $\mathcal{S}_{pod}$ as
$\mathcal{S}_{dis}$ and $f_{pod}$ as $f_{dis}$, and prove the estimates in the
Theorem 4.3. Clearly, $e_{pod}=\left(h^{2}+\Delta t|\ln(\Delta
t)|+\sqrt{\frac{T}{\Delta t}\rho}\right)^{2}$.

###### Remark 4.4

Theorem 4.3 suggests that the optimal smoothing parameter in Eq.(2.7) may
still apply. In particular, if $e_{pod}\leq C\lambda_{n}$ and
${N_{pod}}e_{pod}\leq C\lambda_{n}^{1-1/\alpha}$, we have

$\displaystyle\mathbb{E}\big{[}\|\mathcal{S}f^{*}-\mathcal{S}_{pod}f_{pod}\|_{n}^{2}\big{]}\leq
C\lambda_{n}\|f^{*}\|_{L^{2}(\Omega)}^{2}+\frac{C\sigma^{2}}{n\lambda_{n}^{d/4}},$
(4.30)
$\displaystyle\mathbb{E}\big{[}\|f^{*}-f_{pod}\|_{L^{2}(\Omega)}^{2}\big{]}\leq
C\|f^{*}\|_{L^{2}(\Omega)}^{2}+\frac{C\sigma^{2}}{n\lambda_{n}^{1+d/4}},$
(4.31)

and

$\displaystyle\mathbb{E}\big{[}\|f^{*}-f_{pod}\|_{H^{-1}(\Omega)}^{2}\big{]}\leq
C\lambda^{1/2}_{n}\|f^{*}\|_{L^{2}(\Omega)}^{2}+\frac{C\sigma^{2}}{n\lambda_{n}^{1/2+d/4}}.$
(4.32)

Note that $\rho$ depends on the number of POD basis functions, ${N_{pod}}$,
and quickly decreases to zero when increasing ${N_{pod}}$. Assumptions above
are achievable by taking $h$ and $\Delta t$ to be small and ${N_{pod}}$ to be
small with respect to the DOF of the FEM.

##  5 Numerical examples

In this section, we present numerical examples to investigate the performance
of the proposed POD method in solving parabolic inverse source problems. Two
types of source functions, i.e., letters and circles at different locations,
will be studied. In the example of letters, we verify the convergence
introduced in Theorem 4.3 and apply the iterative method to find optimal
smoothing parameters. In the second example, we mainly study the approximation
property of our POD basis functions.

###  5.1 Recovering Letters

#### General Settings

For the first example, we apply the POD method to recover the source term $f$
in form of discontinuous patterns. To be precise, $f$ is an indicator
function, with the maximum value $1$, of a set in $[0,1]^{2}$ domain whose
shape is decided by a capital letter. For each observation data, we first
apply the FEM with $h=1/32$ and time step $\Delta t=1/32$ to get the solution
of the forward problem with exact source term $f$ at final time $T=1$. We
interpolate the FEM solution to locations of $100\times 100$ evenly
distributed sensor and add a Gaussian unbiased noise with stand derivation
$10^{-3}$ to each observation. It corresponds to a $10\%$ noise level. Initial
value $g=0$ in (2.1).

To construct the POD basis functions, we generate the snapshots source term
decided by letters from $A$ to $O$. The snapshots are generated by solving the
forward problem with each source term via FEM (space mesh size $h=1/32$, time
step $\Delta t=1/32$). Under such settings, $21$ POD basis functions are
constructed to resolve the snapshots of the forward problem with $10^{-4}$
accuracy. In other words, the POD basis functions cover $1-10^{-4}$ of
eigenspace from the snapshots with respect to the $L^{2}$ norm. To solve the
inverse problem, $\lambda$ is set to be $10^{-6}$ and we will discuss how to
find the optimal $\lambda$ later.

#### Computational Cost

In Figure 1, we show some of the letters recovered. The computational cost of
the POD method includes two parts. First, we will construct the POD basis
functions for all the $15$ letters, which will take about $40.3$ seconds on a
desktop computer (3.1GHz, 6-core, i5, Matlab). Then, for each letter, it will
take about $164$ seconds to solve the optimization problem via gradient
descend (relative $L_{2}$ threshold $10^{-5}$) using the POD basis functions.
As a comparison, with the FEM basis functions, it takes $1328.4$ seconds to
recover similar patterns with the same configuration during optimization. In
the case of recovering one single letter, the POD method has increased the
speed by $6.5\times$, and in the case of recovering all $15$ letters, it
achieves about $8\times$ speed-up.

(a) A (b) C (c) G (d) reference A (e) D (f) O

Figure 1: Recovering letters based on basis constructed from ’A’ to ’O’

We also compare the POD method with the FEM in the case when $h=1/20$. We
construct the POD basis from snapshots of the forward problem with $h=1/20$.
$23$ basis functions are constructed to achieve $10^{-4}$ accuracy. As the DOF
during the optimization procedure stays near, the computational cost with POD
basis functions is similar between different $h$’s (i.e., $47$ seconds in
constructing basis functions and $159$ seconds in solving optimization
problems). While with the FEM basis functions, the computation cost drops near
quadratically against the DOF. It takes $340$ seconds for $h=1/20$ case, which
corresponds to $400$ DOF. In Figure 2 we show the comparison between the FEM
and POD method with different $h$’s. The observation noise is the same under
the same $h$. Generally speaking, the restoration by the POD method is less
affected by the noise than the FEM. This is due to the fact that the noise is
perpendicular to the space of the POD basis functions. For a convincing
outcome, the POD method is significantly faster than the FEM.

(a) POD, $h=1/20$ (b) FEM, $h=1/20$

(c) POD, $h=1/32$ (d) FEM, $h=1/32$

Figure 2: Recovering letter ’A’ with different basis and $h$

#### Error Investigation

According to the convergence analysis studied in Theorem 4.3, there are three
factors that affect the accuracy of our algorithm, i.e., mesh size of FEM
basis, regularization factor, and approximation property of the POD basis
functions. As such, we set up three experiments when recovering the letter
’A’. Under each experiment, we generate $100$ set of observation $y$’s and
solve the optimization with the pre-computed POD basis functions. The standard
derivation of the noise is still $10^{-3}$. It corresponds to a $10\%$ noise
level. We compute the mean of the error in each set of parameters and show
them in Figure 3.

##### Dependence on the mesh size

In Figure 3(a) we verify that given $\lambda$ when the mesh size $h$ of the
FEM basis functions goes down, the error of $Sf$ will first decrease and then
converged to a limit that is greater than zero. In the first stage, the order
of convergence against mesh size $h$ is $2$. The limit in the second stage
depends on $\lambda$ and the dependence is non-monotone. The number of POD
basis functions does not change remarkably (also see in the previous example
that compares the POD method with FEM) under different mesh sizes. Along with
the monotonic decreasing of the error against $h\to 0$, it indicates all sets
of POD basis functions all approximate some fixed continuous function space.

##### Dependence on the regularization factor

To further investigate the limiting error, in 3(b), we keep $h=1/32$ fixed and
use different regularization factor $\lambda$’s. We find that the optimal
$\lambda$ for $h=1/32$ with the $10\%$ noise level is between $10^{-7}$ and
$10^{-8}$. The error is a V-shape function near the optimal $\lambda$.

##### Dependence on the number of POD basis functions

The last experiment is to study how the number of POD basis functions affects
the error. In 3(c), we show the mean square error with $100$ sets of random
generated observation. The error goes down as we increase the number of POD
basis functions. The fitted slope in the loglog plot is $-0.0145$, which
indicates the errors exponentially decay with a coefficient $-0.0145$. The
p-value for the fitted slope is smaller than $0.001$.

(a) Mesh size $h$ of the FEM.

(b) regularization factor $\lambda$

(c) POD basis number $n$, fitted slope: $1.45^{-2}$

Figure 3: Error dependence on different parameters

#### Iterative algorithm for finding the optimal $\lambda$

As analyzed in Eq.(2.7), the optimal $\lambda$ relies on
$\|f^{*}\|^{2}_{L^{2}(\Omega)}$. While in the setting of inverse problems, the
norm of input $f^{*}$ is unknown. A natural choice is to find the optimal
$\lambda$ by using a fixed-point iteration. To be precise, we start from an
initial guess like $\lambda_{n,0}=n^{-4/(d+4)}$ and solve for $f_{0}$ with the
POD basis functions. After solving $f_{j}$ the $j$-th step, we update
$\lambda$ as,

$\displaystyle\lambda_{n,j+1}^{1/2+d/8}=n^{-1/2}\left\|S_{\tau,h}f_{j}-m\right\|_{n}\left\|f_{j}\right\|_{L^{2}(\Omega)}^{-1}.$
(5.1)

We will stop the iteration of $\lambda$ until
$\left\|f_{j}\right\|_{L^{2}(\Omega)}$ converges. As shown in Figure 4 such
iterative algorithms, with either POD or FEM basis, finds the correct scale of
$\lambda$ with a similar pattern. It implies that the POD method captures the
whole solution subspace with given known source function and is robust in
iteration. To be noted, in Eq.(2.7), constant $C$ is not $1$ so the fixed-
point iteration Eq.(5.1) only finds the correct scale of $\lambda$ comparing
with Figure 3(b), in which the error due to $\lambda$ is computed with the
knowledge of the ground truth solution $u(\cdot,T)$.

(a) by iterating on FEM basis (b) by iteration on POD basis (c) $\lambda_{n}$
and $L_{2}$ error after each iteration

Figure 4: Finding optimal $\lambda$

###  5.2 Recovering Moving Circles

In addition to letter examples, we also applied the POD method to recover a
ball pattern that moves on some given trajectories in two-dimensional space.
Similar to recovering letters, we construct the POD basis functions from
snapshots that are generated by solving forward problems with classic FEM,
where the mesh size $h=1/32$ and time step $\Delta t=1/32$. $\lambda$ is set
to be $10^{-6}$ and $\sigma=2.75\times 10^{-3}$ which corresponds to a $10\%$
noise level.

#### Moving along a horizontal line

In such case, the source term $f$ is given by indicator function of set
$\Omega_{s}=\\{x\in[0,1]^{2}|(x_{1}-s)^{2}+(x_{2}-0.5)^{2}\leq 0.1^{2}\\}$.
When constructing the POD basis functions, we use $17$ values of $s$ that are
equally distributed between interval $[0.15,0.85]$ including two terminal
points. $14$ basis functions are extracted for $10^{-4}$ accuracy calculated
by eigenvalues. In 5 we use the constructed basis on detecting circles with
the same size. We find that the performance of the POD method is still
acceptable even the centers of the circle do not coincide with the ones when
we construct the POD basis functions.

#### Moving along a ring

In addition, as shown in 6 we recover the same pattern whose center are on a
circle. To be specific, the source term $f$ is given by indicator function of
set
$\Omega_{\theta}=\\{x\in[0,1]^{2}|(x_{1}-0.5-\cos(\theta)/4)^{2}+(x_{2}-0.5-\sin(\theta)/4)^{2}\leq
0.1^{2}\\}$. $33$ basis are extracted for $10^{-4}$ accuracy. Even with a
higher DOF, the POD method still yields $10\times$ acceleration comparing with
FEM ($60$ seconds vs $615$ seconds).

(a) $s=0.3$ (b) $s=0.4$ (c) $s=0.5$ (d) $s=0.6$ (e) $s=0.7$

Figure 5: Recovering Moving Circles with centers point at different $x$, red
circle: ground truth

(a) $\theta=0$ (b) $\theta=\frac{\pi}{4}$ (c) $\theta=\frac{\pi}{2}$ (d)
$\theta=\frac{3\pi}{4}$

(e) $\theta=\pi$ (f) $\theta=\frac{5\pi}{4}$ (g) $\theta=\frac{3\pi}{2}$ (h)
$\theta=\frac{7\pi}{4}$

Figure 6: Recovering Moving Circles with centers point at different $x$, red
circle: ground truth

##  6 Conclusion

We have developed a data-driven and model-based approach for solving parabolic
inverse source problems. The key idea is to exploit (model-based) and
construct (data-based) the intrinsic approximate low-dimensional structure of
the underlying parabolic PDEs that consists of two components – a training
component that constructs a small number of POD basis functions to achieve
significant dimension reduction in the solution space, and a fast algorithm
that computes the optimization problem in inverse source problem by using the
constructed POD basis functions. Hence, we achieve an effective data and
model-based approach for the inverse source problems and overcome the typical
computational bottleneck of FEM in solving such problems.

Under a weak assumption on the regularity of the solution, we provide the
convergence analysis of our POD algorithm in solving the forward parabolic
PDEs and thus obtain the error estimate of the POD algorithm for the parabolic
inverse source problems. Finally, we carried out numerical experiments to
demonstrate the accuracy and efficiency of the proposed method. Through
numerical results, we found that our POD algorithm yields as good
approximations as the reference solution obtained by the FEM. But our
algorithm can be much cheaper. We expect an even better performance of
efficiency can be obtained in 3D problems, which is one of our future topics.
We also studied other issues of the POD algorithm, such as the dependence of
the error on the mesh size, on the regularization parameter in the least-
squares regularized minimization problems, and on the number of POD basis
functions.

## Acknowledgement

The research of W. Zhang is supported by the National Natural Science
Foundation of China No. 11901282 and the Shenzhen Sci-Tech Fund No.
RCBS20200714114941241. The research of Z. Zhang is supported by Hong Kong RGC
grant projects 17300318 and 17307921, National Natural Science Foundation of
China No. 12171406, Seed Funding Programme for Basic Research (HKU), and Basic
Research Programme (JCYJ20180307151603959) of The Science, Technology and
Innovation Commission of Shenzhen Municipality.

## References

  * [1] V. Akçelik, G. Biros, A. Draganescu, J. Hill, O. Ghattas, and B. Waanders.  Dynamic data-driven inversion for terascale simulations: real-time identification of airborne contaminants.  In SC’05: Proceedings of the 2005 ACM/IEEE Conference on Supercomputing, pages 43–43. IEEE, 2005.
  * [2] A. Alla and M. Falcone.  A time-adaptive POD method for optimal control problems.  IFAC Proceedings Volumes, 46(26):245–250, 2013.
  * [3] P. Benner, S. Gugercin, and K. Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM Review, 57(4):483–531, 2015.
  * [4] P. Benner, S. Gugercin, and K. Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM Review, 57(4):483–531, 2015.
  * [5] G. Berkooz, P. Holmes, and J. Lumley.  The proper orthogonal decomposition in the analysis of turbulent flows.  Annual review of fluid mechanics, 25(1):539–575, 1993.
  * [6] Z. Chen, R. Tuo, and W. Zhang.  Stochastic convergence of a nonconforming finite element method for the thin plate spline smoother for observational data.  SIAM Journal on Numerical Analysis, 56(2):635–659, 2018.
  * [7] Z. Chen, W. Zhang, and J. Zou.  Stochastic convergence of regularized solutions and their finite element approximations to inverse source problems.  arXiv:2104.02352, 2021.
  * [8] Jin Cheng and Jijun Liu.  An inverse source problem for parabolic equations with local measurements.  Applied Mathematics Letters, 103:106213, 2020.
  * [9] Wei Cheng, Ling-Ling Zhao, and Chu-Li Fu.  Source term identification for an axisymmetric inverse heat conduction problem.  Computers & Mathematics with Applications, 59(1):142–148, 2010\.
  * [10] S. D’haeyer, B. Johansson, and M. Slodička.  Reconstruction of a spacewise-dependent heat source in a time-dependent heat diffusion process.  IMA Journal of Applied Mathematics, 79(1):33–53, 2014.
  * [11] A. El Badia, A. El Hajj, M. Jazar, and H. Moustafa.  Lipschitz stability estimates for an inverse source problem in an elliptic equation from interior measurements.  Applicable Analysis, 95(9):1873–1890, 2016.
  * [12] A. El Badia, T. Ha Duong, and F. Moutazaim.  Numerical solution for the identification of source terms from boundary measurements.  Inverse Problems in Engineering, 8(4):345–364, 2000.
  * [13] A. El Badia and T. Nara.  An inverse source problem for Helmholtz’s equation from the Cauchy data with a single wave number.  Inverse Problems, 27(10):105001, 2011.
  * [14] G. Garcia, A. Osses, and M. Tapia.  A heat source reconstruction formula from single internal measurements using a family of null controls.  Journal of Inverse and Ill-posed Problems, 21(6):755–779, 2013\.
  * [15] S. Gorelick, B. Evans, and I. Remson.  Identifying sources of groundwater pollution: an optimization approach.  Water Resources Research, 19(3):779–790, 1983.
  * [16] H. Gu, J. Xin, and Z. Zhang.  Error estimates for a POD method for solving viscous G-equations in incompressible cellular flows.  SIAM Journal on Scientific Computing, 43(1):A636–A662, 2021.
  * [17] J. Hadamard.  Lectures on Cauchy’s problem in linear partial differential equations.  Courier Corporation, 2003.
  * [18] J. Hesthaven, G. Rozza, and B. Stamm.  Certified reduced basis methods for parametrized partial differential equations.  Springer, 2016.
  * [19] P. Holmes, J. Lumley, and G. Berkooz.  Turbulence, coherent structures, dynamical systems and symmetry.  Cambridge University Press, 1998.
  * [20] V. Isakov.  Inverse problems for partial differential equations, volume 127\.  Springer, 2006.
  * [21] V. Isakov, S. Leung, and J. Qian.  A three-dimensional inverse gravimetry problem for ice with snow caps.  Inverse Problems & Imaging, 7(2):523, 2013.
  * [22] T. Johansson and D. Lesnic.  Determination of a spacewise dependent heat source.  Journal of computational and Applied Mathematics, 209(1):66–80, 2007.
  * [23] K. Kunisch and S. Volkwein.  Galerkin proper orthogonal decomposition methods for parabolic problems.  Numerische Mathematik, 90(1):117–148, 2001.
  * [24] K. Kunisch, S. Volkwein, and L. Xie.  HJB-POD-based feedback design for the optimal control of evolution problems.  SIAM Journal on Applied Dynamical Systems, 3(4):701–722, 2004.
  * [25] X. Liu.  Identification of indoor airborne contaminant sources with probability-based inverse modeling methods.  Boulder, CO: PhD dissertation. University of Colorado at Boulder, 2008.
  * [26] X. Liu and Z. Zhai.  Inverse modeling methods for indoor airborne pollutant tracking: literature review and fundamentals.  Indoor air, 17(6):419–438, 2007.
  * [27] P. Nelson and S. Yoon.  Estimation of acoustic source strength by inverse methods: Part i, conditioning of the inverse problem.  Journal of sound and vibration, 233(4):639–664, 2000.
  * [28] G. Nunnari, A. Nucifora, and C. Randieri.  The application of neural techniques to the modelling of time-series of atmospheric pollution data.  Ecological Modelling, 111(2-3):187–205, 1998.
  * [29] A. Quarteroni, A. Manzoni, and F. Negri.  Reduced basis methods for partial differential equations: an introduction, volume 92.  Springer, 2015.
  * [30] M. Renardy and R. Rogers.  An introduction to partial differential equations, volume 13.  Springer Science & Business Media, 2006.
  * [31] L. Sirovich.  Turbulence and the dynamics of coherent structures. I. Coherent structures.  Quarterly of applied mathematics, 45(3):561–571, 1987.
  * [32] L. Sirovich.  Turbulence and the dynamics of coherent structures. I. Coherent structures.  Quarterly of applied mathematics, 45(3):561–571, 1987.
  * [33] V. Thomée.  The finite element method for parabolic problems.  In Mathematical Theory of Finite and Boundary Element Methods, pages 135–218. Springer, 1990.
  * [34] S. Volkwein.  Proper orthogonal decomposition: Theory and reduced-order modelling.  Lecture Notes, University of Konstanz, 4(4), 2013.
  * [35] J. Wong and P. Yuan.  A FE-based algorithm for the inverse natural convection problem.  International Journal for numerical methods in fluids, 68(1):48–82, 2012.