{'fraction_non_alphanumeric': 0.0492557225108833, 'fraction_numerical': 0.015693020643168094, 'mean_word_length': 4.2027397260273975, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical solution of an inverse random source problem for the time fractional diffusion equation via PhaseLift

This paper is concerned with the inverse random source problem for a stochastic time fractional diffusion equation, where the source is assumed to be driven by a Gaussian random field. The direct problem is shown to be well-posed by examining the well-posedness and regularity of the solution for the equivalent stochastic two-point boundary value problem in the frequency domain. For the inverse problem, the Fourier modulus of the diffusion coefficient of the random source is proved to be uniquely determined by the variance of the Fourier transform of the boundary data. As a phase retrieval for the inverse problem, the PhaseLift method with random masks is applied to recover the diffusion coefficient from its Fourier modulus. Numerical experiments are reported to demonstrate the effectiveness of the proposed method.


Introduction
In the past two decades, the differential equations involving fractional-order derivatives, known as the fractional differential equations (FDEs), have received increasing attention in applied disciplines. Such models are able to capture more faithfully the dynamics of anomalous diffusion processes in amorphous materials. Consequently, fundamentally different physics can be obtained. For instance, the FDEs can be used to model some anomalous diffusions in a highly heterogeneous aquifer [1], underground environmental problems [23], the relaxation phenomena in complex viscoelastic materials [22], and non-Markovian diffusion processes with memory [34].
Motivated by the significant applications in scientific and industrial fields, the research of inverse problems has gone tremendous developments in the past decade. Recently, the inverse problems on FDEs have become an active research field. In particular, the inverse problems for the time fractional diffusion equations have been extensively studied mathematically and numerically. Generally, the inverse source problems of FDEs are to determine the time-dependent or the spacedependent source functions by using the space-dependent or time-dependent data [2,4,24,48]. There are also cases where boundary conditions are used as the data [36,47]. Although there has been a lot of research done for the inverse problems of FDEs, there is little work on the stochastic inverse problems for the FDEs.
The inverse random source problems refer to the inverse source problems that involve uncertainties. Due to the randomness and uncertainty, compared to the deterministic counterparts, stochastic inverse problems have more difficulties in addition to the existing obstacle of ill-posedness. For the inverse random source scattering problem, where the wave propagation is governed by the stochastic Helmholtz equation driven by a white noise, it is shown in [15] that the correlation of the random source could be determined uniquely by the correlation of the random wave field. Effective computational methods are developed in [5][6][7][8][28][29][30] for the white noise model, where statistical properties, such as the mean and the variance, of the random source are recovered by using the boundary measurements of the random wave field at multiple frequencies.
There is much less work on the inverse random source problems for the FDEs. Consider the stochastic time fractional diffusion equation where ∂ α t denotes the Caputo fractional derivative with α ∈ (0, 1), B H is the fractional temporal Brownian motion with Hurst parameter H ∈ (0, 1) andḂ H denotes the formal derivative of B H with respect to the time t, D is a bounded domain with the Lipschitz boundary ∂D, and f, g are two deterministic functions supported in D. Given f and g, the direct source problem is to determine u; while the inverse source problem is to determine the unknowns f and g that generate a prescribed u. For this kind of random sources perturbed by a time-dependent random noise, the mild solution of the problem can be obtained by using the Mittag-Leffler functions, and it may lead to the reconstruction formulas between the unknowns and the measurements. In [39], f and |g| are proved to be uniquely determined by the mean and covariance of the final data u(x, T ), respectively, under the conditions that α ∈ ( 1 2 , 1) and H = 1 2 , i.e., the white noise case. It is also pointed out that the inverse problem is not stable as a small variance of the data may lead to a huge error on the reconstructions. The inverse problem for a generalized case with α ∈ (0, 1) and H ∈ (0, 1) is considered in [17], where f and |g| are recovered from the same statistics of the final data u(x, T ). Similar results are obtained about the uniqueness and the stability of the inverse problem. If the random sources are perturbed by a space-dependent noise, the mild solution approach is not available anymore since the spatial noise may not be regular enough to guarantee the well-posedness of the problem. New techniques need to be explored to reconstruct the diffusion coefficient of this kind of random sources.
In this paper, we consider the one-dimensional stochastic time fractional diffusion equation where W x is the spatial Brownian motion,Ẇ x denotes the formal derivative with respect to x and is known as the white noise, F (t) is a deterministic function and denotes the diffusion coefficient of the random source, and the Caputo fractional derivative ∂ α t u, α ∈ (0, 1), is defined by To deal with the initial boundary value problem of the stochastic diffusion equation (1.1), we consider an equivalent stochastic two-point boundary value problem in the frequency domain, and study the well-posedness and regularity of the solution based on the estimate of the Green function. It then leads to the well-posedness of (1.1) by showing that the solution to (1.1) is the inverse Fourier transform of the solution to the equivalent problem in the frequency domain. For the inverse problem, the Fourier modulus |F | of the diffusion coefficient F is proved to be uniquely determined by the variance of the solution to the equivalent problem in the frequency domain. However, the recovery of F from |F | is apparently not unique, since we measure |F | instead ofF and lose information about the phase of F . If we could retrieve the phase of F , then it would be trivial to recover F . This kind of problem, i.e., the reconstruction of a signal from the magnitude of its Fourier transform, is generally known as the phase retrieval [40,41]. It arises in many applications such as diffraction imaging, optics and quantum mechanics, and is usually ill-posed and notoriously difficult to solve. A large amount of methods have been proposed to solve the phase retrieval problem [25]. These approaches can be broadly classified into two categories: utilizing either a priori information about the signal F or additional measurements of the modulus |F |. In this work, we adopt the latter and apply the PhaseLift method with random masks to collect additional measurements of the modulus |F |, which may be used to uniquely determine |F |, which is the modulus of the diffusion coefficient of the random source. Numerical experiments are reported to demonstrate the effectiveness of the proposed method.
The rest of the paper is organized as follows. In Section 2, the time fractional derivative and its properties are introduced. Section 3 is concerned with the well-posedness of the direct problem. The inverse problem is addressed in Section 4, which includes the recovery of |F |, the phase retrieval problem, and its numerical method. Section 5 presents two numerical examples to illustrate the effectiveness of the proposed method. The paper is concluded with some general remarks and directions for the future research in Section 6.

The time fractional derivative
In this section, we discuss the Fourier transform of the Caputo fractional derivative, which serves as a basis to convert the initial boundary value problem of the stochastic time fractional differential equation (1.1) into an equivalent stochastic two-point boundary value problem in the frequency domain.
Then the Fourier transform of the fractional derivatives of v satisfies denotes the Fourier transform of v.
Proof. The result is obvious for the case α = 1 if v ∈ H 1 (R). Next we show the assertion for the case α ∈ (0, 1). Define a causal function k α + (t) with α ∈ (0, 1) by Based on the fact In fact, let U R ⊂ C be a simply connected open set with γ R being a closed curve shown in Figure 1. It is clear to note that the mapping z → e −iωz z −α , z ∈ C defines a holomorphic function in U R . It follows from the Cauchy integral theorem that where I R := {z ∈ C : |z| = R} and I 1/R := {z ∈ C : |z| = 1/R} denote the positively oriented a quarter section of the circle with radius R (the dashed curve in Figure 1) and the negatively oriented a quarter section of the circle with radius 1/R (the solid curve in Figure 1), respectively.
Taking the limit of (2.2) as R → +∞ leads to which completes the proof by noting (2.1).
Remark 2.2. Note that (iω) α is a multi-valued function when α is a fractional number. Throughout the paper, we define where sgn(·) denotes the sign function.

The direct problem
In this section, we discuss the direct problem. The well-posedness of the problem (1.1) and the regularity of its mild solution are investigated by studying an equivalent problem in the frequency domain.
3.1. The direct problem in the frequency domain. Since the function F satisfies F (0) = 0, we denote byF the zero extension of F in (−∞, 0) and byF the Fourier transform ofF . Consider the two-point boundary value problem of the stochastic differential equation In the following, we deduce the Green function and present the well-posedness of (3.1).
3.1.1. Green's function. Let s := (iω) α and denote by g ω (x, y) the Green function of (3.1). We consider two cases where g ω takes two different expressions.
If ω = 0, then g ω satisfies where δ is the Dirac delta function. It is known from solving the second order ordinary differential equation with constant coefficients that g ω (x, y) has the general form where A i and B i , i = 1, 2, are to be determined. Using the boundary conditions and the continuity and jump conditions we may easily obtain from Remark 2.2 that If ω = 0, then g 0 (x, y) solves Similarly, we may solve the above two-point boundary value problem and obtain Lemma 3.1. The Green function g ω given in (3.2) and (3.3) satisfies the estimate where C > 0 is a constant independent of ω.
Proof. If ω = 0, then it follows from a simple calculation that If ω = 0, we get where the function h is defined by .
Clearly, h is a nonnegative function for all k > 0. We claim that h(k) is uniformly bounded for all k > 0. On the one hand, due to the fact there exists a constant C 0 such that h(k) < 1 for all k > C 0 , which shows that h is uniformly bounded on (C 0 , ∞). On the other hand, by noting we get that h is also uniformly bounded on (0, C 0 ] due to its smoothness, which completes the proof of the claim. Using the estimates for h, we have from (3.4) that there exists a constant C > 0 independent of ω such that which completes the proof.

3.1.2.
Well-posedness. Based on the Green function g ω , we are able to show that the two-point boundary value problem of the stochastic differential equation (3.1) admits a unique mild solution.
The following result gives the estimate of the mild solution, which plays an important role in the analysis of the time domain problem.
Then the the stochastic differential equation (3.1) admits a unique mild solution given by Moreover, the solution U satisfies the estimate where C > 0 is a constant.
Proof. The existence and uniqueness of the mild solution can be proved similarly as those in [5]. We only show the proof for the estimate (3.6). By Itô's isometry, Lemma 3.1 and Parseval's identity, we get , which completes the proof.
Hereafter, the notation a b stands for a ≤ Cb, where C is a positive constant whose value is not required but should be clear from the context.
Then the initial boundary value problem of the stochastic differential equation (1.1) admits a unique solution u satisfying x ∈ D, t ∈ R, where F −1 denotes the inverse Fourier transform and U is the mild solution of (3.1). Define To show the existence of the solution to (1.1), we prove that the function u defined above is a solution of (1.1). Note thatũ whereF is the zero extension of F on (−∞, 0) which is defined at the beginning of this section. Hence,ũ is a causal function withũ(x, t) = 0 if t ≤ 0, which implies that Moreover, it follows from Parseval's identity and Theorem 3.
which indicates that the Caputo fractional derivative ofũ with respect to t is well-defined. Taking the inverse Fourier transform with respect to ω on both sides of (3.1), we have where we have used Lemma 2.1 for the causal functionũ and the fact u =ũ| R + . The equation (3.10), together with the initial condition (3.8), indicates that u defined in (3.7) is a solution of (1.1) satisfying the estimate (3.9). The uniqueness of the solution of (1.1) is obtained directly by the equivalence between the time domain problem (1.1) and the frequency domain problem (3.1), as well as the well-posedness of (3.1) given in Theorem 3.2.

The inverse problem
In this section, we address the inverse problem which is to reconstruct the diffusion coefficient F of the source from the measured wave field u(0, t) at the point x = 0 for t > 0.
where we have from (3.2) and (3.3) that if ω = 0 and g 0 (0, y) = y − 1, y ∈ D if ω = 0. For where we have used the fact sin(x) x ≤ 1 for any x ∈ R, and the function l 1 is defined by It then suffices to show that l 1 (k) > 0 for all k > 0. Note that It is easy to check that l 2 (k) > 0 and hence l ′ 1 (k) > 0 for all k > 0. In fact, by noting that Hence the function l 1 is increasing for all k > 0 and satisfies which completes the proof. Proof. It follows from (4.1) and Lemma 4.1 that which implies the uniqueness of the reconstruction.

4.1.2.
Instability. By Theorem 4.2, the inverse problem admits a unique solution but it lacks stability due to the fact that the denominator D |g ω (0, y)| 2 dy goes to zero as ω → ∞, which is stated in the following theorem. Proof. By (4.2), we have Following from the same estimate as that of (3.4), we get which completes the proof.

Phase retrieval.
In this section, we discuss the phase retrieval problem. More precisely, we aim to numerically reconstruct |F | from |F |, where the former is the modulus of the diffusion coefficient F of the source in the time domain and the latter is the modulus of the diffusion coefficientF in the frequency domain.
To introduce the numerical methods for the phase retrieval, we begin with presenting a discrete version of the phase retrieval problem, and then introduce an additional measurement based framework named PhaseLift [11], which is adopted to solve the phase retrieval problem. 4.2.1. Discrete phase retrieval problem. Let x = (x 1 , ..., x N ) ⊤ ∈ C N be a signal of length N , and y = (y 1 , ..., y N ) ⊤ ∈ C N be its N -point discrete Fourier transform (DFT). Denote by f (m) the conjugate of the mth column of the N -point DFT matrix, i.e., Then it is easy to check that y m = f (m) , x , where ·, · is the complex inner product defined by x, y := N n=1 x n y n .
The discrete phase retrieval problem is formulated as follows: find x subject to Pioneered by Gerchberg in [21], earlier approaches for the phase retrieval problem are based on alternating projections and can be reformulated as the following least-squares problem: The algorithm requires oversampling by using an M -point DFT with M > N . It attempts to minimize the above non-convex objective by starting with a random initialization and iteratively imposing the time domain and Fourier magnitude constrains using projections. However, since the projections are taken between a convex set and a non-convex set, the solution usually converges to a local minimum, which leads to the limited recovery ability of the algorithm even in the deterministic setting. Recent frameworks to attack the ill-posedness of the phase retrieval problem can be broadly classified into two categories: (1) developing a modified model with a prior information; (2) taking more magnitude measurements. The former aims to reduce the number of unknowns by assuming some a prior information of the signal, such as the support constraints [14,18,32], positivity and real-valuedness [18,19], or sparsity [26,35,44]. Depending on the applications, the latter can be done in various ways, which include the use of masks [27], optical gratings [31], oblique illuminations [16], or short-time Fourier transform magnitude measurements which utilize overlap between adjacent short-time sections [43,45].

4.2.2.
PhaseLift. Based on the semi-definite programming method, the PhaseLift is an effective approach to solve the phase retrieval problem. It has been shown that the PhaseLift may yield robust solutions to various quadratic-constraints. Since phase retrieval results are in quadratic constraints, the PhaseLift is adopted to handle our inverse problem.
There are two main ingredients in this method: multiple structured illuminations and lifting. For the self-contained purpose, we briefly introduce the two components of this method. The details may be found in [11].
First we comment on the multiple structured illuminations. Let x = (x 1 , ..., x N ) ⊤ ∈ C N be the object of interest, and assume that the illumination schemes, which collect the diffraction pattern of the modulated object {w n x n } n=1,··· ,N , are available, where the pattern w = (w 1 , ..., w N ) ⊤ ∈ C N may be selected by the user. There are several ways to implement the recovery in practice, such as masking, optical grating, and oblique illuminations. It is usually preferred to adopt the approach in which fewer diffraction patterns are required for a stable recovery.
Next we present the lifting. Based on the diffraction pattern w mentioned above, suppose that we have quadratic measurements of the form where a (m) can be chosen based on the diffraction pattern w in the following form: is defined in (4.3). The phase retrieval problem is then transformed to the feasibility problem: The core idea of the PhaseLift is to embed the signal x into a higher dimensional space by using the transformation X = xx * . Next, we lift up and interpret the signal x to this rank-one matrix X. Denoting A (m) := a (m) (a (m) ) * , we get where Tr(·) denotes the trace of a matrix. Let A be the linear operator mapping any positive semi-definite matrices X, denoted by X 0, into A(X) := (Tr(A (1) X), · · · , Tr(A (N ) X)) ⊤ . (4.6) Based on the above notation, the phase retrieval problem is equivalent to find X subject to A(X) = z, X 0, rank(X) = 1, where z = (z 1 , · · · , z N ) ⊤ is defined in (4.4). Consequently, the phase retrieval reduces to finding a rank-one positive semi-definite matrix X which satisfies these affine measurement constraints. Equivalently, the phase retrieval problem can be formulated to minimize rank(X) subject to A(X) = z, X 0. (4.8) The equivalence between (4.7) and (4.8) is obvious since b = A(xx * ) according to (4.5) and it has a rank-one solution. After solving (4.8), we factorize the rank-one solution X as xx * , which then leads to the solution of the phase retrieval problem. Therefore, our inverse problem is equivalent to a rank-minimization problem over an affine slice of the positive semi-definite cone. Furthermore, it is reduced to a problem of low-rank matrix completion or matrix recovery, which is a classical optimization problem that has gained tremendous attention in recent years [12,13].
However, it is known that the rank-minimization problem (4.8) is NP-hard. The trace norm is thus explored as a convex surrogate [9,33] for the rank functional in (4.8), which gives the following semi-definite programming problem: minimize Tr(X) subject to A(X) = z, X is Hermitian positive semi-definite. (4.9) The semi-definite programming problem (4.9) is convex, and there exists a wide choice of numerical solvers including the popular Nesterov accelerated first-order method [38]. As far as the relationship between (4.8) and (4.9) is concerned, it is beyond the scope of this work and we refer to [12,13] for the detailed discussion.

Numerical method.
In this work, all the numerical algorithms are implemented in MAT-LAB by modifying the templates for first-order conic solvers (TFOCS) [10]. The TFOCS is a library of MATLAB files which are designed to facilitate the construction of optimal first-order methods for a variety of convex optimization problems including the semi-definite programming problem (4.9) considered in this paper.
To illustrate how we actually handle (4.9), we next briefly introduce the formulation and implementation of a class of optimal first-order methods to solve the following general convex optimization problem: where g : R N → R is convex and smooth and h : R N → R is convex. Consider a class of first-order methods to solve (4.10) based on iterations with a generalized projection (cf. [10]): where · is a chosen norm and t k is the step size. The global convergence of (4.11) is guaranteed if t k is bounded away from zero and g(x k+1 ) has the following upper bound which can be accomplished by assuming that ∇g satisfies a generalized Lipschitz criterion for any x, y belonging to the domain of φ. Here, · * denotes the dual norm of the norm · defined by Then the bound of (4.12) is assured for any t k ≤ L −1 under the assumption (4.13), which leads to the convergence in O(L/ǫ) iterations for a simple algorithm, known as the forward-backward algorithm or proximal gradient descent [20], based on (4.11).
Optimal or accelerated first-order methods are able to improve the bound of number of iterations to O( L/ǫ), and have been studied by many researchers in the past decades (cf. [37,46]). The TFOCS implements a variety of the optimal first-order variants based on a variation of the method described by Nesterov in [38]. One variant, the Auslender and Teboulle method [3], is described as follows: where x 0 is chosen in the domain of φ, x 0 = x 0 , and θ 0 = 1. Here, θ k is usually referred to as the accelerated parameter.

Numerical Experiments
In this section, we discuss the implementation for solving the direct and inverse random source problems, and report two numerical examples to demonstrate the effectiveness of the proposed methods.

5.1.
Numerical discretization of the direct problem. To avoid the inverse crime, we employ the finite difference method to discretize the initial boundary value problem of the stochastic differential equation Define the partition of the time and spatial intervals with nodes t n = nh t , n = 0, 1, ..., N t , where h t = T /N t and h x = 1/N x , respectively. Let u n i be the numerical approximation to u(x i , t n ). The fractional derivative ∂ α t u at (x i , t n ) is approximated by The second order derivative ∂ xx u at (x i , t n ) is approximated by using the central difference method where a d = b means that the random variables a and b have the same distribution and {ξ i } Nx−1 i=0 is a set of independent and identically distributed standard normal random variables, denoted by ξ i ∼ N (0, 1).
Combining (5.1) with (5.2)-(5.3), we conclude that the numerical scheme has the following compact matrix form:

5.2.
Implementation of the PhaseLift. For the problem (1.1), the modulus of the Fourier coefficient |F | can be reconstructed from (4.1) by utilizing the variance V[U (0, ω)] of the Fourier transform of the field u(0, t). To solve the phase retrieval problem numerically, i.e., to uniquely recover the source f := (F (t 1 ), · · · , F (t Nt )) ⊤ at discrete temporal nodes, more measurements are usually required and the PhaseLift with masks introduced in Section 4.2 is used. Let M 1 and M 2 be two masks with M 1 being the N t × N t identity matrix and M 2 being a random binary mask, i.e., an N t × N t diagonal matrix whose entries are randomly chosen as 0 or 1. Using these two masks to perform the structured illumination with sources M i f , we can obtain in total 2N t intensity measurements A(f f * ), where the operator A is defined in (4.6) with N = 2N t , the vectors {a (m) } 2Nt m=1 involved in the operator A are defined by and f (m) is given in (4.3).
To recover the source f , it then suffices to solve the convex optimization problem introduced in Section 4.3 based on the objective functional where · 2 denotes the Euclidean norm and h is the indicator function defined by Here, b is the quadratic Fourier modulus obtained by the measurement {u n 0 } 2Nt n=1 at x = 0 perturbed by the sources where (U 1 0 , · · · , U 2Nt 0 ) ⊤ is the 2N t -point DFT of the measurement (u 1 0 , · · · , u 2Nt 0 ) ⊤ and ω n = 2πnh t /N t . The Auslender and Teboulle method (4.14) is applied to solve the above problem, which is terminated when the relative residual error of our reconstructed resultf is less than a fixed tolerance, i.e., ||A(ff * ) − b|| 2 ≤ 10 −6 ||b|| 2 . Note that the solutionf = (f 1 , · · · ,f Nt ) ⊤ is unique only up to the global phase, i.e., (|cf 1 | 2 , · · · , |cf Nt | 2 ) ⊤ = (|F (t 1 )| 2 , · · · , |F (t Nt )| 2 ) ⊤ is unique, where c is a complex scalar satisfying |c| = 1. In particular, if F (t) is a real-valued nonnegative function, then the recovery is unique. As a result, the absolute value |F | is recovered in the following two numerical examples.
We set N t = 65 and N x = 100, and use the total number of 1000 sample paths to approximate the variance of the solution. In addition, the data is assumed to be polluted by a uniformly distributed noise with the noise level σ = 0.05. In the first example, we set the source function F (t) = sin(t) exp(−t/6) and T = 4π, and in the second example, the source function is taken to be F (t) = sin(2t) cos(3t) with T = π. We present the results for two different cases: α = 0.4 and α = 0.8 in both of the examples, which are shown in Figures 2 and 3, respectively. We comment that (1) the larger the α is, the more regular of the solution is for the direct source problem; (2) the smaller the α is, the more stable of the solution for the inverse source problem. It can be observed from the numerical experiments that the proposed method is effective and robust to reconstruct the modulus of the source functions.

Conclusion
In this paper, we have studied an inverse random source problem for the stochastic time fractional diffusion equation driven by a spatial Gaussian random field. By examining the equivalent two-point stochastic boundary value problem in the frequency domain, we show that the direct source problem is well-posed and the inverse source problem has a unique solution. The ill-posed nature is revealed for the inverse problem, i.e., the modulus of the Fourier coefficient of the source function decays in an α order of the frequency. The inverse problem is then converted into a phase retrieval problem which is to recover the original signal from its Fourier modulus and is implemented via the PhaseLift with random masks. The numerical results show that the method is effective to reconstruct the modulus of the source functions.
The proposed approach based on the Fourier transform can be naturally extended to solve higher dimensional problems. For more complex cases, there are still several interesting problems to be investigated along this line of research. First, the source can be modeled by more general random processes such as spatial fractional Brownian motions. The analysis and results would remain the same for the direct problem. However, the uniqueness of the inverse problem may not be guaranteed due to the correlated kernel of fractional Brownian motions. Second, the Fourier transform requires more constraints about the functions u and f , while the Laplace transform demands less constraints on these functions. Nevertheless, there is no efficient algorithm to deal with the phase retrieval problem on the Laplace transform. In addition, when it comes to the situation that α ∈ (1, 2), more initial conditions are required to ensure the well-posedness of the model equation, which leads to the fact that the transform may not apply in this case. We hope to report the progress on these problems elsewhere in the future.