Numerical solution of an inverse random source problem for the time fractional diffusion equation via PhaseLift

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 [21], the relaxation phenomena in complex viscoelastic materials [20], and non-Markovian diffusion processes with memory [30].

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 [6, 11, 12]. 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 space-dependent source functions by using the space-dependent or time-dependent data [2, 4, 22, 39]. There are also cases where boundary conditions are used as the data [31, 38]. 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 [14] 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, 26] 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

$$\begin{equation*}\begin{cases}_{t}^{\alpha }u\left(x,t\right)-{\Delta}u\left(x,t\right)=f\left(x\right)h\left(t\right)+g\left(x\right){\dot {B}}^{H}\left(t\right),\quad \hfill & \left(x,t\right)\in D{\times}\left(0,T\right),\hfill \\ u\left(x,t\right)=0,\quad \hfill & \left(x,t\right)\in \partial D{\times}\left[0,T\right],\hfill \\ u\left(x,0\right)=0,\quad \hfill & x\in \bar{D}\hfill \end{cases}\end{equation*}$$

where ${\partial }_{t}^{\alpha }$ denotes the Caputo fractional derivative with α ∈ (0, 1), B^H is the fractional temporal Brownian motion with Hurst parameter H ∈ (0, 1) and ${\dot {B}}^{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 [33], 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 $\alpha \in \left(\frac{1}{2},1\right)$ and $H=\frac{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 [16], 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. A related inverse random source problem on the stochastic fractional diffusion equation for the time-dependent noise may be found in [27]. 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

$$\begin{equation}\begin{cases}_{t}^{\alpha }u\left(x,t\right)-{\partial }_{xx}u\left(x,t\right)=F\left(t\right)\dot {{W}_{x}},\quad \hfill & \left(x,t\right)\in \left(0,1\right){\times}{\mathbb{R}}_{+},\hfill \\ u\left(x,0\right)=0,\quad \hfill & x\in \left[0,1\right],\hfill \\ {\partial }_{x}u\left(0,t\right)=0,\quad u\left(1,t\right)=0,\quad \hfill & t\in {\mathbb{R}}_{+},\hfill \end{cases}\end{equation} \tag{ 1.1 }$$

where F is a deterministic function satisfying F(0) = 0 and denotes the diffusion coefficient of the random source, and the Caputo fractional derivative ${\partial }_{t}^{\alpha }u$, α ∈ (0, 1), is defined by

$$\begin{equation*}{\partial }_{t}^{\alpha }u\left(x,t\right)=\frac{1}{{\Gamma}\left(1-\alpha \right)}{\int }_{0}^{t}{\partial }_{s}u\left(x,s\right){\left(t-s\right)}^{-\alpha }\mathrm{d}s\end{equation*}$$

with ${\Gamma}\left(\alpha \right)={\int }_{0}^{\infty }{\mathrm{e}}^{-s}{s}^{\alpha -1}\mathrm{d}s$ being the Gamma function. Here, W_x stands for the spatial Brownian motion satisfying $\mathbb{E}\left[{W}_{x}{W}_{y}\right]=x\wedge y$ for any x, y ∈ (0, 1), and ${\dot {W}}_{x}$ denotes the formal derivative of W_x , which should be interpreted as a distribution and is known as the white noise.

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 $\vert \hat{F}\vert$ 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 $\vert \hat{F}\vert$ is apparently not unique, since we measure $\vert \hat{F}\vert$ instead of $\hat{F}$ 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 [34]. 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 [23]. These approaches can be broadly classified into two categories: utilizing either a priori information about the signal F or additional measurements of the modulus $\vert \hat{F}\vert$. In this work, we adopt the latter and apply the phaselift method with random masks to collect additional measurements of the modulus $\vert \hat{F}\vert$, 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 $\vert \hat{F}\vert$, the phase retrieval problem, and its numerical method. Section 5 presents several 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.

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.

Lemma 2.1. Let v be a causal function, i.e., v(t) = 0 if t ⩽ 0, whose fractional derivative ${\partial }_{t}^{\alpha }v$ is well-defined in ${L}^{2}\left(\mathbb{R}\right)$. Then the Fourier transform of the fractional derivatives of v satisfies

$$\begin{equation*}\mathcal{F}\left[{\partial }_{t}^{\alpha }v\right]\left(\omega \right)={\left(\mathrm{i}\omega \right)}^{\alpha }\hat{v}\left(\omega \right)\quad \forall \enspace \alpha \in \left(0,1\right],\end{equation*}$$

where

$$\begin{equation*}\hat{v}\left(\omega \right)=\mathcal{F}\left[v\right]\left(\omega \right){:=}{\int }_{\mathbb{R}}{\mathrm{e}}^{-\mathrm{i}\omega t}v\left(t\right)\enspace \mathrm{d}t\end{equation*}$$

denotes the Fourier transform of v.

Proof. The result is obvious for the case α = 1 if $v\in {H}^{1}\left(\mathbb{R}\right)$. Next we show the assertion for the case α ∈ (0, 1). Define a causal function ${k}_{+}^{\alpha }\left(t\right)$ with α ∈ (0, 1) by

$$\begin{equation*}{k}_{+}^{\alpha }\left(t\right){:=}\begin{cases}\frac{1}{{\Gamma}\left(1-\alpha \right)}{t}^{-\alpha },\quad \hfill & t{ >}0,\hfill \\ 0,\quad \hfill & t{\leqslant}0,\hfill \end{cases}\end{equation*}$$

which satisfies

$$\begin{equation*}{\partial }_{t}^{\alpha }v\left(t\right)=\frac{1}{{\Gamma}\left(1-\alpha \right)}{\int }_{0}^{t}{\partial }_{s}v\left(s\right){\left(t-s\right)}^{-\alpha }\enspace \mathrm{d}s=\left({\partial }_{t}v\;{\ast}\;{k}_{+}^{\alpha }\right)\left(t\right).\end{equation*}$$

Based on the fact

$$\begin{equation*}\mathcal{F}\left[{\partial }_{t}^{\alpha }v\right]\left(\omega \right)=\mathcal{F}\left[{\partial }_{t}v\right]\left(\omega \right)\mathcal{F}\left[{k}_{+}^{\alpha }\right]\left(\omega \right)=\left(\mathrm{i}\omega \right)\hat{v}\left(\omega \right)\mathcal{F}\left[{k}_{+}^{\alpha }\right]\left(\omega \right),\end{equation*}$$

it suffices to show that the Fourier transform of ${k}_{+}^{\alpha }$ admits the following form (cf [35, Sec. 2.9.2]):

$$\begin{equation}\mathcal{F}\left[{k}_{+}^{\alpha }\right]\left(\omega \right)={\int }_{0}^{\infty }\frac{{\mathrm{e}}^{-\mathrm{i}\omega t}{t}^{-\alpha }}{{\Gamma}\left(1-\alpha \right)}\mathrm{d}t={\left(\mathrm{i}\omega \right)}^{\alpha -1}.\end{equation} \tag{ 2.1 }$$

In fact, let ${U}_{R}\subset \mathbb{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

$$\begin{equation*}z{\mapsto}{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha },\quad z\in \mathbb{C}\end{equation*}$$

defines a holomorphic function in U_R . It follows from the Cauchy integral theorem that

$$\begin{equation}{\int }_{{\gamma }_{R}}{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha }\mathrm{d}z=\left[{\int }_{1/R}^{R}+{\int }_{{I}_{R}}+{\int }_{iR}^{i/R}+{\int }_{{I}_{1/R}}\right]{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha }\enspace \mathrm{d}z=0,\end{equation} \tag{ 2.2 }$$

where ${I}_{R}{:=}\left\{z\in \mathbb{C}:\vert z\vert =R\right\}$ and ${I}_{1/R}{:=}\left\{z\in \mathbb{C}:\vert z\vert =1/R\right\}$ 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.

Clearly, we have

$$\begin{equation*}\underset{R\to \infty }{\mathrm{lim}}\left[{\int }_{{I}_{R}}+{\int }_{{I}_{1/R}}\right]{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha }\enspace \mathrm{d}z=0.\end{equation*}$$

Taking the limit of (2.2) as R → +∞ leads to

$$\begin{equation*}{\int }_{0}^{+\infty }{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha }\enspace \mathrm{d}z+{\int }_{+\mathrm{i}\infty }^{0}{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha }\enspace \mathrm{d}z=0.\end{equation*}$$

Let iωz = s. A simple calculation yields

$$\begin{align*}\hfill {\int }_{0}^{+\infty }{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha }\enspace \mathrm{d}z& ={\int }_{0}^{+\mathrm{i}\infty }{\mathrm{e}}^{-\mathrm{i}\omega z}{z}^{-\alpha }\enspace \mathrm{d}z\hfill \\ \hfill & ={\left(\mathrm{i}\omega \right)}^{\alpha -1}{\int }_{0}^{\infty }{\mathrm{e}}^{-s}{s}^{-\alpha }\enspace \mathrm{d}s={\left(\mathrm{i}\omega \right)}^{\alpha -1}{\Gamma}\left(1-\alpha \right),\hfill \end{align*}$$

which completes the proof by noting (2.1). □

Zoom In Zoom Out Reset image size

Remark 2.2. Note that (iω)^α is a multi-valued function when α is a fractional number. Throughout the paper, we define

$$\begin{equation*}{\left(\mathrm{i}\omega \right)}^{\alpha }{:=}\begin{cases}\vert \omega {\vert }^{\alpha }\enspace \mathrm{exp}\left(\frac{\mathrm{i}\pi \alpha }{2}\enspace \mathrm{sgn}\left(\omega \right)\right),\quad \hfill & \omega \ne 0,\hfill \\ 0,\quad \hfill & \omega =0,\hfill \end{cases}\end{equation*}$$

where sgn  (⋅) denotes the sign function.

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 by $\tilde {F}$ the zero extension of F in (−∞, 0) and by $\hat{F}$ the Fourier transform of $\tilde {F}$. Consider the two-point boundary value problem of the stochastic differential equation

$$\begin{equation}\begin{cases}_{xx}U\left(x,\omega \right)-{\left(\mathrm{i}\omega \right)}^{\alpha }U\left(x,\omega \right)=-\hat{F}\left(\omega \right){\dot {W}}_{x},\quad \hfill & x\in D,\enspace \omega \in \mathbb{R},\hfill \\ {\partial }_{x}U\left(0,\omega \right)=0,\quad U\left(1,\omega \right)=0,\quad \hfill & \omega \in \mathbb{R}.\hfill \end{cases}\end{equation} \tag{ 3.1 }$$

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

$$\begin{equation*}\begin{cases}_{xx}{g}_{\omega }\left(x,y\right)-s{g}_{\omega }\left(x,y\right)=\delta \left(x-y\right),\quad \hfill & x,y\in D,\hfill \\ {\partial }_{x}{g}_{\omega }\left(0,y\right)=0,\quad {g}_{\omega }\left(1,y\right)=0,\quad \hfill & y\in D,\hfill \end{cases}\end{equation*}$$

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

$$\begin{equation*}{g}_{\omega }\left(x,y\right)=\begin{cases}_{1}\left(y\right)\enspace {\mathrm{e}}^{-\sqrt{s}x}+{B}_{1}\left(y\right)\enspace {\mathrm{e}}^{\sqrt{s}x},\hfill & x{< }y,\hfill \\ {A}_{2}\left(y\right)\enspace {\mathrm{e}}^{-\sqrt{s}x}+{B}_{2}\left(y\right)\enspace {\mathrm{e}}^{\sqrt{s}x},\hfill & x{ >}y,\hfill \end{cases}\end{equation*}$$

where A_i and B_i , i = 1, 2, are to be determined. Using the boundary conditions

$$\begin{equation*}{\partial }_{x}{g}_{\omega }\left(0,y\right)=0,\qquad {g}_{\omega }\left(1,y\right)=0,\end{equation*}$$

and the continuity and jump conditions

$$\begin{equation*}\underset{x\to {y}^{+}}{\mathrm{lim}}{g}_{\omega }\left(x,y\right)-\underset{x\to {y}^{-}}{\mathrm{lim}}{g}_{\omega }\left(x,y\right)=0,\qquad \underset{x\to {y}^{+}}{\mathrm{lim}}{\partial }_{x}{g}_{\omega }\left(x,y\right)-\underset{x\to {y}^{-}}{\mathrm{lim}}{\partial }_{x}{g}_{\omega }\left(x,y\right)=1.\end{equation*}$$

we may easily obtain from Remark 2.2 that

$$\begin{equation}{g}_{\omega }\left(x,y\right)=\frac{{\mathrm{e}}^{\sqrt{s}\left(x+y\right)}+{\mathrm{e}}^{\sqrt{s}\vert x-y\vert }-{\mathrm{e}}^{\sqrt{s}\left(2-x-y\right)}-{\mathrm{e}}^{\sqrt{s}\left(2-\vert x-y\vert \right)}}{2\sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right)},\quad x,y\in D,\end{equation} \tag{ 3.2 }$$

where

$$\begin{equation*}\sqrt{s}={\left(\mathrm{i}\omega \right)}^{\frac{\alpha }{2}}=\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{exp}\left(\frac{\mathrm{i}\pi \alpha }{4}\enspace \mathrm{sgn}\left(\omega \right)\right),\qquad \mathfrak{R}\left[\sqrt{s}\right]=\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right){ >}0.\end{equation*}$$

Here the notations $\mathfrak{R}\left[\cdot \right]$ and $\mathfrak{S}\left[\cdot \right]$ stand for the real and imaginary parts of a complex value, respectively.

If ω = 0, then g₀(x, y) solves

$$\begin{equation*}\begin{cases}_{xx}{g}_{0}\left(x,y\right)=\delta \left(x-y\right),\quad \hfill & x,y\in D,\hfill \\ {\partial }_{x}{g}_{0}\left(0,y\right)=0,\quad {g}_{0}\left(1,y\right)=0,\quad \hfill & y\in D.\hfill \end{cases}\end{equation*}$$

Similarly, we may solve the above two-point boundary value problem and obtain

$$\begin{equation}{g}_{0}\left(x,y\right)=\mathrm{max}\left\{x,y\right\}-1,\quad x,y\in D.\end{equation} \tag{ 3.3 }$$

Lemma 3.1. The Green function g_ω given in (3.2) and (3.3) satisfies the estimate

$$\begin{equation*}\begin{cases}{g}_{0}{{\Vert}}_{{L}^{2}\left(D{\times}D\right)}=\frac{1}{6}\hfill & \text{if}\enspace \omega =0,\hfill \\ {\Vert}{g}_{\omega }{{\Vert}}_{{L}^{2}\left(D{\times}D\right)}{\leqslant}C\vert \omega {\vert }^{-\alpha }\hfill & \text{if}\enspace \omega \ne 0,\hfill \end{cases}\end{equation*}$$

where C > 0 is a constant independent of ω.

Proof. If ω = 0, then it follows from a simple calculation that

$$\begin{align*}\hfill {\Vert}{g}_{0}{{\Vert}}_{{L}^{2}\left(D{\times}D\right)}^{2}& ={\int }_{0}^{1}\left[{\int }_{0}^{x}{\left(x-1\right)}^{2}\enspace \mathrm{d}y+{\int }_{x}^{1}{\left(y-1\right)}^{2}\enspace \mathrm{d}y\right]\enspace \mathrm{d}x\hfill \\ \hfill & ={\int }_{0}^{1}\left[x{\left(1-x\right)}^{2}+\frac{{\left(1-x\right)}^{3}}{3}\right]\enspace \mathrm{d}x=\frac{1}{6}.\hfill \end{align*}$$

If ω ≠ 0, we get

$$\begin{align}\hfill & {\Vert}{g}_{\omega }{{\Vert}}_{{L}^{2}\left(D{\times}D\right)}^{2}={\int }_{D}{\int }_{D}\vert {g}_{\omega }\left(x,y\right){\vert }^{2}dxdy\hfill \\ \hfill & \quad =\frac{1}{\vert 2\sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}\hfill \\ \hfill & \qquad {\times}{\int }_{0}^{1}{\int }_{0}^{1}{\left\vert {\mathrm{e}}^{\sqrt{s}\left(x+y\right)}+{\mathrm{e}}^{\sqrt{s}\vert x-y\vert }-{\mathrm{e}}^{\sqrt{s}\left(2-x-y\right)}-{\mathrm{e}}^{\sqrt{s}\left(2-\vert x-y\vert \right)}\right\vert }^{2}\enspace \mathrm{d}x\enspace \mathrm{d}y\hfill \\ \hfill & \quad {\leqslant}\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}{\int }_{0}^{1}{\int }_{0}^{1}\enspace {\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]\left(x+y\right)}+{\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]\vert x-y\vert }\hfill \\ \hfill & \qquad +{\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]\left(2-x-y\right)}+{\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]\left(2-\vert x-y\vert \right)}\mathrm{d}x\mathrm{d}y\hfill \\ \hfill & \quad =\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}\frac{{\mathrm{e}}^{4\enspace \mathfrak{R}\left[\sqrt{s}\right]}-1}{\mathfrak{R}\left[\sqrt{s}\right]}\hfill \\ \hfill & \quad =\vert \omega {\vert }^{-\alpha }\frac{1}{\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right)}\hfill \\ \hfill & \qquad {\times}\frac{{\mathrm{e}}^{4\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right)}-1}{\left[{\mathrm{e}}^{4\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right)}+1+2\enspace {\mathrm{e}}^{2\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right)}\enspace \mathrm{cos}\left(2\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{sin}\left(\frac{\pi \alpha }{4}\right)\right)\right]}\hfill \\ \hfill & \quad =\hspace{-2pt}:\;\vert \omega {\vert }^{-\alpha }h\left(\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right)\right),\hfill \end{align} \tag{ 3.4 }$$

where the function h is defined by

$$\begin{equation*}h\left(k\right)=\frac{{\mathrm{e}}^{4k}-1}{k\left[{\mathrm{e}}^{4k}+1+2\enspace {\mathrm{e}}^{2k}\enspace \mathrm{cos}\left(2k\mathrm{tan}\left(\frac{\pi \alpha }{4}\right)\right)\right]}.\end{equation*}$$

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

$$\begin{equation*}\underset{k\to \infty }{\mathrm{lim}}\;h\left(k\right)=\underset{k\to \infty }{\mathrm{lim}}\;\frac{1-{\mathrm{e}}^{-4k}}{k\left[1+{\mathrm{e}}^{-4k}+2\enspace {\mathrm{e}}^{-2k}\enspace \mathrm{cos}\left(2k\mathrm{tan}\left(\frac{\pi \alpha }{4}\right)\right)\right]}=0,\end{equation*}$$

there exists a constant C₀ such that h(k) < 1 for all k > C₀, which shows that h is uniformly bounded on (C₀, ∞). On the other hand, by noting

$$\begin{equation*}\underset{k\to {0}^{+}}{\mathrm{lim}}h\left(k\right)=1,\end{equation*}$$

we get that h is also uniformly bounded on (0, C₀] 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

$$\begin{equation*}{\Vert}{g}_{\omega }{{\Vert}}_{{L}^{2}\left(D{\times}D\right)}^{2}{\leqslant}h\left(\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right)\right)\vert \omega {\vert }^{-\alpha }{\leqslant}C\vert \omega {\vert }^{-\alpha },\end{equation*}$$

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.

Theorem 3.2. Assume that $F\in {H}^{1}\left({\mathbb{R}}_{+}\right)$. Then the stochastic differential equation (3.1) admits a unique mild solution given by

$$\begin{equation}U\left(x,\omega \right)=-\hat{F}\left(\omega \right){\int }_{D}{g}_{\omega }\left(x,y\right)\enspace \mathrm{d}{W}_{y}.\end{equation} \tag{ 3.5 }$$

Moreover, the solution U satisfies the estimate

$$\begin{equation}\mathbb{E}{\Vert}\mathrm{i}\omega U{{\Vert}}_{{L}^{2}\left(\mathbb{R};{L}^{2}\left(D\right)\right)}^{2}{\leqslant}C{\Vert}F{{\Vert}}_{{H}^{1}\left({\mathbb{R}}_{+}\right)}^{2},\end{equation} \tag{ 3.6 }$$

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

$$\begin{align*}\hfill \mathbb{E}{\Vert}\mathrm{i}\omega U{{\Vert}}_{{L}^{2}\left(\mathbb{R};{L}^{2}\left(D\right)\right)}^{2}& ={\int }_{\mathbb{R}}{\int }_{D}\vert \mathrm{i}\omega \hat{F}\left(\omega \right){\vert }^{2}\mathbb{E}{\left\vert {\int }_{D}{g}_{\omega }\left(x,y\right)\mathrm{d}{W}_{y}\right\vert }^{2}\enspace \mathrm{d}x\enspace \mathrm{d}\omega \hfill \\ \hfill & ={\int }_{\mathbb{R}}\vert \mathrm{i}\omega \hat{F}\left(\omega \right){\vert }^{2}{\Vert}{g}_{\omega }{{\Vert}}_{{L}^{2}\left(D{\times}D\right)}^{2}\enspace \mathrm{d}\omega \hfill \\ \hfill & \lesssim {\int }_{\mathbb{R}}\vert \omega {\vert }^{-\alpha }\vert \mathrm{i}\omega \hat{F}\left(\omega \right){\vert }^{2}\enspace \mathrm{d}\omega \hfill \\ \hfill & {\leqslant}{\int }_{\left\{\omega :\vert \omega \vert {\leqslant}1\right\}}\vert \hat{F}\left(\omega \right){\vert }^{2}\enspace \mathrm{d}\omega +{\int }_{\left\{\omega :\vert \omega \vert { >}1\right\}}\vert \mathrm{i}\omega \hat{F}\left(\omega \right){\vert }^{2}\enspace \mathrm{d}\omega \hfill \\ \hfill & {\leqslant}{\Vert}\hat{F}{{\Vert}}_{{L}^{2}\left(\mathbb{R}\right)}^{2}+{\Vert}\mathrm{i}\omega \hat{F}{{\Vert}}_{{L}^{2}\left(\mathbb{R}\right)}^{2}={\Vert}F{{\Vert}}_{{H}^{1}\left({\mathbb{R}}_{+}\right)},\hfill \end{align*}$$

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.

3.2. The direct problem in the time domain

Using the result obtained for the equivalent problem in the frequency domain (3.1), we are now at the position to show the well-posedness of the time domain problem (1.1).

Theorem 3.3. Assume that $F\in {H}^{1}\left({\mathbb{R}}_{+}\right)$. Then the initial boundary value problem of the stochastic differential equation (1.1) admits a unique solution u satisfying

$$\begin{equation*}\mathbb{E}{\Vert}{\partial }_{t}u{{\Vert}}_{{L}^{2}\left(D\right)}^{2}{\leqslant}C{\Vert}F{{\Vert}}_{{H}^{1}\left({\mathbb{R}}_{+}\right)}^{2}\quad \forall \enspace t\in {\mathbb{R}}_{+},\end{equation*}$$

where C > 0 is a constant.

Proof. Let

$$\begin{equation*}\tilde {u}\left(x,t\right)={\mathcal{F}}^{-1}\left[U\left(x,\cdot \right)\right]\left(t\right),\quad x\in D,\enspace t\in \mathbb{R},\end{equation*}$$

where ${\mathcal{F}}^{-1}$ denotes the inverse Fourier transform and U is the mild solution of (3.1). Define

$$\begin{equation}u\left(x,t\right)=\tilde {u}\left(x,t\right){\vert }_{t\in {\mathbb{R}}_{+}}.\end{equation} \tag{ 3.7 }$$

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

$$\begin{equation*}\tilde {u}\left(x,t\right)={\mathcal{F}}^{-1}\left[U\left(x,\cdot \right)\right]\left(t\right)=-{\int }_{-\infty }^{t}\tilde {F}\left(s\right){\mathcal{F}}^{-1}\left[{\int }_{D}{g}_{\omega }\left(x,y\right)\enspace \mathrm{d}{W}_{y}\right]\left(t-s\right)\enspace \mathrm{d}s,\end{equation*}$$

where $\tilde {F}$ is the zero extension of F on (−∞, 0) which is defined at the beginning of this section. Hence, $\tilde {u}$ is a causal function with $\tilde {u}\left(x,t\right)=0$ if t ⩽ 0, which implies that

$$\begin{equation}u\left(x,0\right)=\tilde {u}\left(x,0\right)=0.\end{equation} \tag{ 3.8 }$$

Moreover, it follows from Parseval's identity and Theorem 3.2 that ${\partial }_{t}\tilde {u}\in {L}^{2}\left(\right.{\Omega};{L}^{2}\left(\mathbb{R};{L}^{2}\left(D\right)\right)$ satisfies

$$\begin{equation}\mathbb{E}{\Vert}{\partial }_{t}\tilde {u}{{\Vert}}_{{L}^{2}\left(\mathbb{R};{L}^{2}\left(D\right)\right)}^{2}=\mathbb{E}{\Vert}\mathrm{i}\omega U{{\Vert}}_{{L}^{2}\left(\mathbb{R};{L}^{2}\left(D\right)\right)}^{2}{\leqslant}C{\Vert}F{{\Vert}}_{{H}^{1}\left({\mathbb{R}}_{+}\right)}^{2},\end{equation} \tag{ 3.9 }$$

which indicates that the Caputo fractional derivative of $\tilde {u}$ with respect to t is well-defined.

Taking the inverse Fourier transform with respect to ω on both sides of (3.1), we have

$$\begin{equation}\begin{cases}_{xx}u\left(x,t\right)-{\partial }_{t}^{\alpha }u\left(x,t\right)=-F\left(t\right){\dot {W}}_{x},\quad \hfill & x\in D,\enspace t\in {\mathbb{R}}_{+},\hfill \\ {\partial }_{x}u\left(0,t\right)=0,\quad u\left(1,t\right)=0,\quad \hfill & t\in {\mathbb{R}}_{+},\hfill \end{cases}\end{equation} \tag{ 3.10 }$$

where we have used Lemma 2.1 for the causal function $\tilde {u}$ and the fact $u=\tilde {u}{\vert }_{{\mathbb{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. □

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.

4.1. The reconstruction of $\vert \hat{F}\vert$

First we consider the reconstruction of the modulus of the Fourier coefficients of F, and investigate the uniqueness and the issue of instability of the inverse problem.

4.1.1. Uniqueness

It follows from (3.5) that the mean and variance of the solution U(x, ω) at x = 0 satisfies

$$\begin{equation*}\mathbb{E}\left[U\left(0,\omega \right)\right]=0\end{equation*}$$

and

$$\begin{equation}\mathbb{V}\left[U\left(0,\omega \right)\right]=\mathbb{E}\left[\vert U\left(0,\omega \right){\vert }^{2}\right]=\vert \hat{F}\left(\omega \right){\vert }^{2}{\int }_{D}\vert {g}_{\omega }\left(0,y\right){\vert }^{2}\enspace \mathrm{d}y,\end{equation} \tag{ 4.1 }$$

where we have from (3.2) and (3.3) that

$$\begin{equation}{g}_{\omega }\left(0,y\right)=\frac{{\mathrm{e}}^{\sqrt{s}y}-{\mathrm{e}}^{\sqrt{s}\left(2-y\right)}}{\sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right)},\quad y\in D\end{equation} \tag{ 4.2 }$$

if ω ≠ 0 and

$$\begin{equation*}{g}_{0}\left(0,y\right)=y-1,\quad y\in D\end{equation*}$$

if ω = 0.

Lemma 4.1. For any $\omega \in \mathbb{R}$, it holds

$$\begin{equation*}{\int }_{D}\vert {g}_{\omega }\left(0,y\right){\vert }^{2}\enspace \mathrm{d}y{ >}0.\end{equation*}$$

Proof. For ω = 0, a simple calculation yields

$$\begin{equation*}{\int }_{D}\vert g\left(0,y\right){\vert }^{2}\enspace \mathrm{d}y={\int }_{0}^{1}{\left(y-1\right)}^{2}\enspace \mathrm{d}y=\frac{1}{3}.\end{equation*}$$

For ω ≠ 0, it holds

$$\begin{align*}\hfill {\int }_{D}\vert {g}_{\omega }\left(0,y\right){\vert }^{2}\mathrm{d}y& =\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}{\int }_{0}^{1}{\left\vert {\mathrm{e}}^{\sqrt{s}y}-{\mathrm{e}}^{\sqrt{s}\left(2-y\right)}\right\vert }^{2}\mathrm{d}y\hfill \\ \hfill & =\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}{\int }_{0}^{1}\left[{\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]y}+{\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]\left(2-y\right)}\right.\hfill \\ \hfill & \quad \left.-2\enspace {e}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]}\enspace \mathrm{cos}\left(2\enspace \mathfrak{S}\left[\sqrt{s}\right]\left(1-y\right)\right)\right]\enspace \mathrm{d}y\hfill \\ \hfill & =\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}\left(\frac{{\mathrm{e}}^{4\enspace \mathfrak{R}\left[\sqrt{s}\right]}-1}{2\enspace \mathfrak{R}\left[\sqrt{s}\right]}-\frac{{\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]}\enspace \mathrm{sin}\left(2\enspace \mathfrak{S}\left[\sqrt{s}\right]\right)}{\mathfrak{S}\left[\sqrt{s}\right]}\right)\hfill \\ \hfill & {\geqslant}\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}\left(\frac{{\mathrm{e}}^{4\enspace \mathfrak{R}\left[\sqrt{s}\right]}-1}{2\enspace \mathfrak{R}\left[\sqrt{s}\right]}-2\enspace {\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]}\right)\hfill \\ \hfill & =\hspace{-2pt}:\;\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}{l}_{1}\left(2\enspace \mathfrak{R}\left[\sqrt{s}\right]\right),\hfill \end{align*}$$

where we have used the fact $\frac{\mathrm{sin}\left(x\right)}{x}{\leqslant}1$ for any $x\in \mathbb{R}$, and the function l₁ is defined by

$$\begin{equation*}{l}_{1}\left(k\right)=\frac{{\mathrm{e}}^{2k}-1}{k}-2\enspace {\mathrm{e}}^{k},\quad k{ >}0.\end{equation*}$$

It then suffices to show that l₁(k) > 0 for all k > 0. Note that

$$\begin{equation*}{l}_{1}^{\prime }\left(k\right)=\frac{\left(2k-1\right){\mathrm{e}}^{2k}-2{k}^{2}\enspace {\mathrm{e}}^{k}+1}{{k}^{2}}=\hspace{-2pt}:\;\frac{{l}_{2}\left(k\right)}{{k}^{2}},\end{equation*}$$

where l₂(k) = (2k − 1)e^2k − 2k²e^k + 1. It is easy to check that l₂(k) > 0 and hence l₁'(k) > 0 for all k > 0. In fact, by noting that

$$\begin{equation*}{l}_{2}^{\prime }\left(k\right)=\left(4k\enspace {\mathrm{e}}^{k}-2{k}^{2}-4k\right){\mathrm{e}}^{k}{ >}\left(4k\left(1+k\right)-2{k}^{2}-4k\right){\mathrm{e}}^{k}{ >}0\quad \forall \enspace k{ >}0,\end{equation*}$$

we get

$$\begin{equation*}{l}_{2}\left(k\right){ >}\underset{k\to {0}^{+}}{\mathrm{lim}}{l}_{2}\left(k\right)=0.\end{equation*}$$

Hence the function l₁ is increasing for all k > 0 and satisfies

$$\begin{equation*}{l}_{1}\left(k\right){ >}\underset{k\to {0}^{+}}{\mathrm{lim}}{l}_{1}\left(k\right)=\underset{k\to {0}^{+}}{\mathrm{lim}}\left[\frac{{\mathrm{e}}^{2k}-1}{k}-2\enspace {\mathrm{e}}^{k}\right]=0,\end{equation*}$$

which completes the proof. □

Theorem 4.2. Assume that $F\in {H}^{1}\left({\mathbb{R}}_{+}\right)$. Then the modulus $\vert \hat{F}\left(\omega \right)\vert$ can be uniquely determined by the data $\mathbb{V}\left[U\left(0,\omega \right)\right]$.

Proof. It follows from (4.1) and Lemma 4.1 that

$$\begin{equation*}\vert \hat{F}\left(\omega \right)\vert ={\left(\frac{\mathbb{V}\left[U\left(0,\omega \right)\right]}{{\int }_{D}\vert {g}_{\omega }\left(0,y\right){\vert }^{2}\mathrm{d}y}\right)}^{\frac{1}{2}}\quad \forall \enspace \omega \in \mathbb{R},\end{equation*}$$

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)|²dy goes to zero as ω → ∞, which is stated in the following theorem.

Theorem 4.3. For any fixed ω ≠ 0, it holds

$$\begin{equation*}{\int }_{D}\vert {g}_{\omega }\left(0,y\right){\vert }^{2}\mathrm{d}y\lesssim \vert \omega {\vert }^{-\alpha }.\end{equation*}$$

Proof. By (4.2), we have

$$\begin{align*}\hfill {\int }_{D}\vert {g}_{\omega }\left(0,y\right){\vert }^{2}\mathrm{d}y& =\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}{\int }_{0}^{1}{\left\vert {\mathrm{e}}^{\sqrt{s}y}-{\mathrm{e}}^{\sqrt{s}\left(2-y\right)}\right\vert }^{2}\mathrm{d}y\hfill \\ \hfill & {\leqslant}\frac{2}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}{\int }_{0}^{1}\left({\left\vert {\mathrm{e}}^{\sqrt{s}y}\right\vert }^{2}+{\left\vert {\mathrm{e}}^{\sqrt{s}\left(2-y\right)}\right\vert }^{2}\right)\mathrm{d}y\hfill \\ \hfill & =\frac{2}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}{\int }_{0}^{1}\left({\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]y}+{\mathrm{e}}^{2\enspace \mathfrak{R}\left[\sqrt{s}\right]\left(2-y\right)}\right)\mathrm{d}y\hfill \\ \hfill & =\frac{1}{\vert \sqrt{s}\left(1+{\mathrm{e}}^{2\sqrt{s}}\right){\vert }^{2}}\frac{{\mathrm{e}}^{4\enspace \mathfrak{R}\left[\sqrt{s}\right]}-1}{\mathfrak{R}\left[\sqrt{s}\right]}.\hfill \end{align*}$$

Following from the same estimate as that of (3.4), we get

$$\begin{equation*}{\int }_{D}\vert {g}_{\omega }\left(0,y\right){\vert }^{2}\mathrm{d}y{\leqslant}\vert \omega {\vert }^{-\alpha }h\left(\vert \omega {\vert }^{\frac{\alpha }{2}}\enspace \mathrm{cos}\left(\frac{\pi \alpha }{4}\right)\right)\lesssim \vert \omega {\vert }^{-\alpha },\end{equation*}$$

which completes the proof. □

4.2. Phase retrieval

In this section, we discuss the phase retrieval problem. More precisely, we aim to numerically reconstruct |F| from $\vert \hat{F}\vert$, 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 coefficient $\hat{F}$ 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 [8], which is adopted to solve the phase retrieval problem.

4.2.1. Discrete phase retrieval problem

Let $\boldsymbol{x}={\left({x}_{1},\dots ,{x}_{N}\right)}^{\top }\in {\mathbb{C}}^{N}$ be a signal of length N, and $\boldsymbol{y}={\left({y}_{1},\dots ,{y}_{N}\right)}^{\top }\in {\mathbb{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.,

$$\begin{equation}{\boldsymbol{f}}^{\left(m\right)}={\left({f}_{1}^{\left(m\right)},\dots ,{f}_{N}^{\left(m\right)}\right)}^{\top }{:=}{\left(1,{\mathrm{e}}^{\mathrm{i}\frac{2\pi \left(m-1\right)}{N}},\dots ,{\mathrm{e}}^{\mathrm{i}\frac{2\pi \left(m-1\right)\left(N-1\right)}{N}}\right)}^{\top }.\end{equation} \tag{ 4.3 }$$

Then it is easy to check that y_m = ⟨ f ^(m), x ⟩, where ⟨⋅, ⋅⟩ is the complex inner product defined by

$$\begin{equation*}\langle \boldsymbol{x},\boldsymbol{y}\rangle {:=}\sum\limits _{n=1}^{N}\bar{{x}_{n}}{y}_{n}.\end{equation*}$$

The discrete phase retrieval problem is formulated as follows:

$$\begin{equation}\begin{aligned}\hfill \text{find}\qquad & \boldsymbol{x}\hfill \\ \hfill \text{subject}\;\text{to}\qquad & {z}_{m}{:=}\vert \langle {\boldsymbol{f}}^{\left(m\right)},\boldsymbol{x}\rangle {\vert }^{2},\quad m=1,\dots ,N.\hfill \end{aligned}\end{equation} \tag{ 4.4 }$$

Pioneered by Gerchberg in [19], earlier approaches for the phase retrieval problem are based on alternating projections and can be reformulated as the following least-squares problem:

$$\begin{equation*}\underset{\boldsymbol{x}}{\mathrm{min}}\sum\limits _{m=1}^{M}{\left({z}_{m}-\vert \langle {\boldsymbol{f}}^{\left(m\right)},\boldsymbol{x}\rangle {\vert }^{2}\right)}^{2}.\end{equation*}$$

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 [13], positivity and real-valuedness [17], or sparsity [24]. Depending on the applications, the latter can be done in various ways, which include the use of masks [25], optical gratings [28], oblique illuminations [15], or short-time Fourier transform magnitude measurements which utilize overlap between adjacent short-time sections [36].

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 [8].

First we comment on the multiple structured illuminations. Let $\boldsymbol{x}={\left({x}_{1},\dots ,{x}_{N}\right)}^{\top }\in {\mathbb{C}}^{N}$ be the object of interest, and assume that the illumination schemes, which collect the diffraction pattern of the modulated object ${\left\{{w}_{n}{x}_{n}\right\}}_{n=1,\dots ,N}$, are available, where the pattern $\boldsymbol{w}={\left({w}_{1},\dots ,{w}_{N}\right)}^{\top }\in {\mathbb{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

$$\begin{equation*}\boldsymbol{b}={\left({b}_{1},\dots ,{b}_{N}\right)}^{\top }{:=}{\left(\vert \langle {\boldsymbol{a}}^{\left(1\right)},\boldsymbol{x}\rangle {\vert }^{2},\dots ,\vert \langle {\boldsymbol{a}}^{\left(N\right)},\boldsymbol{x}\rangle {\vert }^{2}\right)}^{\top },\end{equation*}$$

where a ^(m) can be chosen based on the diffraction pattern w in the following form:

$$\begin{equation*}{\boldsymbol{a}}^{\left(m\right)}={\left({f}_{1}^{\left(m\right)}{w}_{1},\dots ,{f}_{N}^{\left(m\right)}{w}_{N}\right)}^{\top },\end{equation*}$$

where ${\boldsymbol{f}}^{\left(m\right)}={\left({f}_{1}^{\left(m\right)},\dots ,{f}_{N}^{\left(m\right)}\right)}^{\top }$ is defined in (4.3). The phase retrieval problem is then transformed to the feasibility problem:

$$\begin{equation*}\begin{aligned}\hfill \text{find}\qquad & \boldsymbol{x}\hfill \\ \hfill \text{subject}\;\text{to}\qquad & {b}_{m}{:=}\vert \langle {\boldsymbol{a}}^{\left(m\right)},\boldsymbol{x}\rangle {\vert }^{2},\quad m=1,\dots ,N.\hfill \end{aligned}\end{equation*}$$

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 ${\boldsymbol{A}}^{\left(m\right)}{:=}{\boldsymbol{a}}^{\left(m\right)}{\left({\boldsymbol{a}}^{\left(m\right)}\right)}^{{\ast}}$, we get

$$\begin{equation}{b}_{m}=\vert \langle {\boldsymbol{a}}^{\left(m\right)},\boldsymbol{x}\rangle {\vert }^{2}={\boldsymbol{x}}^{{\ast}}{\boldsymbol{a}}^{\left(m\right)}{\left({\boldsymbol{a}}^{\left(m\right)}\right)}^{{\ast}}\boldsymbol{x}=\mathrm{Tr}\left({\boldsymbol{A}}^{\left(m\right)}{\boldsymbol{x}\boldsymbol{x}}^{{\ast}}\right),\end{equation} \tag{ 4.5 }$$

where Tr(⋅) denotes the trace of a matrix. Let $\mathcal{A}$ be the linear operator mapping any positive semi-definite matrices X , denoted by X ⪰ 0, into

$$\begin{equation}\mathcal{A}\left(\boldsymbol{X}\right){:=}{\left(\mathrm{Tr}\left({\boldsymbol{A}}^{\left(1\right)}\boldsymbol{X}\right),\dots ,\mathrm{Tr}\left({\boldsymbol{A}}^{\left(N\right)}\boldsymbol{X}\right)\right)}^{\top }.\end{equation} \tag{ 4.6 }$$

Based on the above notation, the phase retrieval problem is equivalent to

$$\begin{equation}\begin{aligned}\hfill \text{find}\qquad & \boldsymbol{X}\hfill \\ \hfill \text{subject}\;\text{to}\qquad & \mathcal{A}\left(\boldsymbol{X}\right)=\boldsymbol{z},\hfill \\ \hfill & \boldsymbol{X}{\succeq}0,\hfill \\ \hfill & \text{rank}\left(\boldsymbol{X}\right)=1,\hfill \end{aligned}\end{equation} \tag{ 4.7 }$$

where $\boldsymbol{z}={\left({z}_{1},\dots ,{z}_{N}\right)}^{\top }$ 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

$$\begin{equation}\begin{aligned}\hfill \text{minimize}\qquad & \text{rank}\left(\boldsymbol{X}\right)\hfill \\ \hfill \text{subject}\;\text{to}\qquad & \mathcal{A}\left(\boldsymbol{X}\right)=\boldsymbol{z},\hfill \\ \hfill & \boldsymbol{X}{\succeq}0.\hfill \end{aligned}\end{equation} \tag{ 4.8 }$$

The equivalence between (4.7) and (4.8) is obvious since $\boldsymbol{b}=\mathcal{A}\left({\boldsymbol{x}\boldsymbol{x}}^{{\ast}}\right)$ 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 [9, 10].

However, it is known that the rank-minimization problem (4.8) is NP-hard. The trace norm is thus explored as a convex surrogate [29] for the rank functional in (4.8), which gives the following semi-definite programming problem:

$$\begin{equation}\begin{aligned}\hfill \text{minimize}\qquad & \mathrm{Tr}\left(\boldsymbol{X}\right)\hfill \\ \hfill \text{subject}\;\text{to}\qquad & \mathcal{A}\left(\boldsymbol{X}\right)=\boldsymbol{z},\hfill \\ \hfill & \boldsymbol{X}\;\text{is}\;\text{Hermitian}\;\text{positive}\;\text{semi-definite}.\hfill \end{aligned}\end{equation} \tag{ 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 [32]. 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 [9, 10] for the detailed discussion.

4.3. Numerical method

In this work, all the numerical algorithms are implemented in MATLAB by modifying the templates for first-order conic solvers (TFOCS) [7]. 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:

$$\begin{equation}\text{minimize}\quad \phi \left(\boldsymbol{x}\right){:=}g\left(\boldsymbol{x}\right)+h\left(\boldsymbol{x}\right),\end{equation} \tag{ 4.10 }$$

where $g:{\mathbb{R}}^{N}\to \mathbb{R}$ is convex and smooth and $h:{\mathbb{R}}^{N}\to \mathbb{R}$ is convex.

Consider a class of first-order methods to solve (4.10) based on iterations with a generalized projection (cf [7]):

$$\begin{equation}{\boldsymbol{x}}_{k+1}=\underset{\boldsymbol{x}}{\mathrm{arg min}}\left[g\left({\boldsymbol{x}}_{k}\right)+\langle \nabla g\left({\boldsymbol{x}}_{k}\right),\boldsymbol{x}-{\boldsymbol{x}}_{k}\rangle +h\left(\boldsymbol{x}\right)+\frac{1}{2{t}_{k}}{\Vert}\boldsymbol{x}-{\boldsymbol{x}}_{k}{{\Vert}}^{2}\right],\end{equation} \tag{ 4.11 }$$

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

$$\begin{equation}g\left({\boldsymbol{x}}_{k+1}\right){\leqslant}g\left({\boldsymbol{x}}_{k}\right)+\langle \nabla g\left({\boldsymbol{x}}_{k}\right),{\boldsymbol{x}}_{k+1}-{\boldsymbol{x}}_{k}\rangle +\frac{1}{2{t}_{k}}{\Vert}{\boldsymbol{x}}_{k+1}-{\boldsymbol{x}}_{k}{{\Vert}}^{2},\end{equation} \tag{ 4.12 }$$

which can be accomplished by assuming that ∇g satisfies a generalized Lipschitz criterion

$$\begin{equation}{\Vert}\nabla g\left(\boldsymbol{x}\right)-\nabla g\left(\boldsymbol{y}\right){{\Vert}}_{{\ast}}{\leqslant}L{\Vert}\boldsymbol{x}-\boldsymbol{y}{\Vert}\end{equation} \tag{ 4.13 }$$

for any x , y belonging to the domain of ϕ. Here, ||⋅||_* denotes the dual norm of the norm ||⋅|| defined by

$$\begin{equation*}{\Vert}g{{\Vert}}_{{\ast}}=\mathrm{sup}\enspace \left\{\langle h,g\rangle :{\Vert}h{\Vert}{\leqslant}1\right\}.\end{equation*}$$

Then the bound of (4.12) is assured for any t_k ⩽ L⁻¹ under the assumption (4.13), which leads to the convergence

$$\begin{equation*}\left\vert \phi \left({\boldsymbol{x}}_{k}\right)-\underset{\boldsymbol{x}}{\mathrm{inf}}\phi \left(\boldsymbol{x}\right)\right\vert {\leqslant}{\epsilon}\end{equation*}$$

in O(L/) iterations for a simple algorithm, known as the forward–backward algorithm or proximal gradient descent [18], based on (4.11).

Optimal or accelerated first-order methods are able to improve the bound of number of iterations to $O\left(\sqrt{L/{\epsilon}}\right)$, and have been studied by many researchers in the past decades (e.g., [37]). The TFOCS implements a variety of the optimal first-order variants based on a variation of the method described by Nesterov in [32]. One variant, the Auslender and Teboulle method [3], is described as follows:

$$\begin{equation}\begin{aligned}\hfill & {\boldsymbol{y}}_{k}=\left(1-{\theta }_{k}\right){\boldsymbol{x}}_{k}+{\theta }_{k}{\bar{\boldsymbol{x}}}_{k},\hfill \\ \hfill & {\bar{\boldsymbol{x}}}_{k+1}=\underset{\boldsymbol{x}}{\mathrm{arg min}}\enspace \langle \nabla g\left({\boldsymbol{y}}_{k}\right),\boldsymbol{x}\rangle +\frac{1}{2}{\theta }_{k}L{\Vert}\boldsymbol{x}-{\bar{\boldsymbol{x}}}_{k}{{\Vert}}^{2}+h\left(\boldsymbol{x}\right),\hfill \\ \hfill & {\boldsymbol{x}}_{k+1}=\left(1-{\theta }_{k}\right){\boldsymbol{x}}_{k}+{\theta }_{k}{\bar{\boldsymbol{x}}}_{k}+1,\hfill \\ \hfill & {\theta }_{k+1}=\frac{2}{1+\sqrt{1+4/{\theta }_{k}^{2}}},\hfill \end{aligned}\end{equation} \tag{ 4.14 }$$

where x ₀ is chosen in the domain of ϕ, ${\bar{\boldsymbol{x}}}_{0}={\boldsymbol{x}}_{0}$, and θ₀ = 1. Here, θ_k is usually referred to as the accelerated parameter.

In this section, we discuss the implementation for solving the direct and inverse random source problems, and report several 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

$$\begin{equation*}\begin{cases}_{t}^{\alpha }u\left(x,t\right)-{\partial }_{xx}u\left(x,t\right)=F\left(t\right)\dot {{W}_{x}},\quad \hfill \\ u\left(x,0\right)=0,\quad \hfill \\ {\partial }_{x}u\left(0,t\right)=0,\quad u\left(1,t\right)=0\quad \hfill \end{cases}\end{equation*}$$

over the interval (0, 1) × (0, T) for some T > 0.

Define the partition of the time and spatial intervals with nodes

$$\begin{equation*}{t}_{n}=n{h}_{t},\enspace \quad n=0,1,\dots ,{N}_{t},\qquad {x}_{i}=i{h}_{x},\quad \enspace i=0,1,\dots ,{N}_{x},\end{equation*}$$

where h_t = T/N_t and h_x = 1/N_x , respectively. Let ${u}_{i}^{n}$ be the numerical approximation to u(x_i , t_n ). The fractional derivative ${\partial }_{t}^{\alpha }u$ at (x_i , t_n ) is approximated by

$$\begin{align*}\hfill {\partial }_{t}^{\alpha }u\left({x}_{i},{t}_{n}\right)& =\frac{1}{{\Gamma}\left(1-\alpha \right)}{\int }_{0}^{{t}_{n}}\frac{\partial u\left({x}_{i},s\right)}{\partial s}\frac{1}{{\left({t}_{n}-s\right)}^{\alpha }}\mathrm{d}s\hfill \\ \hfill & =\frac{1}{{\Gamma}\left(1-\alpha \right)}\sum\limits _{j=1}^{n}{\int }_{{t}_{j-1}}^{{t}_{j}}\frac{\partial u\left({x}_{i},s\right)}{\partial s}\frac{1}{{\left({t}_{n}-s\right)}^{\alpha }}\mathrm{d}s\hfill \\ \hfill & \approx \frac{1}{{\Gamma}\left(1-\alpha \right)}\sum\limits _{j=1}^{n}\frac{{u}_{i}^{j}-{u}_{i}^{j-1}}{{h}_{t}}{\int }_{{t}_{j-1}}^{{t}_{j}}{\left({t}_{n}-s\right)}^{-\alpha }\mathrm{d}s\hfill \\ \hfill & =\frac{1}{{\Gamma}\left(1-\alpha \right)}\frac{1}{{h}_{t}^{\alpha }}\frac{1}{1-\alpha }\sum\limits _{j=1}^{n}\left({u}_{i}^{j}-{u}_{i}^{j-1}\right)\left[{\left(n-j+1\right)}^{1-\alpha }-{\left(n-j\right)}^{1-\alpha }\right]\hfill \\ \hfill & =\frac{1}{{\Gamma}\left(2-\alpha \right)}\frac{1}{{h}_{t}^{\alpha }}\left[{u}_{i}^{n}+\sum\limits _{j=1}^{n-1}{u}_{i}^{j}\left({\left(n-j+1\right)}^{1-\alpha }-2{\left(n-j\right)}^{1-\alpha }\right.\right.\hfill \\ \hfill & \quad \left.\left.+{\left(n-j-1\right)}^{1-\alpha }\right)-{u}_{i}^{0}\left({n}^{1-\alpha }-{\left(n-1\right)}^{1-\alpha }\right)\right].\hfill \end{align*}$$

The second order derivative ∂_xx u at (x_i , t_n ) is approximated by using the central difference method

$$\begin{equation*}{\partial }_{xx}u\left({x}_{i},{t}_{n}\right)\approx \frac{{u}_{i+1}^{n}-2{u}_{i}^{n}+{u}_{i-1}^{n}}{{h}_{x}^{2}}.\end{equation*}$$

The white noise ${\dot {W}}_{x}$ at x = x_i is approximated by the increment [W(x_i+1) − W(x_i )]/h_x of the Brownian motion W satisfying

$$\begin{equation*}\frac{W\left({x}_{i+1}\right)-W\left({x}_{i}\right)}{{h}_{x}}d{=}\frac{{\xi }_{i}}{\sqrt{{h}_{x}}},\end{equation*}$$

where $a\;d{=}\;b$ means that the random variables a and b have the same distribution and ${\left\{{\xi }_{i}\right\}}_{i=0}^{{N}_{x}-1}$ is a set of independent and identically distributed standard normal random variables, denoted by ξ_i ∼ N(0, 1).

Using the above approximations, we obtain the following implicit scheme:

$$\begin{align*}\hfill & \frac{1}{{\Gamma}\left(2-\alpha \right){h}_{t}^{\alpha }}\left[{u}_{i}^{n}+\sum\limits _{j=1}^{n-1}{u}_{i}^{j}\left({\left(n-j+1\right)}^{1-\alpha }-2{\left(n-j\right)}^{1-\alpha }+{\left(n-j-1\right)}^{1-\alpha }\right)\right.\hfill \\ \hfill & \qquad \left.-{u}_{i}^{0}\left({n}^{1-\alpha }-{\left(n-1\right)}^{1-\alpha }\right)\right]-\frac{{u}_{i+1}^{n}-2{u}_{i}^{n}+{u}_{i-1}^{n}}{{h}_{x}^{2}}\hfill \\ \hfill & \qquad =F\left({t}_{n}\right)\frac{W\left({x}_{i+1}\right)-W\left({x}_{i}\right)}{{h}_{x}}.\hfill \end{align*}$$

By denoting $\sigma =\frac{{h}_{x}^{2}}{{\Gamma}\left(2-\alpha \right){h}_{t}^{\alpha }}$, the above scheme can be rewritten as

$$\begin{align}\hfill & -{u}_{i+1}^{n}+\left(\sigma +2\right){u}_{i}^{n}-{u}_{i-1}^{n}\hfill \\ \hfill & \quad ={h}_{x}F\left({t}_{n}\right)\left[W\left({x}_{i+1}\right)-W\left({x}_{i}\right)\right]-\sigma {u}_{i}^{0}\left({n}^{1-\alpha }-{\left(n-1\right)}^{1-\alpha }\right)\hfill \\ \hfill & \quad -\sigma \sum\limits _{j=1}^{n-1}{u}_{i}^{j}\left({\left(n-j+1\right)}^{1-\alpha }-2{\left(n-j\right)}^{1-\alpha }+{\left(n-j-1\right)}^{1-\alpha }\right)\hfill \\ \hfill & \quad ={h}_{x}F\left({t}_{n}\right)\left[W\left({x}_{i+1}\right)-W\left({x}_{i}\right)\right]\hfill \\ \hfill & \quad -\sigma \sum\limits _{j=1}^{n-1}{u}_{i}^{j}\left({\left(n-j+1\right)}^{1-\alpha }-2{\left(n-j\right)}^{1-\alpha }+{\left(n-j-1\right)}^{1-\alpha }\right)\hfill \\ \hfill & \quad =\hspace{-2.5pt}:\;{G}_{i}^{n}\hfill \end{align} \tag{ 5.1 }$$

for n = 1, ⋯, N_t and i = 0, 1, ⋯, N_x − 1, where we have used the initial condition

$$\begin{equation*}{u}_{i}^{0}=0\quad \forall i=0,\dots ,{N}_{x}.\end{equation*}$$

Moreover, according to the boundary condition ∂_x u(0, t) = 0 in (1.1) and noting that

$$\begin{equation*}{\partial }_{x}u\left(0,t\right)=\underset{{h}_{x}\to 0}{\mathrm{lim}}\frac{u\left({h}_{x},t\right)-u\left(-{h}_{x},t\right)}{2{h}_{x}}=0,\end{equation*}$$

we define

$$\begin{equation*}{u}_{-1}^{n}{:=}{u}_{1}^{n}\quad \forall \enspace n=0,\dots ,{N}_{t}.\end{equation*}$$

Hence, when i = 0, (5.1) leads to

$$\begin{equation}-2{u}_{1}^{n}+\left(\sigma +2\right){u}_{0}^{n}={G}_{0}^{n}\quad \forall \enspace n=1,\dots ,{N}_{t}.\end{equation} \tag{ 5.2 }$$

By the boundary condition u(1, t) = 0 in (1.1), the following boundary condition is taken into account to the numerical scheme:

$$\begin{equation*}{u}_{{N}_{x}}^{n}=0\quad \forall \enspace n=0,\dots ,{N}_{t},\end{equation*}$$

and thus (5.1), when i = N_x − 1, turns to be

$$\begin{equation}\left(\sigma +2\right){u}_{{N}_{x}-1}^{n}-{u}_{{N}_{x}-2}^{n}={G}_{{N}_{x}-1}^{n}\quad \forall \enspace n=1,\dots ,{N}_{t}.\end{equation} \tag{ 5.3 }$$

Combining (5.1) with (5.2) and (5.3), we conclude that the numerical scheme has the following compact matrix form:

$$\begin{equation*}\left[\begin{matrix}\hfill \sigma +2\hfill & \hfill -2\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill \sigma +2\hfill & \hfill -1\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill \sigma +2\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill {\vdots}\hfill & \hfill \dots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill -1\hfill & \hfill \sigma +2\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill \sigma +2\hfill \end{matrix}\right]\left[\begin{matrix}\hfill {u}_{0}^{n}\hfill \\ \hfill {u}_{1}^{n}\hfill \\ \hfill {u}_{2}^{n}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {u}_{{N}_{x}-1}^{n}\hfill \end{matrix}\right]=\left[\begin{matrix}\hfill {G}_{0}^{n}\hfill \\ \hfill {G}_{1}^{n}\hfill \\ \hfill {G}_{2}^{n}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {G}_{{N}_{x}-1}^{n}\hfill \end{matrix}\right]\quad \forall \enspace n=1,\dots ,{N}_{t}.\end{equation*}$$

5.2. Implementation of the phaselift

For the problem (1.1), the modulus of the Fourier coefficient $\vert \hat{F}\vert$ can be reconstructed from (4.1) by utilizing the variance $\mathbb{V}\left[U\left(0,\omega \right)\right]$ of the Fourier transform of the field u(0, t). To solve the phase retrieval problem numerically, i.e., to uniquely recover the source $\boldsymbol{f}{:=}{\left(F\left({t}_{1}\right),\dots ,F\left({t}_{{N}_{t}}\right)\right)}^{\top }$ at discrete temporal nodes, more measurements are usually required and the phaselift with masks introduced in section 4.2 is used.

Let M ₁ and M ₂ be two masks with M ₁ being the N_t × N_t identity matrix and M ₂ 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 $\mathcal{A}\left(\boldsymbol{f}{\boldsymbol{f}}^{{\ast}}\right)$, where the operator $\mathcal{A}$ is defined in (4.6) with N = 2N_t , the vectors ${\left\{{\boldsymbol{a}}^{\left(m\right)}\right\}}_{m=1}^{2{N}_{t}}$ involved in the operator $\mathcal{A}$ are defined by

$$\begin{equation*}{\boldsymbol{a}}^{\left(m\right)}=\begin{cases}{\boldsymbol{M}}_{1}{\boldsymbol{f}}^{\left(m\right)},\hfill & m=1,\dots ,{N}_{t},\hfill \\ {\boldsymbol{M}}_{2}{\boldsymbol{f}}^{\left(m-{N}_{t}\right)},\hfill & m={N}_{t}+1,\dots ,2{N}_{t},\hfill \end{cases}\end{equation*}$$

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

$$\begin{equation*}\phi \left(\boldsymbol{f}\right)=g\left(\boldsymbol{f}\right)+h\left(\boldsymbol{f}\right)=\frac{1}{2}{\Vert}\boldsymbol{b}-\mathcal{A}\left(\boldsymbol{f}{\boldsymbol{f}}^{{\ast}}\right){{\Vert}}_{2}^{2}+\lambda \enspace \mathrm{Tr}\left(\boldsymbol{f}{\boldsymbol{f}}^{{\ast}}\right)+h\left(\boldsymbol{f}\right),\end{equation*}$$

where ||⋅||₂ , λ > 0 is a regularization parameter, and h is the indicator function defined by

$$\begin{equation*}h\left(\boldsymbol{f}\right)=\begin{cases}0\hfill & \text{if}\enspace \boldsymbol{f}{\boldsymbol{f}}^{{\ast}}\;\text{is}\;\text{Hermitian}\;\text{positive}\;\text{semi-definite,}\hfill \\ \infty \hfill & \text{else}.\hfill \end{cases}\end{equation*}$$

Here, b is the quadratic Fourier modulus obtained by the measurement ${\left\{{u}_{0}^{n}\right\}}_{n=1}^{2{N}_{t}}$ at x = 0 perturbed by the sources ${\left[{\left({\boldsymbol{M}}_{1}\boldsymbol{f}\right)}^{\top },{\left({\boldsymbol{M}}_{2}\boldsymbol{f}\right)}^{\top }\right]}^{\top }$ with masks. More precisely, $\boldsymbol{b}={\left({b}_{1},\dots ,{b}_{2{N}_{t}}\right)}^{\top }$ has entries

$$\begin{equation*}{b}_{n}=\frac{\mathbb{V}\left[{U}_{0}^{n}\right]}{{\int }_{0}^{1}\vert {g}_{{\omega }_{n}}\left(0,y\right){\vert }^{2}\mathrm{d}y},\quad n=1,\dots ,2{N}_{t},\end{equation*}$$

where ${\left({U}_{0}^{1},\dots ,{U}_{0}^{2{N}_{t}}\right)}^{\top }$ is the 2N_t -point DFT of the measurement ${\left({u}_{0}^{1},\dots ,{u}_{0}^{2{N}_{t}}\right)}^{\top }$ 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 result $\tilde {\boldsymbol{f}}$ is less than a fixed tolerance, i.e., ${\Vert}\mathcal{A}\left(\tilde {\boldsymbol{f}}{\tilde {\boldsymbol{f}}}^{{\ast}}\right)-\boldsymbol{b}{{\Vert}}_{2}{\leqslant}1{0}^{-6}{\Vert}\boldsymbol{b}{{\Vert}}_{2}$. Note that the solution $\tilde {\boldsymbol{f}}={\left({\tilde {f}}_{\hspace{-1.5pt}1},\dots ,{\tilde {f}}_{{\hspace{-1.5pt}N}_{t}}\right)}^{\top }$ is unique only up to the global phase, i.e., ${\left(\vert c{\tilde {f}}_{\hspace{-1.5pt}1}{\vert }^{2},\dots ,\vert c{\tilde {f}}_{{\hspace{-1.5pt}N}_{t}}{\vert }^{2}\right)}^{\top }={\left(\vert F\left({t}_{1}\right){\vert }^{2},\dots ,\vert F\left({t}_{{N}_{t}}\right){\vert }^{2}\right)}^{\top }$ 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 three 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π. In the second example, the source function is taken to be F(t) = sin(2t) cos(3t) with T = π. In the third example, a discontinuous function F is considered:

$$\begin{equation*}F\left(t\right)=\begin{cases}0,\quad \hfill & t\in \left[0,T/5\right),\hfill\\ 2,\quad \hfill & t\in \left[T/5,2T/5\right),\hfill \\ 0.5,\quad \hfill & t\in \left[2T/5,3T/5\right),\hfill \\ 1.5,\quad \hfill & t\in \left[3T/5,4T/5\right),\hfill \\ 0,\quad \hfill & t\in \left[4T/5,T\right],\hfill \end{cases}\end{equation*}$$

where T = π. We present the results for two different cases: α = 0.4 and α = 0.8 in all of the examples. The reconstructed solutions are plotted against the exact solution in figures 2–4, for the three examples, respectively. In the third example, the results are given for three different noise levels σ = 0, σ = 0.03, and σ = 0.05 in order to test the stability of the method. It can be observed from the numerical experiments that the proposed method is effective and robust to reconstruct the modulus of the source functions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

The research of P Li is supported in part by the NSF Grant No. DMS-1912704. The research of X Xu was partially supported by NSFC 11621101, 12071430 and the Fundamental Research Funds for the Central Universities.

No new data were created or analysed in this study.