{'fraction_non_alphanumeric': 0.03740133363669776, 'fraction_numerical': 0.017484292139944343, 'mean_word_length': 4.414897649734647, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving linear Bayesian inverse problems using a fractional total variation-Gaussian (FTG) prior and transport map

The Bayesian inference is widely used in many scientific and engineering problems, especially in the linear inverse problems in infinite-dimensional setting where the unknowns are functions. In such problems, choosing an appropriate prior distribution is an important task. In particular, when the function to infer has much detail information, such as many sharp jumps, corners, and the discontinuous and nonsmooth oscillation, the so-called total variation-Gaussian (TG) prior is proposed in function space to address it. However, the TG prior is easy to lead the blocky (staircase) effect in numerical results. In this work, we present a fractional order-TG (FTG) hybrid prior to deal with such problems, where the fractional order total variation (FTV) term is used to capture the detail information of the unknowns and simultaneously uses the Gaussian measure to ensure that it results in a well-defined posterior measure. For the numerical implementations of linear inverse problems in function spaces, we also propose an efficient independence sampler based on a transport map, which uses a proposal distribution derived from a diagonal map, and the acceptance probability associated to the proposal is independent of discretization dimensionality. And in order to take full advantage of the transport map, the hierarchical Bayesian framework is applied to flexibly determine the regularization parameter. Finally we provide some numerical examples to demonstrate the performance of the FTG prior and the efficiency and robustness of the proposed independence sampler method.


Introduction
The Bayesian inference methods [10,14] have been popular among solving inverse problems, where the estimation results uncertainties can be quantified by learning the statistical information such as moments, confidence intervals, and the marginalizing. A typical Bayesian inference is learning from the measurement data by incorporating some prior information to yield the posterior distribution and has different features compared with classical deterministic regularization methods. As we know, the most practical inverse problems are often highly ill-posed due to the limited and noisy measurement data and the prior distribution plays a significant role on the performance of the Bayesian inference. Recently the Bayesian inverse problems in the infinite dimensional setting have been extensively studied and [7,32] privode a Bayesian framework for the inverse problems in function spaces where the unknowns are of infinite dimension. The Gaussian measures are widely used as the prior distributions for the infinite dimensional Bayesian inverse problems. In fact such a choice has many advantages about theories and computation in the infinite dimensional Bayesian inference [7].
In our work, we focus on the linear inverse problem in infinite-dimensional setting. Let X be separable Hilbert space, equipped with the Borel σ-algebra, and A : X → R n a linear forward model (unknows-toobservation). We wish to solve the inverse problem of finding the unknown function u ∈ X (in this paper we shall restrict ourselves to the situation where u is a real-valued function defined in R m , m = 1, 2) from measurement y ∈ R n , which is usually generated by y = Au + η, (1.1) where the noise η is assumed to be a n-dimensional zero-mean Gaussian random variable with covariance matrix Σ η . In Bayesian inversion, the unknown model input and the measurement data are usually regarded as random variables. The posterior distribution can be in a closed form only in some special cases. For example, if the prior and noise are Gaussian, the posterior distribution is Gaussian for the linear forward model. However the Gaussian prior distributions is usually applied to recover the unknowns with smooth property and in many practical problems, especially in the image reconstructions, are not suitable for modeling such inversion functions with smooth, sharp jumps or discontinuities [43]. In order to model such functions, the total variation (TV) regularization [30] was proposed in the deterministic inverse problem context. The TV prior for the Bayesian inference is applied to a range of inverse problems [1,18]. However, with the discretization dimension increasing, the posterior distribution based on TV prior may not converge to a well-defined infinite-dimensional measure [18]. And the TV prior distributions differ significantly from the Gaussian measures, which many analysis techniques and numerical implementation methods based on Gaussian priors can not be simply and directly extended to non-Gaussian priors. Recently, [43] has proposed a TV-Gaussian (TG) prior to deal with this issue and the corresponding numerical method is given. The main idea of the TG prior is that the TV term is effective for identifying the sharp jumps and the Gaussian distribution is mainly used to ensure the posterior measure is well-defined in the function space.
In fact, although the total variational regularization can catch the sharp jumps of the unknowns, the reconstruction results are easy to have the blocky effect [38]. And the Riemann-Liouville fractional order derivative based total variational regularization constructed in [45] has a better result than that of the total variational regularization in addressing the image reconstruction problems. Inspired by the success of this nonlocal regularization, in this paper, we propose a hybrid FTV-Gaussian (FTG) prior that we shall replace the total variational term in the hybrid TG prior with a fractional order total variational (FTV) term. Namely, it uses a FTV term to capture the detail information in the unknowns and the Gaussian distribution makes sure that the posterior distribution is a well-defined probability measure in the function space. The FTG prior can effectively reduces blocky effects and have better performance than the TG for the detail information in the unknowns, especially for recovering textures in images. More details are present in Section 4 about the numerical examples.
Practically we need to explore the posterior distribution in the Bayesian inference. Markov chain Monte Carlo (MCMC) simulations are flexible and widely used to draw samples from the complex posterior probability distribution in Bayesian inferences [34,27,10]. As we all know, the standard MCMC algorithms, such as Metropolis-Hastings algorithm, can become arbitrarily slow as the discretization mesh of the unknown function is refined [4,28], and in this case the algorithm is said to be dimension-dependent. [4] presented a dimension-independent MCMC algorithm that is obtained by constructing a preconditioned Crank-Nicolson (pCN) discretization of a stochastic partial differential equation that preserves the reference measure. The sampling efficiency of MCMC algorithms for the infinite-dimensional Bayesian inverse problems can be further improved by incorporating the data information in the proposal design that can reflect the local or global geometry of the target distribution, like doing in finite dimensional counterparts, such as the stochastic Newton MCMC [20], the dimension-independent likelihood-informed MCMC [5] and the adaptive independence sampler [8]. The proposal distribution for the independence sampler is represented as a mixture of a finite number of specially parametrized Gaussian measures in [8]. [24] constructed a more effective and adaptive proposal based transport maps to accelerate MCMC and in [25], a multifidelity preconditioner is proposed to increase the efficiency of MCMC sampling. This multifidelity approach exploits low-fidelity models to construct a proposal distribution that approximates the posterior distribution via a transport map, then uses this proposal distribution to perform MCMC sampling of the original (high-fidelity) posterior distribution. However, in this work, we propose a proposal distribution based on a diagonal transport map to accelerate the independence sampler for the linear inverse problem in infinite-dimensional setting. And the proposed independence sampler based on a diagonal map is a dimension-independent method, which the acceptance probability associated to the proposal derived from the transport map is independent of discretization dimensionality. In order to make full use of the transport map, the hierarchical Bayesian framework [11] is also applied here, where the regularization parameter can be flexibly determined. Loosely speaking, this approach has two stages: firstly we construct a diagonal transport map that can approximately push forward the reference measure to the posterior measure (not low-fidelity models in [25]) through solving a optimization problem and at the same time the regularization parameter can be determined, and secondly the posterior measure is explored by an independence sampler using a proposal distribution derived from the diagonal map. Our numerical results in Section 4.2, 4.1 illustrate the high efficiency of the diagonal map-based independence sampler algorithm and its robustness with respect to various parameters, such as the noise level and the parameters in hyper-prior.
In summary, the main contributions of this paper are two-fold: (1) We propose a FTG hybrid prior to handle unknowns with some complex oscillation features and textural information that cannot be well modeled by TV-Gaussian hybrid prior. In particular, it is able to effectively reduces blocky effects. (2) We also provide an efficient independence sampler with a proposal distribution based on a diagonal transport map, which is more efficient than standard pCN. (3) By using FTG hybrid prior, transport map and hierarchical Bayesian which can choose regularization parameters flexibly in combination, we significantly improve the accuracy and efficiency of Bayesian inference for several linear inverse problems of varying dimension, involving the deconvolution problem, inverse source problem, limited computed tomography and image denoising. The rest of the paper is organized as follows. In Section 2, we introduce the fundamentals of Bayesian inference and the fractional order total variation Gaussian (FTG) prior. Section 3 first describes the hierarchical modeling, and then introduces the transport maps in the context of Bayesian inverse problems, explains how the diagonal transport maps can be constructed from reference measure and discusses a numerical solution for the optimization problem. And the independence sampler using a proposal distribution derived from a diagonal map is described in this section. Section 4 presents a range of numerical examples for the linear inverse problems from the one-dimensional deconvolution problem and inverse source identification problem to the two-dimensional limited computed tomography (CT) reconstruction in medical imaging and the image denoising. The paper ends with some conclusions in Section 5.

The FTG priors 2.1 The Bayesian framework and hybrid priors
We first give a brief introduction to the Bayesian inference framework for the infinite-dimensional linear inverse problems. According to (1.1), the likelihood function which is law of y conditional on u is L y (u) ∝ exp(−Φ y (u)), (2.1) where the data-misfit function where f Σ := Σ − 1 2 f 2 2 , for any positive symmetric matrix Σ, and the · 2 is the Euclidean norm. We then choose a prior probabilistic measure of u denoted by µ pr which is related to the posterior measure of u,denoted by µ y , through the Radon-Nikodym (R-N) derivative [32] dµ y dµ pr where Z = X exp(−Φ y (u))dµ pr (u) is a normalization constant, which can be interpreted as the Bayes' rule in the infinite-dimensional setting. The most popular prior µ pr is chosen as the Gaussian measure µ 0 , i.e. µ pr = µ 0 = N (0, C 0 ), a zero mean and covariance operator C 0 Gaussian measure on X. Note that the C 0 is symmetric positive and of trace class. Howerer, to better reconstruct the detail informations of unknows, the hybrid prior [43] is necessary which takes the form dµ pr dµ 0 (u) ∝ exp(−J(u)), (2.4) where J(u) is the additional prior (or the regularization) information on u. As a result, the R-N derivative of µ y with respect to µ 0 is In this case, it is well-defined on X for the posterior measure µ y associated with the linear forward model satisfied with certain conditions [7,32]. And in [43], a TV-Gaussian (TG) prior is presented that total variation (TV) term is used to deal with sharp jumps of the function, and the Gaussian measure is used as a reference measure so that it results in a well-defined posterior measure in the function space. However, in this paper we introduce the fractional order TG (FTG) prior to better address unknown function in the linear forward, especially in the image reconstruction.

The FTG prior
Here we provide the formulation of the FTG prior and we start with briefly reviewing some basic knowledge on fractional derivatives [31,16].
Definition 2.1 (Riemann-Liouville (RL) fractional derivative). Let n − 1 < α ≤ n, n ∈ R + , and [a, b] ⊂ R, then Riemann-Liouville fractional derivative of a function f are defined as follows: the left derivative where the Γ (·) is the standard Gamma function and the corresponding Riesz-RL (center) derivative is defined as Then the fractional order α total variation of a function ϕ is given by and the corresponding norm: where p is any positive integer, and Ω ⊆ R m . In [44], it has been proved that the fractional sobolev space W α p (Ω) is the Banach space, and has the following embedding relation that is, for any u ∈ W α 1 (Ω), there is a constant C > 0 such that In addition, [44] has showed that for any u ∈ W α 1 (Ω), the fractional order total variation T V α (u) = Ω |∇ α u|dx. And for any 1 ≤ p < ∞ fractional order Sobolev space W α p (Ω) is separable [40,41]. According to the definition of fractional Sobolev space, the norm of W α 2 (Ω) space is induced by the inner product below We can obtain that W α 2 (Ω) is the separable Hilbert space. Thus we choose the X = W α 2 (Ω) and consider the embedding relation (2.6). Then we can get the following α−order FTV regularization term or FTG prior where λ is a regularization parameter, and u T V α = Ω |∇ α u|dx. [40,41] have show that the equation (2.4) is a well-behaved prior under certain assumptions on the FTG prior J, and the posterior µ y is a well-defined probability measure on the X and it is also Lipschitz in the data y. Moreover, under some additional assumptions, the posterior measure can be well approximated by a measure defined in a finite-dimensional space.

Diagonal map-based independence sampler
As we know, the commonly used approach to exploring the posterior is the Markov chain Monte Carlo [9], which can generate a stream of samples to estimate some statistical information of complex target probability distributions. However, samples generated by MCMC are necessarily correlated among successive samples, which means smaller effective sample sizes and larger errors in posterior estimates. Thus the algorithm usually needs large number of samples to make accurate estimates, in particular, in the infinite-dimensional setting, which leads to bring a heavy computation burden. In this section, in order to increase the efficiency of independence sampler, we construct a fixed proposal distribution that approximates the posterior distribution at hand through a diagonal transport map that pushes forward a simple reference distribution to the target probability distribution. And, simultaneously the regularization parameters for the independence sampler can be determined automatically via (3.10) after constructing a transport map at hand. Moreover, we describe that diagonal map based independence sampler is well-defined in the function space.
We start by briefly reviewing the independence sampler MCMC. To this end, we define the measures on the product space X × X, where the µ is given a proposal distribution. When the ν is absolute continuous with respect to ν T [33,7], we then can define the acceptance probability . It can be seen that the acceptance probability of the independence sampler is well-defined if and only if the ν is absolute continuous with respect to ν T , which obviously need that µ y and µ are equivalent each other. And it suffices to require that µ and µ 0 are equivalent since µ y and µ 0 are equivalent. The dimensionindepandent pCN algorithm and independence sampler can be obtained by simply choosing the Gaussian prior as proposal distribution which is equivalent to the µ 0 and invariant [7].
The standard pCN algorithm used the proposal as follow where ω ∼ N (0, C 0 ) and the β plays the role of the locality parameter, and then v ∼ N ( 1 − β 2 u, β 2 C 0 ). And the associated acceptance probability is where the λ is regularization parameter. And the Algorithm 1 describes the standard pCN approach, which is used to compare the sampling with our independence sampler. However, this simply proposed distribution according to the prior works well only when the data and likelihood is not too informative. When the posterior is far away from the prior, and the Φ y varies extremely depending on where it is evaluated, the independence sampler will not work well. Next we introduce an efficient proposal measure derived from a diagonal map for the independence sampler.

Algorithm 1: The standard pCN
Input: Date-misfit function Φ y (u), prior J(u; λ) and µ 0 , regularization parameter λ, total number of steps K. Output: pCN samples of the target distribution, Calculate the acceptance probability a(u (i) , v) via the (3.2) ;

Construction of transport maps
Before describing the transport maps, we first introduce the hierarchical Bayesian formulation in order to make the most of it, which we can flexibly determine the regularization parameter.

Hierarchical Bayesian formulation
The probability density function corresponding to the target measure µ y is the posterior π y , and the probability density function corresponding to Gaussian measure µ 0 is denoted by π 0 . Then according to (2.5), we can obtain that We rewrite u ∈ R d to represent the unknowns and let Σ η = σ 2 I and the covariance operator C 0 ∈ R d×d , where σ is the standard deviation and I is a n-dimensional identity matrix. Then the posterior density π y can be obtained as The posterior density π y provides the complete distribution of u relying on the observations y. Following the same steps as in [6], we can compute the MAP point with the FTG prior, u M AP := arg max u π y (u), which is equivalent to the following minimization problem where the regularization parameter λ plays a critical role in classical inverse problems. It is essential to choose a suitable λ for Bayesian inverse problems. However, as we know, it is nontrivial task to ascertain the regularization parameter in almost inverse problems. Thanks to hierarchical Bayesian modeling [11,9], we can overcome the difficulty flexibly. The idea is to let the data y determine the parameters in the hope of effectively diminishing the effect of the initial (prior) assumptions of their values on the inverse solution. Then, the unknown function u and the regularization parameter λ can be identified at the same time. In hierarchical Bayesian framework, λ (a hyper-parameter) can be regarded as a random variable. If we choose Gamma distribution G(λ; k, ϑ) as the hyper-prior [10] for λ, i.e., where the positive constants k and ϑ is called shape parameter and rate parameter respectively. Then the posterior density can be written as and the MAP estimate for the posterior density (3.7) can be easily obtained by minimizing the following functional One salient feature of the functional is that the regularization parameter λ can be computed by data-driven procedure. And the minimization problem of (3.8) can be solved by the iteratively reweighted approach [19], which is usually applied in compressed sensing. However, we focus on the relationship between λ and u from the (3.8) rather than its minimizer. Now taking the partial derivative of R(u, λ) with respect to λ, we can obtain its second order partial derivative is ∂ 2 R(u,λ) After construction of a transport map, the above equation can be used to calculate the regularization parameter for our independence sampler.

Optimal transport
In [22], the measure-preserving transport maps constructed via the solution of an optimization problem between continuous probability measures is first applied to Bayesian inference. [22] presented a variational approach to the construction of transport maps explicitly that pushes forward the prior measure to the posterior measure, which entirely avoid Markov chain simulation. The transport approach of Benjamin et al. [25] and Parno et al. [24] instead follows a precondition MCMC sampling. [24] used transport map to obtain a proposal distributions that can more effectively explore the target density and is adapted as the MCMC sampling proceeds. And [25] proposed a multifidelity approach that the low-fidelity model is used to construct a transport map and the high-fidelity posterior distribution is explored using a non-Gaussian proposal distribution derived from the transport map. Below, we follow [21,25] to denote some notions of transport maps. We will refer to the posterior measure µ y and Gaussian measure µ ref as the target and reference measures on R d , respectively. A transport map T : R d → R d is a deterministic coupling that pushes forward µ ref to µ y , satisfying In other words, µ y (B) = µ ref (T −1 (B)) for any Borel set B ∈ R d . Both the reference measure µ ref and the target measure µ y are absolutely continuous with respect to the Lebesgue measure on R d that assure the existence of transport maps satisfying (3.11). Of course, there may be infinitely many such transport maps between the reference measure and the target measure. One way of guaranteeing the uniqueness of map is to introduce a transport cost function and minimizes it simultaneously satisfying the constraint (3.11). This minimization problem is called the Monge problem [35,37,36], and its solution is the optimal transport map. If the cost is taken to be a quadratic form in [2,3], the optimal transport map is exactly the Knothe-Rosenblatt rearrangement [2,3,29]. However, in this paper we directly assume that the transport map is a triangular diffeomorphism T such that ∇T 0 (i.e. monotone increasing) as in [24,22,21,25], instead of being particularly concerned with the optimality aspect of the transport.
As noted above, the lower triangular maps take the form where u ∈ R d and T k : R k → R (k = 1, 2, . . . , d) is kth component function of the transport map T . In this setting, the T k of the map depends only on the first i input variables, and it holds that det ∇T > 0 (see [21] for more details). Since the reference and the target measures are absolutely continuous, existence and uniqueness of such a lower triangular transport map (i.e. Knothe-Rosenblatt rearrangement) are guaranteed [2,3,29]. To obtain a useful approximation of the transport map, we will define a map-induced density π ref (u) and minimize the distance between the density π ref (u) of Gaussian measure µ ref and this mapinduced density. The next subsections describes the setup of this optimization problem.

Optimization problems
To set up the optimization problem, following [21], the density form of (3.11) can be written as T π ref = 1 τ π y , where τ is the normalization constant of posterior density and the transport map T only move mass of u but not including the λ. Now consider the pushforward of the reference density under the map T , and it is defined as where det ∇T −1 (u) denotes the determinant of the Jacobian of the inverse of the map at u. Then we can obtain the map-induced density where T is an approximation of a transport map T , which will be obtained via numerical optimization. If the reference density and the map-induced density are equal, i.e. π ref = π ref , the T can exactly satisfy T µ ref = µ y . Thus we can minimize a distance between π ref and π ref to obtain the T . In this paper, we use the Kullback-Leibler (KL) divergence to measure the distance between distributions as in [24,25]. Let the D KL (·||·) and E π ref [·] denote the KL divergence and integration with respect to the reference measure, respectively. Then a minimizer of the optimization problem: is a valid approximation of a transport map [22], where the constraint ∇ T 0 suffices to enforce mmonotonicity of a triangular map and T is some space of smooth lower triangular functions from R d to R d . If T is rich enough, we will obtain D KL (π ref || π ref ) = 0. Then the solution T of this optimization problem will satisfy (3.11) [21,22].
Furthermore, form the (3.14), the objective function of this optimization problem can be written as where the E π ref [log π ref + log τ ] is independent of the map T and thus a constant that can be ignored for the purposes of optimization.

Diagonal approximation of triangular map
To obtain the numerical solution of this optimization problem (3.15), the infinite-dimensional function space T must be replaced with a finite-dimensional subspace T ⊂ T . In [25], each component function of the approximation map T k is parameterized with the integrated-squared ansatz, which can enforce the monotonicity constraints explicitly and capture nonlinear dependencies in the target measure. However, this parameterization of the map is computationally expensive in the infinite dimensional setting. In this paper, the T k is parameterized by expanding it in a diagonal basis of univariate polynomials for the linear inverse problems.
Let each component of the map be written as T k (a k ; x), k = 1, 2, . . . , d, where a k ∈ R q+1 is a column vector coefficients. Then, we can express each component of the transport map T ∈ T as where x k is the kth component of x, and P q (a k ; x k ) is a q degree univariate polynomial with respect to x k , which forces T to be lower triangular. It is easy to see that the number of coefficients is proportional to that of the mode parameters, and as a result the diagonal approximation of transport map T k can be rapidly calculated by numerical optimization of (3.15). In fact, the diagonal parameterization (3.17) is the same to that using the multi-index sets J D k with p = q in [21].

Numerical optimization
There is need to approximate the expectation with respect to reference measure in the objective of (3.15) in the process of numerical optimization. We approximate the expectation E π ref [·] by its sample-average approximation (SAA) [17], i.e., a Monte Carlo estimator with M independent samples, denoted by can be calculated more accurately, as the cardinality of the sample set grows. Then, the coefficients (a 1 , a 2 , . . . , a d ) is denoted by F and we can obtain the optimization problem where we have simplified the monotonicity constraint ∇ T 0 by using the fact that ∇ T is lower triangular and according to the (3.10) we can explicitly obtain the λ in π y , ) via the SAA. In our work, we use the alternating direction algorithm to solve this optimization problem, which is summarized in Algorithm 2. Then, the above numerical optimization problem (3.18), in fact, just require to optimize the parameter F in T .
If the unknowns are only endowed with a Gaussian prior µ 0 for the linear forward model inverse problems, then the posterior is the Gaussian, i.e. u ∼ N (µ p , Σ p ), where where A ∈ R n×d is the discretization of the forward model A, and the transport map is linear and available in closed form: where Z 0 C 0 Z T 0 = Σ p and z 0 = µ p . The initial point of the optimization problem (3.18) is chosen as this linear map (3.20), i.e. F 0 = [z 0 Z 0 ] T for the alternating direction method in Algorithm 2. In fact, the optimization problem can become unconstrained depending on the parameterization of the transport map (e.g. integrated exponential parametrization or squared-integrated parameterization, see [21,25] for more details) that the constraint can be automatically satisfied.

Linear diagonal map-based independence sampler
We expect that the diagonal transport map approximately pushforward the reference onto the posterior, it is reasonable to consider an independence sampler for the diagonal transport map together with the reference Algorithm 2: Alternating direction method Input: Target density π y (x, λ), Gamma distribution parameters k and ϑ, samples Compute the parameter λ k+1 by Update the coefficients F k+1 by solving the optimization problem (3.18) for fixed λ k+1 ; 5 end 6 return Direct transport map T .
distribution serves as a proposal distribution. And the independence sampler with the proposal that is close to target space through a transport map can greatly reduce integrated autocorrelation time during sampling, which can clearly improve sampling efficiency.
In this paper, the approximation space for the map T is taken by the first-order polynomials (linear), i.e. q = 1 in (3.17), which is able to catch the main information of posterior and be rapidly calculated for the numerical optimization. Once we have a valid approximation of a transport map T between reference measure and target measure at hand, then the regularization parameters is determined by the (3.19). Moreover, the proposed distribution is set as the T µ ref , which is also a Gaussian measure, dented by N (m 1 , C 1 ), because of the linear diagonal map, and the µ ref can be taken by any Gaussian measures in the sampling process. Obviously, the T µ ref and µ 0 are equivalent (see [8] for more details). Further, if we also take the Gaussian measure N (m 1 , C 1 ) as our Gaussian prior µ 0 in the sampling process, then the proposal distribution T µ ref is apparently reversible, and hence invariant, with respect to this µ 0 [7]. Then, we can get for the independence sampler. In this work, the complete linear diagonal map-based independence sampler is summarized in Algorithm 3. And in Algorithm 3, we set the initial state u (1) = E π ref [ T (x)] that belongs to approximation target space, which can accelerate the convergence of the Markov chains. Although the regularization parameter λ is determined in advance, as will be demonstrated with our numerical results in Section 4, the additional FTG prior still has an outstanding performance in terms of handling the complex structure and detail information of model inputs.

Numerical examples
In this section, we present some examples and applications to demonstrates the performance of our FTG prior and the linear diagonal map-based independence sampler. These applications and examples for the linear inverse problems range from the one-dimensional deconvolution problem and inverse source identification problem to the two-dimensional limited computed tomography reconstruction in medical imaging and the image denoising.
We use the Grünwald formula to discretize the Riemann-Liouville fractional derivative. If 0 < α ≤ 1, i.e. n = 1, the Riemann-Liouville fractional derivative is approximated by the following standard Grünwald formula Pushforward proposed stater onto target, v = T (r); 6 Calculate the acceptance probability a(u (i) , v) using the (3.1): and if 1 < α ≤ 2, i.e. n = 2, the shifted Grünwald formula is used where l = 1, 2, . . . , d − 1, and α In addition, if no special explanation is given, we choose the following exponential covariance for the Gaussian prior µ 0 , where γ and ν are the parameters.
To obtain the transport map, in this paper, we use the interior point algorithm to solve the (3.18) for fixed λ in Algorithm 2. If no special explanation is given, the numerical optimization problem is performed with MATLAB's fmincon optimizer, where the step tolerance (StepTolerance) is set to 10 −6 ; the Specify-ConstraintGradient and SpecifyOdjectiveGradient are set to true (we provide the gradients for the solver in the linear inverse problems, see the [22] for more detail). And the reference distribution µ ref is set the same as the Gaussian distribution µ 0 in solving the numerical optimization process, and the chain was run for 1 × 10 5 steps for our diagonal map-based independent sampler in the Algorithm 3 and then all sample is used to calculate the posterior mean.
We calculate the average relative error (RelErr) to measure the difference between the reconstructed result (the posterior mean) and the target. For the posterior mean x and the target x 0 , RelErr is defined as This two-dimensional version of a model that occurs in optical imaging. In this application, g stands for the blurred image intensity, f represents light source intensity and the kernel A characterizes blurring effects that occur during image formation. In this example, we consider the one-dimensional version of above model, i.e. Ω = [l, r], and its kernel is where the ξ and δ are positive parameters. Obviously, we can use the standard numerical quadrature to get the accurate approximation of Af = g. Then the midpoint quadrature is applied to discretize the equation (4.5) and a discrete linear system Af = g is obtained, where the A = (a ij ) d×d and where h = (r − l)/d. And the observations are obtained by the y = Af + η, where the η is the Gaussian white noise.

Set up of inverse problems
In order to illustrate get the numerical results, we take ξ = 1/(δ √ 2π), δ = 0.02 and Ω = [0, 1]. We consider the light source and set d = 120. The measurement data are obtained by y = Af + η and η is assumed to be the Gaussian noise with zero mean and standard deviation 4.84 × 10 −2 , 9.7 × 10 −3 , 4.8 × 10 −3 , which corresponds to 5%, 1%, 0.05% noise, respectively, with respect to the maximum norm of the output Af . Note that the measurement data y is computed from a twice finer grid. The Gaussian prior µ 0 is taken by zero mean and set γ = 0.016, ν = 0.0003 and the shape and rate parameter of Gamma distribution are set to k = 2×10 3 , ϑ = 1, respectively. We then construct a linear diagonal transport map from the Gaussian distribution µ 0 to the posterior (3.7) using the alternating direction method in Algorithm 2. We use M = 1000 samples of the Gaussian distribution µ 0 to approximate the expected value via the SAA in the objective function (see Section 3.1.5). The transport map T and the regularization parameter λ are then used as precondition for the independent sampler as shown in Algorithm 3. And the proposal µ ref is the Gaussian distribution with zeros mean and the standard deviation 2 × 10 −3 in Algorithm 3.

Result
We first compare the performance of FTG prior to that of TG prior in two cases. For 0 < α ≤ 1, the posterior mean and absolute error with TG prior and with FTG prior are plotted in Figure 4.1 and the relative error RelErr are listed in Table 1. From this table and figure, in general, we can conclude that the reconstruction results using FTG prior gradually converges to that using TG prior if the fractional order α → 1 − . However, from the Table 1, our FTG prior with α = 0.95, 0.99 have better numerical results compared with the TG prior. And the reconstructed results with FTG prior for α = 0.9, 0.95 outperforms that with TG prior and can catch some details such as the corner points, seeing the absolute error curve in Figure 4.1. Therefore if the reconstructed target has much this details information, as discussed the example CT reconstruction in Section 4.3, the FTG prior with the α → 1 − can obtain better recovery result, compared with the TG prior. In addition, we can find that the FTG with α = 0.2, 0.5 have similar the RelErr value in Table 1. It seems that our Gaussian prior µ 0 is proper for this example in the process of solving numerical optimization. This is because the Gaussian prior µ 0 will play main role in the numerical results for the a small fractional order FTG prior.  The posterior mean and absolute error with FTG prior for 1 < α ≤ 2 are plotted in Figure 4.2 and the relative error RelErr are listed in Table 2. As we can see, from the Table 2, the TG prior yield a lower RelErr values compared to the FTG results. However, in the Figure 4.2, our FTG prior with α = 1.01, 1.05, 1.1 can eliminate well the staircase effect in term of a smooth part in the reconstruction target f (x) and also catch the piecewise constant structure. Thus, if the reconstructed target has much oscillation information, our FTG prior with α → 1 + can not only outperforms the TG, but also maintain an acceptable relative error.  Next, we study the influence of noise level on the numerical results. The posterior mean and posterior standard deviation for different noise level 5%, 1%, 0.5%, using the linear diagonal map-based independence sampler for FTG prior with α = 0.95 and TG prior, are plotted in Figure 4.3. Comparing the numerical results, it is observed that the posterior mean is able to approximate the exact f (x) well as the noise level decreases and at the same time, the posterior standard derivation gradually becomes smaller regardless of the TG or FTG prior. Finally, we study the sensitivity of the inversion results with respect to the parameter (k, θ) of Gamma distribution and list the corresponding numerical results in Table 3. We find that when fixing one of the parameter, changing the other parameter had little effect on the results. And according to the relationship (3.10), the parameter (k, θ) are depended on the value of the u T V α in some extent.

Problem setup
We consider the following initial-boundary value problem for the non-homogeneous heat equation.
where the ∆ is the Laplace operator, Ω is a bounded domain of R m , T > 0, and ψ and the source term f are given functions.
Our task is to determine the heat source on the usual initial-boundary conditions with the assistance of the final temperature data. In fact, this is inversely determined and usually ill-posed. In this paper the inverse problem of determining f (x) from the knowledge of ψ(x) and the final temperature measurement v(x, T) = φ(x), x ∈ Ω. (4.11) As shown in [13,12], when the data φ and ψ satisfy suitable conditions, the above linear inverse problem has a unique solution. Analogy to [42], we use the finite difference method (FDM) to solve the one-dimensional heat equation (4.8)−(4.10), and the equation (4.8) is discretized by using Crank-Nicolson method, i.e. using the forward-difference for the time derivative and a weighted combination of backward-difference and forwarddifference approximations for the remainder of the equation. Let Ω be the interval (0, r), and ∆t is the time step size and use the notation, v n (x, t n ) where t n = n∆t, n = 1, 2, . . . , N, t 0 = 0 and t N = T, then we can get the time-discrete equation where 0 ≤ w ≤ 1 is the weights. We use the uniform grid and let ∆x be the space step size, x j = j∆x, j = 0, 1, 2, . . . , d + 1 be a set of discrete points that defines the grid. The operator ∆ is discretized by using the second-order central difference and apply the initial conditions, and then we can obtain the fully discrete equation . . , f (x d ) T , and D + , D − are the discretization matrices of the 1/∆t − w∆, 1/∆t − (1 − w)∆, respectively. Both sides of the equation (4.13) multiply the inverse matrix of D + and recurse, which can yield where + and then we can get the matrix equation b = Hf . Thus, the solution of heat equation can be obtained by solve the matrix equation and at the same time we can get the observations through y = Hf + D N V 0 + η where the D N V 0 + H· is the discretization of our linear forward model A and the η is the noise.

Set up of inverse problems
In order to illustrate the numerical results, we take T = 1, r = 12, w = 0.5 and the initial temperature In this example, we consider the heat source: 10 ≤ x < 11.25, 0 otherwise, and set d = 150, N = 120. The measurement data are obtained by y = Hf + D N V 0 + η and η is assumed to be the Gaussian noise with zero mean and standard deviation 6.185 × 10 −3 , 6.185 × 10 −4 , 6.185 × 10 −5 , which corresponds to 1%, 0.1%, 0.01% noise, respectively, with respect to the maximum norm of the output Hf + D N V 0 . Note that the measurement data y is computed from a twice finer grid. The Gaussian prior µ 0 is taken by zero mean and set γ = 0.03, ν = 0.0009 and the shape and rate parameter about Gamma distribution are set to k = 1 × 10 4 , ϑ = 1, respectively. Then using the alternating direction method, we construct a linear diagonal transport map between the Gaussian distribution µ 0 and the posterior (3.7). Analogy to above deconvolution problem, we use M = 1000 samples of the Gaussian distribution µ 0 to approximate the expected value via the SAA in the objective function (see Section 3.1.5). And then the regularization parameter λ can be determined by the formula (3.10), which together with the diagonal map T is used as the precondition for the independence sampler as shown in Algorithm 3. The proposal µ ref is the Gaussian distribution with zero mean and the standard deviation 4 × 10 −3 in the sampling process.

Result
At first, we consider the robustness of the diagonal map-based independence sampler Algorithm 3 with respect to the noise level. The posterior mean and posterior standard deviation for different noise level 1%, 0.1%, 0.01%, using our independence sampler with FTG for α = 1.1 and TG prior, are plotted in Figure  4.4. Comparing the numerical results, it is observed that the posterior mean is more and more consistent with the exact f (x) as the noise level becomes smaller and simultaneously, the posterior standard derivation presents a stable trend and becomes smaller gradually whether it is for TG or FTG. In the following numerical results for this example, we fix the noise level as 0.1%. Similar to the above deconvolution example, below we compare the reconstruction performance in two cases. One situation is for the FTG prior with 1 < α ≤ 2. The posterior mean and absolute error using FTG with α = 1.05, 1.1, 1.5 and TG prior are plotted in Figure 4.5 and the relative error RelErr are listed in Table 4. As we can see from this tables, the RelErr values using our FTG prior for 1 < α ≤ 2 are consistent with the TG prior and are getting smaller and smaller with the fractional order α approaching to 1 + . And the FTG prior with α = 1.05, 1.1 has better performs than the TG prior in term of the complex oscillation structure of the heat source from the Figure 4.5. In this figure and the corresponding absolute error curve, our FTG prior with α = 1.05, 1.1 can eliminate well the staircase effect whether it is sharp points or a smooth part in the reconstruction target f (x) and at the same time also catch the piecewise constant part in it, compared with the numerical result using the TG prior. With the fractional order α → 1 + , FTG prior can not only reconstruct the complex oscillation structures of target well, but also maintain a relative error similar to that of TG prior, cf. the Table 4 and Figure 4.5. Therefore, obviously, the FTG prior with some appropriate α in (1,2] can obtain better recovery result than that using the TG prior, when the reconstructed target has much complex oscillation information, as discussed the image denoising example in Section 4.4.
Another situation is for the FTG prior with 0 < α ≤ 1. The posterior mean and absolute error using FTG prior with α = 0.5, 0.9, 0.95, 0.99 are plotted in Figure 4.6 and the relative error RelErr are listed in Table 5. From this table and figure, we conclude that the reconstruction results using FTG prior gradually converges to that using TG prior with the α → 1 − . In Table 5, we can find that the FTG with α = 0.2 has lower RelErr value than that with α = 0.5. It seems that the selected Gaussian prior µ 0 is appropriate for   this example in the process of solving the numerical optimization, which are consistent with the above the example. The Gaussian prior µ 0 will play main role in the numerical results for the a small fractional order FTG prior. And from this table, the RelErr values generated by the FTG prior are approaching to that of TG prior as the fractional order α tends to 1 − . However, from absolute error curve in the Figure 4.6, the reconstructed results with FTG prior for α = 0.99 outperforms that with TG prior in some details such as the corner points.   We now compare the efficiency performance of the linear diagonal map-based independence sampler in Algorithm 3 and standard pCN in Algorithm 1. In this example, the stepsize β is chosen so that the resulting acceptance probability is in the range 20% − 30% and all other parameters of pCN take the same values as our diagonal map-based independence sampler method. And we draw 10 6 samples from the posterior with the first 0.5 × 10 6 samples discarded as burn-in period for the standard pCN.
First, in Figure 4.7 we show the trace plots of the two methods for the unknown function f (x) at x = 5.9603 using the FTG with α = 0.9, 1, 1 and TG prior. High mixing rate in the plots indicate that successive iterations are highly independent and that the series of iterations have converged. The trace plots, in this respect, indicate that our diagonal map-based independence sampler method achieves a much faster mixing rate than the standard pCN. Next we compute the autocorrelation functions (ACF) of various quantities with the samples drawn by the two methods to further compare the efficiency. We consider the points at x = 0.0795, 5.9603, 11.9205 and plot the ACF results in Figure 4.8 using the FTG with α = 0.9, 1, 1 and TG prior. From this figure, the ACF for all three points of the chains generated by the diagonal-map based independence sampler decay faster than that corresponding to the standard pCN, which also suggests that the our method achieves a significantly better performance. The average acceptance rate of the independence sampler is about 77% and that of the standard pCN is about 31%.
Another common measure of the sampling efficiency of MCMC is the effective sample size (ESS) [15]. The ESS is computed by where K is the total sample size and is the integrated autocorrelation time, and an estimate of the number of effectively independent draws in the chain can be obtained by the ESS. We compute the ESS of the unknown function f (x) at each grid point through our diagonal map-based independence sampler method and show the results in Figure 4.9. The results show that our independence sampler algorithm on average produces significantly more effective independent samples, which agrees with faster decay of the ACF of this algorithm in Figure 4.8.
indep. sampler with a diagonal map standard pCN FTG with α = 0.9

Limited computed tomography reconstruction
We consider a classical inverse problem of X-ray computed tomography (CT), where X-rays travel from sources to detectors passing an object of interest. The intensities from multiple sources are measured at the detectors, the goal is to reconstruct the density of the object.

Problem setup
In this case, the forward operator is characterized by the Radon transform [26,23], which can model the attenuation of the x-ray when traveling from the source to the detector through the target, Shepp-Logan phantom, (cf. where the dx ⊥ is the Lebesgue measure along the line L and the f is called an attenuation coefficient function [23]. Then the result is the logarithm of the ratio of the intensity between the source, I 0 , and the detector, I d , through the Beer's law where the y is called sinogram, which is also called the measurement (1.1) and can take the form where η ∼ N (0, σ 2 I) is assumed to be the measurement errors and n is sum of the number of the x-ray in all discrete angles. If the domain of f is discretized into a grid and the density (i.e. the f ) is assumed to be constant within each grid, then the line integrals are approximated as where d is the total number of grid cell and the a ij is the length of the intersection between line L i and the jth grid cell. Then (4.17) and (4.18) yield a linear model y = Af + η, where the matrix A is defined by A = (a ij ) n×d .

Set up of inverse problems
In this case, the target square Shepp-Logan phantom f (x), defined on Ω = [−1, 1] 2 , is divided uniformly into 64 × 64 grid cell (or pixels), seeing the top left in Figure 4.11, i.e. d = 4096. We consider the measurement data y = Af +η and η adds Gaussian noise with zero mean and standard deviation 0.0115, which corresponds to 0.07% noise with respect to the maximum norm of the output Af . Note that the measurement data y is computed from a twice finer (128 × 128) grid Shepp-Logan phantom (at same measurement angles but finer arrangement of X-rays) and interpolated to lower resolution, while we use a 64 × 64 grid for inference. The Gaussian prior µ 0 is taken by zero mean and covariance matrix C 0 = 10 −5 I and the FTG prior with the fractional order α = 0.5, 0.8, 0.9, 1.1, 1.2, 1.5. The shape and rate parameter of Gamma distribution are taken by k = 2.55 × 10 6 , ϑ = 1, respectively. We then construct a transport map from the Gaussian distribution µ 0 to the posterior (3.7). The numerical optimization problem (3.18) is performed with MATLAB's fmincon optimizer, where the step tolerance (StepTolerance) is set to 10 −3 ; the SpecifyConstraintGradient and SpecifyOdjectiveGradient are set to true and we use M = 4096 samples of the Gaussian distribution µ 0 to approximate the expected value through the SAA in the objective function (see Section 3.1.5). The transport map and the regularization parameter via the (3.19) are then used to precondition for our independent sampler as shown in Algorithm 3. And the proposal µ ref is the Gaussian distribution with zero mean and the standard deviation 10 −5 in the Algorithm 3. We compare the posterior mean using our independent sampler with FTG to TG prior and the result of filtered back-projection (FBP) [23] inversion technique. The structural similarity (SSIM) [39], which is good at measuring quality of the reconstructed images in terms of image structure, is also adopted here.

Result
Top left in Figure 4.11 shows the target Shepp-Logan phantom 64×64 pixel and the area on the image where the red line passes are used to make a line plot in It can be seen that the reconstructed images with FTG prior for α = 1.1, 1.2, 1.5, in Figure 4.11 and the RelErr and SSIM values in Table 6, are consistent with the TG prior and outperforms results through the FBP reconstruction method. And the FTG prior with α = 1.1, 1.2, 1.5 can eliminate the blocky effect and be able to well remove under-sampling artifacts while preserving high resolution information, cf. the bottom row in Figure 4.12 and Figure 4.11. From the Table 6, the FTG prior with α = 0.5, 0.8, 0.9 can yield the lower RelErr value compared with TG and the FBP. Because the phantom of Shepp-Logan has much details information, and for the FTG prior with α → 1 the reconstruction results can better preserve this small details such as textural information, cf. the middle row in Generally speaking, the Shepp-Logan phantom image is of piecewise constant gray level, which is the case where the TG prior is most applicable. However, our FTG prior with α → 1 still has a good performance for the reconstruction result.

Image denoising
Image denoising is a simple yet heavily addressed problem in image processing. The noise in digital images can be caused by the failure or poor performance of image sensors, or failure of the data transmission process. The denoising process deals with the removal of noise from the noisy images considered.

Problem setup
The denoising problem is often modeled by a degradation model, which, in this paper, is given by where the f represents the original image, y stands for a degraded/noisy image (given data) and the random variable η is assumed a Gaussian white noise with the known standard deviation. For a given degraded image y, the image denoising is to reconstruct the original noise-free image f . In fact, the degradation model (4.19) is equivalent to the model (1.1) for the A as an identity operator.

Set up the inverse problrm
In this example, the original image−cameraman f (x), defined on Ω = [−1, 1] 2 , is choosing 128×128 pixel, i.e. the dimension d = 1.6384×10 4 . The noisy image y is taken as the original image corrupted by Gaussian noise with zero mean and standard deviation 0.03, which corresponds to 3% noise with respect to the maximum norm of the original image, seeing the top right of Figure 4.13. The Gaussian prior µ 0 is taken to be the standard Gaussian distribution and the FTG prior with the fractional order α = 0.5, 0.9, 1.1, 1.5, 1.8. The shape and rate parameter about Gamma distribution are set to k = 3 × 10 4 , ϑ = 1, respectively. Then, using the alternating direction method we construct a transport map between the standard Gaussian distribution µ 0 and the posterior (3.7). The numerical optimization problem (3.18), analogy to above CT reconstruction, is performed with MATLAB's fmincon optimizer, where the the step tolerance (StepTolerance) is set to 10 −3 ; the SpecifyConstraintGradient and SpecifyOdjectiveGradient are set to true and we use M = 2000 samples of the standard Gaussian distribution µ 0 to approximate the expected value via the SAA in the objective function (see Section 3.1.5). The transport map T and the regularization parameter λ via the formula (3.19) are then used as precondition for the linear diagonal map-based independence sampler. And the proposal µ ref is taken to be a Gaussian distribution with zero mean and the standard deviation 0.01 in the sampling process. The Peak signal-to-noise ratio (PSNR), which does well in measuring quality of the reconstructed images in terms of image denoising, is also adopted in this case.

Result
Top left in Figure 4.13 shows the original image, the Cameraman, 128 × 128 pixel and the area on the image where the red line passes are used to make a line plot in Figure 4.14. The middle and bottom row in Figure  4.11 are our reconstructed images with the FTG prior for the fractional order α = 0.5, 0.9, 1.1, 1.5, 1.8 and TG using 3% noise with respect to the maximum norm of the original image. Figure 4.14 shows the reconstruction images for the lines plot in term of the red line areas in original image (see top left in Figure  4.13).
It can be seen that the reconstructed images with FTG prior in Figure 4.11 and the PSNR and SSIM values in Table 7 are in agreement with the TG prior. Of course, the TG prior yields a better result in term of the PSNR value from the Table 7. Because the original image has some edges. However, from Table 7, the FTG prior with α = 1.1, 1.5, 1.8 yield the higher SSIM value compared with the TG prior. This is because the original image also has some texture generating the nonconstant gray level in the domain, which leads to significant blocky effect for the reconstructed result with TG prior, cf. the top left in Figure 4.14, while the FTG prior can eliminate well the blocky effect caused by TG and better preserve some small details such as textural information and corner points. Thus, for the target image with some small details, our FTG with some appropriate α in (1,2] prior not only maintains a low relative error, which is consistent with that of TG, but also has a better reconstructed result in terms of the structure.   Figure 4.13). The corresponding reconstruction results using the FTG with the fractional order α = 0.5, 0.9, 1.1, 1.5, 1.8 and TG prior.

Conclusions
In this work, we have presented a FTG prior for infinite-dimensional Bayesian inverse problems. We use the FTV term to improve the ability to capture the detail information and use the Gaussian reference measure to ensure that it results in a well defined posterior measure. And the hierarchical Bayesian framework is also applied here, where the regularization parameter can be flexibly determined. Moreover, we also propose an efficient diagonal map-based independence sampler for the linear inverse problems in infinitedimensional setting. This sampler has two stages: firstly we construct a diagonal map that can approximately pushforwards the reference measure to the posterior measure and secondly the posterior measure is explored by an independence sampler using a proposal distribution derived from the diagonal map. Finally, we provide some numerical examples to demonstrate the performance of the FTG prior and the efficiency and robustness of the proposed independence sampler method. We find that the FTG prior has better performance than the TG for the detail information in the unknowns, especially for recovering textures of image. A natural extension of the present work is to use the FTG prior in other applications such as the reconstruction of rough surfaces. And the diagonal map-based preconditioner for independence sampler can be applied to the nonlinear infinite-dimensional inverse problems.