A Low-Rank Multigrid Method for the Stochastic Steady-State Diffusion ProblemThis work was supported by the U.S. Department of Energy Office of Advanced Scientific Computing Research, Applied Mathematics program under award DE-SC0009301 and by the U.S. National Science Foundation under grant DMS1418754.

# A Low-Rank Multigrid Method for the Stochastic Steady-State Diffusion Problem††thanks: This work was supported by the U.S. Department of Energy Office of Advanced Scientific Computing Research, Applied Mathematics program under award DE-SC0009301 and by the U.S. National Science Foundation under grant DMS1418754.

Howard C. Elman Department of Computer Science and Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742 ().    Tengfei Su Applied Mathematics & Statistics, and Scientific Computation Program, University of Maryland, College Park, MD 20742 ().
###### Abstract

We study a multigrid method for solving large linear systems of equations with tensor product structure. Such systems are obtained from stochastic finite element discretization of stochastic partial differential equations such as the steady-state diffusion problem with random coefficients. When the variance in the problem is not too large, the solution can be well approximated by a low-rank object. In the proposed multigrid algorithm, the matrix iterates are truncated to low rank to reduce memory requirements and computational effort. The method is proved convergent with an analytic error bound. Numerical experiments show its effectiveness in solving the Galerkin systems compared to the original multigrid solver, especially when the number of degrees of freedom associated with the spatial discretization is large.

s
\newsiamremark

remarkRemark \headersLow-Rank Multigrid for Stochastic Diffusion Problem H. C. Elman and T. Su

tochastic finite element method, multigrid, low-rank approximation

{AMS}

35R60, 60H15, 60H35, 65F10, 65N30, 65N55

## 1 Introduction

Stochastic partial differential equations (SPDEs) arise from physical applications where the parameters of the problem are subject to uncertainty. Discretization of SPDEs gives rise to large linear systems of equations which are computationally expensive to solve. These systems are in general sparse and structured. In particular, the coefficient matrix can often be expressed as a sum of tensor products of smaller matrices [6, 13, 14]. For such systems it is natural to use an iterative solver where the coefficient matrix is never explicitly formed and matrix-vector products are computed efficiently. One way to further reduce costs is to construct low-rank approximations to the desired solution. The iterates are truncated so that the solution method handles only low-rank objects in each iteration. This idea has been used to reduce the costs of iterative solution algorithms based on Krylov subspaces. For example, a low-rank conjugate gradient method was given in [9], and low-rank generalized minimal residual methods have been studied in [2, 10].

In this study, we propose a low-rank multigrid method for solving the Galerkin systems. We consider a steady-state diffusion equation with random diffusion coefficient as model problem, and we use the stochastic finite element method (SFEM, see [1, 7]) for the discretization of the problem. The resulting Galerkin system has tensor product structure and moreover, quantities used in the computation, such as the solution sought, can be expressed in matrix format. It has been shown that such systems admit low-rank approximate solutions [3, 9]. In our proposed multigrid solver, the matrix iterates are truncated to have low rank in each iteration. We derive an analytic bound for the error of the solution and show the convergence of the algorithm. We demonstrate using benchmark problems that the low-rank multigrid solver is often more efficient than a solver that does not use truncation, and that it is especially advantageous in reducing computing time for large-scale problems.

An outline of the paper is as follows. In Section 2 we state the problem and briefly review the stochastic finite element method and the multigrid solver for the stochastic Galerkin system from which the new technique is derived. In Section 3 we discuss the idea of low-rank approximation and introduce the multigrid solver with low-rank truncation. A convergence analysis of the low-rank multigrid solver is also given in this section. The results of numerical experiments are shown in Section 4 to test the performance of the algorithm, and some conclusions are drawn in the last section.

## 2 Model problem

Consider the stochastic steady-state diffusion equation with homogeneous Dirichlet boundary conditions

 {−∇⋅(c(x,ω)∇u(x,ω))=f(x)in D×Ω,u(x,ω)=0on ∂D×Ω. (2.1)

Here is a spatial domain and is a sample space with -algebra and probability measure . The diffusion coefficient is a random field. We consider the case where the source term is deterministic. The stochastic Galerkin formulation of Eq. 2.1 uses a weak formulation: find satisfying

 ∫Ω∫Dc(x,ω)∇u(x,ω)⋅∇v(x,ω)%dxdP=∫Ω∫Df(x)v(x,ω)dxdP (2.2)

for all . The problem is well posed if is bounded and strictly positive, i.e.,

 0

so that the Lax-Milgram lemma establishes existence and uniqueness of the weak solution.

We will assume that the stochastic coefficient is represented as a truncated Karhunen-Love (KL) expansion [11, 12], in terms of a finite collection of uncorrelated random variables :

 c(x,ω)≈c0(x)+m∑l=1√λlcl(x)ξl(ω) (2.3)

where is the mean function, is the th eigenpair of the covariance function , and the eigenvalues are assumed to be in non-increasing order. In Section 4 we will further assume these random variables are independent and identically distributed. Let be the joint density function and be the joint image of . The weak form of Eq. 2.1 is then given as follows: find s.t.

 ∫Γρ(ξ)∫Dc(x,ξ)∇u(x,ξ)⋅∇v(x,ξ)%dxdξ=∫Γρ(ξ)∫Df(x)v(x,ξ)dxdξ (2.4)

for all .

### 2.1 Stochastic finite element method

We briefly review the stochastic finite element method as described in [1, 7]. This method approximates the weak solution of Eq. 2.1 in a finite-dimensional subspace

 Whp=Sh⊗Tp=span{ϕ(x)ψ(ξ)∣ϕ(x)∈Sh,ψ(ξ)∈Tp}, (2.5)

where and are finite-dimensional subspaces of and . We will use quadrilateral elements and piecewise bilinear basis functions for the discretization of the physical space , and generalized polynomial chaos [17] for the stochastic basis functions . The latter are -dimensional orthogonal polynomials whose total degree doesn’t exceed . The orthogonality relation means

 ∫Γψr(ξ)ψs(ξ)ρ(ξ)dξ=δrs∫Γψ2r(ξ)ρ(ξ)dξ.

For instance, Legendre polynomials are used if the random variables have uniform distribution with zero mean and unit variance. The number of degrees of freedom in is

 Nξ=(m+p)!m!p!.

Given the subspace, now one can write the SFEM solution as a linear combination of the basis functions,

 uhp(x,ξ)=Nx∑j=1Nξ∑s=1ujsϕj(x)ψs(ξ), (2.6)

where is the dimension of the subspace . Substituting Eqs. 2.6 and 2.3 into Eq. 2.4, and taking the test function as any basis function results in the Galerkin system: find , s.t.

 Au=f. (2.7)

The coefficient matrix can be represented in tensor product notation [14],

 A=G0⊗K0+m∑l=1Gl⊗Kl, (2.8)

where are the stiffness matrices and correspond to the stochastic part, with entries

 G0(r,s) =∫Γψr(ξ)ψs(ξ)ρ(ξ)dξ,K0(i,j)=∫Dc0(x)∇ϕi(x)∇ϕj(x)dx, (2.9) Gl(r,s) =∫Γξlψr(ξ)ψs(ξ)ρ(ξ)dξ,Kl(i,j)=∫D√λlcl(x)∇ϕi(x)∇ϕj(x)dx,

. The right-hand side can be written as a tensor product of two vectors:

 f=g0⊗f0, (2.10)

where

 g0(r) =∫Γψr(ξ)ρ(ξ)dξ,r=1,…,Nξ, (2.11) f0(i) =∫Df(x)ϕi(x)dx,i=1,…,Nx.

Note that in the Galerkin system Eq. 2.7, the matrix is symmetric and positive definite. It is also blockwise sparse (see Fig. 2.1) due to the orthogonality of . The size of the linear system is in general very large (). For such a system it is suitable to use an iterative solver. Multigrid methods are among the most effective iterative solvers for the solution of discretized elliptic PDEs, capable of achieving convergence rates that are independent of the mesh size, with computational work growing only linearly with the problem size [8, 15].

### 2.2 Multigrid

In this subsection we discuss a geometric multigrid solver proposed in [4] for the solution of the stochastic Galerkin system Eq. 2.7. For this method, the mesh size varies for different grid levels, while the polynomial degree is held constant, i.e., the fine grid space and coarse grid space are defined as

 Whp=Sh⊗Tp,W2h,p=S2h⊗Tp, (2.12)

respectively. Then the prolongation and restriction operators are of the form

 P=I⊗P,R=I⊗PT, (2.13)

where is the same prolongation matrix as in the deterministic case. On the coarse grid we only need to construct matrices , and

 A2h=G0⊗K2h0+m∑l=1Gl⊗K2hl. (2.14)

The matrices are the same for all grid levels.

LABEL:alg:mg describes the complete multigrid method. In each iteration, we apply one multigrid cycle (Vcycle) for the residual equation

 Ac(i)=r(i)=f−Au(i) (2.15)

and update the solution and residual . The Vcycle function is called recursively. On the coarsest grid level () we form matrix and solve the linear system directly. The system is of order since where is a very small number on the coarsest grid. The smoothing function (Smooth) is based on a matrix splitting and stationary iteration

 us+1=us+Q−1(f−Aus), (2.16)

which we assume is convergent, i.e., the spectral radius The algorithm is run until the specified relative tolerance or maximum number of iterations is reached. It is shown in [4] that for , the convergence rate of this algorithm is independent of the mesh size , the number of random variables , and the polynomial degree .

\@float

algocf[htbp]     \end@float

## 3 Low-rank approximation

In this section we consider a technique designed to reduce computational effort, in terms of both time and memory use, using low-rank methods. We begin with the observation that the solution vector of the Galerkin system Eq. 2.7

 u=[u11,u21,…,uNx1,…,u1Nξ,u2Nξ,…,uNxNξ]T∈RNxNξ

can be restructured as a matrix

 U=mat(u)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝u11u12⋯u1Nξu21u22⋯u2Nξ⋮⋮⋱⋮uNx1uNx2⋯uNxNξ⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠∈RNx×Nξ. (3.1)

Then (2.7) is equivalent to a system in matrix format,

 A(U)=F, (3.2)

where

 A(U)=K0UGT0+m∑l=1KlUGTl, (3.3) F=mat(f)=mat(g0⊗f0)=f0gT0.

It has been shown in [3, 9] that the “matricized” version of the solution can be well approximated by a low-rank matrix when is large. Evidence of this can be seen in Fig. 3.1, which shows the singular values of the exact solution for the benchmark problem discussed in Section 4. In particular, the singular values decay exponentially, and low-rank approximate solutions can be obtained by dropping terms from the singular value decomposition corresponding to small singular values.

Now we use low-rank approximation in the multigrid solver for Eq. 3.2. Let be the th iterate, expressed in matricized format111In the sequel, we use and interchangeably to represent the equivalent vectorized or matricized quantities., and suppose is represented as the outer product of two rank- matrices, i.e., , where , . This factored form is convenient for implementation and can be readily used in basic matrix operations. For instance, the sum of two matrices gives

 V(i)1W(i)T1+V(i)2W(i)T2=[V(i)1,V(i)2][W(i)1,W(i)2]T. (3.4)

Similarly, can also be written as an outer product of two matrices:

 A(V(i)W(i)T) =(K0V(i))(G0W(i))T+m∑l=1(KlV(i))(GlW(i))T (3.5) = [K0V(i),K1V(i),…,KmV(i)][G0W(i),G1W(i),…,GmW(i)]T.

If are used to represent iterates in the multigrid solver and , then both memory and computational (matrix-vector products) costs can be reduced, from to . Note, however, that the ranks of the iterates may grow due to matrix additions. For example, in Eq. 3.5 the rank may increase from to in the worst case. A way to prevent this from happening, and also to keep costs low, is to truncate the iterates and force their ranks to remain low.

### 3.1 Low-rank truncation

Our truncation strategy is derived using an idea from [9]. Assume , , , and is truncated to rank with , , and . First, compute the QR factorization for both and ,

 ~V=Q~VR~V,~W=Q~WR~W, so ~X=Q~VR~VRT~WQT~W. (3.6)

The matrices and are of size . Next, compute a singular value decomposition (SVD) of the small matrix :

 R~VRT~W=^Vdiag(σ1,…,σ~k)^WT (3.7)

where are the singular values in descending order. We can truncate to a rank- matrix where is specified using either a relative criterion for singular values,

 √σ2k+1+⋯+σ2~k≤ϵrel% √σ21+⋯+σ2~k (3.8)

or an absolute one,

 k=max{k∣σk≥ϵabs}. (3.9)

Then the truncated matrices can be written in MATLAB notation as

 V=Q~V^V(:,1:k),W=Q~W^W(:,1:k)diag(σ1,…,σk). (3.10)

Note that the low-rank matrices obtained from Eq. 3.8 and Eq. 3.9 satisfy

 ∥X−~X∥F≤ϵrel∥~X∥F (3.11)

and

 ∥X−~X∥F≤ϵabs√~k−k, (3.12)

respectively. The right-hand side of Eq. 3.12 is bounded by since in general . The total cost of this computation is . In the case where becomes larger than , we compute instead a direct SVD for , which requires a matrix-matrix product to compute and an SVD, with smaller total cost .

### 3.2 Low-rank multigrid

The multigrid solver with low-rank truncation is given in LABEL:alg:low-rank. It uses truncation operators and , which are defined using a relative and an absolute criterion, respectively. In each iteration, one multigrid cycle (Vcycle) is applied to the residual equation. Since the overall magnitudes of the singular values of the correction matrix decrease as converges to the exact solution (see Fig. 3.2 for example), it is suitable to use a relative truncation tolerance inside the Vcycle function. In the smoothing function (Smooth), the iterate is truncated after each smoothing step using a relative criterion

 ∥Trel1(U)−U∥F≤ϵrel∥Fh−Ah(Uh0)∥F (3.13)

where , , and are arguments of the Vcycle function, and is the residual at the beginning of each V-cycle. In Line LABEL:lst:line:13, the residual is truncated via a more stringent relative criterion

 ∥Trel2(Rh)−Rh∥F≤ϵrelh∥Fh−Ah(Uh0)∥F (3.14)

where is the mesh size. In the main while loop, an absolute truncation criterion Eq. 3.9 with tolerance is used and all the singular values of below are dropped. The algorithm is terminated either when the largest singular value of the residual matrix is smaller than or when the multigrid solution reaches the specified accuracy.

\@float

algocf[htbp]     \end@float

Note that the post-smoothing is not explicitly required in LABEL:alg:low-rank and LABEL:alg:mg, and we include it just for sake of completeness. Also, in LABEL:alg:low-rank, if the smoothing operator has the form , then for any matrix with a low-rank factorization , application of the smoothing operator gives

 S(X)=S(VWT)=(S2V)(S1W)T, (3.15)

so that the result is again the outer product of two matrices of the same low rank. The prolongation and restriction operators Eq. 2.13 are implemented in a similar manner. Thus, the smoothing and grid-transfer operators do not affect the ranks of matricized quantities in LABEL:alg:low-rank.

### 3.3 Convergence analysis

In order to show that LABEL:alg:low-rank is convergent, we need to know how truncation affects the contraction of error. Consider the case of a two-grid algorithm for the linear system , where the coarse-grid solve is exact and no post-smoothing is done. Let be the coefficient matrix on the coarse grid, let be the error associated with , and let be the residual. It is shown in [4] that if no truncation is done, the error after a two-grid cycle becomes

 e(i+1)notrunc=(A−1−P¯A−1R)A(I−Q−1A)νe(i), (3.16)

and

 ∥e(i+1)notrunc∥A≤Cη(ν)∥e(i)∥A, (3.17)

where is the number of pre-smoothing steps, is a constant, and as . The proof consists of establishing the smoothing property

 ∥A(I−Q−1A)νy∥2≤η(ν)∥y∥A,∀y∈RNxNξ, (3.18)

and the approximation property

 ∥(A−1−P¯A−1R)y∥A≤C∥y∥2,∀y∈RNxNξ, (3.19)

and applying these bounds to Eq. 3.16.

Now we derive an error bound for LABEL:alg:low-rank. The result is presented in two steps. First, we consider the Vcycle function only; the following lemma shows the effect of the relative truncations defined in Eqs. 3.14 and 3.13.

###### Lemma 1

Let and let be the associated error. Then

 ∥e(i+1)∥A≤C1(ν)∥e(i)∥A, (3.20)

where, for small enough and large enough , independent of the mesh size .

###### Proof

For , let be the quantity computed after application of the smoothing operator at step before truncation, and let be the modification obtained from truncation by of Eq. 3.13. For example,

 ~u(i)1=u(i)+Q−1(f−Au(i)),u(i)1=Trel1(~u(i)1). (3.21)

Denote the associated error as . From Eq. 3.13, we have

 (3.22)

Similarly, after smoothing steps,

 e(i)ν =(I−Q−1A)νe(i)+Δ(i)ν (3.23) =(I−Q−1A)νe(i)+(I−Q−1A)ν−1δ(i)1+⋯+(I−Q−1A)δ(i)ν−1+δ(i)ν,

where

 ∥δ(i)s∥2≤ϵrel∥r(i)∥2,s=1,…,ν. (3.24)

In Line LABEL:lst:line:13 of LABEL:alg:low-rank, the residual is truncated to via Eq. 3.14, so that

 ∥r(i)ν−~r(i)ν∥2≤ϵrelh∥r(i)∥2. (3.25)

Let . Referring to Eqs. 3.23 and 3.16, we can write the error associated with as

 e(i+1) =e(i)ν−P¯A−1Rr(i)ν (3.26) =(I−P¯A−1RA)e(i)ν−P¯A−1Rτ(i) =e(i+1)notrunc+(A−1−P¯A−1R)AΔ(i)ν−P¯A−1Rτ(i) =e(i+1)notrunc+(A−1−P¯A−1R)(AΔ(i)ν+τ(i))−A−1τ(i).

Applying the approximation property Eq. 3.19 gives

 ∥(A−1−P¯A−1R)(AΔ(i)ν+τ(i))∥A≤C(∥AΔ(i)ν∥2+∥τ(i)∥2). (3.27)

Using the fact that for any matrix ,

 supy≠0∥By∥A∥y∥A=supy≠0∥A1/2By∥2∥A1/2y∥2=supz≠0∥A1/2BA−1/2z∥2∥z∥2=∥A1/2BA−1/2∥2, (3.28)

we get

 ∥A(I−Q−1A)ν−sδ(i)s∥2 ≤∥A1/2∥2∥(I−Q−1A)ν−sδ(i)s∥A (3.29) ≤∥A1/2∥2∥A1/2(I−Q−1A)ν−sA−1/2∥2∥δ(i)s∥A ≤ρ(I−Q−1A)ν−s∥A1/2∥22∥δ(i)s∥2

where is the spectral radius. We have used the fact that is a symmetric matrix (assuming is symmetric). Define . Then Eqs. 3.25 and 3.24 imply that

 ∥AΔ(i)ν∥2+∥τ(i)∥2 ≤ϵrel(d1(ν)+h)∥r(i)∥2 (3.30) ≤ϵrel(d1(ν)+h)∥A1/2∥2∥e(i)∥A.

On the other hand,

 ∥A−1τ(i)∥A=(A−1τ(i),τ(i))1/2 ≤∥A−1∥1/22∥τ(i)∥2 (3.31) ≤ϵrelh∥A−1∥1/22∥r(i)∥2 ≤ϵrelh∥A−1∥1/22∥A1/2∥2∥e(i)∥A.

Combining Eqs. 3.31, 3.30, 3.27, 3.26 and 3.17, we conclude that

 ∥e(i+1)∥A≤C1(ν)∥e(i)∥A (3.32)

where

 C1(ν)=Cη(ν)+ϵrel(C(d1(ν)+h)+h∥A−1∥1/22)∥A1/2∥2. (3.33)

Note that , is bounded by a constant, and is of order [14]. Thus, for small enough and large enough , is bounded below 1 independent of .

Next, we adjust this argument by considering the effect of the absolute truncations in the main while loop. In LABEL:alg:low-rank, the Vcycle is used for the residual equation, and the updated solution and residual are truncated to and , respectively, using an absolute truncation criterion as in Eq. 3.9. Thus, at the th iteration (), the residual passed to the Vcycle function is in fact a perturbed residual, i.e.,

 r(i)=~r(i)+β=Ae(i)+β,where ∥β∥2≤√Nξϵabs. (3.34)

It follows that in the first smoothing step,

 ~u(i)1=u(i)+Q−1(f−Au(i)+β),u(i)1=Trel1(~u(i)1), (3.35)

and this introduces an extra term in (see Eq. 3.23),

 Δ(i)ν=(I−Q−1A)ν−1δ(i)1+⋯+(I−Q−1A)δ(i)ν−1+δ(i)ν−(I−Q−1A)ν−1Q−1β. (3.36)

As in the derivation of Eq. 3.29, we have

 ∥A(I−Q−1A)ν−1Q−1β∥2≤ρ(I−Q−1A)ν−1∥A1/2∥22∥Q−1∥2∥β∥2. (3.37)

In the case of a damped Jacobi smoother (see Eq. 4.4), is bounded by a constant. Denote . Also note that . Then Eqs. 3.31 and 3.30 are modified to

 ∥AΔ(i)ν∥2+∥τ(i)∥2 (3.38) ≤ϵrel(d1(ν)+h)∥r(i)∥2+d2(ν)∥β∥2 ≤ϵrel(d1(ν)+h)∥A1/2∥2∥e(i)∥A+(d2(ν)+ϵrel(d1(ν)+h))∥β∥2,

and

 ∥A−1τ(i)∥A≤ϵrelh∥A−1∥1/22∥A1/2∥2∥e(i)∥A+ϵrelh∥A−1∥1/22∥β∥. (3.39)

As we truncate the updated solution , we have

 u(i+1)=~u(i+1)+γ,where ∥γ∥2≤√Nξϵabs. (3.40)

Let

 C2(ν)=Cd2(ν)+ϵrel(C(d1(ν)+h)+h∥A−1∥1/22)+∥A1/2∥2. (3.41)

From Eqs. 3.41, 3.40, 3.39 and 3.38, we conclude with the following theorem:

###### Theorem 3.1

Let denote the error at the th iteration of LABEL:alg:low-rank. Then

 ∥e(i+1)∥A≤C1(ν)∥e(i)∥A+C2(ν)√Nξϵabs, (3.42)

where for large enough and small enough , and is bounded by a constant. Also, Eq. 3.42 implies that

 ∥e(i)∥A≤Ci1(ν)∥e(0)∥A+1−Ci1(ν)1−C1(ν)C2(ν)√Nξϵabs, (3.43)

i.e., the -norm of the error for the low-rank multigrid solution at the th iteration is bounded by . Thus, LABEL:alg:low-rank converges until the -norm of the error becomes as small as .

It can be shown that the same result holds if post-smoothing is used. Also, the convergence of full (recursive) multigrid with these truncation operations can be established following an inductive argument analogous to that in the deterministic case (see, e.g., [5, 8]). Besides, in LABEL:alg:low-rank, the truncation on imposes a stopping criterion, i.e.,

 ∥~r(i+1)∥2 ≤∥~r(i+1)−r(i+1)∥2+∥r(i+1)∥2 (3.44) ≤√Nξϵabs+tol∗r0.

In Section 4 we will vary the value of and see how the low-rank multigrid solver works compared with LABEL:alg:mg where no truncation is done.

###### Remark 1

It is shown in [14] that for Eq. 2.7, with constant mean and standard deviation ,

 ∥A∥2=α(c0+σCmaxp+1m∑l=1√λl∥cl(x)∥∞), (3.45)

where is the maximal root of an orthogonal polynomial of degree , and is a constant independent of , , and . If Legendre polynomials on the interval are used, . Since both and in Theorem 3.1 are related to , the convergence rate of LABEL:alg:low-rank will depend on . However, if the eigenvalues decay fast, this dependence is negligable.

###### Remark 2

If instead a relative truncation is used in the while loop so that

 r(i+1)=~r(i+1)+β=Ae(i+1)+β,where ∥β∥2≤ϵrel∥~r(i+1)∥, (3.46)

then a similar convergence result can be derived, and the algorithm stops when

 ∥~r(i+1)∥2≤tol∗r1−ϵ% rel. (3.47)

However, the relative truncation in general results in a larger rank for , and the improvement in efficiency will be less significant.