Christophe Audouze
GMRES algorithm for nonsymmetric linear random algebraic equations

Computational Engineering and Design Group
School of Engineering Sciences
University of Southampton
Southampton SO17 1BJ
United Kingdom
c.audouze@soton.ac.uk
Prasanth B. Nair

This paper is concerned with Krylov methods for solving linear random algebraic equations of the form

$${A}(\xi)x(\xi)={b}(\xi) \eqno{(1)}$$

where $\xi\in \mathbb{R}^M$ denotes a vector of random variables with specified probability density function $\mathcal{P}(\xi): \mathbb{R}^{M} \mapsto \mathbb{R}^{+}$, $A(\xi): \mathbb{R}^M \mapsto \mathbb{R}^{n\times n}$ is an invertible matrix for all values of $\xi$ drawn from $\mathcal{P}(\xi)$ (a.s. invertible) and $b(\xi): \mathbb{R}^{M} \mapsto \mathbb{R}^{n}$. $x(\xi): \mathbb{R}^M \mapsto \mathbb{R}^{n}$ is the solution vector whose statistics are sought to be computed.

Our focus is on the particular case when $A(\xi)$ and $b(\xi)$ are provided in a polynomial chaos (PC) basis, i.e., $A(\xi) = \sum_{i=0}^{P_{A}} A_{i} \varphi_{i}(\xi)$ and $b(\xi) = \sum_{i=0}^{P_{b}} b_{i} \varphi_{i}(\xi)$, where $\varphi_{i}(\xi)$ are multivariate Hermite, Legendre or general monomial basis functions that are constructed to be orthonormal with respect to $\mathcal{P}(\xi)$ ( $\langle \varphi_{i} \varphi_{j} \rangle = \int \varphi_{i} \varphi_{j} {\rm d}\mathcal{P}=\delta_{ij}$). Parametrized linear random algebraic equations of this form are routinely encountered in numerical solution of stochastic (or randomly parametrized) partial differential equations (SPDEs) [1-3]. A well known method for solving this class of problems is the Ghanem-Spanos projection scheme [1], wherein the solution is approximated as $\widehat{x}_{P}(\xi) = \sum_{i=0}^{P} x_{i} \varphi_{i}(\xi)$, where $x_{i} \in \mathbb{R}^{n}$ are undetermined coefficient vectors. The undetermined coefficient vectors are estimated via Galerkin projection as follows:

$$A(\xi) \sum_{i=0}^{P} x_{i} \varphi_{i}(\xi) - b(\xi)\perp \varphi_{k}(\xi),~~$$for$$~k=0,1,\ldots,P. \eqno{(2)}$$

Using the standard Hilbert space inner product for random vectors, the above conditions lead to a structured deterministic matrix system of equations for the undetermined coefficient vectors [2]. The sparsity structure of these equations can be exploited to design efficient preconditioned Krylov subspace algorithms.

For the case when $A(\xi)$ is symmetric positive definite (SPD), it can be easily shown that the above conditions are equivalent to directly minimizing the $A-$norm error $\vert\vert x^{*}(\xi) - \sum_{i=0}^{P}x_{i}\varphi_{i}(\xi)\vert\vert _{A(\xi)}$, where $x^{*}(\xi)$ is the exact solution of (1). In order to study which error norm is minimized for the case when $A(\xi)$ is nonsymmetric, it is instructive to write the residual error corresponding to the approximation $\widehat{x}_{P}(\xi)$ in the form

$$\vert\vert r(\xi)\vert\vert _2^2 = \sum_{i=0}^{P} r_{i}^Tr_{i} \l... ...ngle + \sum_{i=P+1}^{P_r} r_i^{T} r_{i} \langle \varphi_i^2 \rangle,\eqno{(3)}$$

where $r_{i}$ is the $i$th term in the PC expansion of the residual error $r(\xi) = b(\xi) - A(\xi) \widehat{x}_{P}(\xi)$, i.e., $r_i = b_{i} - \frac{1}{\langle \varphi_i^2 \rangle} \sum_{j=0}^{P_{A}}\sum_{l=0}^{P} \langle \varphi_{i}\varphi_{j}\varphi_{l} \rangle A_{j}x_{l}$. It can be shown that when the Ghanem-Spanos method is applied to solve nonsymmetric linear random algebraic equations, only the first term (i.e., the residual error projected onto the first $P$ basis functions) is minimized [3]. This can potentially lead to numerical instabilities due to increase in the true residual when the number of PC expansion terms ($P$) is increased. One way to circumvent this problem would be to employ a Petrov-Galerkin projection scheme, i.e., orthogonalize the residual with respect to $A(\xi)\varphi_{k}(\xi)$ in (2) - this would ensure minimization of the true residual error norm. However, this condition essentially means we now need to solve the normal equations $A(\xi)^{T}A(\xi) x(\xi) = A(\xi)^{T} b(\xi)$, which makes this approach numerically unattractive.

In the present work, we propose a numerically stable Krylov method based on the GMRES algorithm [4] (which we refer to as the sGMRES algorithm) for minimizing the stochastic residual error norm corresponding to the $P$ term approximation $x_{P}(\xi)$.The sGMRES algorithm combines the ideas developed in [3] in the context of conjugate gradient methods with the standard GMRES algorithm for deterministic nonsymmetric linear algebraic equations. The key idea is to iteratively minimize the residual error norm over a Krylov subspace whose basis are constructed using the Arnoldi procedure. We present a general derivation of the sGMRES algorithm, while pointing out the similarities and differences with the deterministic GMRES algorithm. We also provide a theoretical convergence analysis of the sGMRES algorithm and show how restarting procedures in conjunction with function decomposition schemes can be employed to significantly reduce the computational cost and memory requirements for high-dimensional problems ( $M\sim \mathcal{O}(100)$).

We present numerical results for two different test-cases. As a first validation, we consider a structured random perturbation of a fixed nonsymmetric matrix. We then consider a more realistic situation where the linear equations arise from semi-discretization of the stochastic convection-diffusion equation. Detailed comparison studies are made with the standard Ghanem-Spanos method to illustrate the advantages of the sGMRES algorithm.

Acknowledgements: This research is supported by the United Kingdom Engineering and Physical Sciences Research Council (EPSRC) Grant No. EP/F006802/1.

REFERENCES

[1] R. Ghanem and P. Spanos, ``Stochastic finite elements: A spectral approach'', Springer Verlag, New York, 1991.

[2] C. E. Powell and H. C. Elman, ``Block-diagonal preconditioning for spectral stochastic finite-element systems'', IMA Journal of Numerical Analysis 2009 29(2):350-375.

[3] P. Hakansson, P. B. Nair, ``Conjugate gradient methods for parametrized linear random algebraic equations'', submitted.

[4] Y. Saad, M. H. Schultz, ``GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems'', SIAM Journal on Scientific and Statistical Computing, Vol. 7, Issue 3 (1986), 856-869.