A non-intrusive reduced basis EKI for time-fractional diffusion inverse problems

TL;DR

This paper introduces a non-intrusive reduced basis method combined with ensemble Kalman inversion to efficiently solve time-fractional diffusion inverse problems, significantly reducing computational costs while maintaining accuracy.

Contribution

The paper presents a novel POD-DSRBF surrogate modeling approach that accelerates EKI for TFDIPs by decoupling offline and online computations.

Findings

Achieves significant computational speed-up.

Maintains high accuracy in inverse problem solutions.

Demonstrated effectiveness on nonlinear TFDIPs.

Abstract

In this study, we consider an ensemble Kalman inversion (EKI) for the numerical solution of time-fractional diffusion inverse problems (TFDIPs). Computational challenges in the EKI arise from the need for repeated evaluations of the forward model. We address this challenge by introducing a non-intrusive reduced basis (RB) method for constructing surrogate models to reduce computational cost. In this method, a reduced basis is extracted from a set of full-order snapshots by the proper orthogonal decomposition (POD), and a doubly stochastic radial basis function (DSRBF) is used to learn the projection coefficients. The DSRBF is carried out in the offline stage with a stochastic leave-one-out cross-validation algorithm to select the shape parameter, and the outputs for new parameter values can be obtained rapidly during the online stage. Due to the complete decoupling of the offline and online stages, the proposed non-intrusive RB method -- referred to as POD-DSRBF -- provides a powerful tool to accelerate the EKI approach for TFDIPs. We demonstrate the practical performance of the proposed strategies through two nonlinear time-fractional diffusion inverse problems. The numerical results indicate that the new algorithm can achieve significant computational gains without sacrificing accuracy.

Tables

Table 1. Table 1: The numerical results given by two different methods.

Direct EKI			RB-EKI
$\delta (\%)$	${\overline{\theta }}^1$	${\overline{\theta }}^2$	$\#$ iter	CPU(s)	${\overline{\theta }}^1$	${\overline{\theta }}^2$	$\#$iter	CPU(s)
1	0.2034	0.7017	37	305.11	0.2003	0.7004	42	4.03
3	0.2082	0.6931	32	267.74	0.2130	0.6959	31	2.98
5	0.2426	0.7005	9	76.13	0.2484	0.6985	9	0.87

Table 2. Table 2: The numerical results with four other source locations given by RB-EKI.

$\theta$	${\overline{\theta }}^1$	${\overline{\theta }}^2$	$\#$ iter	CPU(s)
(0.8, 0.4)	0.7969	0.4101	30	2.89
(0.8, 0.2)	0.7572	0.1736	30	2.88
(0.9, 0.2)	0.7691	0.1032	47	4.53
(0.9, 0.1)	0.7765	0.0483	46	4.43

Equations

$$\{\begin{array}{ll}{}^c{D}_t^{\alpha }u(x,t)=\nabla \cdot (\kappa (x)\nabla u(x,t))+f(x,t),& x\in \Omega ,t\in (0,T)\\ \mathcal{B}u(x,t)=g(x,t),& x\in \partial \Omega ,t\in (0,T)\\ u(x,0)=u_0(x),& x\in \Omega ,\end{array}$$

$${}^c{D}_t^{\alpha }u(x,t)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{\partial u(x,s)}{\partial s}\frac{ds}{(t-s{)}^{\alpha }},0<\alpha <1,$$

$$A(u_h(t;\theta ),\theta )=0,$$

$$y_{obs}=g(u_h;\theta )+\xi .$$

$${\bar{\theta }}_n=\frac{1}{N_e}\sum _{j=1}^{N_e}{\theta }_n^{(j)},{\bar{\omega }}_n=\frac{1}{N_e}\sum _{j=1}^{N_e}f({\theta }_n^{(j)}),$$

$$C_n^{\theta \omega }=\frac{1}{N_e-1}\sum _{j=1}^{N_e}({\theta }_n^{(j)}-{\bar{\theta }}_n)(f({\theta }_n^{(j)})-{\bar{\omega }}_n{)}^T,$$

$$C_n^{\omega \omega }=\frac{1}{N_e-1}\sum _{j=1}^{N_e}(f({\theta }_n^{(j)})-{\bar{\omega }}_n)(f({\theta }_n^{(j)})-{\bar{\omega }}_n{)}^T.$$

$${\theta }_{n+1}^{(j)}={\theta }_n^{(j)}+C_n^{\theta \omega }(C_n^{\omega \omega }+{\gamma }_n\Gamma {)}^{-1}(y_{n+1}^{(j)}-{\omega }_n^{(j)}),$$

$$y_{n+1}^{(j)}=y_{obs}+{\xi }_{n+1}^{(j)},{\xi }_{n+1}^{(j)}\sim N(0,\Gamma ),$$

$${\gamma }_n^N\Vert {\Gamma }^{1/2}(C_n^{\omega \omega }+{\gamma }_n^N\Gamma {)}^{-1}(y^{(j)}-{\bar{\omega }}_n)\Vert \ge \rho \Vert {\Gamma }^{-1/2}(y^{(j)}-{\bar{\omega }}_n)\Vert .$$

$$\Vert {\Gamma }^{-1/2}(y_{obs}-{\bar{\omega }}_n)\Vert \le \tau \Vert {\Gamma }^{-1/2}(y_{obs}-f({\theta }^{†}))\Vert ,$$

$$\mathbf{S}={\mathbf{U\Sigma V}}^T,$$

$${\mathbb{V}}_p=\text{mbox}Span\{{\mathbf{u}}_1,{\mathbf{u}}_2,\cdots ,{\mathbf{u}}_p\}.$$

$$\underset{{\mathbf{U}}_p\in {\mathbb{R}}^{n_h\times p}}{\min }\Vert \mathbf{S}-{\mathbf{U}}_p{\mathbf{U}}_p^T\mathbf{S}{\Vert }_F^2=\sum _{k=p+1}^Q{\sigma }_k^2,$$

$$\frac{{\sum }_{i=1}^p{\sigma }_i^2}{{\sum }_{i=1}^Q{\sigma }_i^2}>{ϵ}_{pod},$$

$$u_p(t;\theta )=\sum _{k=1}^p{\mathbf{u}}_k{a}_k(t;\theta ),$$

$$\displaystyle\begin{array}[]{lr}\mathcal{D}=\Big{\{}\{(t,\theta),\mathbf{U}^{T}u_{h}(t;\theta)\}:t\in\mathcal{T}_{tr},\theta\in\Theta_{tr}\Big{\}},\\ \mathcal{T}_{tr}=\{t_{1},\cdots,t_{{N^{tr}_{t}}}\}\subset\mathcal{T},\\ \Gamma_{tr}=\{\theta_{1},\cdots,\theta_{{N^{tr}_{\theta}}}\}\subset\Theta.\\ \end{array}$$

$${\mathcal{U}}_{Z,\mathcal{E}}:=\text{Span}\{\phi ({\epsilon }_j\Vert \cdot -z_j{\Vert }_{{\ell }^2}|z_j\in Z,{\epsilon }_j\sim \mathcal{E}\}$$

$$u(x)=\sum _{z\in Z}{c}_j{\phi }_j(x):=\sum _{z\in Z}{c}_j\phi ({\epsilon }_j\Vert x-z_j{\Vert }_{{\ell }^2})$$

$$u_{X,Z,\epsilon }:=\arg \underset{u\in {\mathcal{U}}_{Z,\mathcal{E}}}{\inf }\sum _{x\in X}|u-f{|}^2.$$

$$A_{\epsilon }(Z,Z)\mathbf{c}=f(Z),$$

$${\epsilon }^{\ast }=\arg \underset{\epsilon >0}{\min }\Vert e(\epsilon ){\Vert }_{{\ell }^2}.$$

$$[A_{\epsilon }(Z,Z){]}_{ij}=\phi ({\epsilon }_j(z_i-z_j))\text{mbox}and[f(Z){]}_j=f(z_j)$$

$${\mathbf{Q}}_k=[a_k(t_i;{\theta }_j){]}_{i,j},1\le i\le N_t^{tr},1\le j\le N_{\theta }^{tr},$$

$${\mathbf{Q}}_k\approx {\tilde{\mathbf{Q}}}_k=\sum _{l=1}^{q_k}{\lambda }_l^k{\Psi }_l^k({\Phi }_l^k{)}^T.$$

$$\displaystyle\begin{array}[]{lr}t\mapsto\tilde{\Psi}_{l}^{k}(t),\,\text{ trained from}\{(t_{i},\Psi_{l}^{k}(t_{i})):i=1,\cdots,N^{tr}_{t}\},\\ \theta\mapsto\tilde{\Phi}^{k}_{l}(\theta),\,\text{trained from}\{(\theta_{i},\Phi_{l}^{k}(\theta_{i})):i=1,\cdots,N^{tr}_{\theta}\},\\ \end{array}$$

$$a_k(t;\theta )\approx {\tilde{a}}_k(t;\theta )=\sum _{l=1}^{q_k}{\lambda }_l^k{\tilde{\Psi }}_l^k({\tilde{\Phi }}_l^k{)}^T,(t,\theta )\in \mathcal{T}\times \Theta .$$

$$\tilde{a}(t,\theta )=\Phi (t,\theta ;\epsilon ),$$

$${\theta }_0^{(j)}\sim \pi (\theta ),y^{(j)}=y_{obs}+{\xi }^{(j)},{\xi }^{(j)}\sim N(0,\Gamma ),j=1,\cdots ,N_e.$$

$${\omega }_n^{(j)}={\tilde{f}}_p({\theta }_n^{(j)}),j=1,\cdots ,N_e$$

Full text

1 Introduction

Consider the following time-fractional diffusion equation:

[TABLE]

where $\Omega$ is a bounded domain with boundary $\partial\Omega$ and $T>0$ is a fixed value, $u$ represents the state variable, $k(x)$ is the diffusion coefficient and the right hand side $f$ denotes source terms. Here $\mathcal{B}$ is a boundary operator, ${}^{c}D^{\alpha}_{t}\,(0<\alpha<1)$ denotes the Caputo fractional derivative of order $\alpha$ with respect to $t$ and it is defined by

[TABLE]

where $\Gamma(\cdot)$ is the Gamma function. The system (1) has been widely used to model anomalous subdiffusion in which the mean square variance grows slower than that in a Gaussian process, and has found a number of important practical applications [1]. The forward problems for time-fractional diffusion equation (TFDE) have been extensively studied, see [2, 3, 4] and references therein. However, in practice there are many inputs unknown in the time fractional diffusion model, such as the fractional derivatives, the heat source and the diffusion coefficient and so on. These are time fractional diffusion inverse problems (TFDIPs). In general, the TFDIPs are much more difficult to solve than the direct problem for TFDE due to the fact that the solution does not satisfy the general conditions of well-posedness [5]. There are, however, many works on both theoretical and numerical aspects over the last few years. We refer to the papers [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16] for a rather incomplete list, which range from recovering structural parameters of TFDEs to source terms and initial states.

In this paper, we consider the ensemble Kalman-based inversion method (EKI) [17, 18, 19] to solve TFDIPs. The EKI is a derivate-free approach that lies at the interface between the optimization and Bayesian approaches [20, 21]. Since it only requires the evaluation of the forward operator but not its derivative, this approach has been successfully used to solve PDE-constrained parameter identification problems [18]. A potential drawback in the application of EKI to inverse problems is it requires a sufficiently large ensemble size to guarantee reliable estimations. This is undesirable to solve the TFDIPs as each ensemble realization requires a solution of the forward model of TFDE. It is well known that the simulation of fractional-order systems is computationally demanding due to their long-range history dependence. To address this challenge, surrogate modeling has been extensively developed for decades, see [23, 24, 25, 26, 27]. The key idea is to replace the high-fidelity full order system with a carefully constructed computationally cheap approximated model, to reduce the memory needs and CPU time. Specifically, the surrogate model is often constructed offline and subsequently used online when running the EKI algorithm.

Projection-type reduced order models (ROM) are one possible realization of this idea. Here, the governing equations of the forward model are projected onto a limited number of reduced basis (RB) functions, which can be calculated via the proper orthogonal decomposition [28] or with greedy algorithm [29], to obtain an inexpensive reduced order model. During the online stage, for a new parameter value, the reduced model is constructed as a linear combination of the pre-computed RB functions, where the expansion coefficients are computed by projecting the full-order equation onto the RB space. However, such scheme is problem-dependent and of an intrusive nature. Additionally, it causes computational inefficiency for complex nonlinear problems. To tackle these challenges, non-intrusive methods [30, 31] have been developed to construct a surrogate of the reduced coefficient so that the reduced coefficient can be recovered without requiring a projection of the full-order model. In this method, the data-fit surrogate models are constructed using interpolation and regression techniques to directly learn the map between the input parameters and the POD expansion coefficients. There is also substantial past work on ROM for time-dependent problems, e.g. the non-intrusive frameworks based on radial basis function interpolation [32, 33], and the Gaussian process regression [34]. One can refer to [35] for a comprehensive survey.

Motivated by the recent developments in reduced order modeling, we propose a new non-intrusive RB method, namely POD-DSRBF, for constructing surrogate models to accelerate the EKI for TFDIPs. After extracting the RB functions from a set of snapshots by POD, a learning approach using a doubly stochastic radial basis function (DSRBF) is used to establish a mapping from parameter values to projection coefficients onto the RB space. Equipped with a stochastic leave-one-out cross-validation (LOOCV) parameter selection for the shape parameter of the RBF, the efficiency of the DSRBF can be further enhanced. We shall discuss the basic idea and the efficient implementation of the algorithms. To the best of our knowledge, this is the first investigation of the EKI combining with the non-intrusive RB method to solve TFDIPs. Numerical results indicate that the RB-EKI approach shows good performance in both accuracy and efficiency for TFDIPs.

The structure of the paper is as follows. In the next section, we review the formulation of time fractional diffusion inverse problems and the solution approach via EKI. In section 3, we introduce the POD-DSRBF approach to construct accurate surrogate models. In section 4, we use two different time fractional diffusion inverse problems to demonstrate the accuracy and efficiency of our method. Finally, we conclude the paper with Section 5.

2 Background and problem formulation

In this section, we first give a brief overview of the time fractional diffusion inverse problems (TFDIPs). Then we will introduce the EKI to solve TFDIPs.

2.1 Mathematical formulation

Unlike integer order differential equations, the analytical solutions of time fractional diffusion equation are usually not available. We consider a numerical discretization of the system (1), described by

[TABLE]

where $u_{h}(t;\theta):\mathcal{T}\times\Theta\rightarrow\mathbb{R}^{n_{h}}$ is the discrete solution, $n_{h}$ is the dimension of the finite-dimensional discretization of the spatial domain, and $\theta\in\Theta\subset\mathbb{R}^{d}$ is the $d$-dimensional parameter vector. The discrete operator $A$ denotes a numerical approximation, e.g., by the finite element or finite difference method. Typically, the $n_{h}$ is large in order to resolve the details of the system with high accuracy. The goal of TFDIPs is to estimate the unknown parameters $\theta$ from noisy observations of the states $u_{h}$ given by

[TABLE]

Here $g$ is a discretized observation operator mapping from the states and parameters to the observables, $\xi\in\mathbb{R}^{m}$ is the measurement error. We assume that the error $\xi$ is a Gaussian random vector with mean zero and covariance matrix $\Gamma\in\mathbb{R}^{m\times m}$, i.e., $\xi\sim N(0,\Gamma)$. The system model (2) and observation model (3) together define a forward model $y=f(\theta)$ that maps the unknown parameter to the observable data.

2.2 Ensemble Kalman inversion

Ensemble Kalman inversion (EKI) is a recently proposed inversion methodology that lies at the interface between the optimization and Bayesian approaches [17, 20]. It applies the ensemble Kalman filter (EnKF) [22] to the inverse problem setting by introducing a trivial dynamics for the unknown parameters. In this work, we follow closely the framework in [18].

To formulate the TFDIPs in a Bayesian framework, we model the parameter $\theta$ as a random variable, endow it with a prior distribution $\pi(\theta)$. Assume we derive $N_{e}$ initial ensemble $\theta_{0}^{(j)}\,(j\in\{1,\cdots,N_{e}\})$ from the prior $\pi(\theta)$. The algorithm works by iteratively updating an ensemble of candidate solutions $\{\theta^{(j)}_{n}\}^{N_{e}}_{j=1}$ from iteration index $n$ to $n+1$. Define the ensemble mean

[TABLE]

and convariances

[TABLE]

[TABLE]

Then the EKI update formulae are given by

[TABLE]

where

[TABLE]

and the regularization parameter $\gamma_{n}$ is choosing by

[TABLE]

As an iterative regularization method, the iterative EKI is terminated according to the following discrepancy principle

[TABLE]

where $\tau$ is a constant and $\theta^{{\dagger}}$ denotes the truth properties.

Noted that the total cost of an $N_{e}$ size ensemble of EKI is approximately $N_{e}J$ forward model evaluations where $J$ is the total number of iterations. Therefore, each iteration could involve over $N_{e}$ TFDE solvers even for relatively small $J$. Hence, the cost of EKI for TFDIPs is prohibited as it requires a relatively large ensemble size to guarantee its performance. It is thus natural to construct a surrogate of the forward model before the data are available. In the next section, we will focus on the reduced basis method, which is widely used in applied mathematics and engineering.

3 The reduced basis EKI method

3.1 The proper orthogonal decomposition and the reduced basis space

The goal of model order reduction is to reduce the high dimensional model (2), in the sense of reducing the number of degrees of freedom, with the expectation of computational savings in terms of both storage and CPU times. One of the convenient tools for model order reduction is the reduced basis (RB) method. It seeks the approximate solution to (2) in a reduced space spanned by a set of parameter-independent RB functions.

To generate an RB space, we consider the set of $Q=n_{t}n_{\theta}$ snapshots, $\{u_{h}(t_{i};\theta_{j})|i=1,\cdots,n_{t},j=1,\cdots,n_{\theta}\}$, which are snapshots at $n_{t}$ different time instances $t_{1},\cdots,t_{n_{t}}\in\mathcal{T}$ and $n_{\theta}$ different inputs $\theta_{1},\cdots,\theta_{n_{\theta}}\subset\Theta$. Define the snapshot matrix $\mathbf{S}\in\mathbb{R}^{n_{h}\times Q}$, which collects the snapshots $u_{h}(t_{i};\theta_{j})$ as its columns. Thus, each row in the snapshot matrix corresponds to a spatial location and each column corresponds to a snapshot.

With this snapshot matrix, a (thin) singular value decomposition (SVD) is then performed to obtain the reduced space:

[TABLE]

where $\mathbf{U}\in\mathbb{R}^{n_{h}\times Q}$ and $\mathbf{V}\in\mathbb{R}^{Q\times Q}$ are the left and right singular vectors of $\mathbf{S}$ respectively. The singular values $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{Q}\geq 0$ of $\mathbf{S}$ give the diagonal matrix $\mathbf{\Sigma}=\mbox{diag}\{\sigma_{1},\sigma_{2},\cdots,\sigma_{Q}\}\in\mathbb{R}^{Q\times Q}$. The POD basis $\mathbf{U}_{p}=[\mathbf{u}_{1},\mathbf{u}_{2},\cdots,\mathbf{u}_{p}]\in\mathbb{R}^{n_{h}\times p}$ is defined as the $p$ left singular vectors of $\mathbf{S}$ that correspond to the $p$ largest singular values. Then the RB space $\mathbb{V}_{p}$ is defined as

[TABLE]

Among all orthonormal bases of size $p$, the POD basis minimizes the error of snapshot reconstruction

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm. Thus, the singular values provide a guidance for choosing the size of the POD basis. A typical approach is to choose $p$ such that

[TABLE]

where $\epsilon_{pod}$ is a user-specified tolerance.

Once the RB space is available, the field $u_{h}$ can be approximated by a linear expansion in the POD basis:

[TABLE]

where $a_{k}(t;\theta)$ denote the POD expansion coefficients. One can evaluate the POD coefficient for a given parameter value $\theta$ by interpolation [36] or the Galerkin procedure [28]. In this study, we propose a learning approach, namely doubly stochastic radial basis function (DSRBF), to calculate the POD expansion coefficients. Note that, given a snapshot $u_{h}(t;\theta)$, we can compute its representation in the POD basis via the coefficients $a_{k}(t;\theta)=\mathbf{u}_{k}^{T}u_{h}(t;\theta),\,k=1,\cdots,p$, where we have used that the POD basis vectors are orthonormal. Our learning task is now transformed into learning a model for the POD coefficients $a_{k}(t;\theta)$ from a training data

[TABLE]

3.2 Doubly stochastic radial basis function

In this subsection, we briefly describe the doubly stochastic radial basis function (DSRBF) method. We refer interested readers to [37] for more details of the development of the DSRBF strategy.

In DSRBF, we first pick a translation-invariant radial kernel $K=\varphi(\|\boldsymbol{\cdot}-\boldsymbol{\cdot}\|_{\ell^{2}}):\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ with a function $\varphi:\mathbb{R}\to\mathbb{R}$ known as the radial basis function (RBF). Commonly used RBF, include the Gaussian $\varphi(r)=\exp(-r^{2})$, the multiquadrics (MQ) $\varphi(r)=(1+r^{2})^{\beta/2}$ with $\beta=1$ and inverse multiquadrics (IMQ) with $\beta=-1$. Then, we can define the finite-dimensional doubly stochastic trial space in the form of

[TABLE]

for some quasi-uniform trial centers $Z$ and some stochastic shape parameters following the probability distribution $\mathcal{E}$.

For learning problems using DSRBF, one seeks an approximant $u\in\mathcal{U}_{Z,\mathcal{E}}$ of an unknown function $f$ from a collection of examples, called training data. Recall that any trial function is a linear combination of the basis used in defining (12) and is in the form of

[TABLE]

with $\varepsilon_{j}\sim\mathcal{E}$ and for some coefficients $\mathbf{c}=[c_{1},\ldots,c_{n_{Z}}]^{T}\in\mathbb{R}^{n_{Z}}$. To have a well-posed fully-discretized problem, we make an observation $\varepsilon=\{\varepsilon_{j}\}_{j=1}^{n_{Z}}\sim\mathcal{E}^{n_{Z}}$ and define the least-squares numerical solution

[TABLE]

Using the training data $\{(z_{i},y_{i}):z_{i}\in Z,y_{i}=f(z_{i}),i=1,\ldots,n_{Z}\}$, we can yield the square unsymmetric system

[TABLE]

where $\mathbf{c}$ is the vector of unknown coefficients, and $f(Z)$ is the data vector. The matrix $A_{\varepsilon}(Z,Z)$ is an $n_{Z}\times n_{Z}$ coefficient matrix with entries $[A_{\varepsilon}(Z,Z)]_{ij}=\varphi(\varepsilon_{j}\|z_{i}-z_{j}\|_{\ell^{2}})$ for $z_{i},z_{j}\in Z$. Therefore, all unknown coefficients $\{c_{j}\}$ can be easily obtained by a standard matrix solver. Once $\{c_{j}\}$ have been determined, the approximation solution $u_{X,Z,\varepsilon}$ can be evaluated from (13).

For efficiency, one should choose an ‘optimal’ shape parameter $\varepsilon_{j}\sim\mathcal{E}$. In this work, we use a stochastic leave-one-out cross validation (LOOCV) proposed in [37] to select this parameter. The stochastic LOOCV method basically consists of the following two steps:

• (1)

For the linear equation $A_{\varepsilon}(Z,Z)\mathbf{c}=f(Z)$, whose matrix depends on a single shape parameter $\varepsilon>0$, we collect a set of observations $\widetilde{\varepsilon}=[\widetilde{\varepsilon}_{1},\ldots,\widetilde{\varepsilon}_{N_{\text{obs}}}]^{T}$ of the ‘optimal’ shape parameter.

• (2)

Generate a random vector $[\varepsilon_{j}]_{j=1}^{n_{Z}}\sim\mathcal{E}^{n_{Z}}=[\chi^{2}(\text{mean}(\widetilde{\varepsilon}))]^{n_{Z}}$.

The key ingredient of stochastic LOOCV is how to select the ‘optimal’ shape parameter $\widetilde{\varepsilon}_{j}$. To this end, we define a stochastic LOOCV cost vector $e(\varepsilon)=[e_{1},\ldots,e_{n_{Z}}]^{T}$ using Algorithm 1. Then the stochastic-LOOCV optimal shape parameter can be defined as

[TABLE]

For $j=1$ to some ${N_{\text{obs}}}>0$, we apply (15) to estimate its optimal shape parameter by some minimization algorithm. Then, store result as $\widetilde{\varepsilon}_{j}$.

We now present the DSRBF strategy. The DSRBF method consists of the following steps which are outlined here:

• •

Generate a random vector $[\varepsilon_{j}]_{j=1}^{n_{Z}}\sim\mathcal{E}^{n_{Z}}=[\chi^{2}(\text{mean}(\widetilde{\varepsilon}))]^{n_{Z}}$ using stochastic LOOCV.

• •

Construct the matrix system in (14), where

[TABLE]

for $1\leq i,j\leq n_{Z}$.

• •

Seek for numerical solution in the form of (13) by solving (14) for the expansion coefficient $\mathbf{c}\in\mathbb{R}^{n_{Z}}$.

3.3 RB-based EKI algorithm

In this subsection, we will use the DSRBF combining with a tensor decomposition [34] to learn the POD expansion coefficients. For the $k$th entry $a_{k}=\mathbf{u}_{k}^{T}u_{h}$ of the projection coefficients, the training data can be written in a matrix as

[TABLE]

as a result of the tensor grid between the time and the parameter locations in the training data.

We employ the SVD again to decompose the data of an expansion coefficient into several time- and parameter-modes

[TABLE]

Here $\Psi_{l}^{k}$ and $\Phi^{k}_{l}$ are the $l$th discrete time- and parameter-modes for the $k$th projection coefficient, respectively, $\lambda^{k}_{l}$ is the $l$-th singular value for the same coefficient, and $q_{k}\,(1\leq q_{k}\leq p)$ is the corresponding truncation number.

DSRBFs are now trained to approximate the corresponding continuous modes as

[TABLE]

where $\tilde{\Psi}_{l}^{k}(t)$ and $\tilde{\Phi}^{k}_{l}(\theta)$ are the $l$th continuous time and parameter DSRBF modes for the $k$th expansion coefficient, respectively. Then, a continuous learning function $\tilde{a}_{k}(t,\theta)$ for the $k$th expansion coefficient $a_{k}$ can be recovered as

[TABLE]

The resulted DSRBF approximation of the POD coefficient vector $a(t;\theta)$ can be also represented as follows:

[TABLE]

where $\Phi(t,\theta;\varepsilon)$ is the trained DSRBF model for the POD coefficient $a(t;\theta)$ and $\varepsilon$ is the shape parameter of the DSRBF. The complete POD-DSRBF algorithm is displayed in Algorithm 2.

Once the DSRBF is trained, we can predict the reduced coefficient $\tilde{a}(t;\theta^{*})$ for a new given parameter $\theta^{*}$, during the online stage. The corresponding reduced solution is given by

[TABLE]

where $\mathbf{U}_{p}$ is the POD basis set. Then the associated model outputs are

[TABLE]

If $p\ll n_{h}$, the dimension of the unknown state in (18) and (19) are greatly reduced compared with that of the original full model (2) and (3). Thus, (18) and (19) define a reduced-order model $y_{p}=\tilde{f}_{p}(\theta)$ that maps the parameter $\theta$ to an approximation of the model outputs $y_{p}$.

It is clear that after obtaining the reduced-order model, we can then replace the forward model $f$ by its approximation $\tilde{f}_{p}$, and obtain the RB-based EKI algorithm. The pseudo-code of the RB-based EKI algorithm is presented in Algorithm 3. Notice that the computational cost of generating $N_{e}$ samples using $\tilde{f}_{p}$ in prediction step requires nothing but sampling of the DSRBF model of (17) with $N_{e}$ samples of $\theta$. This cost is minimal because it does not require any simulations of the forward model.

4 Numerical Examples

In order to illustrate the accuracy and efficiency of the RB-EKI approach for solving the time fractional diffusion inverse problems, in this section, we present numerical experiments with two type different inverse problems. The first example, adapted from [38, 39], considers a heat source inversion problem. The second example is the problem of inferring the spatially-varying diffusion coefficient [27, 40]. For numerical examples, we set $\Omega=[0,1]^{2}$ and $T=1$. In all our tests, unless otherwise specified, we shall use the following parameters $\rho=1/\tau=0.7,{N_{\text{rv}}}=15,{N_{\text{obs}}}=10$. We employ the MQ function in the DSRBF and the MATLAB function fminbnd to solve (15). All the computations are performed using MATLAB 2015a on an Intel-i5 desktop computer.

4.1 Example 1: 2D heat source inversion

In this example, we consider the following model

[TABLE]

The aim is to determine the source location $\theta=(\theta^{1},\theta^{2})$ from noisy measurements of the $u$-field at a finite set of locations and times. The prior distributions on $\theta^{i}$ are independent and uniform, i.e., $\theta^{i}\sim U(0,1)$. Below, unless stated otherwise, the observations are generated by selecting the values of the states at a uniform $3\times 3$ sensor network. At each sensor location, three measurements are taken at time $t=\{0.25,0.75,1\}$, which corresponds to a total of 27 measurements. For any given values of $\theta$, the Eq. (23) is solved by using a finite difference/ spectral approximations ([3]) with time step $\Delta t=0.01$ and polynomial degree $P=6$. In order not to commit an ‘inverse crime’, we generate the data by solving the forward problem at a much higher resolution than that used in the inversion, i.e., with $P=10$. In the examples below, unless otherwise specified, the Caputo fractional derivative of order $\alpha$ is 0.5.

4.1.1 Solution to the forward problem

In this subsection, we assess the accuracy of the reduced order models for the forward model. We construct a reduced order model using POD-DSRBF whose solution captures the output of the forward model over the support of the prior distribution. It should be noted that at this stage, measurement data do not enter the inference procedure. Therefore, we can construct the surrogate model off-line and save it for future use.

To generate the snapshot set, we sample 100 uniformly distributed parameter location of $\Theta$. Among $100$ time steps, $n_{t}=50$ ones are included in the database, and this results in an ensemble of $5000$ snapshots. The domain $\Omega$ is discretized to $21\times 21$ regular grid. The decay of the singular values of the snapshot matrix indicates the approximation power of the POD basis to the information contained in the snapshots. The first 30 singular values are shown in Figure 1. The singular values decay exponentially, and this indicates that the information contained in the snapshots can be well captured by a few POD basis functions. Through the POD, a set of $p=6$ orthogonal bases are extracted as the RB functions with a tolerance $\epsilon_{pod}=99.99\%$, shown in Figure 2.

To validate the accuracy of the POD-DSRBF, the RB solutions are recovered at 3 time points for $M=400$ different values of $\theta$. We define the relative error $\epsilon_{a}$ and the projection error $\epsilon_{p}$ as

[TABLE]

[TABLE]

respectively. We also compute the coefficient learning error $\epsilon_{c}$

[TABLE]

Figure 3 illustrates the evolution of the approximation error $\epsilon_{a}$ and the projection error $\epsilon_{p}$ with respect to the number of POD basis, $p$, obtained by DSRBF and the classical RBF. We can observe that both DSRBF and RBF enjoy fast decay with respect to the number of $p$ and resemble the projection error. As the number of $p$ increases, classical RBF saturates quickly and while DSRBF can continue to decrease and saturates at a lower level. The saturation is because the coefficient learning error is dominant over the projection error when the number of $p$ is large. When more training data is available, the saturation level can be further reduced. Figure 4 shows the approximation error $\epsilon_{a}$ and the coefficient learning error $\epsilon_{c}$ with respect to the size of training data for a fixed number of $p=20$. As we increase the number of training set, both DSRBF and RBF get more accurate results. Moreover, the DSRBF outperform the classical constant LOOCV-optimal shape parameter formulations. Figure 5 shows the CPU times for finding LOOCV-optimal shape parameters and the stochastic-LOOCV algorithm respectively. It can be seen that DSRBF took less than one-tenth of the time of the classical LOOCV algorithm. This indicates that the DSRBF can effectively improve the accuracy of the results.

4.1.2 Solution to the inverse problem

In this section, it will be shown how the RB model can be used for accelerating the EKI for solving the model inverse problem, i.e., Example 1. To better present the results, we shall perform the following two types of approaches: the conventional EKI or direct EKI based on the full model evaluation and RB-EKI based on the POD-DSRBF surrogate model.

In order to assess the performance of both methods, we define the relative error and data misfit as

[TABLE]

and

[TABLE]

Here, $\bar{\theta}$ are $\theta^{{\dagger}}$ are the ensemble mean approximates and exact solutions respectively. For simplicity, we consider a diagonal measurement error covariance $\Gamma=\delta^{2}I$. We first set the exact parameter $\theta=(0.2,0.7)$. For our EKI methods, we choose $N_{e}=100$ ensemble members.

The exact and numerical results for the source location $\theta$, obtained using full model and RB model with various levels of noise added into the data, namely, $\delta={1\%,3\%,5\%}$, are presented in Figure 6. The corresponding relative error $e_{\theta}$ at each iteration is displayed on the left side of Figure 7. From these figures, we can see that though the distance between the initial ensemble mean and the truth is large, the approximation converges towards the exact solution as the iteration progresses. In particular, the final ensemble mean approximates the truth very well, especially for small noise deviation. The data misfit along with a line marking the stopping point defined by (7) is shown on the right side of Figure 7. That is, once the data misfit reaches the dashed line, the convergence criterion defined by (7) is considered satisfied. From Figures 6 and 7, we can see that with the stopping criterion (7), the EKI can yield an accurate and stable approximation. In addition, Figure 7 clearly shows that both the relative error $e_{\theta}$ and the data misfit $E_{\theta}$ decrease as the iteration progresses. Moreover, the smaller the deviation in the data is, the more iterations are needed, which means that the cost of obtaining a good approximation from the initial ensemble to the truth increases. This is possibly due to the support of the induced posterior distribution shrinks when the standard deviation in the data decreases [41].

Table 1 summarize the numerical results with various levels of noise $\delta$ in the data obtained using the direct full model and POD-DSRBF model, respectively. In the tables, $\bar{\theta}$ denote the final iteration reconstruction arising from EKI. The ensemble mean of $\theta$ is in excellent agreement with the exact solution with up to $5\%$ noise in the data. However, the estimate becomes less accurate with an increased noise level, see also Figure 7. The numerical results obtained by the two approaches are nearly the same, which is consistent with the discussion from the last section that the POD-DSRBF model is accurate enough to approximate the full model. The foregoing numerical verifications indicate that the EKI with the full model simulation and POD-DSRBF model could yield practically identical numerical results for the inverse heat source problem. However, the computing time required by RB-EKI is only a small fraction of that by direct full model simulation. To see that, consider the case with $\delta=3\%$, the RB-EKI takes $2.98s$ to complete the simulation, whereas direct EKI takes $267.74s$. Therefore, the speedup of the RB-EKI approach over the direct EKI approach is dramatic.

Besides the great computational savings, another advantage of using RB-EKI is that it is reusable when there are new measurements coming after constructing the POD-DSRBF model. Table 2 shows the results for four other source locations with $3\%$ noise in the data. The POD-DSRBF model is the same as before. It can be seen that the method can always infer the exact value without performing any additional direct forward simulation as in the direct EKI. This demonstrates the flexibility of the RB-EKI approach for solving the time fractional heat source identification problems.

4.1.3 Possibility of recovering the fractional order $\alpha$

In this section, we consider to recover the fractional order $\alpha$ in the model (23) with the source location simultaneously. In this case, the unknown parameter is $\theta=(\theta^{1},\theta^{2},\alpha)\in\mathbb{R}^{3}$. For the results presented below, a uniform prior was adopted for $\theta$ on the parameter space $[0,1]^{3}$. To generate data for inversion, the ‘true’ parameter is choose to be $(\theta^{1}=0.25,\theta^{1}=0.75,\alpha=0.8)$. In the following tests, the measurements were taken at three different times $t=\{0.25,0.75,1\}$ using 25 locations on a uniform $5\times 5$ grid covering the domain $\Omega$, thus leading to $75$ measurements. To construct the POD-DSRBF model, we sample the parameter space $[0,1]^{3}$ at 200 uniformly distributed random points, and add the full-order solutions at 50 time points $t\in\mathcal{T}_{tr}=\{0.02,0.04,\cdots,1\}$ into the database $\mathcal{D}$. These samples will be used as snapshots and training data. In this example, we use $p=10$ POD basis functions to construct the RB space, and the truncation number of SVDs for all expansion coefficients are set as $q_{k}=10$.

The numerical results by RB-EKI and the direct EKI with $5\%$ noise in the data are shown in Figure 8. The approximation converges towards the exact solution as the iteration progresses. In particular, the final ensemble mean is in excellent agreement with the exact one for both RB-EKI and direct EKI. The result by RB-EKI is practically identical to that by the direct EKI, but its computational expense is only a small fraction of that of direct EKI. The CPU time of evaluating the conventional EKI is about $256.86s$, while the CPU time of RB-EKI is about $2.82s$.

4.2 Example 2: Estimating the diffusion coefficient

In this example, we illustrate the RB-EKI approach on the nonlinear inverse problem of estimating the diffusion coefficient. Consider the following two dimensional time-fractional PDEs

[TABLE]

The goal is to determine the diffusion coefficient $\kappa(x)$ from noisy measurements of the $u$-field at a finite set of locations and times. We assume the log-diffusivity field $\log\kappa(x)$ is endowed with a Gaussian process prior, with mean zero and an isotropic covariance kernel:

[TABLE]

for which we choose variance $\sigma^{2}=1$ and a length scale $l=0.2$. This prior allows the field to be easily parameterized with a Karhunen-Loeve expansion:

[TABLE]

where $\lambda_{i}$ and $\phi_{i}(x)$ are the eigenvalues and eigenfunctions, respectively, of the integral operator on $[0,1]^{2}$ defined by the kernel $K$, and the parameter $\theta^{i}$ are endowed with independent standard normal priors. These parameters then become the targets of inference. In the numerical simulation, we truncate the Karhunen-Loeve expansion at $d=9$ modes that preserve $90\%$ energy of the prior distribution. The true solution is directly drawn from the prior distribution. The measurement sensors of $u$ are using 49 locations on uniform $7\times 7$ grid covering the domain $\Omega$. Similar to Example 1, at each sensor location, three measurements are taken at time $t=\{0.25,0.75,1\}$. The observational errors are taken to be additive and Gaussian:

[TABLE]

with $\xi_{j}\sim N(0,\delta^{2})$.

To construct the POD-DSRBF model, we sample the parameter space at 500 random points, and the full-order database is constructed in the same way as Example 1. All full-order samples will be used as snapshots and training data. It is observed that $p=10$ POD basis functions already portray accurately the solution manifold, with the total energy $\epsilon_{pod}=99.99\%$. The true diffusion coefficient and the initial ensemble mean of the EKI are shown in Figure 9. The final iteration reconstruction by RB-EKI for Example 2 with various levels of noise in the data is shown in Figure 10. With up to $5\%$ noise in the data, the ensemble mean is still in good agreement with the exact one. Taking into consideration the ill-posedness of the problems, the results presented here are quite satisfactory. Direct EKI simulation is not conducted in this example since it is expected to much more expensive.

In order to study the convergence of the results, we also plot the relative errors and data misfits of RB-EKI at each iteration. The numerical results together with a horizontal line marking the stopping points defined by the (7) are displayed in Figure 11. The accuracy error $e_{\theta}$ of the RB-EKI exhibits the typical ‘semiconvergence’ phenomenon: the approximation solution converges towards the exact solution up to a certain iteration number, and beyond this point, it deviates from the exact solution. Therefore, an appropriate stopping rule is crucial to obtain an accurate and stable approximation solution. The right side of Figure 11 shows that the data misfits curve is relatively flat over a wide range of iteration numbers, and thus the rule (7) requires an accurate estimation of the amount of data noise, in order that the RB-EKI could yield a good approximation. However, this is not always available in practical situations, and heuristic approaches, such as the generalized cross-validation and the L-curve criterion, may be employed. We leave this to our next exploration.

5 Summary

In this paper we developed an efficient non-intrusive RB ensemble Kalman inversion approach to solving the time fractional diffusion inverse problems. In the framework of this method, a reduced order model is constructed by POD-DSRBF. The identification of the optimal shape parameter of the DSRBF is performed during the offline stage through a stochastic LOOCV algorithm. Then the prediction steps of the EKI are calculated from a large number of realizations generated by the reduced order model with virtually no additional computational cost. Thus, the computational cost can be significantly reduced. With the accuracy and the efficiency validated by numerical examples, the proposed RB-EKI algorithm is shown to be a powerful tool for solving time fractional diffusion inverse problems.

Ongoing works will extend the techniques to more complex fractional order model inverse problems, e.g. time-space fractional diffusion equation. Note that the reduced models under investigation are global in that they are valid over the prior parametric domain, which underlies the challenge of the high dimensionality of the parameter space as well as the intrinsic properties of the forward model. Strategies for building multi-fidelity reduced models adaptively would be very interesting.