Solvers and precondtioners based on Gauss-Seidel and Jacobi algorithms for non-symmetric stochastic Galerkin system of equations

To account for uncertainties in model, data and parameters, physical phenomena are often modeled as partial differential equations (PDEs) with random coefficients. While Monte Carlo sampling methods remain the standard for providing a probabilistic characterization to the solution of these problems, there is growing need for approaches that address the computational challenge associated with statistical sampling of large scale computational models. More recently, methods based on stochastic Galerkin approximations [8, 3] have been increasingly explored as an alternative to Monte Carlo for solving these problems because of their computational and analytical properties, as well as for their ability to provide an approximation to the deterministic mapping that exists between random input parameters and random solution.

Stochastic finite element methods based on intrusive stochastic Galerkin methods ([14, 12, 4, 10]) and non-intrusive stochastic collocation methods ([16, 19, 11, 2]) have gained popularity in recent years. Here the intrusive and the non-intrusive methods are defined as methods which require significant changes to software already available for solving deterministic PDEs and methods that can use this legacy software without modifications, respectively. Both methods exploit solution regularity to improve on the convergence rates of Monte Carlo methods. The first approach translates the stochastic PDE into a coupled set of deterministic PDEs while the second samples the stochastic PDE at a predetermined set of collocation points, resulting in a set of uncoupled deterministic PDEs. The solution at these collocation points is then used to interpolate the solution in the entire random input domain. Extending legacy software to support stochastic collocation methods is simpler than supporting stochastic Galerkin methods (SGMs). Moreover, intrusive SGMs require specialized linear solvers. However, the resulting set of PDEs in the stochastic Galerkin system is much smaller in number than that in the collocation method. For a canonical random diffusion problem, it is shown [5] that SGM using iterative Krylov-based linear solvers and mean-based preconditioning [13] is more efficient than the non-intrusive sparse grid collocation method.

While the stochastic Galerkin method is often considered to be a fully intrusive method, there are in fact a variety of solver approaches for the stochastic Galerkin method with varying level of intrusiveness. In the present work, the less intrusive Gauss-Seidel and Jacobi mean-based algorithms are explored to solve stochastic Galerkin system of equations. We consider these methods to be less intrusive than the Krylov-based methods as they allow reuse of existing deterministic solvers. Moreover preconditioning techniques for Krylov-based methods based on Gauss-Seidel and Jacobi ideas are also explored and compared to traditional mean-based preconditioning. Various iterative methods and preconditioners based on matrix splitting methods to solve stochastic Galerkin system of equations are compared in [14]. The $C_{i}$-splitting methods discussed in section 3.1.3 in [14] are similar to the Gauss-Seidel and Jacobi algorithms discussed in this paper. In [14], it was reported that the mean-based preconditioner was more efficient in terms of computational time than the other $C_{i}$-splitting methods considered. A symmetric Gauss-Seidel preconditioner was used in [14] which involves one forward and one backward Gauss-Seidel iteration. In the present paper, we use only one forward Gauss-Seidel iteration for preconditioning and hence the computational cost of this preconditioner is half that of the symmetric Gauss-Seidel preconditioner. Although, these solvers and preconditioners can be used to solve any stochastic Galerkin system of equations, we focus here on solving non-symmetric stochastic Galerkin system of equations. Gauss-Seidel and Jacobi algorithms are carefully formulated such that the number of matrix-vector products are minimized. Unlike the mean-based preconditioner, which has information from the mean stiffness matrix only, preconditioners based on Gauss-Seidel and Jacobi algorithms make use of higher order information of the stochastic stiffness matrix. Another preconditioner that makes use of higher order information is the Kronecker product preconditioner defined in [17]. This preconditioner has the Kronecker product structure of the Galerkin stiffness matrix and is constructed to be close to the stochastic stiffness matrix. All these techniques are then compared by solving a non-symmetric test problem with random diffusion coefficient. These comparisons demonstrate a trade-off in computational cost versus intrusiveness with the Krylov-based methods using an approximate Gauss-Seidel or Jacobi mean preconditioner being the most efficient. Main contribution in this paper is to adapt Gauss-Seidel and Jacobi algorithms such that the legacy software can be reused to solve stochastic partial differential equations and to reduce the computational cost of Krylov based iterative methods by using these algorithms as preconditioners in which the number of matrix-vector products is reduced by minimizing the duplicated computations.

This paper is organized as follows. Two models of the input random field are developed in section 3 which dictate very different behavior for the stochastic solution methods considered next. Section 4 describes the stochastic Galerkin method. Various solver and preconditioning methods for solving the stochastic Galerkin system of equations are introduced in section 5. In section 6, a test problem governed by the advection-diffusion equation with a random diffusion coefficient is formulated. In section 7, numerical experiments are carried out to compare the efficiency of various solver and preconditioning methods that have been introduced in section 5. Finally section 8 provides the concluding remarks.

Let $D$ be an open subset of $\mathbb{R}^{n}$ and $(\Omega,\Sigma,P)$ be a complete probability space with sample space $\Omega$, $\sigma$-algebra $\Sigma$ and probability measure $P$. We are interested in studying the following stochastic partial differential equation: find a random field, $u(x,\omega):D\times\Omega\rightarrow\mathbb{R}$ such that the following holds $P$-almost surely ($P$-a.s.):

subject to the boundary condition

where $\mathcal{L}$ is a differential operator and $\mathcal{B}$ is a boundary operator.

Uncertainty in stochastic PDEs often arises by treating coefficients in the differential operator $\mathcal{L}(x,\omega;u)$ as random fields. For computational purposes, each stochastic coefficient $a(x,\omega)$ must be discretized in both spatial and stochastic domains. To this end, it is often approximated with a truncated series expansion that separates the spatial variable $x$ from the stochastic variable $\omega$ resulting in a representation by a finite number of random variables. In the present problem, two cases of random field models are considered. In the first case, the random field is assumed to be charcaterized by a truncated Karhunen-Loève expansion obtained by specifying the covariance function and a product uniform measure on the associated random variables. In the second case, the random field is assumed to have a log-normal distribution, that is $a(x,\omega)=\exp{(g(x,\omega))}$ where $g(x,\omega)$ is a Gaussian random field, and is approximated by a truncated polynomial chaos expansion.

Let $H^{1}_{0}(D)$ be the subspace of the Sobolev space $H^{1}(D)$ that vanishes on the boundary $\partial D$ and is equipped with the norm $\|u\|_{H^{1}_{0}(D)}=[\int_{D}|\nabla u|^{2}dx]^{\frac{1}{2}}$. Problem (1) can then be written in the following equivalent variational form [9]: find $u\in H^{1}_{0}(D)\otimes L_{2}(\Omega)$ such that

where $b(u,v)$ is a continuous and coercive bilinear form and $l(v)$ is a continuous bounded linear functional. In the stochastic Galerkin method, we seek the solution of the variational problem (10) in a tensor product space $X_{h}\otimes Y_{p}$, where, $X_{h}\subset H^{1}_{0}(D)$ is finite dimensional space of continuous polynomials corresponding to the spatial discretization of $D$ and $Y_{p}\subset L_{2}(\Omega)$ is the space of random variables spanned by polynomial chaos [8] of order up to $p$. Then the finite dimensional approximation $u_{X_{h}Y_{p}}(x,\omega)$ of the exact solution $u(x,\omega)$ on the tensor product space $X_{h}\otimes Y_{p}$ is given as the solution to

In equation (11) the input random field $a(x,\omega)$ in the bilinear form $b(u_{X_{h}Y_{p}},v)$ can be approximated using either a K-L expansion or a polynomial chaos expansion depending on a choice of model for the random field. The finite dimensional approximation $u_{X_{h}Y_{p}}(x,\omega)$ is represented with truncated polynomial chaos expansion, where the multidimensional polynomial chaos are orthogonal with respect to the probability measure of the underlying random variables. The resulting set of coupled PDEs are then discretized using standard techniques such as the finite element or finite difference methods. In the present problem, the set of coupled PDEs are discretized using the finite element method and the resulting system of linear equations can be written as

where $f_{k}=E\{f(x,\boldsymbol{\xi})\psi_{k}\}$, $c_{ijk}=E\{\xi_{i}\psi_{j}\psi_{k}\}$ and $\hat{P}=M$ when $a(x,\omega)$ is approximated by a truncated K-L expansion, or $c_{ijk}=E\{\psi_{i}\psi_{j}\psi_{k}\}$ and $\hat{P}=\hat{N}_{\xi}$ when $a(x,\omega)$ is approximated by a polynomial chaos expansion. Here $\{K_{i}\in\mathbb{R}^{N_{x}\times N_{x}}\}_{i=0}^{\hat{P}}$ are the polynomial chaos coefficients of the stiffness matrix (section (3.4) of [13]) and $\{u_{j}\in\mathbb{R}^{N_{x}}\}_{j=0}^{N_{\xi}}$ are the polynomial chaos coefficients of the discrete solution vector

Equation (12) can be written in the form of a global stochastic stiffness matrix of size $((N_{\xi}+1)\times N_{x})$ by $((N_{\xi}+1)\times N_{x})$ as

where $K^{j,k}=\sum_{i=0}^{\hat{P}}c_{ijk}K_{i}$. We will denote this system as $\bar{K}\bar{u}=\bar{f}$. In practice it is prohibitive to assemble and store the global stochastic stiffness matrix in this form, rather each block of the stochastic stiffness matrix can be computed from the $\{K_{i}\}$ when needed. The stochastic stiffness matrix is block sparse, which means that some of the off-diagonal blocks are zero matrices, because of the fact that the $c_{ijk}$ defined above vanishes for certain combination of $i,j$ and $k$. Unlike the off-diagonal blocks, the diagonal blocks have contributions from the mean striffness matrix $K_{0}$ and thus dominate over the off-diagonals. These properties of block sparsity and diagonal dominance are exploited to develop solvers and preconditioners based on Jacobi and Gauss-Seidel algorithms.

In this section, various solver techniques and preconditioning methods for solving the linear algebraic equations arising from stochastic Galerkin discretizations (12) are described. The solver methods discussed are: relaxation methods, namely, a Jacobi mean method and a Gauss-Seidel mean method, and Krylov-based iterative methods [15]. Also various stochastic preconditioners used to accelerate convergence of the Krylov methods are discussed, including mean-based [13], Gauss-Seidel mean, approximate Gauss-Seidel mean, approximate Jacobi mean and Kronecker product [17] preconditioners. The relaxation schemes can be viewed as fixed point iterations on a preconditioned system [15]. In Jacobi and Gauss-Seidel methods, mean splitting is used rather than traditional diagonal block splitting as it allows use of the same mean matrix $K_{0}$ for all inner deterministic solves (and thus reuse of the preconditioner $P_{0}\approx K_{0}$).

In this work, a stochastic steady state advection-diffusion equation [20] with Dirichlet boundary conditions is used as a test problem for various solution methods. Assume $a(x,\omega):D\times\Omega\rightarrow\mathbb{R}$ to be a random field that is bounded and strictly positive, that is,

For the steady-state advection-diffusion equation, we wish to compute a random field $u(x,\omega):D\times\Omega\rightarrow\mathbb{R}$, such that the following holds $P$-almost surely ($P$-a.s.):

where $\vec{w}=[w_{x},w_{y}]^{T}$ is the advection direction, $w_{x},w_{y}>0$.

Problem (20) can then be written in the following equivalent variational form [9]: find $u\in H^{1}_{0}(D)\otimes L_{2}(\Omega)$ such that

where $b(u,v)$ is the continuous and coercive (from assumption (19)) bilinear form given by

and $l(v)$ is the continuous bounded linear functional given by

From the Lax-Milgram lemma, equation (22) has unique a solution in $H^{1}_{0}(D)\otimes L_{2}(\Omega)$.

To compare the performance of different solvers and preconditioners discussed above, a 2-D stochastic diffusion equation and a 2-D stochastic advection-diffusion equation presented in section 6 are solved using the stochastic Galerkin method described in section 4. In the above two problems, the diffusion coefficient is modeled as both a random field discretized using a truncated K-L expansion with uniform random variables (section 3.1) and a log-normal random field discretized using a truncated polynomial chaos expansion (section 3.2). In the first case, the orthogonal polynomials used in the stochastic Galerkin method are tensor products of 1-D Legendre polynomials and in the second case tensor products of Hermite polynomials are used. For simplicity a constant unit force $f(x,\omega)=1$ is used as the right-hand-side in equation (20). The spatial dimensions are discretized using standard finite element mesh with linear quadrilateral elements. In advection-diffusion equation, the parameter, $\vec{w}=[1,1]^{T}$. The corresponding stochastic Galerkin linear system is constructed using the Stokhos and Epetra packages in Trilinos. For the Jacobi solver, Gauss-Seidel solver and Gauss-Seidel preconditioner the non-symmetric linear systems obtained from the discretization of stochastic advection-diffusion equation are solved via multi-grid preconditioned GMRES provided by the AztecOO and ML packages in Trilinos.

For the numerical comparisons, we chose to discretize the domain $D=[-0.5,0.5]\times[-0.5,0.5]$ into a $64\times 64$ grid resulting in a total number of nodes, $N_{x}=4096$. In figure (1) the solution time for the stochastic Galerkin method, scaled by the deterministic solution time at the mean of the random field, is compared for different mesh sizes ($32\times 32$, $64\times 64$, $96\times 96$ and $128\times 128$) for the diffusion problem demonstrating that the solution time does not depend strongly on the mesh size (as is to be expected for the multi-grid preconditioner). The scaled solution time for these solvers and preconditioning techniques as a function of the standard deviation of the input random field, stochastic dimension, and polynomial order are then tabulated in Tables 1–6. In the tables, the number of Krylov iterations for the aforementioned preconditioners and iterations for Gauss-Seidel and Jacobi solvers are provided in parentheses. In the tables, MB, AGS, AJ, GS and KP are the mean-based, approximate Gauss-Seidel, approximate Jacobi, Gauss-Seidel and Kronecker-product preconditioners respectively. GS in the solution methods refers to the Gauss-Seidel mean algorithm (Algorithm-2). “Jacobi” refers to the Jacobi mean algorithm (Algorithm-1). The solution tolerance for all of the stochastic Galerkin solvers is $1e^{-12}$. For the Gauss-Seidel and Jacobi solvers, the inner solver tolerance is $3e^{-13}$. All the computations are performed using a single core of an 8 core, Intel Xeon machine with 2.66 GHz and 16GB Memory.

Figures (2) and (3) show the plots of relative residual error vs iteration count for the stochastic Galerkin system with stochastic dimension 4 and polynomial order 4 and standard deviation 0.1. It can be observed that the matrix-free Krylov solver with the Gauss-Seidel preconditioner takes the least number of iterations in case of uniform random field and Gauss-Seidel solver in case of log-normal random field, whereas the Jacobi solver takes highest number of iterations in both cases for a given tolerance. However in terms of solution time, the matrix-free Krylov solver with the approximate Gauss-Seidel preconditioner is the most efficient compared to all other stochastic Galerkin solvers. Comparison in terms of iteration count alone is misleading to evaluate the computational cost because, in each iteration the cost of preconditioner as observed in the case of Gauss-Seidel preconditioner could be very high resulting in higher computational cost even with small number of iterations. Hence we also compare the solution time for all preconditioners and solvers.

In the tables, “Div” means diverged. In the case of diffusion coefficient modeled with uniform random random variables, with small variance ($\sigma=0.1$), it can be observed from Tables 1 and 2 that more intrusive Krylov-based stochastic Galerkin solvers are more efficient than less intrusive Gauss-Seidel and Jacobi solvers. Moreover the approximate Gauss-Seidel and Jacobi preconditioners are a significant improvement over the traditional mean-based approach. However as the variance of the random field increases, we see from Table 3 that the Gauss-Seidel and Jacobi solvers suffer considerably, whereas the the Krylov-based approaches (excluding the Gauss-Seidel preconditioner) still perform quite well. This is not unexpected, as the operator becomes more indefinite as the variance increases. Since the efficiency of the preconditioners and solvers based on Gauss-Seidel and Jacobi algorithms depends heavily on the efficiency of the deterministic solver, a good preconditioner for the mean stiffness matrix will significantly improve the above mentioned stochastic Galerkin preconditioners and solvers.

In the case of the log-normal random field, we can see from Tables 4 and 5 that Gauss-Seidel and Jacobi solvers have not performed well in terms of solution time. The Gauss-Seidel solver is comparable to the Krylov solver with the AGS preconditioner only in case of log-normal random field and at higher polynomial chaos order. For higher variance of the random field, we see from Table 6 that the Jacobi solver diverges. This problem might be addressed to some extent by using the true diagonal matrix $K^{k,k}=\sum_{i=0}^{M}c_{ikk}K_{i}$ from global stochastic stiffness matrix as the left-hand-side in the Jacobi solver and preconditioner instead of the mean matrix $K_{0}$. From the Krylov iteration count provided in the parentheses of the tables, we can observe that the preconditioners and solvers based on Gauss-Seidel and Jacobi are robust with respect to stochastic dimension and polynomial order, however they are not robust with respect to the variance of the input random field. The numerical results clearly show that the approximate Gauss-Seidel preconditioner performs betters than the mean-based preconditioner. It was reported in [14] that the mean-based preconditioner was more efficient in terms of computational time than the Gauss-Seidel preconditioner. It could be due to the use of symmetric Gauss-Seidel preconditioner and duplicated matrix-vector products in the Gauss-Seidel preconditioner. In our work, we use non-symmetric Gauss-Seidel preconditioner which is half expensive as that of its symmetric version and we also minimize the duplicated matrix-vector products in Gauss-Seidel algorithm (Algorithm-2).

In this work, various preconditioners for Krylov-based methods and solver methods based on Gauss-Seidel and Jacobi method are introduced. Results are compared with Krylov based methods (GMRES) with mean-based preconditioning. The less intrusive Gauss-Seidel/Jacobi approaches did not perform well. However these solvers can be used when legacy software has to be used to solve SPDEs with the stochastic Galerkin descritization. The use of approximate Gauss-Seidel or Jacobi preconditioners yields a significant improvement over traditional mean-based preconditioning. The Kronecker product preconditioner proposed in [17] was in between mean-based preconditioner and approximate Gauss-Seidel preconditioner. Block and recycled Krylov methods present additional promising alternatives for improve the efficiency of the Jacobi and Gauss-Seidel solvers.

Acknowledgements

This work was partially supported by the US Department of Energy through the NNSA Advanced Scientific Computing and Office of Science Advanced Scientific Computing Research programs.