# Efficient Stochastic Galerkin Spectral Methods for Optimal Control Problems Constrained by Fractional PDEs with Uncertain Inputs

Authors: Shengzhu Shi; Dazhi Zhang; Jiamin Li
Journal: Journal of Sensors
 (2022)

---

## Abstract

This paper is devoted to designing fast solvers and efficient preconditioners for the optimal control problems (OCPs) constrained by stochastic fractional elliptic equations. We first prove the existence and uniqueness of the stochastic optimal control solution and then derive the stochastic optimality system. For the numerical approximation, we use the stochastic Galerkin spectral methods, which apply the stochastic Galerkin method for the discretization of random variables and employ the spectral-Galerkin approach for the approximation of spatial variables. To solve the large coupled saddle-point system resulted from discretization, we adopt the most commonly used MINRES method and a more effective PPCG method in low-rank matrix iteration format. Specially, we develop the efficient preconditioners based on the matrix decomposition method and the virtual variable method. We also study the eigenvalue distribution of the corresponding preconditioned matrix. Besides, for the approximation of the discretized state equation, we use the mean-based approximation and the Ullmann approximation to handle different values of the variance of the random inputs. Finally, we present numerical experiments to demonstrate the effectiveness of our solvers and preconditioners.


---

## Body

## 1. Introduction

Optimal control theory is an important field of applied mathematics. It is usually of interest on how to determine forces or boundary conditions such that the result of engineering or physical PDE systems arrives at some desired state. Recently, investigations have demonstrated that the PDEs involving fractional order tend to provide a better description for modeling complex process in various fields, e.g., in mechanics [1], biology [2], finance [3], and image processing [4]. Because of their high accuracy and wild applications, the fractional PDEs have been a new and active research area in scientific study. We also see a growing interest on their related optimal control problems (OCPs) [5–9]. Furthermore, the considered control systems are modeled with measurement data [10, 11] or random variable inputs (see, for instance, groundwater flows simulation [12]). These factors may cause uncertainties in coefficients, boundary conditions, forcing terms, and the geometry of the domain [13]. Therefore, the deterministic modeling sometimes cannot satisfy the real industrial demand. The study on uncertainty quantification (UQ) has attracted more and more attention nowadays [13–16].In this paper, we shall focus on the fast solvers and the efficient preconditioners for the optimal control problems constrained by the fractional model investigated in [17]. To provide more reliable predictions, we take into account the uncertainties in coefficients. Usually, these uncertainties are described as random processes. We make use of the finite-dimensional noise assumption [18] and approximate the uncertainties by a finite number of random variables. Then, the stochastic OCPs are transferred to parametric OCPs. The existence and uniqueness are consequently derived by Lion’s argument [19] and the well-posedness of the solution of the fractional PDEs [17]. The stochastic optimality system is obtained via a variational inequality method. Especially, we employ the formal Lagrangian method to provide a more intuitive way to construct the adjoint system [20]. The uncertain inputs and nonlocal property of the fractional operator bring the challenges in the design of the efficient numerical algorithms and the corresponding fast solvers.To address the uncertain inputs, we mention some works in UQ study. With the introduction of the generalized polynomial chaos [14], there are mainly two approaches for the discretization of the stochastic part: the stochastic collocation method (SCM) [16] and the stochastic Galerkin method (SGM) [15]. The SCM is more popular since it leads to a set of uncoupled deterministic systems and enjoys the advantage to use the existing spatial algorithms. SCM and its related research have received fruitful results on stochastic OCPs [10, 11, 21]. However, as mentioned in [22], when considering the stochastic OCPs, both spatial and stochastic dimensions are involved in the formulation of the cost functionals. Hence, the cost functional may also include the standard derivation of the state function, which is of our concern. The SCM and its related methods lose their advantage in this case since they no longer lead to uncoupled deterministic systems because of the standard derivation term. The SGM appears to be more competitive for this case; however, the finite-dimensional noise assumption usually leads to a moderate or high dimension in stochastic space and hence brings huge computational cost when solving the discretized optimality system. The fast solvers and preconditioning technique play the key role of the SGM implementation.In addition, due to the nonlocal property of the fractional operator, we employ the efficient spectral approach for the discretization of the spatial part [17]. Although Mao and Shen [17] have developed efficient numerical algorithms for the fractional PDEs, how to establish the fast solvers for the discretized optimality system, which is a saddle-point system of large scale, remains to be studied.Since the preconditioning technique plays a crucial rule in accelerating the convergence of the iterative algorithms, we pay our attention to the construction of the efficient preconditioners. To handle the discretization of the spatial part, we apply the matrix decomposition method [17, 23] and propose a new virtual variable method, which is suited for the discretized optimality system. For the stochastic part, we consider the mean-based preconditioner (MBP) and the Ullmann preconditioner (UP). We develop the preconditioners for the most commonly used MINRES iteration method and a more efficient PPCG iteration method [24]. The fast solvers for the preconditioners and the eigenvalue distribution of the preconditioned matrix are studied. In addition, to reduce the computational cost, the iterative methods are presented in the low-rank matrix format [25–27] when the variance of the uncertainty inputs are small.The rest of the paper is organized as follows. Section2 is devoted to the study of the stochastic OCPs and the stochastic optimality system. We present the numerical algorithm in Section 3. We introduce the MINRES and PPCG methods in the low-rank matrix format and propose the corresponding efficient preconditioners in Sections 4 and 5, respectively. In Section 6, we present some numerical experiments to demonstrate the effectiveness of our solvers. Finally, in Section 7, we draw some conclusions.


## 2. The Stochastic Optimal Control Problems and the Optimality Systems

In this section, we shall consider the fractional OCPs with uncertain inputs. We first demonstrate the well-posedness of the stochastic fractional elliptic PDEs. Then, we prove the existence and uniqueness of the corresponding OCPs. Besides these, we derive the optimality systems which will be used in the numerical approximation later.


### 2.1. The Stochastic Fractional Elliptic PDEs

LetΩ,F,P denote a complete probability space, where Ω is a sample space, F⊂2Ω is a σ-algebra of events, and P:F⟶0,1,PΩ=1 is an appropriate probability measure. Let I=−1,1 and D=I×I. Let a:D×Ω⟶ℝ represent a random field, which means that ax,· is a real-valued random variable defined on Ω,F,P for each x∈D.Let−1∂xs and x∂1s be the left-sided and right-sided Riemann-Liouville fractional derivatives of order s, respectively [28]. Consider the following fractional elliptic PDEs with fractional derivatives from both sides [17]:
(1)Ls,γux,ω=gx,ω,x,ω∈D×Ω,ux,ω=0,x,ω∈∂D×Ω,where1<s,γ<2, x=x,y∈D, g∈L2D⊗L2Ω is a random field representing the distributed control, and Ls,γ is a fractional elliptic operator with uncertain inputs given by
(2)Ls,γux,ω=x∂1s/2d+x,ω−1∂xs/2u+−1∂xs/2d−x,ωx∂1s/2u+y∂1γ/2e+x,ω−1∂yγ/2u+−1∂yγ/2e−x,ωy∂1γ/2u.To ensure the well-posedness of the solution for (1), we suppose that the random diffusion coefficients d+, d−, e+, and e− are bounded and uniformly coercive; i.e., there exist constants 0<d¯<d¯<∞ and 0<e¯<e¯<∞ such that
(3)d¯≤d+x,ω,d−x,ω≤d¯,e¯≤e+x,ω,e−x,ω≤e¯,a.e.x,ω∈D×Ω.In this paper, we assume that the coefficients adopt the finite-dimensional noise assumption. That is, they can be represented by a finite number of prescribed centered, normalized, and mutually uncorrelated random variablesξ=ξmm=1M, where M∈ℕ and ξm=ξmω:Ω⟶Γm⊂ℝ. The preceding assumptions enable the transformation from the stochastic system with random event ω to the deterministic parametric system with variable ξ. Thus, for any random field vx,ω in this paper, we assume that
(4)vx,ω=vx,ξω=v0x+∑m=1Mvmxξmω≔∑m=0Mvmxξm,ξ0=1.We also suppose that the random variablesξmm=1M are independent and ξm is characterized by a probability density function ρm:Γm⟶0,1. Let ρξ=∏m=1Mρmξm be the joint probability density of ξ and Γ=∏m=1MΓm⊂ℝM be the support of such probability density. Then, given the joint probability density of ρξ and a deterministic Banach space BD, we define the associated stochastic space BρD=BD⊗Lρ2Γ equipped with the norm vBρD≔∫Γv·,ξBD2ρξdξ1/2<∞. And, for any u,v∈Lρ2D=L2D⊗Lρ2Γ, the inner product is given by u,vD,ρ=∫Γ∫Dux,ξvx,ξρξdxdξ.In addition, lets>0 and denote by HsI the usual Sobolev space equipped with the norm ·HsI. For s,γ>0, define the Hilbert space
(5)Hs,γD≔v∈L2D:v·,yHsI∈L2Iandvx,·HγI∈L2I,endowed with the norm:(6)vHs,γD≔v·,yHsIL2I2+vx,·HγIL2I21/2.In addition, fors,γ>1/2, zero trace spaces H0s,γD can be defined as the closure of C0∞D with respect to the Hs,γD norm.Notice that(7)∂−1xsu,vI=∂−1xs/2u,x∂1s/2vI,∀u,v∈H0s/2I,and(8)∂x1su,vI=∂x1s/2u,−1∂xs/2vI,∀u,v∈H0s/2I,hold for1<s<2, where ·,·I is the duality between H0s/2I and its dual space H0s/2I′ [17]. Define the weak deterministic formulation of the stochastic equation (1) as follows: for a given g∈H0s/2,γ/2′ρD, find u∈H0s/2,γ/2ρD, such that
(9)Au,v=Gv,∀v∈H0s/2,γ/2ρD,whereAu,v:H0s/2,γ/2ρD×H0s/2,γ/2ρD⟶ℝ is the bilinear form defined as
(10)Au,v≔d+x,ξ−1∂xs/2u,−1∂xs/2vD,ρ+d−x,ξx∂1s/2u,x∂1s/2vD,ρ+e+x,ξ−1∂yγ/2u,−1∂yγ/2vD,ρ+e−x,ξy∂1γ/2u,y∂1γ/2vD,ρ,andGv:H0s/2,γ/2ρ′D⟶ℝ is the functional defined as Gv=g,vD,ρ.We obtain the following theorem.Theorem 1.
For1<s,γ<2 and g∈H0s/2,γ/2ρ′D, problem (9) admits a unique solution satisfying
(11)uH0s/2,γ/2ρD≤CgH0s/2,γ/2ρ′D.Proof.
By the assumption (3) and Lemma 3 in [17], one can prove that A·,· is continuous and coercive in H0s/2,γ/2ρD×H0s/2,γ/2ρD, which implies this theorem in view of the well-known Lax-Milgram lemma (for more information, we refer to [17], which consider the case of the deterministic fractional equation).


### 2.2. The Stochastic Optimal Control Problems

The goal of the stochastic OCPs is to find the appropriate distributed control functiong such that the associated state function u of the fractional system (1) is the best attainable approximation to the prescribed target function ud by some measure.In this paper, we consider the case where the cost functional not only involves the expectation of the distance between the state function and the target function but also the standard deviation of the state function. We model this in the following way: for any givenud∈Lρ2D, find u,g∈U=H0s/2,γ/2ρD⊗Lρ2D such that
(12)minu,g∈UJu,g=α2u−udLρ2D2+β2stduL2D2+λ2gLρ2D2,subject to the state equation (1).The first term in the cost functional (12) measures the expectation of the distance between the state function u and the target function ud. The second term penalizes the standard deviation stdu which is given by stdux,·=Eux,·−Eux,·21/2. The positive constants α and β can be adjusted to measure the importance of the first two terms. The last term, called the Tikhonov regularization term, is added to avoid the ill-posedness. Usually, the Tikhonov parameter λ is a value around λ=10−2, but it can be showed later that our algorithms also work for smaller λ.Theorem 2.
There exists a unique optimal solutionu^,g^∈U to problems (12) and (1).Proof.
Define the fractional control-to-state operatorS:Lρ2D⟶Lρ2D, g⟶u, where u is the associated weak solution of system (1) with the given distributed control function g. Then, the cost functional (12) reduces to the following quadratic optimization problem:
(13)ming∈Uad=Lρ2DJg=α2Sg−udLρ2D2+β2stdSgL2D2+λ2gLρ2D2.
By virtual of Theorem1, H0s/2,γ/2ρD↪Lρ2D and S is a linear and continuous mapping from Lρ2D into Lρ2D, J is continuous and strictly convex. Therefore, J is weakly lower semicontinuous. Furthermore, since λ>0, one can restrict the control problem on a closed, convex, and bounded subset of Lρ2D (see [20], Theorem 2.16). Finally, with the above conditions, we derive the existence and uniqueness of the optimal solution u^,g^ with the minimizing sequence argument (see [20], Theorem 2.14).


### 2.3. The Stochastic Optimality Systems

We devote this subsection to deriving the optimality systems of our model problems. The optimality systems, which are also called the first-order necessary optimality conditions, are essential since they provide a better description for the solution of the OCPs and pave a way for the construction of the numerical algorithms. According to the Lagrange method, we derive the following theorem demonstrating the stochastic optimality systems.Theorem 3.
The controlg, together with the optimal state u and the adjoint state p, is optimal for problems (12) and (1) if and only if the triple g,u,p satisfies the optimality system
(14)Au,v=g,vD,ρ,∀v∈H0s/2,γ/2ρD,Ap,q=αud−α+βu+βEu,qD,ρ,∀q∈H0s/2,γ/2ρD,λg−p,hD,ρ=0,∀h∈Lρ2D.Proof.
Let us start with recalling the condition that the optimal solution has to satisfy. Ifg^ is a solution of problems (12) and (1), then it solves the variational inequality
(15)J′g∧g−g∧≥0,∀g∈Uad=Lρ2D, where J′g^ is the Gâteaux derivative of the reduced cost functional (13) with respect to g at g^. Then, with some calculations [20, 22], it appears that
(16)J′g^g−g^=αSg^−ud+βSg^−ESg^,S∗g−g^D,ρ+λg^,g−g^D,ρ=−S∗αud−α+βSg^+βESg^+λg^,g−g^D,ρ,where S∗ is the adjoint operator of the fractional control-to-state operator S satisfying
(17)Su,vD,ρ=u,S∗vD,ρ,∀u,v∈Lρ2D.
In the sequel, let us determine the adjoint operatorS∗. Recalling (9) and observing that
(18)Ls,γu,vD,ρ=Au,v=u,Ls,γvD,ρ,∀u,v∈H0s/2,γ/2ρD,
we define the following adjoint equation:(19)Ls,γpx,ω=vx,ω,x,ω∈D×Ω,px,ω=0,x,ω∈∂D×Ω.
Sinceu=Sg is the weak solution of system (1) with the distributed control g, one obtains that
(20)Sg,vD,ρ=Sg,Ls,γpD,ρ=u,Ls,γpD,ρ=Au,p=g,pD,ρ,
holds for allg,v∈Lρ2D. Thus, from (17) and (20), one derives the adjoint operator S∗:Lρ2D⟶Lρ2D and S∗v=p, where p∈H0s/2,γ/2ρD is the associate weak solution of the adjoint system (19) with right-hand side input v. In addition, in view of Uad=Lρ2D, one can rewrite the variational inequality (15) as
(21)−p+λg^,hD,ρ=0,∀h∈Lρ2D,where p is the associate weak solution of the adjoint system (19) with right-hand side input v=αud−α+βSg^+βESg^. The identity (21) is also referred as the gradient equation. On the other hand, it can be verified that the above conditions are also sufficient due to the convexity of the reduced cost functional (13).
In conclusion, the above method for deriving the optimality systems starts with the variational inequality, invoking the first-order derivative of the cost functional and the concept of the adjoint operator. It remains for us on how to construct the adjoint operatorS∗ in a more intuitive way. For this purpose, we refer to the formal Lagrange method. Let p be the Lagrange multiplier. Define a Lagrangian function associated with the cost functional (12) and the weak formulation of the state equation (9) as
(22)JLu,g;p=Ju,g+Au,p−Gp.
By setting the Gâteaux derivative of the Lagrangian function (22) with respect to the control variable u equal to zero, one gets the adjoint equation (19) with v=αud−α+βu+βEu. Notice that the procedure above is similar to the computing part of the exact Lagrange method. However, the theory analysis of the latter one is more difficult and complicated. We refer to [20] and the references therein for more detail descriptions about the formal Lagrange method and the exact Lagrange method.


## 2.1. The Stochastic Fractional Elliptic PDEs

LetΩ,F,P denote a complete probability space, where Ω is a sample space, F⊂2Ω is a σ-algebra of events, and P:F⟶0,1,PΩ=1 is an appropriate probability measure. Let I=−1,1 and D=I×I. Let a:D×Ω⟶ℝ represent a random field, which means that ax,· is a real-valued random variable defined on Ω,F,P for each x∈D.Let−1∂xs and x∂1s be the left-sided and right-sided Riemann-Liouville fractional derivatives of order s, respectively [28]. Consider the following fractional elliptic PDEs with fractional derivatives from both sides [17]:
(1)Ls,γux,ω=gx,ω,x,ω∈D×Ω,ux,ω=0,x,ω∈∂D×Ω,where1<s,γ<2, x=x,y∈D, g∈L2D⊗L2Ω is a random field representing the distributed control, and Ls,γ is a fractional elliptic operator with uncertain inputs given by
(2)Ls,γux,ω=x∂1s/2d+x,ω−1∂xs/2u+−1∂xs/2d−x,ωx∂1s/2u+y∂1γ/2e+x,ω−1∂yγ/2u+−1∂yγ/2e−x,ωy∂1γ/2u.To ensure the well-posedness of the solution for (1), we suppose that the random diffusion coefficients d+, d−, e+, and e− are bounded and uniformly coercive; i.e., there exist constants 0<d¯<d¯<∞ and 0<e¯<e¯<∞ such that
(3)d¯≤d+x,ω,d−x,ω≤d¯,e¯≤e+x,ω,e−x,ω≤e¯,a.e.x,ω∈D×Ω.In this paper, we assume that the coefficients adopt the finite-dimensional noise assumption. That is, they can be represented by a finite number of prescribed centered, normalized, and mutually uncorrelated random variablesξ=ξmm=1M, where M∈ℕ and ξm=ξmω:Ω⟶Γm⊂ℝ. The preceding assumptions enable the transformation from the stochastic system with random event ω to the deterministic parametric system with variable ξ. Thus, for any random field vx,ω in this paper, we assume that
(4)vx,ω=vx,ξω=v0x+∑m=1Mvmxξmω≔∑m=0Mvmxξm,ξ0=1.We also suppose that the random variablesξmm=1M are independent and ξm is characterized by a probability density function ρm:Γm⟶0,1. Let ρξ=∏m=1Mρmξm be the joint probability density of ξ and Γ=∏m=1MΓm⊂ℝM be the support of such probability density. Then, given the joint probability density of ρξ and a deterministic Banach space BD, we define the associated stochastic space BρD=BD⊗Lρ2Γ equipped with the norm vBρD≔∫Γv·,ξBD2ρξdξ1/2<∞. And, for any u,v∈Lρ2D=L2D⊗Lρ2Γ, the inner product is given by u,vD,ρ=∫Γ∫Dux,ξvx,ξρξdxdξ.In addition, lets>0 and denote by HsI the usual Sobolev space equipped with the norm ·HsI. For s,γ>0, define the Hilbert space
(5)Hs,γD≔v∈L2D:v·,yHsI∈L2Iandvx,·HγI∈L2I,endowed with the norm:(6)vHs,γD≔v·,yHsIL2I2+vx,·HγIL2I21/2.In addition, fors,γ>1/2, zero trace spaces H0s,γD can be defined as the closure of C0∞D with respect to the Hs,γD norm.Notice that(7)∂−1xsu,vI=∂−1xs/2u,x∂1s/2vI,∀u,v∈H0s/2I,and(8)∂x1su,vI=∂x1s/2u,−1∂xs/2vI,∀u,v∈H0s/2I,hold for1<s<2, where ·,·I is the duality between H0s/2I and its dual space H0s/2I′ [17]. Define the weak deterministic formulation of the stochastic equation (1) as follows: for a given g∈H0s/2,γ/2′ρD, find u∈H0s/2,γ/2ρD, such that
(9)Au,v=Gv,∀v∈H0s/2,γ/2ρD,whereAu,v:H0s/2,γ/2ρD×H0s/2,γ/2ρD⟶ℝ is the bilinear form defined as
(10)Au,v≔d+x,ξ−1∂xs/2u,−1∂xs/2vD,ρ+d−x,ξx∂1s/2u,x∂1s/2vD,ρ+e+x,ξ−1∂yγ/2u,−1∂yγ/2vD,ρ+e−x,ξy∂1γ/2u,y∂1γ/2vD,ρ,andGv:H0s/2,γ/2ρ′D⟶ℝ is the functional defined as Gv=g,vD,ρ.We obtain the following theorem.Theorem 1.
For1<s,γ<2 and g∈H0s/2,γ/2ρ′D, problem (9) admits a unique solution satisfying
(11)uH0s/2,γ/2ρD≤CgH0s/2,γ/2ρ′D.Proof.
By the assumption (3) and Lemma 3 in [17], one can prove that A·,· is continuous and coercive in H0s/2,γ/2ρD×H0s/2,γ/2ρD, which implies this theorem in view of the well-known Lax-Milgram lemma (for more information, we refer to [17], which consider the case of the deterministic fractional equation).


## 2.2. The Stochastic Optimal Control Problems

The goal of the stochastic OCPs is to find the appropriate distributed control functiong such that the associated state function u of the fractional system (1) is the best attainable approximation to the prescribed target function ud by some measure.In this paper, we consider the case where the cost functional not only involves the expectation of the distance between the state function and the target function but also the standard deviation of the state function. We model this in the following way: for any givenud∈Lρ2D, find u,g∈U=H0s/2,γ/2ρD⊗Lρ2D such that
(12)minu,g∈UJu,g=α2u−udLρ2D2+β2stduL2D2+λ2gLρ2D2,subject to the state equation (1).The first term in the cost functional (12) measures the expectation of the distance between the state function u and the target function ud. The second term penalizes the standard deviation stdu which is given by stdux,·=Eux,·−Eux,·21/2. The positive constants α and β can be adjusted to measure the importance of the first two terms. The last term, called the Tikhonov regularization term, is added to avoid the ill-posedness. Usually, the Tikhonov parameter λ is a value around λ=10−2, but it can be showed later that our algorithms also work for smaller λ.Theorem 2.
There exists a unique optimal solutionu^,g^∈U to problems (12) and (1).Proof.
Define the fractional control-to-state operatorS:Lρ2D⟶Lρ2D, g⟶u, where u is the associated weak solution of system (1) with the given distributed control function g. Then, the cost functional (12) reduces to the following quadratic optimization problem:
(13)ming∈Uad=Lρ2DJg=α2Sg−udLρ2D2+β2stdSgL2D2+λ2gLρ2D2.
By virtual of Theorem1, H0s/2,γ/2ρD↪Lρ2D and S is a linear and continuous mapping from Lρ2D into Lρ2D, J is continuous and strictly convex. Therefore, J is weakly lower semicontinuous. Furthermore, since λ>0, one can restrict the control problem on a closed, convex, and bounded subset of Lρ2D (see [20], Theorem 2.16). Finally, with the above conditions, we derive the existence and uniqueness of the optimal solution u^,g^ with the minimizing sequence argument (see [20], Theorem 2.14).


## 2.3. The Stochastic Optimality Systems

We devote this subsection to deriving the optimality systems of our model problems. The optimality systems, which are also called the first-order necessary optimality conditions, are essential since they provide a better description for the solution of the OCPs and pave a way for the construction of the numerical algorithms. According to the Lagrange method, we derive the following theorem demonstrating the stochastic optimality systems.Theorem 3.
The controlg, together with the optimal state u and the adjoint state p, is optimal for problems (12) and (1) if and only if the triple g,u,p satisfies the optimality system
(14)Au,v=g,vD,ρ,∀v∈H0s/2,γ/2ρD,Ap,q=αud−α+βu+βEu,qD,ρ,∀q∈H0s/2,γ/2ρD,λg−p,hD,ρ=0,∀h∈Lρ2D.Proof.
Let us start with recalling the condition that the optimal solution has to satisfy. Ifg^ is a solution of problems (12) and (1), then it solves the variational inequality
(15)J′g∧g−g∧≥0,∀g∈Uad=Lρ2D, where J′g^ is the Gâteaux derivative of the reduced cost functional (13) with respect to g at g^. Then, with some calculations [20, 22], it appears that
(16)J′g^g−g^=αSg^−ud+βSg^−ESg^,S∗g−g^D,ρ+λg^,g−g^D,ρ=−S∗αud−α+βSg^+βESg^+λg^,g−g^D,ρ,where S∗ is the adjoint operator of the fractional control-to-state operator S satisfying
(17)Su,vD,ρ=u,S∗vD,ρ,∀u,v∈Lρ2D.
In the sequel, let us determine the adjoint operatorS∗. Recalling (9) and observing that
(18)Ls,γu,vD,ρ=Au,v=u,Ls,γvD,ρ,∀u,v∈H0s/2,γ/2ρD,
we define the following adjoint equation:(19)Ls,γpx,ω=vx,ω,x,ω∈D×Ω,px,ω=0,x,ω∈∂D×Ω.
Sinceu=Sg is the weak solution of system (1) with the distributed control g, one obtains that
(20)Sg,vD,ρ=Sg,Ls,γpD,ρ=u,Ls,γpD,ρ=Au,p=g,pD,ρ,
holds for allg,v∈Lρ2D. Thus, from (17) and (20), one derives the adjoint operator S∗:Lρ2D⟶Lρ2D and S∗v=p, where p∈H0s/2,γ/2ρD is the associate weak solution of the adjoint system (19) with right-hand side input v. In addition, in view of Uad=Lρ2D, one can rewrite the variational inequality (15) as
(21)−p+λg^,hD,ρ=0,∀h∈Lρ2D,where p is the associate weak solution of the adjoint system (19) with right-hand side input v=αud−α+βSg^+βESg^. The identity (21) is also referred as the gradient equation. On the other hand, it can be verified that the above conditions are also sufficient due to the convexity of the reduced cost functional (13).
In conclusion, the above method for deriving the optimality systems starts with the variational inequality, invoking the first-order derivative of the cost functional and the concept of the adjoint operator. It remains for us on how to construct the adjoint operatorS∗ in a more intuitive way. For this purpose, we refer to the formal Lagrange method. Let p be the Lagrange multiplier. Define a Lagrangian function associated with the cost functional (12) and the weak formulation of the state equation (9) as
(22)JLu,g;p=Ju,g+Au,p−Gp.
By setting the Gâteaux derivative of the Lagrangian function (22) with respect to the control variable u equal to zero, one gets the adjoint equation (19) with v=αud−α+βu+βEu. Notice that the procedure above is similar to the computing part of the exact Lagrange method. However, the theory analysis of the latter one is more difficult and complicated. We refer to [20] and the references therein for more detail descriptions about the formal Lagrange method and the exact Lagrange method.


## 3. Numerical Approximation

By virtue of Theorem3, one can derive the optimal solution (g,u,p) by solving the stochastic optimality system (14). In this section, we consider the numerical approximation to (14) which involves both deterministic part in the physical domain D and stochastic part in the probability domain Γ.


### 3.1. Spatial Discretization and the Jacobi Polynomials

For the spatial part, due to the nonlocal property of the fractional derivatives, the main challenge in developing numerical scheme is how to form and compute mass matrix and stiffness matrix efficiently. To this end, we adopt the spectral approach and take the advantage of the properties of the Jacobi polynomials. LetPN be the space of polynomials of degree less than or equal to N and denote XN=v∈PN:v±1=0. For two-dimensional spatial case, take the same number of modes in each direction and define XN=XN×XN. It can be verified that [23, 29]
(23)XN=spanϕkjx=ϕkxϕjy:k,j=0,1,⋯,N−2,where ϕk·=Lk·−Lk+2· and Lk· is the Legendre polynomial of degree k.Then, the spatial mass matrix can be written asMx=My⊗Mx, where
(24)Mxlk=∫Iϕk−1xϕl−1xdx,My=Mx,k,l=1,⋯,N−1.It can be shown that, with the given basisϕk, Mx and My are symmetric penta-diagonal and can be solved easily [23].In what follows, we shall exploit the structure of the stiffness matrixSx.LetSt,La, St,Ra, and Ma be the matrixes with respect to the coefficient function a and the fractional order t, which elements are defined as
(25)St,Lalk=∫Iaz∂−1zt/2ϕk−1z∂−1zt/2ϕl−1zdz,St,Ralk=∫Iaz∂z1t/2ϕk−1z∂z1t/2ϕl−1zdz,Malk=∫Iazϕk−1zϕl−1zdz,where k,l=1,⋯,N−1, z=x, or z=y.Then, if the coefficients adopt the one-dimensional finite chaos assumption(26)d±x,y,ξ=d±x,ξ=∑m=0Mdm±xξm≔∑m=0Mdm,x±ξm,e±x,y,ξ=e±y,ξ=∑m=0Mem±yξm≔∑m=0Mem,y±ξm,we obtain themth spatial stiffness matrix
(27)Sxm=My⊗Ss,Ldm,x++Ss,Rdm,x−+Sγ,Lem,y++Sγ,Rem,y−⊗Mx.More generally, if the coefficients adopt the two-dimensional finite chaos assumption, we assume that the covariance kernel is separable. According to [30], the two-dimensional eigen-problem of the separable covariance kernel can be expressed as the products of two corresponding one-dimensional eigen-problem. Thus, we assume that
(28)d±x,y,ξ=∑m=0Mdm±xd~m±yξm≔∑m=0Mdm,x±dm,y±ξm,e±x,y,ξ=∑m=0Mem±xe~m±yξm≔∑m=0Mem,x±em,y±ξm.And consequently, we obtain themth spatial stiffness matrix
(29)Sxm=Mdm,y+⊗s,LSdm,x++Mdm,y−⊗s,RSdm,x−+γ,LSem,y+⊗Mem,x++γ,RSem,y−⊗Mem,x−.Next, we recall how to form (25) efficiently. Let Jns,γzs,γ>−1 be the nth Jacobi polynomial with respect to the weight function ωs,γz≔1−zs1+zγ [29]. Then, by virtue of Lemma 4 in [17], one can rewrite St,Lalk and St,Ralk as
(30)St,Lalk=∫I1+z2−tazϕ^k−1zϕ^l−1zdz,St,Ralk=∫I1−z2−tazϕ~k−1zϕ~l−1zdz,where
(31)ϕ^kz=Γk+2Γk−t/2+2Jkt/2,1−t/2z−Γk+3Γk−t/2+3Jk+1t/2,1−t/2z,ϕ~kz=Γk+2Γk−t/2+2Jk1−t/2,t/2z+Γk+3Γk−t/2+3Jk+11−t/2,t/2z.It can be shown that (30) and the rest part of (25) can be efficiently formed by Jacobi-Gauss quadrature formula [29].


### 3.2. Stochastic Discretization and the Generalized Polynomial Chaos Expansion

For the approximation in stochastic space, we use the stochastic Galerkin method. Let us recall the definition of generalized polynomial chaos [14]. Let ψqξ:Γ⟶ℝ be the multivariate polynomials represented by the products of orthonormal univariate polynomials
(32)ψqξ≔∏m=1Mφqmξm,q=q1,⋯,qM,where φqm is the orthogonal polynomial which is determined by the distribution of the given data. Define the total degree (TD) set by
(33)Mp=q∈ℕM:q=∑m=1Mqm≤p.Then, we have the stochastic approximation space(34)YM,p=ψqξ:q∈Mp∈Lρ2Γ,which is referred as the set of basis functions for the generalized polynomial chaos expansion (generalized PCE). The number of the generalized PCE basis functions isP=dimYM,p=dimMp=M+p!/M!p!, where p is the highest degree of the orthogonal polynomial φqm. Then, one can use generalized PCE to approximate the stochastic part of the random field.


### 3.3. Stochastic Galerkin Spectral Methods (SGSMs)

In a nutshell, SGSMs employ the spectral approach for spatial approximation and the generalized PCE for stochastic approximation. The SGSMs for (14) read: find uNP,pNP,gNP∈XN⊗YM,p such that
(35)AuNP,vnp=gNP,vnpD,ρ,∀vnp∈XN⊗YM,p,ApNP,qnp=αINud−α+βuNP+βEuNP,qnpD,ρ,∀qnp∈XN⊗YM,p,λgNP−pNP,hnpD,ρ=0,∀hnp∈XN⊗YM,p,whereIN is the interpolation operator associated with tensor product of Legendre-Gauss-Lobatto points in spatial part. Since XN⊗YM,p∈H0s/2,γ/2ρD, the existence and uniqueness of the discretized optimality system (35) can be established as in Theorems 2 and 3. Then, with the above spaces XN and YM,p, we expand the numerical solution uNP,pNP,gNP as
(36)uNP=∑n,punpϕnψp,pNP=∑n,ppnpϕnψp,gNP=∑n,pgnpϕnψp,where ϕn=ϕkjx with n=N−1∗j−1+k, 1≤k,j≤N−1 and ψp=ψpξ. Let U, P, and G be the matrix of size J×P with elements unp, pnp, and gnp, respectively.Setu=vecU, p=vecP, and g=vecG. Here, by x=vecX, we denote the column vector of size m·n obtained by stacking X∈ℝm×n atop one another from left to right, and we define the operator mat· such that X=matx. In addition, we denote the right-hand side as
(37)Udli=αINud,ϕlψiD,ρ,ud=vecUd.Then, subsisting (36) into (35) and taking vnp=qnp=hnp=ϕnxψpξ, we present (35) in tensor product notation
(38)λM0−M0Mα,βS−MS0gup=0ud0,where M=Mξ⊗Mx, Mα,β=Mξα,β⊗Mx, and S=∑m=0MMξm⊗Sxm. Here, ⊗ is the Kronecker product, Sxm is the mth spatial stiffness matrix, and Mx is the spatial mass matrix which are given in Section 3.1. Mξ, Mξm, and Mξα,β are referred to the stochastic matrixes defined via
(39)Mξmij=∫Γξmψiξψjξρξdξ,0≤m≤M,Mξ=Mξ0,Mξα,β=α+βMξ−βM¯ξ,M¯ξ=diag∫Γψ02ξρξdξ,0,⋯,0.By virtue of the orthogonality and the three term formulation, the matrixMξ and Mξα,β are diagonal and Mξm are symmetric and sparse.


## 3.1. Spatial Discretization and the Jacobi Polynomials

For the spatial part, due to the nonlocal property of the fractional derivatives, the main challenge in developing numerical scheme is how to form and compute mass matrix and stiffness matrix efficiently. To this end, we adopt the spectral approach and take the advantage of the properties of the Jacobi polynomials. LetPN be the space of polynomials of degree less than or equal to N and denote XN=v∈PN:v±1=0. For two-dimensional spatial case, take the same number of modes in each direction and define XN=XN×XN. It can be verified that [23, 29]
(23)XN=spanϕkjx=ϕkxϕjy:k,j=0,1,⋯,N−2,where ϕk·=Lk·−Lk+2· and Lk· is the Legendre polynomial of degree k.Then, the spatial mass matrix can be written asMx=My⊗Mx, where
(24)Mxlk=∫Iϕk−1xϕl−1xdx,My=Mx,k,l=1,⋯,N−1.It can be shown that, with the given basisϕk, Mx and My are symmetric penta-diagonal and can be solved easily [23].In what follows, we shall exploit the structure of the stiffness matrixSx.LetSt,La, St,Ra, and Ma be the matrixes with respect to the coefficient function a and the fractional order t, which elements are defined as
(25)St,Lalk=∫Iaz∂−1zt/2ϕk−1z∂−1zt/2ϕl−1zdz,St,Ralk=∫Iaz∂z1t/2ϕk−1z∂z1t/2ϕl−1zdz,Malk=∫Iazϕk−1zϕl−1zdz,where k,l=1,⋯,N−1, z=x, or z=y.Then, if the coefficients adopt the one-dimensional finite chaos assumption(26)d±x,y,ξ=d±x,ξ=∑m=0Mdm±xξm≔∑m=0Mdm,x±ξm,e±x,y,ξ=e±y,ξ=∑m=0Mem±yξm≔∑m=0Mem,y±ξm,we obtain themth spatial stiffness matrix
(27)Sxm=My⊗Ss,Ldm,x++Ss,Rdm,x−+Sγ,Lem,y++Sγ,Rem,y−⊗Mx.More generally, if the coefficients adopt the two-dimensional finite chaos assumption, we assume that the covariance kernel is separable. According to [30], the two-dimensional eigen-problem of the separable covariance kernel can be expressed as the products of two corresponding one-dimensional eigen-problem. Thus, we assume that
(28)d±x,y,ξ=∑m=0Mdm±xd~m±yξm≔∑m=0Mdm,x±dm,y±ξm,e±x,y,ξ=∑m=0Mem±xe~m±yξm≔∑m=0Mem,x±em,y±ξm.And consequently, we obtain themth spatial stiffness matrix
(29)Sxm=Mdm,y+⊗s,LSdm,x++Mdm,y−⊗s,RSdm,x−+γ,LSem,y+⊗Mem,x++γ,RSem,y−⊗Mem,x−.Next, we recall how to form (25) efficiently. Let Jns,γzs,γ>−1 be the nth Jacobi polynomial with respect to the weight function ωs,γz≔1−zs1+zγ [29]. Then, by virtue of Lemma 4 in [17], one can rewrite St,Lalk and St,Ralk as
(30)St,Lalk=∫I1+z2−tazϕ^k−1zϕ^l−1zdz,St,Ralk=∫I1−z2−tazϕ~k−1zϕ~l−1zdz,where
(31)ϕ^kz=Γk+2Γk−t/2+2Jkt/2,1−t/2z−Γk+3Γk−t/2+3Jk+1t/2,1−t/2z,ϕ~kz=Γk+2Γk−t/2+2Jk1−t/2,t/2z+Γk+3Γk−t/2+3Jk+11−t/2,t/2z.It can be shown that (30) and the rest part of (25) can be efficiently formed by Jacobi-Gauss quadrature formula [29].


## 3.2. Stochastic Discretization and the Generalized Polynomial Chaos Expansion

For the approximation in stochastic space, we use the stochastic Galerkin method. Let us recall the definition of generalized polynomial chaos [14]. Let ψqξ:Γ⟶ℝ be the multivariate polynomials represented by the products of orthonormal univariate polynomials
(32)ψqξ≔∏m=1Mφqmξm,q=q1,⋯,qM,where φqm is the orthogonal polynomial which is determined by the distribution of the given data. Define the total degree (TD) set by
(33)Mp=q∈ℕM:q=∑m=1Mqm≤p.Then, we have the stochastic approximation space(34)YM,p=ψqξ:q∈Mp∈Lρ2Γ,which is referred as the set of basis functions for the generalized polynomial chaos expansion (generalized PCE). The number of the generalized PCE basis functions isP=dimYM,p=dimMp=M+p!/M!p!, where p is the highest degree of the orthogonal polynomial φqm. Then, one can use generalized PCE to approximate the stochastic part of the random field.


## 3.3. Stochastic Galerkin Spectral Methods (SGSMs)

In a nutshell, SGSMs employ the spectral approach for spatial approximation and the generalized PCE for stochastic approximation. The SGSMs for (14) read: find uNP,pNP,gNP∈XN⊗YM,p such that
(35)AuNP,vnp=gNP,vnpD,ρ,∀vnp∈XN⊗YM,p,ApNP,qnp=αINud−α+βuNP+βEuNP,qnpD,ρ,∀qnp∈XN⊗YM,p,λgNP−pNP,hnpD,ρ=0,∀hnp∈XN⊗YM,p,whereIN is the interpolation operator associated with tensor product of Legendre-Gauss-Lobatto points in spatial part. Since XN⊗YM,p∈H0s/2,γ/2ρD, the existence and uniqueness of the discretized optimality system (35) can be established as in Theorems 2 and 3. Then, with the above spaces XN and YM,p, we expand the numerical solution uNP,pNP,gNP as
(36)uNP=∑n,punpϕnψp,pNP=∑n,ppnpϕnψp,gNP=∑n,pgnpϕnψp,where ϕn=ϕkjx with n=N−1∗j−1+k, 1≤k,j≤N−1 and ψp=ψpξ. Let U, P, and G be the matrix of size J×P with elements unp, pnp, and gnp, respectively.Setu=vecU, p=vecP, and g=vecG. Here, by x=vecX, we denote the column vector of size m·n obtained by stacking X∈ℝm×n atop one another from left to right, and we define the operator mat· such that X=matx. In addition, we denote the right-hand side as
(37)Udli=αINud,ϕlψiD,ρ,ud=vecUd.Then, subsisting (36) into (35) and taking vnp=qnp=hnp=ϕnxψpξ, we present (35) in tensor product notation
(38)λM0−M0Mα,βS−MS0gup=0ud0,where M=Mξ⊗Mx, Mα,β=Mξα,β⊗Mx, and S=∑m=0MMξm⊗Sxm. Here, ⊗ is the Kronecker product, Sxm is the mth spatial stiffness matrix, and Mx is the spatial mass matrix which are given in Section 3.1. Mξ, Mξm, and Mξα,β are referred to the stochastic matrixes defined via
(39)Mξmij=∫Γξmψiξψjξρξdξ,0≤m≤M,Mξ=Mξ0,Mξα,β=α+βMξ−βM¯ξ,M¯ξ=diag∫Γψ02ξρξdξ,0,⋯,0.By virtue of the orthogonality and the three term formulation, the matrixMξ and Mξα,β are diagonal and Mξm are symmetric and sparse.


## 4. The Preconditioned MINRES Method in Low-Rank Matrix Format

Noticing that the coefficient matrix of (38) is of a saddle-point form [24], we denote it by A and rewrite it in saddle-point structure
(40)A=ABTB0,where A=λM,0;0,Mα,β is a symmetric and positive definite matrix and B=−M,S is a full row rank matrix. The matrix (40) is generally symmetric but indefinite. Therefore, the preconditioned MINRES method arises naturally to solve the system (38) with the preconditioner of the following form
(41)PM=A00BA−1BT.The detailed algorithm of MINRES can be found in [31]. Next, we introduce the MINRES method in low-rank matrix format and then consider the preconditioning technique when using the stochastic Galerkin method.


### 4.1. The Low-Rank Representations

The stochastic Galerkin approximation always results in solving large coupled system. Therefore, we proceed the above iterative methods in matrix format, and we do not form the matrixes and the preconditioners explicitly. When the variance of the input is small compared with their mean, we employ the low-rank (LR) technique in the iterative method [25]. The LR technique is an effective and popular technique to reduce the storage and computational cost when solving the stochastic problems with the stochastic Galerkin method [26, 27, 32, 33].With right-hand side of (38) in the low-rank matrix format, all operations of the algorithm can be implemented efficiently in low-rank format (see, for instance, [25] and the references therein). More precisely, the matrix-vector products M·x, Mα,β·x, and S·x are performed in matrix-function format. For X=matx=XLXR with XL∈ℝJ×r and XR∈ℝr×P (r≪J,P) of low-rank format, one has
(42)MX=MxXMξ≔MxXLXRMξ,Mα,βX=MxXMξα,β≔MxXLXRMξα,β,SX=∑m=0MSxmXMξm≔∑m=0MSxmXLXRMξm.Therefore, define(43)YL,YR=MXL,XR:YL=MxXL,YR=XRMξ,YL,YR=Mα,βXL,XR:YL=MxXL,YR=XRMξα,β,Y=SXL,XR:Y=∑m=0MSxmXLXRMξm,where Mx is the spatial mass matrix operators and Sxm are the spatial stiffness matrix operators.As mentioned above, to reduce the computational and storage cost, neitherMx nor Sxm are formed explicitly. Let xLp be the pth column of XL,XLp=matxLp. We set the spatial mass matrix operator as
(44)MxXL=vecMxXLpMyp=1r∈ℝJ×r.We denote themth spatial stiffness matrix operator as
(45)SxmXL=vecs,LSdm,x+XLpMdm,y++s,RSdm,x−XLpMdm,y−+Mem,x+XLpγ,LSem,y++Mem,x−XLpγ,RSem,y−p=1r∈ℝJ×r.As for the truncation operatorT, we adopt the most standard QR and SVD (singular value decomposition) truncation method [25]. The other truncation operator options can be found in [33, 34]. Furthermore, for low-rank matrix X=XLXR,Y=YLYR∈ℝJ×P where XL∈ℝJ×rX, XR∈ℝrX×P, YL∈ℝJ×rY, and YR∈ℝrY×P, we compute the low-rank matrix inner product X,YLR as [25]
(46)X,YLR=vecXTvecY=traceXTY=traceXRTXLTYLYR=traceYRXRTXLTYL=vecXRYRTTvecXLTYL.With the above low-rank representations and by virtue of [34], we obtain the LR-MINRES matrix iterative method.


### 4.2. Preconditioning for the MINRES Method

Now, we consider the preconditioning for the MINRES method. Rewrite (41) as
(47)PM=λM000Mα,β000E,where E=1/λM+SMα,β−1S.The spectrum analysis shows that the MINRES method will terminate in at most three iterations with preconditionerPM. However, it is expansive to solve the above preconditioners in real practical implement. One needs to find suitable approximations which not only lead to less iteration count but also can be solved cheaply. Notice that the matrixes M and Mα,β can be solved cheaply. Hence, the difficulty arises from how to approximate the (3,3) block E in (47).Moreover, different from the deterministic OCPs, one challenge in constructing the preconditioner for the system derived from the stochastic Galerkin method is how to approximate the stiffness matrixS=∑m=0MMξm⊗Sxm effectively.LetS^ represent some approximation for S. Enlightened by the deterministic problem [24, 35], one can construct approximations
(48)PE1=S^Mα,β−1S^or(49)PE2=1α+βS^+α+βλMM−1S^+α+βλMforE. The above approximations can be easily solved once the fast solver for S^ or S^+α+β/λM is developed.We note here thatPE1 is the most commonly used approximation and achieves good behavior when λ is not too small. However, the spectrum analysis in [24] shows that the eigenvalues of PE1 depend to a large extent on the Tikhonov parameter λ; thus, the iteration count grows quickly as λ decreases. To solve this, [35] proposed a more robust approach PE2 which can effectively reduce the affection of the Tikhonov parameter λ and was adopted widely in solving stochastic OCPs [27].In view of the spectral approach in spatial space and the matrix decomposition method [17, 23], we construct a new, more accurate, and efficient approximation
(50)PE3=1α+βS^M−1S^+1λM,with respect for the most commonly used mean-based preconditioner (MBP) and Ullmann preconditioner (UP) forS^. We will develop the detailed fast solver algorithms for (50) with respect to MBP and UP in Subsections 4.3 and 4.4, respectively.Here, we present the spectral bounds of the preconditioned system coupled withPM, where E=PE3.Theorem 4.
LetA be preconditioned by PM defined by (47) with E=PE3 defined by (50). Let μ be an eigenvalue of PM−1A. If S^ satisfies S∧−1S−I≤η<1, then μ satisfies one of the following:
(51)μ=1,−cMα+βα1+η≤μ≤121−1+4cM1−η2,or12+1cM1−η≤μ≤1+cMα+βα1+η,where cM=λmaxM/λminM is the condition number of M.Proof.
LetA=ABTB0 be preconditioned by P=A00E. Similar to the process of Proposition 3.2 in [24], the eigenvalues of P−1A are either
(52)1,121−1+4tor121+1+4t,where t is the eigenvalue of E−1BA−1BT.
To obtain the bounds of the eigenvalues ofE−1BA−1BT where E=PE3, we consider the Raleigh quotient
(53)Rx=xTSMα,β−1S+1/λMxxT1/α+βS^M−1S^+1/λMx,
for nonzero vectorx of suitable size. Since S∧−1S−I≤η<1 implies 1−η2≤xTS2x/xTS∧2x≤1+η2, it can be proved that
(54)Rx≤maxα+βαxTSM−1SxxTS^M−1S^x,1≤maxα+βαλmaxM−1xTS2xλminM−1xTS∧2x,1≤cMα+βα1+η2,Rx≥minxTSM−1SxxTS^M−1S^x,1≥1cM1−η2,where cM=λmaxM/λminM.
Finally, using (52) and the above estimates yields the spectral bounds (51).


### 4.3. Fast Solver forPE3 Coupled with Mean-Based Preconditioner

Consider the mean-based preconditioner (MBP) which is given by(55)S^≔SM=Mξ⊗Sx0.Then, with the MBP and the low-rank inputF=FLFR with FL∈ℝJ×r, FR∈ℝr×P, we rewrite X=PE3−1F as
(56)1α+βSx0Mx−1Sx0+1λMxXLXRMξ=FLFR,or equivalently(57)1α+βSx0Mx−1Sx0+1λMxxLp=fLp,p=1,⋯,r,XR=FRMξ−1∈ℝr×P,wherexLp represents the pth column of matrix XL∈ℝJ×r. Since Mξ is diagonal, one can easily decompose the large J×P×J×P coupled system into P (r with low-rank notation) uncoupled small J×J spatial systems.Recall (30) and notice that the spatial stiffness matrices Sx0 and Sy0 are full even when the coefficients are constants; therefore, it is expensive to solve the first r equations of (57) in a direct version. [17, 29] provide an efficient separation of variables method—matrix decomposition method to solve Mx and Sx0 efficiently. The total cost of this method is only ON3.To this end, we introduce the variableyLp=Mx−1Sx0xLp and rewrite the equality w.r.t. xLp of (57) as the following system in a matrix format
(58)1α+βSx0YLpMy+MxYLpSy0+1λMxXLpMy=FLp,MxYLpMy−Sx0XLpMy+MxXLpSy0=0,p=1,⋯,r,where XLp=matxLp, YLp=matyLp, and FLp=matfLp. Next, following the process of full matrix decomposition method, we consider the generalized eigenvalue problem of Mx and Sx0 (resp. My and Sy0) and denote
(59)MxEx=Sx0ExΛx,MyEy=Sy0EyΛy,where Λxresp.Λy is the diagonal matrix with diagonal entries λixi=1N−1 (resp. λjyj=1N−1) and the columns of Ex (resp. Ey) are the corresponding eigenvectors. Set XLp=ExX~LpEyT and YLp=ExY~LpEyT. Multiplying the left of system (58) by Sx0Ex−1, right by Sy0Ey−T and by virtue of (59), we arrive at
(60)1α+βY~LpΛy+ΛxY~Lp+1λΛxX~LpΛy=F~Lp,ΛxY~LpΛy−X~LpΛy+ΛxX~Lp=0,p=1,⋯,r,where F~Lp=Sx0Ex−1FLpSy0Ey−T. Then, it follows from (60) that
(61)X~Lpij=λα+βλixλjyF~Lpijλλix+λjy2+α+βλixλjy2,p=1,⋯,r,i,j=1,⋯,N−1.Finally, we haveXLp=ExX~LpEyT and xLp=vecXLp. We note here that, though we denote the new variable yLp, it does not need to be computed and do not add any computational cost. We define this strategy as virtual variable method, which is very useful in solving the preconditioner of such structure. This technique will be repeatedly used in the remainder of this section.In summary, the full matrix decomposition and virtual variable method for solving(62)1α+βSMM−1SM+1λMx=fwith MBP (SM=Mξ⊗Sx0) in low-rank format consists of the following steps:
(1)
LetF=matf and decompose it into low-rank format F=FLFR(2)
ComputeXR=FRMξ−1(3)
LetxLp be the pth column of XL∈ℝJ×r. For p=1,⋯,r, compute xLp by the following steps:
(i)
Compute the eigenvalues and eigenvectors of the generalized eigenvalue problems (59)(ii)
ComputeF~Lp=Sx0Ex−1FLpSy0Ey−T(iii)
ObtainX~Lp from (61)(iv)
SetXLp=ExX~LpEyT and xLp=vecXLp(4)
SetX=XLXR and x=vecX


### 4.4. Fast Solver forPE3 Coupled with Ullmann Preconditioner

MBP is best suited for most uncertain problems; however, its performance deteriorates as the variance of the random inputs increases. This is mainly because the off-diagonal blocks ofS become more significant. To this end, [36] developed the Ullmann preconditioner (UP) which is constructed by all components of S.
(63)S^≔SU=∑m=0MtraceSxmTSx0traceSx0TSx0Mξm⊗Sx0≔MξU⊗Sx0.It has been proved thatSU is sparse, symmetric, positive definite and MξU minimizes S−M^ξ⊗Sx0F for all M^ξ∈ℝP×P [26].Notice that a traditional method to computetm=traceSxmTSx0 is
(64)tm=traceSxmTSx0=vecSxmT·vecSx0.However, when solving the problem involving the fractional order, to achieve the desired result, the discretization parameterJ is usually very large. Therefore, the computational cost of forming SU by the traditional method might be very expensive.To solve this, recall the structure of the spatial stiffness matrixSxm given in (27) for one-dimensional case and (29) for two-dimensional case and the fact that
(65)traceA⊗B=traceA⊗traceB.For one-dimensional case (27), setting Sxm=s,LSdm,x++s,RSdm,x− and Sym=γ,LSem,y++γ,RSem,y−, we have
(66)tm=traceSxmTSx0=traceMyTMytraceSxmTSx0+traceMyTSy0traceSxmTMx+traceSymTMytraceMxTSx0+traceSymTSy0traceMxTMx=mTmamTa0+bmTb0+mTb0amTm+bmTmmTa0,where m=vecMx=vecMy, am=vecSxm, and bm=vecSym. For two-dimensional case, one can compute tm in a similar way, which we omit the detailed computation here. In this way, we can obtain SU with much less computational cost compared with the traditional method (64).Next, we give the fast solver forPE3 coupled with the approximation SU. Letting x=vecX, y=M−1SUx, and Y=maty, we rewrite X=PE3−1F as the following system in matrix format
(67)1α+βSx0YMξU+1λMxXMξ=F,MxYMξ−Sx0XMξU=0.Different from the MBP case, we cannot decompose the large coupled system (67) into P uncoupled spatial systems directly. To solve this, we recall the partial matrix decomposition method [29]. Consider the generalized eigenvalue problem of MξU and Mξ and denote
(68)MξUEU=MξEUΛU.LetX=X~EUT and Y=Y~EUT and multiply the right of (67) by MξEU−T. It can be verified that (67) is equivalent to
(69)1α+βSx0Y~ΛU+1λMxX~=F~,MxY~−Sx0X~ΛU=0,where F~=FMξEU−T. More precisely, letting x~p (resp. f~p) be the pth column of X~ (resp. F~) and λUp be the pth diagonal entry of the diagonal matrix ΛU, we rewrite (69) as
(70)λUp21α+βSx0Mx−1Sx0x~p+1λMxx~p=f~p,p=1,⋯,P.It is observed that the structure of (70) is similar to that of the first r equations in (57). Then, letting X~p=matx~p and following the process of the full matrix decomposition method introduced above, we deduce that
(71)X~~pij=λα+βλixλjyF~~pijλUp2λλix+λjy2+α+βλixλjy2,p=1,⋯,P,i,j=1,⋯,N−1,where F~~p=Sx0Ex−1F~pSy0Ey−T and the parameters λix, λjy, Ex, Ey are given in (59). Finally, we obtain X~p=ExX~~pEyT and X=X~EUT.In general, the UP can reduce the iteration counts and thus does a better job than MBP when handling the system with large variance inputs.In summary, the full matrix decomposition and virtual variable method for solving(72)1α+βSUM−1SU+1λMx=f,with UPSU=MξU⊗Sx0 consists of the following steps:
(1)
Compute the eigenvalues and eigenvectors of the generalized eigenvalue problems (68)(2)
LetF=matf and compute F~=FMξEU−T(3)
Letx~p be the pth column of X~∈ℝJ×P. For p=1,⋯,P, compute x~p by the following steps:
(i)
Compute the eigenvalues and eigenvectors of the generalized eigenvalue problems (59)(ii)
Letf~p be the pth column of F~, F~p=matf~p and compute F~~p=Sx0Ex−1F~pSy0Ey−T(iii)
ObtainX~~p from (71)(iv)
SetX~p=ExX~~pEyT and x~p=vecX~p(4)
SetX=X~EUT and x=vecX


## 4.1. The Low-Rank Representations

The stochastic Galerkin approximation always results in solving large coupled system. Therefore, we proceed the above iterative methods in matrix format, and we do not form the matrixes and the preconditioners explicitly. When the variance of the input is small compared with their mean, we employ the low-rank (LR) technique in the iterative method [25]. The LR technique is an effective and popular technique to reduce the storage and computational cost when solving the stochastic problems with the stochastic Galerkin method [26, 27, 32, 33].With right-hand side of (38) in the low-rank matrix format, all operations of the algorithm can be implemented efficiently in low-rank format (see, for instance, [25] and the references therein). More precisely, the matrix-vector products M·x, Mα,β·x, and S·x are performed in matrix-function format. For X=matx=XLXR with XL∈ℝJ×r and XR∈ℝr×P (r≪J,P) of low-rank format, one has
(42)MX=MxXMξ≔MxXLXRMξ,Mα,βX=MxXMξα,β≔MxXLXRMξα,β,SX=∑m=0MSxmXMξm≔∑m=0MSxmXLXRMξm.Therefore, define(43)YL,YR=MXL,XR:YL=MxXL,YR=XRMξ,YL,YR=Mα,βXL,XR:YL=MxXL,YR=XRMξα,β,Y=SXL,XR:Y=∑m=0MSxmXLXRMξm,where Mx is the spatial mass matrix operators and Sxm are the spatial stiffness matrix operators.As mentioned above, to reduce the computational and storage cost, neitherMx nor Sxm are formed explicitly. Let xLp be the pth column of XL,XLp=matxLp. We set the spatial mass matrix operator as
(44)MxXL=vecMxXLpMyp=1r∈ℝJ×r.We denote themth spatial stiffness matrix operator as
(45)SxmXL=vecs,LSdm,x+XLpMdm,y++s,RSdm,x−XLpMdm,y−+Mem,x+XLpγ,LSem,y++Mem,x−XLpγ,RSem,y−p=1r∈ℝJ×r.As for the truncation operatorT, we adopt the most standard QR and SVD (singular value decomposition) truncation method [25]. The other truncation operator options can be found in [33, 34]. Furthermore, for low-rank matrix X=XLXR,Y=YLYR∈ℝJ×P where XL∈ℝJ×rX, XR∈ℝrX×P, YL∈ℝJ×rY, and YR∈ℝrY×P, we compute the low-rank matrix inner product X,YLR as [25]
(46)X,YLR=vecXTvecY=traceXTY=traceXRTXLTYLYR=traceYRXRTXLTYL=vecXRYRTTvecXLTYL.With the above low-rank representations and by virtue of [34], we obtain the LR-MINRES matrix iterative method.


## 4.2. Preconditioning for the MINRES Method

Now, we consider the preconditioning for the MINRES method. Rewrite (41) as
(47)PM=λM000Mα,β000E,where E=1/λM+SMα,β−1S.The spectrum analysis shows that the MINRES method will terminate in at most three iterations with preconditionerPM. However, it is expansive to solve the above preconditioners in real practical implement. One needs to find suitable approximations which not only lead to less iteration count but also can be solved cheaply. Notice that the matrixes M and Mα,β can be solved cheaply. Hence, the difficulty arises from how to approximate the (3,3) block E in (47).Moreover, different from the deterministic OCPs, one challenge in constructing the preconditioner for the system derived from the stochastic Galerkin method is how to approximate the stiffness matrixS=∑m=0MMξm⊗Sxm effectively.LetS^ represent some approximation for S. Enlightened by the deterministic problem [24, 35], one can construct approximations
(48)PE1=S^Mα,β−1S^or(49)PE2=1α+βS^+α+βλMM−1S^+α+βλMforE. The above approximations can be easily solved once the fast solver for S^ or S^+α+β/λM is developed.We note here thatPE1 is the most commonly used approximation and achieves good behavior when λ is not too small. However, the spectrum analysis in [24] shows that the eigenvalues of PE1 depend to a large extent on the Tikhonov parameter λ; thus, the iteration count grows quickly as λ decreases. To solve this, [35] proposed a more robust approach PE2 which can effectively reduce the affection of the Tikhonov parameter λ and was adopted widely in solving stochastic OCPs [27].In view of the spectral approach in spatial space and the matrix decomposition method [17, 23], we construct a new, more accurate, and efficient approximation
(50)PE3=1α+βS^M−1S^+1λM,with respect for the most commonly used mean-based preconditioner (MBP) and Ullmann preconditioner (UP) forS^. We will develop the detailed fast solver algorithms for (50) with respect to MBP and UP in Subsections 4.3 and 4.4, respectively.Here, we present the spectral bounds of the preconditioned system coupled withPM, where E=PE3.Theorem 4.
LetA be preconditioned by PM defined by (47) with E=PE3 defined by (50). Let μ be an eigenvalue of PM−1A. If S^ satisfies S∧−1S−I≤η<1, then μ satisfies one of the following:
(51)μ=1,−cMα+βα1+η≤μ≤121−1+4cM1−η2,or12+1cM1−η≤μ≤1+cMα+βα1+η,where cM=λmaxM/λminM is the condition number of M.Proof.
LetA=ABTB0 be preconditioned by P=A00E. Similar to the process of Proposition 3.2 in [24], the eigenvalues of P−1A are either
(52)1,121−1+4tor121+1+4t,where t is the eigenvalue of E−1BA−1BT.
To obtain the bounds of the eigenvalues ofE−1BA−1BT where E=PE3, we consider the Raleigh quotient
(53)Rx=xTSMα,β−1S+1/λMxxT1/α+βS^M−1S^+1/λMx,
for nonzero vectorx of suitable size. Since S∧−1S−I≤η<1 implies 1−η2≤xTS2x/xTS∧2x≤1+η2, it can be proved that
(54)Rx≤maxα+βαxTSM−1SxxTS^M−1S^x,1≤maxα+βαλmaxM−1xTS2xλminM−1xTS∧2x,1≤cMα+βα1+η2,Rx≥minxTSM−1SxxTS^M−1S^x,1≥1cM1−η2,where cM=λmaxM/λminM.
Finally, using (52) and the above estimates yields the spectral bounds (51).


## 4.3. Fast Solver forPE3 Coupled with Mean-Based Preconditioner

Consider the mean-based preconditioner (MBP) which is given by(55)S^≔SM=Mξ⊗Sx0.Then, with the MBP and the low-rank inputF=FLFR with FL∈ℝJ×r, FR∈ℝr×P, we rewrite X=PE3−1F as
(56)1α+βSx0Mx−1Sx0+1λMxXLXRMξ=FLFR,or equivalently(57)1α+βSx0Mx−1Sx0+1λMxxLp=fLp,p=1,⋯,r,XR=FRMξ−1∈ℝr×P,wherexLp represents the pth column of matrix XL∈ℝJ×r. Since Mξ is diagonal, one can easily decompose the large J×P×J×P coupled system into P (r with low-rank notation) uncoupled small J×J spatial systems.Recall (30) and notice that the spatial stiffness matrices Sx0 and Sy0 are full even when the coefficients are constants; therefore, it is expensive to solve the first r equations of (57) in a direct version. [17, 29] provide an efficient separation of variables method—matrix decomposition method to solve Mx and Sx0 efficiently. The total cost of this method is only ON3.To this end, we introduce the variableyLp=Mx−1Sx0xLp and rewrite the equality w.r.t. xLp of (57) as the following system in a matrix format
(58)1α+βSx0YLpMy+MxYLpSy0+1λMxXLpMy=FLp,MxYLpMy−Sx0XLpMy+MxXLpSy0=0,p=1,⋯,r,where XLp=matxLp, YLp=matyLp, and FLp=matfLp. Next, following the process of full matrix decomposition method, we consider the generalized eigenvalue problem of Mx and Sx0 (resp. My and Sy0) and denote
(59)MxEx=Sx0ExΛx,MyEy=Sy0EyΛy,where Λxresp.Λy is the diagonal matrix with diagonal entries λixi=1N−1 (resp. λjyj=1N−1) and the columns of Ex (resp. Ey) are the corresponding eigenvectors. Set XLp=ExX~LpEyT and YLp=ExY~LpEyT. Multiplying the left of system (58) by Sx0Ex−1, right by Sy0Ey−T and by virtue of (59), we arrive at
(60)1α+βY~LpΛy+ΛxY~Lp+1λΛxX~LpΛy=F~Lp,ΛxY~LpΛy−X~LpΛy+ΛxX~Lp=0,p=1,⋯,r,where F~Lp=Sx0Ex−1FLpSy0Ey−T. Then, it follows from (60) that
(61)X~Lpij=λα+βλixλjyF~Lpijλλix+λjy2+α+βλixλjy2,p=1,⋯,r,i,j=1,⋯,N−1.Finally, we haveXLp=ExX~LpEyT and xLp=vecXLp. We note here that, though we denote the new variable yLp, it does not need to be computed and do not add any computational cost. We define this strategy as virtual variable method, which is very useful in solving the preconditioner of such structure. This technique will be repeatedly used in the remainder of this section.In summary, the full matrix decomposition and virtual variable method for solving(62)1α+βSMM−1SM+1λMx=fwith MBP (SM=Mξ⊗Sx0) in low-rank format consists of the following steps:
(1)
LetF=matf and decompose it into low-rank format F=FLFR(2)
ComputeXR=FRMξ−1(3)
LetxLp be the pth column of XL∈ℝJ×r. For p=1,⋯,r, compute xLp by the following steps:
(i)
Compute the eigenvalues and eigenvectors of the generalized eigenvalue problems (59)(ii)
ComputeF~Lp=Sx0Ex−1FLpSy0Ey−T(iii)
ObtainX~Lp from (61)(iv)
SetXLp=ExX~LpEyT and xLp=vecXLp(4)
SetX=XLXR and x=vecX


## 4.4. Fast Solver forPE3 Coupled with Ullmann Preconditioner

MBP is best suited for most uncertain problems; however, its performance deteriorates as the variance of the random inputs increases. This is mainly because the off-diagonal blocks ofS become more significant. To this end, [36] developed the Ullmann preconditioner (UP) which is constructed by all components of S.
(63)S^≔SU=∑m=0MtraceSxmTSx0traceSx0TSx0Mξm⊗Sx0≔MξU⊗Sx0.It has been proved thatSU is sparse, symmetric, positive definite and MξU minimizes S−M^ξ⊗Sx0F for all M^ξ∈ℝP×P [26].Notice that a traditional method to computetm=traceSxmTSx0 is
(64)tm=traceSxmTSx0=vecSxmT·vecSx0.However, when solving the problem involving the fractional order, to achieve the desired result, the discretization parameterJ is usually very large. Therefore, the computational cost of forming SU by the traditional method might be very expensive.To solve this, recall the structure of the spatial stiffness matrixSxm given in (27) for one-dimensional case and (29) for two-dimensional case and the fact that
(65)traceA⊗B=traceA⊗traceB.For one-dimensional case (27), setting Sxm=s,LSdm,x++s,RSdm,x− and Sym=γ,LSem,y++γ,RSem,y−, we have
(66)tm=traceSxmTSx0=traceMyTMytraceSxmTSx0+traceMyTSy0traceSxmTMx+traceSymTMytraceMxTSx0+traceSymTSy0traceMxTMx=mTmamTa0+bmTb0+mTb0amTm+bmTmmTa0,where m=vecMx=vecMy, am=vecSxm, and bm=vecSym. For two-dimensional case, one can compute tm in a similar way, which we omit the detailed computation here. In this way, we can obtain SU with much less computational cost compared with the traditional method (64).Next, we give the fast solver forPE3 coupled with the approximation SU. Letting x=vecX, y=M−1SUx, and Y=maty, we rewrite X=PE3−1F as the following system in matrix format
(67)1α+βSx0YMξU+1λMxXMξ=F,MxYMξ−Sx0XMξU=0.Different from the MBP case, we cannot decompose the large coupled system (67) into P uncoupled spatial systems directly. To solve this, we recall the partial matrix decomposition method [29]. Consider the generalized eigenvalue problem of MξU and Mξ and denote
(68)MξUEU=MξEUΛU.LetX=X~EUT and Y=Y~EUT and multiply the right of (67) by MξEU−T. It can be verified that (67) is equivalent to
(69)1α+βSx0Y~ΛU+1λMxX~=F~,MxY~−Sx0X~ΛU=0,where F~=FMξEU−T. More precisely, letting x~p (resp. f~p) be the pth column of X~ (resp. F~) and λUp be the pth diagonal entry of the diagonal matrix ΛU, we rewrite (69) as
(70)λUp21α+βSx0Mx−1Sx0x~p+1λMxx~p=f~p,p=1,⋯,P.It is observed that the structure of (70) is similar to that of the first r equations in (57). Then, letting X~p=matx~p and following the process of the full matrix decomposition method introduced above, we deduce that
(71)X~~pij=λα+βλixλjyF~~pijλUp2λλix+λjy2+α+βλixλjy2,p=1,⋯,P,i,j=1,⋯,N−1,where F~~p=Sx0Ex−1F~pSy0Ey−T and the parameters λix, λjy, Ex, Ey are given in (59). Finally, we obtain X~p=ExX~~pEyT and X=X~EUT.In general, the UP can reduce the iteration counts and thus does a better job than MBP when handling the system with large variance inputs.In summary, the full matrix decomposition and virtual variable method for solving(72)1α+βSUM−1SU+1λMx=f,with UPSU=MξU⊗Sx0 consists of the following steps:
(1)
Compute the eigenvalues and eigenvectors of the generalized eigenvalue problems (68)(2)
LetF=matf and compute F~=FMξEU−T(3)
Letx~p be the pth column of X~∈ℝJ×P. For p=1,⋯,P, compute x~p by the following steps:
(i)
Compute the eigenvalues and eigenvectors of the generalized eigenvalue problems (59)(ii)
Letf~p be the pth column of F~, F~p=matf~p and compute F~~p=Sx0Ex−1F~pSy0Ey−T(iii)
ObtainX~~p from (71)(iv)
SetX~p=ExX~~pEyT and x~p=vecX~p(4)
SetX=X~EUT and x=vecX


## 5. The PPCG Method in Low-Rank Matrix Format

Enlightened by [24, 37], we find another more efficient symmetric Krylov subspace iterative method—the preconditioned projected conjugate gradient (PPCG) method whose preconditioner takes the form
(73)PP=00−M0FS−MS0,where F=λSM−1S+Mα,β.In this section, we first introduce the PPCG method in low-rank matrix format and then consider the preconditioning technique.


### 5.1. The Algorithm of PPCG Method in Low-Rank Matrix Format

We present the PPCG method in the low-rank matrix format in Algorithm1.Algorithm 1: LR-PPCG matrix iterative method.
Intput: Matrix functions M, Mα,β, S and P−1. Right-hand side Ud in low-rank format. Truncation operator T w.r.t. relative accuracy εrel.Output: Matrixes G,U.R2+=−UdGL+,GR+,UL+,UR+,R2,R2,L,R2,R=P−10,0,R2+Q1,L=−GL+, Q1,R=GR+Q2,L=−UL+, Q2,R=UR+γ0=R2,U+LRα0=γ0/λQ1,MQ1LR+Q2,Mα,βQ2LRGL=α0Q1,L, GR=Q1,RUL=α0Q2,L, UR=Q2,Rr=Mα,βUL,UR+λSGL,GR−UdF/UdFk=0whiler>toldoR1,L+,R1,R+=αkλMQ1,L,Q1,RR2+=R2+αkMα,βQ2,L,Q2,RGL+,GR+,UL+,UR+,R2,R2,L,R2,R=P−1R1,L+,R1,R+,R2+δ=R1+,G+LR+R2+,U+/γkQ1=−GL+GR++δQ1, Q1,L,Q1,R=TQ1Q2=−UL+UR++δQ2, Q2,L,Q2,R=TQ2γk+1=R2,U+LRαk+1=γk+1/λQ1,MQ1LR+Q2,Mα,βQ2LRG=G+αk+1Q1, GL,GR=TGU=U+αk+1Q2, UL,UR=TUr=Mα,βUL,UR+λSGL,GR−UdF/UdFk=k+1end whileAlgorithm 2: Algorithm for GL+,GR+,UL+,UR+,R2,R2,L,R2,R=P−1R1,L+,R1,R+,R2+ in Algorithm 1.
Intput: Matrix functions M, M−1, S and PF−1. Input variables R1,L+∈ℝJ×r, R1,R+∈ℝr×P and R2+∈ℝJ×R.Output: Matrixes GL+, GR+, UL+, UR+, R2, R2,L and R2,R.VL,VR=M−1−R1,L+,R2,L+H=SVL,VRR2=R2+−H, R2,L,R2,R=TR2UL+,UR+=PF−1R2,L,R2,RW=SUL+,UR+, WL,WR=TWGL+,GR+=M−1WL,WRIn Algorithm1, the matrix functions M, Mα,β, and S and the low-rank matrix inner-product X,YLR are defined the same as in Subsection 4.1. PF is the approximation for the (2,2) block of (73), which will be discussed in the next subsection.


### 5.2. Preconditioning for the PPCG Method

Notice that the PPCG method will terminate in only one iteration with preconditionerPP defined by (73). However, similar to the case of the MINRES iteration, it is not cheap to solve the (2,2) block F of (73) directly. Usually, the following two approximations
(74)PF1=λS^M−1S^,PF2=λS^+α+βλMM−1S^+α+βλM,are used for approximatingF. Again, by virtue of the matrix decomposition method, we propose a new approximation for F as
(75)PF3=λS^M−1S^+α+βM,where S^ is the mean-based preconditioner or Ullmann preconditioner for the global stiffness matrix S.Next, we give the spectral bounds of the preconditioned system coupled withPP, where F=PF3.Theorem 5.
LetA be preconditioned by PP defined by (73) with F=PF3 defined by (75). Let μ be an eigenvalue of PP−1A. If S^ satisfies S∧−1S−I≤η<1, then μ satisfies
(76)minαα+β,1cM1−η2≤μ≤cM1+η2where cM=λmaxM/λminM.Proof.
By virtue of Theorem 4.2 in [24], the eigenvalue of PP−1A with F=PF3 satisfies either
(77)μ=1orλSM−1S+Mα,βx=μλS^M−1S^+α+βMx.
To obtain the bounds ofμ, we consider the Raleigh quotient
(78)Rx=xTλSM−1S+Mα+βxxTλS^M−1S^+α+βMx.
Similar to the proof of Theorem4, we get the desired result (76).
The detailed fast solver algorithms for (75) with respect to MBP and UP can be obtained similar to that of Subsections 4.3 and 4.4, respectively. We omit the detailed process here.


## 5.1. The Algorithm of PPCG Method in Low-Rank Matrix Format

We present the PPCG method in the low-rank matrix format in Algorithm1.Algorithm 1: LR-PPCG matrix iterative method.
Intput: Matrix functions M, Mα,β, S and P−1. Right-hand side Ud in low-rank format. Truncation operator T w.r.t. relative accuracy εrel.Output: Matrixes G,U.R2+=−UdGL+,GR+,UL+,UR+,R2,R2,L,R2,R=P−10,0,R2+Q1,L=−GL+, Q1,R=GR+Q2,L=−UL+, Q2,R=UR+γ0=R2,U+LRα0=γ0/λQ1,MQ1LR+Q2,Mα,βQ2LRGL=α0Q1,L, GR=Q1,RUL=α0Q2,L, UR=Q2,Rr=Mα,βUL,UR+λSGL,GR−UdF/UdFk=0whiler>toldoR1,L+,R1,R+=αkλMQ1,L,Q1,RR2+=R2+αkMα,βQ2,L,Q2,RGL+,GR+,UL+,UR+,R2,R2,L,R2,R=P−1R1,L+,R1,R+,R2+δ=R1+,G+LR+R2+,U+/γkQ1=−GL+GR++δQ1, Q1,L,Q1,R=TQ1Q2=−UL+UR++δQ2, Q2,L,Q2,R=TQ2γk+1=R2,U+LRαk+1=γk+1/λQ1,MQ1LR+Q2,Mα,βQ2LRG=G+αk+1Q1, GL,GR=TGU=U+αk+1Q2, UL,UR=TUr=Mα,βUL,UR+λSGL,GR−UdF/UdFk=k+1end whileAlgorithm 2: Algorithm for GL+,GR+,UL+,UR+,R2,R2,L,R2,R=P−1R1,L+,R1,R+,R2+ in Algorithm 1.
Intput: Matrix functions M, M−1, S and PF−1. Input variables R1,L+∈ℝJ×r, R1,R+∈ℝr×P and R2+∈ℝJ×R.Output: Matrixes GL+, GR+, UL+, UR+, R2, R2,L and R2,R.VL,VR=M−1−R1,L+,R2,L+H=SVL,VRR2=R2+−H, R2,L,R2,R=TR2UL+,UR+=PF−1R2,L,R2,RW=SUL+,UR+, WL,WR=TWGL+,GR+=M−1WL,WRIn Algorithm1, the matrix functions M, Mα,β, and S and the low-rank matrix inner-product X,YLR are defined the same as in Subsection 4.1. PF is the approximation for the (2,2) block of (73), which will be discussed in the next subsection.


## 5.2. Preconditioning for the PPCG Method

Notice that the PPCG method will terminate in only one iteration with preconditionerPP defined by (73). However, similar to the case of the MINRES iteration, it is not cheap to solve the (2,2) block F of (73) directly. Usually, the following two approximations
(74)PF1=λS^M−1S^,PF2=λS^+α+βλMM−1S^+α+βλM,are used for approximatingF. Again, by virtue of the matrix decomposition method, we propose a new approximation for F as
(75)PF3=λS^M−1S^+α+βM,where S^ is the mean-based preconditioner or Ullmann preconditioner for the global stiffness matrix S.Next, we give the spectral bounds of the preconditioned system coupled withPP, where F=PF3.Theorem 5.
LetA be preconditioned by PP defined by (73) with F=PF3 defined by (75). Let μ be an eigenvalue of PP−1A. If S^ satisfies S∧−1S−I≤η<1, then μ satisfies
(76)minαα+β,1cM1−η2≤μ≤cM1+η2where cM=λmaxM/λminM.Proof.
By virtue of Theorem 4.2 in [24], the eigenvalue of PP−1A with F=PF3 satisfies either
(77)μ=1orλSM−1S+Mα,βx=μλS^M−1S^+α+βMx.
To obtain the bounds ofμ, we consider the Raleigh quotient
(78)Rx=xTλSM−1S+Mα+βxxTλS^M−1S^+α+βMx.
Similar to the proof of Theorem4, we get the desired result (76).
The detailed fast solver algorithms for (75) with respect to MBP and UP can be obtained similar to that of Subsections 4.3 and 4.4, respectively. We omit the detailed process here.


## 6. Numerical Results

In this section, we present some numerical experiments to show the computational efficiency of our algorithms. Without loss of generality, we assume thate±=1 and d+=d−. We assume that the random inputs d+=d− is represented by the truncated KL expansion characterized by the covariance function
(79)C2x1,x2=σ2exp−x1−x2l2−y1−y2l2,∀xi=xi,yi∈D,i=1,2,respectively, whereσ is the standard deviations. According to [30], the eigenpairs λi,ψixk=1∞ of (79) can be expressed as λi=σλk1λj2 and ψix=ψk1xψj2y where λk1,ψk1xk=1∞ and λj2,ψj2yj=1∞ are eigenpairs of
(80)C1x1,x2=exp−x1−x2l,∀x1,x2∈I,withl2=l. The analytical expressions for the eigenpairs λk,ψkxk=1∞ of (80) are given explicitly in [18]. More precisely, the independent random variables ξmm=1M considered here are uniformly distributed over −3,3.We assume that the desired stateud=sinπxsinπy. Unless otherwise stated, we always set the parameters α=β=1, the fractional orders s=γ=1.8, the mean of the random input Ed+=Ed−=1, the order of the TD space p=3, and discretization parameter J=N−12=16384. All matrix iterations herein are terminated when the relative residual error, measured in the Frobenius norm, is reduced to tol=10−6. We note here that the truncation tolerance trunctol should be smaller enough to guarantee the convergence of the iterative method. However, the smaller the trunctol is, the more CPU time will cost. Therefore, the parameter trunctol should be chosen carefully when implementing the iterative method in low-rank format. In our experiments, we set trunctol=10−8.We first investigate the efficiency of the MINRES preconditioner coupled withPE3 and the PPCG preconditioner coupled with PF3 compared with the traditional approaches PE1, PE2, PF1, and PF2 with different values of the Tikhonov parameter λ. More precisely, we consider the system (38) with two-dimensional (2D) random inputs characterized by (79), and we assume that σ=0.01 and l2=2. Noticing that the behavior of the above preconditioners does not affect by the choice of S^, we employ the simplest MBP S^=SM.Tables1 and 2 report the total iteration counts, the total CPU times of the simulations in standard format and low-rank format, respectively. We observe that in both formats, all preconditioners behave well when λ=10−2. However, the iteration count increases significantly as λ decreases when using the preconditioner PE1 and PF1. The rest preconditioners behave robust. They maintain a relatively equal iteration count with the different value of λ. Moreover, the preconditioners PE3 and PF3 we construct tend to cost less iteration count and perform better than the traditional options. This is because the effectiveness of the iterative solver depends largely on how well the approximation can be. We also did further experiments with different values of J, M, and p, and we obtain similar observations. Specially, the results show that the iteration count keeps the same for different parameters J, M, and p with preconditioners PE3 and PF3, which coincides with the spectrum analysis Theorems 4 and 5.Table 1
Comparison of iterations and CPU times (in seconds) to solve system (38) with given random inputs using MINRES/PPCG method with PE1, PE2, PE3/PF1, PF2, PF3, and SM for J=16384, l2=2, M=13, P=560, σ=0.01, and λ∈10−2,10−3,10−4,10−5.
Standard formatMINRES+SMPPCG+SMPE1PE2PE3PF1PF2PF3λ=10−2# iterCPU time21305.4029426.9615216.56781.73892.98340.42λ=10−3# iterCPU time39565.9831455.3317244.2914154.23892.50450.44λ=10−4# iterCPU time911318.0629425.0519272.8425266.75892.26560.06λ=10−5# iterCPU time2273290.5029425.9519273.3862645.23782.55560.31Table 2
Comparison of iterations and CPU times (in seconds) to solve system (38) with given random inputs using MINRES/PPCG method in low-rank format with PE1, PE2, PE3/PF1, PF2, PF3, and SM for J=16384, l2=2, M=13, P=560, σ=0.01, and λ∈10−2,10−3,10−4,10−5.
Low-rank formatMINRES+SMPPCG+SMPE1PE2PE3PF1PF2PF3λ=10−2# iterCPU time21224.9329339.5815152.72752.02866.31320.80λ=10−3# iterCPU time39424.3631360.8717171.6114104.40866.12428.49λ=10−4# iterCPU time931021.9629332.7019194.2626186.70866.23535.30λ=10−5# iterCPU time2332490.1129325.2619178.4563446.27752.82534.81Next, we demonstrate how the standard derivationσ and the number of the truncated term M influence the behavior of the solvers. Notice that the truncated term M is usually based on the speed of decay of the eigenvalues since ∑m=1∞λm=D. We consider the covariance function C2x1,x2 with different length l2. We report the number of the stochastic discretization parameter P and the parameter ∑m=1Mλm/D in Table 3 with different variant tuples l2,M. We use the MINRES preconditioner PE3 and the PPCG preconditioner PF3 coupled with the approximation SM.Table 3
The number of the stochastic discretization parameterP and the parameter ∑m=1Mλm/D with variant tuples l2,M.
l2,M(1,15)(2,13)(2,12)(3,9)(4,6)P81656045522084∑m=1Mλm/D0.81240.90340.89750.92050.9203It is first indeed evident from Tables4 and 5 that the low-rank approach is more efficient for large-scale simulations and small variance problems. This can be explained by the fact that the efficiency of the low-rank format method mainly depends on the ranks of the numerical solutions. The smaller the rank (compared with P) is, the better the performance is. The ranks of G and U are reported in the tables denoted by rankg/u. It is second observed from Tables 4 and 5 that the iteration count also increases as the increase of σ. This is caused by that the MBP SM cannot approximate S well when σ is large, which means the off-diagonal blocks of S become more significant.Table 4
Comparison of iterations and CPU times (in seconds) to solve system (38) with given random inputs and the ranks of truncated solutions (gNP/uNP) using MINRES method in standard format and low-rank format with PE3 and SM for J=16384, σ∈0.01,0.04,0.09,0.16, and variant tuples l2,M,P.
l2,M,PMINRES+PE3+SM (standard/low-rank)(1,15, 816)(2,12,455)(3,9,220)(4,6,84)σ=0.01# iterCPU timeRank(g/u)15/1515/1515/1515/15387.12/268.39159.15/115.6557.91/45.0115.53/13.0785/9959/6845/4824/29σ=0.04# iterCPU timeRank(g/u)21/21535.80/433.87157/17321/21223.14/186.27109/12121/1981.09/64.2181/9019/1919.70/18.8447/51σ=0.09# iterCPU timeRank(g/u)27/27687.16/612.64211/23227/27287.02/261.78145/15927/27104.28/100.06107/11525/2525.84/26.6958/59σ=0.16# iterCPU timeRank(g/u)35/35887.52/867.63275/36235/35372.22/365.40180/21735/35135.07/137.40126/13635/3536.17/38.8265/66Table 5
Comparison of iterations and CPU times (in seconds) to solve system (38) with given random inputs and the ranks of truncated solutions (gNP/uNP) using PPCG method in standard format and low-rank format with PF3 and SM for J=16384, σ∈0.01,0.04,0.09,0.16, and variant tuples l2,M,P.
l2,M,PPPCG+PF3+SM (standard/low-rank)(1,15, 816)(2,12,455)(3,9,220)(4,6,84)σ=0.01# iterCPU timeRank(g/u)3/371.96/35.39115/1043/329.82/16.0179/713/310.50/6.6258/523/32.77/1.9532/32σ=0.04# iterCPU timeRank(g/u)6/6127.06/98.87168/1566/652.06/41.90117/1105/515.70/13.2292/845/54.12/3.7151/48σ=0.09# iterCPU timeRank(g/u)8/8162.43/147.41229/2048/866.63/62.04161/1418/823.52/23.73116/1088/86.19/6.3260/55σ=0.16# iterCPU timeRank(g/u)12/12234.40/244.09286/25512/1296.63/101.52191/17112/1233.92/37.37135/12211/118.28/9.0067/64Finally, we study the behavior of the preconditionersSM and SU. Notice that the mean we considered here is the constant 1. Table 6 reports the total iteration counts and the total CPU times of solving the system (38) with given random inputs with different σ and different correlation length l2. We can see that SU can reduce the iteration count compared with SM, especially for large variance problems.Table 6
Comparison of iterations and CPU times (in seconds) to solve system (38) with given random inputs using standard MINRES method with PE3 and PPCG method with PF3 and SM, SU for J=16384, σ∈0.1,0.2,0.3,0.4, and l∈1,2,3,4.
MINRES+PE3PPCG+PF3SMSUSMSUσ=0.1l=1,M=15# iterCPU time27682.3827702.928162.888164.70l=2,M=12# iterCPU time27287.4325270.79866.59759.95l=3,M=9# iterCPU time27104.692390.23823.50721.33l=4,M=6# iterCPU time2728.242324.3286.2275.63σ=0.2l=1,M=15# iterCPU time411032.4837961.3415288.4413256.28l=2,M=12# iterCPU time41436.1837400.2415118.5413105.33l=3,M=9# iterCPU time41158.6035137.361541.711234.39l=4,M=6# iterCPU time3940.763132.821410.35118.42σ=0.3l=1,M=15# iterCPU time671687.23591534.7026486.2423439.89l=2,M=12# iterCPU time69734.3957616.9727209.0321164.89l=3,M=9# iterCPU time67258.9853208.152567.721952.99l=4,M=6# iterCPU time6365.794749.692417.291712.59σ=0.4l=1,M=15# iterCPU time1353400.251152990.2654995.1345840.50l=2,M=12# iterCPU time1411500.561091178.2555414.5043329.64l=3,M=9# iterCPU time127491.6293364.8251135.533698.02l=4,M=6# iterCPU time111115.567781.404431.112921.11


## 7. Conclusions

In this paper, we derived the theory analysis and efficient numerical algorithms for the stochastic optimal control problems constrained by fractional elliptic equations. We presented the proof of the uniqueness and existence of the optimal control problems, and we deduced the stochastic optimality system. We applied the stochastic Galerkin spectral methods (SGSMs) to the numerical approximation.To solve the large coupled system resulted from the SGSMs, we adopt the MINRES and PPCG methods in both standard and low-rank format. Specially, to handle the system with different value of the variance, we consider the mean-based preconditioner (MBP) and Ullmann preconditioner (UP) to approximate the stiffness matrixS. We can see that the MBP achieves a better behavior when the standard derivation σ is small and UP can reduce the iteration count when the standard derivation σ is large.More importantly, with spectral approach in spatial space, we developed more accurate and efficient preconditionersPE3 and PF3 with respect to the above two approximations (MBP and UP) by the matrix decomposition method and the virtual variable method for both MINRES method and PPCG method. These preconditioners are proved to achieve better performance than the traditional approach.Finally, let us mention that the MINRES method can be applied to the more general case only if the derived system is symmetric, and it is the most commonly used method when handling the large saddle-point system. However, we can see that the PPCG method is more suitable for our problems. The PPCG method can reduce the iteration count and CPU time effectively compared with the MINRES method.



---

*Source: 6369492-2022-07-22.xml*