\documentclass{report}
\usepackage{amsmath,amssymb}
\setlength{\parindent}{0mm}
\setlength{\parskip}{1em}
\begin{document}
\begin{center}
\rule{6in}{1pt} \
{\large Christophe Audouze \\
{\bf 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
\\
{\tt c.audouze@soton.ac.uk}\\
Prasanth B. Nair\end{center}

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)$ ({\em 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 para\-metrized) 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),~~\text{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 $||x^{*}(\xi) -
\sum_{i=0}^{P}x_{i}\varphi_{i}(\xi)||_{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
$$
||r(\xi)||_2^2 = \sum_{i=0}^{P} r_{i}^Tr_{i} \langle \varphi_i^2\rangle +
\sum_{i=P+1}^{P_r} r_i^{T} r_{i} \langle \varphi_i^2
\rangle,\label{true_r} \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 {\em 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.

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

\begin{center}
{\bf REFERENCES}
\end{center}

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

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

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

\noindent
[4] Y. Saad, M. H. Schultz, {\it{``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.


\end{document}
