Asymptotic convergence of spectral inverse iterations for stochastic eigenvalue problems

TL;DR

This paper analyzes the asymptotic convergence of spectral inverse iteration algorithms for computing eigenpairs of elliptic operators with random coefficients, combining stochastic and spatial approximations.

Contribution

It provides an error analysis for the convergence of spectral inverse iteration and subspace iteration in stochastic eigenvalue problems, supported by numerical experiments.

Findings

Spectral inverse iteration converges asymptotically to the smallest eigenvalue and eigenvector.

Numerical experiments confirm the theoretical convergence results.

The algorithms effectively handle eigenvalue crossings in stochastic parameter spaces.

Abstract

We consider and analyze applying a spectral inverse iteration algorithm and its subspace iteration variant for computing eigenpairs of an elliptic operator with random coefficients. With these iterative algorithms the solution is sought from a finite dimensional space formed as the tensor product of the approximation space for the underlying stochastic function space, and the approximation space for the underlying spatial function space. Sparse polynomial approximation is employed to obtain the first one, while classical finite elements are employed to obtain the latter. An error analysis is presented for the asymptotic convergence of the spectral inverse iteration to the smallest eigenvalue and the associated eigenvector of the problem. A series of detailed numerical experiments supports the conclusions of this analysis. Numerical experiments are also presented for the spectral subspace iteration, and convergence of the algorithm is observed in an example case, where the eigenvalues cross within the parameter space. The outputs of both algorithms are verified by comparing to solutions obtained by a sparse stochastic collocation method.

Equations

$$||u-u_{h,\mathcal{A},k}|{|}_{L_{\nu }^2(\Gamma )\otimes L^2(D)}\lesssim h^{1+l}+(\#{\mathcal{A}}_{ϵ}{)}^{-r}+{\lambda }_{1/2}^k$$

$$||\mu -{\mu }_{h,\mathcal{A},k}|{|}_{L_{\nu }^2(\Gamma )}\lesssim h^{2l}+(\#{\mathcal{A}}_{ϵ}{)}^{-r}+{\lambda }_{1/2}^k,$$

$$\mathbb{E}[v]={\int }_{\Omega }v(\omega )d\mathcal{P}(w)$$

$$\mathcal{P}\left(\omega \in \Omega :a_{\min }\le \underset{x\in D}{\mathrm{essinf}}a(x,\omega )\le \underset{x\in D}{\mathrm{esssup}}a(x,\omega )\le a_{\max }\right)=1.$$

$$\left\{\begin{array}[]{ll}-\nabla\cdot(a(\cdot,\omega)\nabla u(\cdot,\omega))=\mu(\omega)u(\cdot,\omega)&\textrm{ in }D\\ u(\cdot,\omega)=0&\textrm{ on }\partial D\end{array}\right.$$

$$a(x,\omega )=a_0(x)+\sum _{m=1}^{\infty }{a}_m(x)y_m(\omega ),$$

$$\underset{x\in D}{\mathrm{essinf}}{a}_0(x)>\sum _{m=1}^{\infty }||a_m|{|}_{L^{\infty }(D)}$$

$$\sum _{m=1}^{\infty }||a_m|{|}_{L^{\infty }(D)}^{p_0}<\infty$$

$$\sum _{m=1}^{\infty }||a_m|{|}_{W^{s,\infty }(D)}^{p_s}<\infty$$

$$||a_m|{|}_{L^{\infty }(D)}\le C{m}^{-\varsigma },\varsigma >1,m=1,2,\dots$$

$$(A(y)v)(x):=-\nabla \cdot (a(x,y)\nabla v(x)),x\in D,y\in \Gamma ,$$

$$a(x,y)=a_0(x)+\sum _{m=1}^{\infty }{a}_m(x)y_m.$$

$$A(y)u(y)=\mu (y)u(y)\forall y\in \Gamma .$$

$$b(y;u(\cdot ,y),v)=\mu (y)⟨u(\cdot ,y),v{⟩}_{L^2(D)}\forall v\in H_0^1(D),$$

$$b(y;u,v):={\int }_{D}a(\cdot ,y)\nabla u\cdot \nabla vdx.$$

$$E(\tau ):=\{z\in {\mathbb{C}}^{\infty }|\mathrm{dist}(z_m,[-1,1])\le {\tau }_m\}.$$

$${\tau }_m=\min \left\{C_1||a_m|{|}_{L^{\infty }(D)}^{p_0-1},C_2||a_m|{|}_{W^{s,\infty }(D)}^{p_s-1}\right\},m=1,2,\dots$$

$$\underset{v_h\in V_h}{\inf }||v-v_h|{|}_{L^2(D)}\le C{h}^{l+1}||v|{|}_{H^{l+1}(D)}$$

$$\underset{v_h\in V_h}{\inf }||v-v_h|{|}_{H_0^1(D)}\le C{h}^l||v|{|}_{H^{l+1}(D)}$$

$$b(y;u_h(\cdot ,y),v_h)={\mu }_h(y)⟨u_h(\cdot ,y),v_h{⟩}_{L^2(D)}\forall v_h\in V_h,$$

$$|\mu (y)-{\mu }_h(y)|\le C{h}^{2l},$$

$$||u(\cdot ,y)-u_h(\cdot ,y)|{|}_{L^2(D)}\le C{h}^{1+l}||u(\cdot ,y)|{|}_{H^{1+l}(D)}.$$

$$\left({\mathbf{K}}^{(0)}+\sum _{m=1}^{\infty }{\mathbf{K}}^{(m)}{y}_m\right){\mathbf{u}}_h(y)={\mu }_h(y)\mathbf{M}{\mathbf{u}}_h(y)\forall y\in \Gamma ,$$

$${\mathbf{K}}_{ij}^{(m)}={\int }_D{a}_m\nabla {\phi }_i\cdot \nabla {\phi }_{j}dx,m=0,1,\dots$$

$${\mathbf{M}}_{ij}={\int }_D{\phi }_i{\phi }_{j}dx.$$

$$\mathbb{E}[v]={\int }_{\Gamma }v(y)d\nu (y)$$

$$({\mathbb{N}}_0^{\infty }{)}_c:=\{\alpha \in {\mathbb{N}}_0^{\infty }|\#\mathrm{supp}(\alpha )<\infty \},$$

$${\Lambda }_{\alpha }(y):=\prod _{m\in \mathrm{supp}\alpha }{L}_{{\alpha }_m}(y_m),$$

$$v(y)=\sum _{\alpha \in ({\mathbb{N}}_0^{\infty }{)}_c}{v}_{\alpha }{\Lambda }_{\alpha }(y)$$

$$c_{\alpha \beta \gamma }$$

Full text

Keywords: Eigenvalue problems, error estimation, finite elements, sparse tensor approximation

AMS subject classifications: 65C20, 65N12, 65N15, 65N25, 65N30

1 Introduction

During the recent years numerical solution of stochastic partial differential equations (sPDE) has attracted a lot of attention and become a well-established field. However, the field of stochastic eigenvalue problems (sEVP) and their numerical solution is still in its infancy. It is natural that, after the source problem, more effort is put on addressing the eigenvalue problem.

A few different algorithms have recently been suggested for computing approximate eigenpairs of sEVPs. As with sPDEs, the solution methods are typically divided into intrusive and non-intrusive ones. A benchmark for non-intrusive methods is the sparse collocation algorithm suggested and thoroughly analyzed by Andreev and Schwab in [1]. An attempt towards a Galerkin-based (intrusive) method was made by Verhoosel et al. in [20], though this method omits uniform normalization of the eigenmodes. Very recently Meidani and Ghanem proposed a spectral power iteration, in which the eigenmodes are normalized using a quadrature rule over the parameter space [16]. The algorithm has been further developed and studied by Sousedík and Elman in [19]. However, neither of the papers present a comprehensive error analysis for the method.

Inspired by the original method of Meidani and Ghanem we have suggested a purely Galerkin-based spectral inverse iteration, in which normalization of the eigenmodes is achieved via solution of a simple nonlinear system [10]. This method, and its generalization to a spectral subspace iteration, is the focus of the current paper. Although the algorithms in [16] and [19] differ from ours in the way normalization is performed, the basic principles are still the same and hence our results on convergence should apply to these methods as well.

In this work we consider computing eigenpairs of an elliptic operator with random coefficients. We assume a physical domain $D\subset\mathbb{R}^{d}$ and, in order to capture the random dimension of the system, a parameter domain $\Gamma\subset\mathbb{R}^{\infty}$ with associated measure $\nu$. One may think of a parametrization that arises from Karhunen-Loève representations of the random coefficients in the system, for instance. Discretization in space is achieved by standard FEM and associated with a discretization parameter $h$, whereas discretization in the random dimension is achieved using collections of certain multivariate polynomials. These collections are represented by multi-index sets $\mathcal{A}_{\epsilon}$ of increasing cardinality $\#\mathcal{A}_{\epsilon}$ as $\epsilon\to 0$.

In the current paper we present a step-by-step analysis that leads to the main result: the asymptotic convergence of the spectral inverse iteration towards the exact eigenpair $(\mu,u)$. In this context the eigenpair of interest is the ground state, i.e., the smallest eigenvalue and the associated eigenfunction of the system. However, analogously to the classical inverse iteration, the computation of other eigenpairs may be possible by using a suitably chosen shift parameter $\lambda\in\mathbb{R}$. We show that under sufficient assumptions the iterates of the algorithm $(\mu_{k},u_{k})$ for $k=1,2,\ldots$ obey

[TABLE]

and

[TABLE]

where $l\in\mathbb{N}$ is the degree of polynomials used in the spatial discretization and $r>0$ depends on the properties of the region to which the solution, as a function of the parameter vector, admits a complex-analytic extension. The quantity $\lambda_{1/2}$ reflects the gap between the two smallest eigenvalues of the system and should be less than one.

The first term in the formulas (1) and (2) is justified by standard theory for Galerkin approximation of eigenvalue problems, a simple consequence of which we have recapped in Theorem 1. The second term can be deduced from Theorem 2, which bounds the Galerkin approximation errors by residuals of certain polynomial approximations of the solution. Using best $P$-term polynomial approximations, we see that these residuals are ultimately expected to decay at an algebraic rate $r>0$, see [15] and [6]. Finally, the third term follows from Theorem 3, which states that asymptotically the iterates of the spectral inverse iteration converge to a fixed point in geometric fashion. Here the analogy to classical inverse iteration is evident. Each of these three important steps that comprise the main result is separately verified through detailed numerical examples.

A variant of our algorithm for spectral subspace iteration is also presented. No analysis of this algorithm is given, but the numerical experiments support the conclusion that it convergences towards the exact subspace of interest, and that the rate of convergence is analogous to what we would expect from classical theory. This is despite the fact that the individual eigenmodes, as defined by the pointwise order of magnitude of the eigenvalues, are not continuous functions over the parameter space due to an eigenvalue crossing. To the authors’ knowledge such a scenario has not yet been considered in the scientific literature.

The rest of the paper is organized as follows. Our model problem and its fundamental properties are assessed in Sections 2 and 3. A detailed review of the discretization of the spatial and stochastic approximation spaces is given in Section 4. Analysis of the spectral inverse iteration, supported by thorough numerical experiments, is given in Section 5. Finally, the algorithm of spectral subspace iteration and numerical experiments of its convergence are presented in Section 6.

2 Problem statement

In this work we consider eigenvalue problems of elliptic operators with random coefficients. It is assumed that the random coefficients admit a parametrization with respect to countably many independent and bounded random variables. As a model problem we consider the eigenvalue problem of a diffusion operator with a random diffusion coefficient. It will be evident, however, that our methods and analysis in fact cover a much broader class of problems.

2.1 Model problem

Let $(\Omega,\mathcal{F},\mathcal{P})$ be a probability space, $\Omega$ being the set of outcomes, $\mathcal{F}$ a $\sigma$-algebra of events, and $\mathcal{P}$ a probability measure defined on $\Omega$. We denote by $L^{2}_{\mathcal{P}}(\Omega)$ the space of square integrable random variables on $\Omega$ and define for $v\in L^{2}_{\mathcal{P}}(\Omega)$ the expected value

[TABLE]

and variance $\mathrm{Var}[v]=\mathbb{E}[(v-\mathbb{E}[v])^{2}]$.

Let $D\subset\mathbb{R}^{d}$ be a bounded convex domain with a sufficiently smooth boundary and assume a diffusion coefficient $a\!:D\times\Omega\to\mathbb{R}$ that is a random field on $D$. The diffusion coefficient is assumed to be strictly uniformly positive and uniformly bounded, i.e., for some positive constants $a_{\min}$ and $a_{\max}$ it holds that

[TABLE]

We now formulate the model problem as: find functions $\mu\!:\Omega\to\mathbb{R}$ and $u\!:D\times\Omega\to\mathbb{R}$ such that the equations

[TABLE]

hold $\mathcal{P}$-almost surely. In order to make the solutions physically meaningful we also impose a normalization condition $||u(\cdot,\omega)||_{L^{2}(D)}=1$ that should hold $\mathcal{P}$-almost surely.

2.2 Parametrization of the random input

We make the assumption that the input random field admits a representation of the form

[TABLE]

where $\{y_{m}\}_{m=1}^{\infty}$ are mutually independent and bounded random variables. For simplicity, we assume here that each $y_{m}$ is uniformly distributed. Thus, after possible rescaling, the dependence on $\omega$ is now parametrized by the vector $y=(y_{1},y_{2},\ldots)\in\Gamma:=[-1,1]^{\infty}$. We denote by $\nu$ the underlying uniform product probability measure and by $L^{2}_{\nu}(\Gamma)$ the corresponding weighted $L^{2}$-space.

The usual convention is that the parametrization (5) results from a Karhunen-Loève expansion, which gives $a(x,\omega)$ as a linear combination of the eigenfunctions of the associated covariance operator. The distinguishing feature of the Karhunen-Loève expansion compared to other linear expansions is that it minimizes the mean square truncation error [8].

It is easy to see that $a_{0}\in L^{\infty}(D)$ and

[TABLE]

are sufficient conditions to ensure the assumption (3). In order to ensure analyticity of the eigenpair $(\mu,u)$ with respect to the parameter vector $y=(y_{1},y_{2},\ldots)$ we assume that

[TABLE]

for some $p_{0}\in(0,1)$ and that for a certain level of smoothness $s\in\mathbb{N}$ we have $a_{0}\in W^{s,\infty}(D)$ and

[TABLE]

for some $p_{s}\in(0,1)$. In particular, we consider the interesting case of algebraic

[TABLE]

decay of the coefficients in the series (5).

2.3 Parametric eigenvalue problem and its variational formulation

With the diffusion coefficient given by (5), the model problem (4) becomes an eigenvalue problem of the operator

[TABLE]

where

[TABLE]

Thus, we obtain the parametric eigenvalue problem: find $\mu\!:\Gamma\to\mathbb{R}$ and $u\!:\Gamma\to H_{0}^{1}(D)$ such that

[TABLE]

We denote by $\sigma(A(y))$ the set of eigenvalues of $A(y)$ for $y\in\Gamma$.

For any fixed $y\in\Gamma$ the problem (9) reduces to a single deterministic eigenvalue problem. In variational form this is given by: find $\mu(y)\in\mathbb{R}$ and $u(\cdot,y)\in H_{0}^{1}(D)$ such that

[TABLE]

where

[TABLE]

Under assumption (6) the bilinear form $b(y;u,v)$ is continuous and elliptic. Thus, as in [10] and [1], we deduce that the problem (10) admits a countable number of real eigenvalues and corresponding eigenfunctions that form an orthogonal basis of $L^{2}(D)$.

3 Analyticity of eigenmodes

A key issue in the analysis of parametric eigenvalue problems is that eigenvalues may cross within the parameter space. Here we first disregard this possibility and recap the main results from [1] for simple eigenvalues that are sufficiently well separated from the rest of the spectrum. After this we briefly comment on the case of possibly clustered eigenvalues and associated invariant subspaces.

3.1 Eigenmodes of simple eigenvalues

Let us first restrict our analysis to eigenvalues that are strictly nondegenerate. We call an eigenvalue $\mu$ of problem (9) strictly nondegenerate if

• (i)

$\mu(y)$ is simple as an eigenvalue of $A(y)$ for all $y\in\Gamma$ and

• (ii)

the minimum spectral gap $\inf_{y\in\Gamma}\operatorname*{dist}(\mu(y),\sigma(A(y))\backslash\{\mu(y)\})$ is positive.

In the case of strictly nondegenerate eigenvalues, the eigenpair $(\mu,u)$ is in fact analytic with respect to the parameter vector $y$.

Proposition 1**.**

Consider a strictly nondegenerate eigenvalue $\mu$ of the problem (9) and the corresponding eigenfunction $u$ normalized so that $||u(y)||_{L^{2}(D)}=1$ for all $y\in\Gamma$. For $s\in\mathbb{N}$ assume that $a_{0}\in W^{s,\infty}(D)$ and the assumptions (6) – (8) hold for some $p_{0},p_{s}\in(0,1)$. Given $\tau=(\tau_{1},\tau_{2},\ldots)\in\mathbb{R}_{+}^{\infty}$ define

[TABLE]

Then there exists $C_{1}>0$ independent of $m$ such that with $C_{2}>0$ arbitrary and $\tau$ given by

[TABLE]

the eigenpair $(\mu,u)$ can be extended to a jointly complex-analytic function on $E(\tau)$ with values in $\mathbb{C}\times(H^{s+1}(D)\cap H_{0}^{1}(D))$.

Proof.

This is analogous to Corollary 2 of Theorem 4 in [1]. ∎

It is well known that for elliptic operators on a connected domain $D$ the smallest eigenvalue is simple [11]. Thus, Proposition 1 may at least be applied for the smallest eigenvalue of problem (9).

3.2 Finite dimensional subspaces

Let us now consider invariant subspaces for which the corresponding cluster of eigenvalues is sufficiently well separated from the rest of the spectrum. Assume a cluster $\mathcal{M}(y)=\{\mu_{q}(y)\}_{q=1}^{Q}$ of eigenvalues of (9) so that

• (i)

each $\mu_{q}(y)$ is of finite multiplicity as an eigenvalue of $A(y)$ for all $y\in\Gamma$ and

• (ii)

the minimum spectral gap $\inf_{y\in\Gamma}\operatorname*{dist}(\mathcal{M}(y),\sigma(A(y))\backslash\mathcal{M}(y))$ is positive.

It is in general difficult to consider the analyticity of each of the eigenmodes separately. However, we might still expect the associated invariant subspace to be analytic as a function of $y$. More precisely, let $\{u_{q}(y)\}_{q=1}^{Q^{\prime}}$ be a maximal collection of linearly independent eigenfunctions corresponding to the eigenvalues $\mathcal{M}(y)$ for all $y\in\Gamma$. It is not completely unreasonable to assume that $\operatorname*{span}\{u_{q}(y)\}_{q=1}^{Q^{\prime}}$ is analytic, in a suitable sense, as a function of the parameter vector $y$. This assumption is the basis of our Algorithm of spectral subspace iteration given in Section 6. For more information on the regularity of perturbed eigenvalues see [13] and [14].

4 Stochastic finite elements

Proposition 1, under sufficient assumptions, guarantees the existence of an analytic eigenpair for problem (9). It now makes sense to look for the eigenvalue in the space $L^{2}_{\nu}(\Gamma)$ and the eigenfunction in the space $L^{2}_{\nu}(\Gamma)\otimes H^{1}_{0}(D)$. The space $H^{1}_{0}(D)$ may be discretized by means of the traditional finite element method. For the discretization of $L^{2}_{\nu}(\Gamma)$, we follow the usual convention in stochastic Galerkin methods and construct a basis of orthogonal polynomials of the input random variables. Orthogonal polynomials for various probability distributions exist and the use of these as the approximation basis has been observed to yield optimal rates of convergence [18], [21]. Here we consider uniformly distributed random variables which lead to the choice of tensorized Legendre polynomials.

4.1 Galerkin discretization in space

Let $V_{h}\subset H_{0}^{1}(D)$ denote a finite dimensional approximation space associated with the discretization parameter $h>0$. We assume approximation estimates

[TABLE]

and

[TABLE]

that are standard for piecewise polynomials of degree $l$.

Fix $y\in\Gamma$ and let $(\mu_{h},u_{h})$ be the solution to the variational equation

[TABLE]

where $b(y;\cdot,\cdot)$ is as in (10). Then we have the following bounds for the discretization error.

Theorem 1**.**

Assume (11) and (12). For $y\in\Gamma$ let $\mu(y)$ be a simple eigenvalue of (10) and $\mu_{h}(y)$ an eigenvalue of (13) such that $\lim_{h\to 0}\mu_{h}(y)=\mu(y)$. Let $u(\cdot,y)\in H^{1+l}(D)$ and $u_{h}(\cdot,y)\in V_{h}$ denote the associated eigenfunctions normalized in $L^{2}(D)$. Then there exists $C>0$ such that

[TABLE]

and

[TABLE]

as $h\to 0$.

Proof.

This follows from the theory of Galerkin approximation for variational eigenvalue problems. See Section 8 in [3] and Section 9 in [7]. ∎

Let $V_{h}=\operatorname*{span}\{\varphi_{i}\}_{i\in J}$ where $J:=\{1,2,\ldots,N\}$. Then (13) can be written as a parametric matrix eigenvalue problem: find $\mu_{h}\!:\Gamma\to\mathbb{R}$ and $\mathbf{u}_{h}\!:\Gamma\to\mathbb{R}^{N}$ such that

[TABLE]

where $u_{h}(x,y)=\sum_{i\in J}\varphi_{i}(x)(\mathbf{u}_{h})_{i}(y)$. The coefficient matrices are given by

[TABLE]

and

[TABLE]

For each fixed $y\in\Gamma$ the problem (16) reduces to a positive-definite generalized matrix eigenvalue problem.

4.2 Legendre chaos

Recall that $y=(y_{1},y_{2},\ldots)\in\Gamma$ is a vector of mutually independent uniform random variables and $\nu$ is the underlying constant product probablity measure. Now

[TABLE]

whenever the integral is finite. We define $(\mathbb{N}_{0}^{\infty})_{c}$ to be the set of all multi-indices with finite support, i.e.,

[TABLE]

where $\mathop{\mathrm{supp}}(\alpha)=\{m\in\mathbb{N}\ |\ \alpha_{m}\not=0\}$. Given a multi-index $\alpha\in(\mathbb{N}_{0}^{\infty})_{c}$ we now define the multivariate Legendre polynomial

[TABLE]

where $L_{p}(x)$ denotes the univariate Legendre polynomial of degree $p$. We will assume the normalization $\mathbb{E}[\Lambda_{\alpha}^{2}]=1$ for all $\alpha\in(\mathbb{N}_{0}^{\infty})_{c}$.

The system $\{\Lambda_{\alpha}(y)\ |\ \alpha\in(\mathbb{N}_{0}^{\infty})_{c}\}$ forms an orthonormal basis of $L^{2}_{\nu}(\Gamma)$. Therefore, we may write any square integrable random variable $v$ in a series

[TABLE]

with convergence in $L^{2}_{\nu}(\Gamma)$. The expansion coefficients are given by $v_{\alpha}=\mathbb{E}[v\Lambda_{\alpha}]$.

Due to the orthogonality of the Legendre polynomials we have $\mathbb{E}[\Lambda_{\alpha}]=\delta_{\alpha 0}$ and $\mathbb{E}[\Lambda_{\alpha}\Lambda_{\beta}]=\delta_{\alpha\beta}$ for all $\alpha,\beta\in(\mathbb{N}_{0}^{\infty})_{c}$. Moreover, we denote

[TABLE]

We now obtain moment matrices $G^{(m)}\in\mathbb{R}^{P\times P}$ for $m\in\mathbb{N}_{0}$ and $G^{(\alpha)}\in\mathbb{R}^{P\times P}$ for $\alpha\in\mathcal{A}$ by setting $[G^{(m)}]_{\alpha\beta}=c_{m\alpha\beta}$ and $[G^{(\alpha)}]_{\beta\gamma}=c_{\alpha\beta\gamma}$.

4.3 Sparse polynomial approximation in the parameter domain

We fix a finite set $\mathcal{A}\subset(\mathbb{N}_{0}^{\infty})_{c}$ and employ the approximation space $W_{\mathcal{A}}=\operatorname*{span}\{\Lambda_{\alpha}\}_{\alpha\in\mathcal{A}}\subset L^{2}_{\nu}(\Gamma)$. We let $P_{\mathcal{A}}$ and $R_{\mathcal{A}}$ denote the underlying projection and residual operators so that $v\in L^{2}_{\nu}(\Gamma)$ is approximated by

[TABLE]

and the approximation error is given by $R_{\mathcal{A}}(v)=v-P_{\mathcal{A}}(v)$. Since

[TABLE]

where $\mathcal{A}^{c}=\{\alpha\in(\mathbb{N}_{0}^{\infty})_{c}\ |\ \alpha\notin\mathcal{A}\}$, we conclude that the choice of the multi-index set $\mathcal{A}$ ultimately determines the accuracy of our expansion.

Remark 1**.**

The orthogonality of the basis $\{\Lambda_{\alpha}\}_{\alpha\in\mathcal{A}}$ yields easy formulas for the mean and variance:

[TABLE]

We proceed as in [15] and use best $P$-term approximations to prove convergence of the approximation error.

Proposition 2**.**

Let $H$ be a Hilbert space. Assume that $v\!:\Gamma\to H$ admits a complex-analytic extension in the region

[TABLE]

with

[TABLE]

Then for each $\epsilon>0$ and $0<r<\varrho-\frac{1}{2}$ there exists $\mathcal{A}_{\epsilon}\subset(\mathbb{N}_{0}^{\infty})_{c}$ such that

[TABLE]

and as $\epsilon\to 0$ we have

[TABLE]

Proof.

This is essentially Proposition 3.1 in [15], the proof of which can be found in [5]. Here we recapitulate the main ideas.

Let $M\in\mathbb{N}$, $\mathcal{A}\subset\mathbb{N}_{0}^{M}$ and $P=\#\mathcal{A}$. Set $v_{M}(z):=v(z_{1},\ldots,z_{M},0,0,\ldots)$ so that $P_{\mathcal{A}}(v)=P_{\mathcal{A}}(v_{M})$ and

[TABLE]

Define $\eta_{m}=\tau_{m}+\sqrt{1+\tau_{m}^{2}}$. It holds that

[TABLE]

and thus, by Lemmas 4.3. and 4.4 in [2], we obtain

[TABLE]

for any $M\geq C(\varrho)P^{r/(\varrho-1)}$. On the other hand, the proof of Lemma A.3 in [5] now applies to $v_{M}$, and we see that the Legendre coefficients (note the different scaling as compared to [5])

[TABLE]

satisfy

[TABLE]

We set $\vartheta:=(r+1/2)^{-1}$ so that $\varrho^{-1}<\vartheta<2$. Using the geometric series formula we compute

[TABLE]

where we have used the inequality $1+x\leq e^{x}$ for $x\geq 0$. Now assume that $v^{*}=\{v^{*}_{i}\}_{i=1}^{\infty}$ are the Legendre coefficients $\{(v_{M})_{\alpha}\}_{\alpha\in\mathbb{N}_{0}^{M}}\subset H$ of $v_{M}$ ordered by decreasing $H$-norm and that $\mathcal{A}$ represents a selection of the $P$ largest coefficients. Using (4.3) we may now write

[TABLE]

The claim follows from combining (26) and (4.3). ∎

4.4 Stochastic Galerkin approximation of vectors and matrices

We now generalize the concept of sparse polynomial approximation to vector and matrix valued functions. Assume that the dimensions of the approximation spaces $V_{h}$ and $W_{\mathcal{A}}$ are $N$ and $P$ respectively. We denote by $W_{\mathcal{A}}^{N}$ (or $W_{\mathcal{A}}^{N\times N}$) the space of functions $\mathbf{v}\!:\Gamma\to\mathbb{R}^{N}$ (or $\mathbf{A}\!:\Gamma\to\mathbb{R}^{N\times N}$) whose every component is in $W_{\mathcal{A}}$. Whenever $\mathbf{v}\in W_{\mathcal{A}}^{N}$ and $\alpha\in\mathcal{A}$ we set $v_{\alpha i}=(\mathbf{v}_{i})_{\alpha}$ and use $\mathbf{v}_{\alpha}$ to denote the vector of coefficients $\{v_{\alpha i}\}_{i\in J}\in\mathbb{R}^{N}$. Moreover, we associate any $v\in W_{\mathcal{A}}$ with the array of coefficients $\hat{v}:=\{v_{\alpha}\}_{\alpha\in\mathcal{A}}\in\mathbb{R}^{P}$ and similarly any $\mathbf{v}\in W_{\mathcal{A}}^{N}$ with the array of coefficients $\hat{\mathbf{v}}:=\{v_{\alpha i}\}_{\alpha\in\mathcal{A},i\in J}\in\mathbb{R}^{PN}$.

We denote by $\langle\cdot,\cdot\rangle_{\mathbb{R}^{N}_{\mathbf{M}}}$ the inner product on $\mathbb{R}^{N}$ induced by the positive definite matrix $\mathbf{M}$ and by $||\cdot||_{\mathbb{R}^{N}_{\mathbf{M}}}$ the associated norm. Furthermore, we let $||\cdot||_{\mathbb{R}^{P}}$ denote the standard norm on $\mathbb{R}^{P}$ and $||\cdot||_{\mathbb{R}^{P}\otimes\mathbb{R}^{N}_{\mathbf{M}}}$ denote the tensorized norm on $\mathbb{R}^{PN}$ given by

[TABLE]

Remark 2**.**

Observe that if $v\in W_{\mathcal{A}}\otimes V_{h}$ is written as $v(x,y)=\sum_{i\in J}\varphi_{i}(x)\mathbf{v}_{i}(y)$, then

[TABLE]

Let us consider the linear system defined by a parametric matrix $\mathbf{A}\in W_{\mathcal{A}}^{N\times N}$. The Galerkin approximation of this system is: given $\mathbf{f}\in W_{\mathcal{A}}^{N}$ find $\mathbf{v}\in W_{\mathcal{A}}^{N}$ such that

[TABLE]

Using the notation introduced earlier we may write (32) as the fully discrete system: given $\hat{\mathbf{f}}\in\mathbb{R}^{PN}$ find $\hat{\mathbf{v}}\in\mathbb{R}^{PN}$ such that

[TABLE]

where $\mathbf{A}_{\alpha}=\mathbb{E}[\mathbf{A}\Lambda_{\alpha}]\in\mathbb{R}^{N\times N}$. The existence of a solution, i.e. the invertibility of the coefficient matrix, is guaranteed by the following Lemma.

Lemma 1**.**

If $\mathbf{A}\in W_{\mathcal{A}}^{N\times N}$ is a parametric matrix such that $\mathbf{A}(y)$ is positive-definite for every $y\in\Gamma$, then for any $\mathbf{f}\in W_{\mathcal{A}}^{N}$ there exists a unique $\mathbf{v}\in W_{\mathcal{A}}^{N}$ such that

[TABLE]

Furthermore,

[TABLE]

where $\lambda(y)$ is the smallest eigenvalue of $\mathbf{A}(y)$ for each $y\in\Gamma$.

Proof.

Observe that the system (34) is equivalent to the variational form

[TABLE]

The left hand side of (36) is a symmetric and elliptic bilinear form so the existence of a unique solution is guaranteed by the Lax-Milgram Lemma. Moreover, the associated coefficient matrix in (33) is positive definite.

Now let $\tilde{\lambda}\in\mathbb{R}$ be such that $\tilde{\lambda}<\inf_{y\in\Gamma}\lambda(y)$. The matrix $\mathbf{A}(y)-\tilde{\lambda}\mathbf{I}_{N}$, where $\mathbf{I}_{N}$ is an identity matrix, is positive definite for all $y\in\Gamma$. Thereby the eigenvalues of the associated coefficient matrix should be positive. Let $\chi$ be an eigenvalue of (33), i.e., there exists $\mathbf{w}\in W_{\mathcal{A}}^{N}$ such that

[TABLE]

Then

[TABLE]

and we deduce that $\chi>\tilde{\lambda}$. Equation (35) now follows from taking the limit $\tilde{\lambda}\to\inf_{y\in\Gamma}\lambda(y)$. ∎

5 Spectral inverse iteration

In this section we introduce the Algorithm of spectral inverse iteration, analyze its asymptotic convergence, and present numerical examples to support our analysis. The spectral inverse iteration, see [10], can be considered as an extension of the classical inverse iteration to the case of parametric matrix eigenvalue problems. In the spectral version each of the elementary operations is computed in Galerkin sense via projecting to the sparse polynomial basis $W_{\mathcal{A}}$. Optimal convergence of the Algorithm requires that the eigenmode of interest, i.e., the smallest eigenvalue of the parametric matrix, is strictly nondegenerate.

5.1 Algorithm description

Fix a finite set of multi-indices $\mathcal{A}\subset(\mathbb{N}_{0}^{\infty})_{c}$ and let $P=\#\mathcal{A}$. The spectral inverse iteration for the system (16) is now defined in Algorithm 1. One should note that, if the projections in the Algorithm were precise, the Algorithm would correspond to performing classical inverse iteration pointwise over the parameter space $\Gamma$. We expect the Algorithm to converge to an approximation of the eigenvector corresponding to the smallest eigenvalue of the system.

Algorithm 1** (Spectral inverse iteration).**

Fix $tol>0$ and let $\mathbf{u}^{(0)}\in W_{\mathcal{A}}^{N}$ be an initial guess for the eigenvector. For $k=1,2,\ldots$ do

1. (1)

Solve $\mathbf{v}\in W_{\mathcal{A}}^{N}$ from the linear equation

[TABLE] 2. (2)

Solve $s\in W_{\mathcal{A}}$ from the nonlinear equation

[TABLE] 3. (3)

Solve $\mathbf{u}^{(k)}\in W_{\mathcal{A}}^{N}$ from the linear equation

[TABLE] 4. (4)

Stop if $||\mathbf{u}^{(k)}-\mathbf{u}^{(k-1)}||_{L^{2}_{\nu}(\Gamma)\otimes\mathbb{R}^{N}_{\mathbf{M}}}<tol$ and return $\mathbf{u}^{(k)}$ as the approximate eigenvector.

Once the approximate eigenvector $\mathbf{u}^{(k)}\in W_{\mathcal{A}}^{N}$ has been computed, the corresponding eigenvalue $\mu^{(k)}\in W_{\mathcal{A}}$ may be evaluated from the Rayleigh quotient, as in [10], or alternatively from the linear system

[TABLE]

Lemma 1 guarantees the invertibility of the linear system (39) and, assuming that $s(y)>0$ for all $y\in\Gamma$, the invertibility of the systems (41) and (42). The nonlinear system (40) may be solved using for instance Newton’s method.

Remark 3**.**

For the computation of non-extremal eigenmodes, one may proceed as in [10] and replace $\mathbf{K}(y)$ in (39) with $(\mathbf{K}(y)-\lambda\mathbf{M})$, where $\lambda\in\mathbb{R}$ is a suitably chosen parameter. In this case we expect the algorithm to converge to an eigenpair for which the eigenvalue is close to $\lambda$. Note, however, that now the existence of a unique solution to (39) is not necessarily guaranteed by Lemma 1.

We try to write Algorithm 1 in a computationally more convenient form. The projections in the algorithm can be computed explicitly using the notation introduced in Section 4. It is easy to verify that equations (39) – (41) become

[TABLE]

respectively. Given $\hat{s}=\{s_{\alpha}\}_{\alpha\in\mathcal{A}}\in\mathbb{R}^{P}$ we define matrices

[TABLE]

where $M(\mathcal{A}):=\max\{m\in\mathbb{N}\ |\ \exists\alpha\in\mathcal{A}\textrm{ s.t. }\alpha\not=0\}$ and $I_{P}\in\mathbb{R}^{P\times P}$ and $\mathbf{I}_{N}\in\mathbb{R}^{N\times N}$ are identity matrices. We also define the nonlinear function $F\!:\mathbb{R}^{P}\times\mathbb{R}^{PN}\to\mathbb{R}^{P}$ via

[TABLE]

and let $F^{s}\!:\mathbb{R}^{P}\times\mathbb{R}^{P}\to\mathbb{R}^{P}$ and $F^{v}\!:\mathbb{R}^{PN}\times\mathbb{R}^{PN}\to\mathbb{R}^{P}$ denote the associated bilinear forms given by $F^{s}_{\alpha}(\hat{s},\hat{t}):=\hat{s}\cdot G^{(\alpha)}\hat{t}$ and $F^{v}_{\alpha}(\hat{\mathbf{v}},\hat{\mathbf{w}}):=\hat{\mathbf{v}}\cdot(G^{(\alpha)}\otimes\mathbf{M})\hat{\mathbf{w}}$. Now Algorithm 1 may be rewritten in the following form.

Algorithm 2** (Spectral inverse iteration in tensor form).**

Fix $tol>0$ and let $\hat{\mathbf{u}}^{(0)}=\{u_{\alpha i}^{(0)}\}_{\alpha\in\mathcal{A},i\in J}\in\mathbb{R}^{PN}$ be an initial guess for the eigenvector. For $k=1,2,\ldots$ do

1. (1)

Solve $\hat{\mathbf{v}}=\{v_{\alpha i}\}_{\alpha\in\mathcal{A},i\in J}\in\mathbb{R}^{PN}$ from the linear system

[TABLE] 2. (2)

Solve $\hat{s}=\{s_{\alpha}\}_{\alpha\in\mathcal{A}}\in\mathbb{R}^{P}$ from the nonlinear system

[TABLE]

with the initial guess $s_{\alpha}=||\hat{\mathbf{v}}||_{\mathbb{R}^{P}\otimes\mathbb{R}^{N}_{\mathbf{M}}}\delta_{\alpha 0}$ for $\alpha\in\mathcal{A}$. 3. (3)

Solve $\hat{\mathbf{u}}^{(k)}=\{u_{\alpha i}^{(k)}\}_{\alpha\in\mathcal{A},i\in J}\in\mathbb{R}^{PN}$ from the linear system

[TABLE] 4. (4)

Stop if $||\hat{\mathbf{u}}^{(k)}-\hat{\mathbf{u}}^{(k-1)}||_{\mathbb{R}^{P}\otimes\mathbb{R}^{N}_{\mathbf{M}}}<tol$ and return $\hat{\mathbf{u}}^{(k)}$ as the approximate eigenvector.

The approximate eigenvalue $\hat{\mu}^{(k)}\in\mathbb{R}^{P}$ may now be solved from the equation

[TABLE]

where $\hat{e}_{1}=\{\delta_{\alpha 0}\}_{\alpha\in\mathcal{A}}\in\mathbb{R}^{P}$.

Remark 4**.**

In [10] Newton’s method with the initial guess $s_{\alpha}=||\mathbf{v}_{\alpha}||_{\mathbb{R}^{N}_{\mathbf{M}}}$ was suggested for the system of equations (47). Here the initial guess is somewhat different and corresponds to $s_{0}=||\mathbf{v}||_{L^{2}_{\nu}(\Gamma)\otimes\mathbb{R}^{N}_{\mathbf{M}}}$ (and $s_{\alpha}=0$ for $\alpha\not=0$).

In general it is not guaranteed that the Newton iteration for the system (47) converges to a solution. The following proposition will give some insight to the conditions under which this happens to be the case.

Proposition 3**.**

Fix $\hat{\mathbf{v}}\in\mathbb{R}^{PN}$ and let $\hat{s}^{(0)}=\{s^{(0)}_{\alpha}\}_{\alpha\in\mathcal{A}}\in\mathbb{R}^{P}$ be given by $s_{\alpha}^{(0)}=||\hat{\mathbf{v}}||_{\mathbb{R}^{P}\otimes\mathbb{R}^{N}_{\mathbf{M}}}\delta_{\alpha 0}$ for $\alpha\in\mathcal{A}$. Assume that there is a norm $||\cdot||_{*}$ on $\mathbb{R}^{P}$ and $r>0$ such that

[TABLE]

for all $\hat{s},\hat{t}$ in $B(\hat{s}^{(0)},r):=\{\hat{s}\in\mathbb{R}^{P}\ |\ ||\hat{s}-\hat{s}^{(0)}||_{*}\leq r\}$. If

[TABLE]

then the Newton method for $F(\cdot,\hat{\mathbf{v}})=0$ with the initial guess $\hat{s}^{(0)}$ converges to a unique solution in $B(\hat{s}^{(0)},r)$.

Proof.

This is a direct application of the Newton-Kantorovich theorem for the equation $F(\cdot,\hat{\mathbf{v}})=0$, see [12] (Theorem 6, 1.XVIII). Note that the first derivative (Jacobian) of $F(\cdot,\hat{\mathbf{v}})$ at $\hat{s}^{(0)}$ is $2||\hat{\mathbf{v}}||_{\mathbb{R}^{P}\otimes\mathbb{R}^{N}_{\mathbf{M}}}I_{P}$ and the second derivative is represented by the tensor of coefficients $2c_{\alpha\beta\gamma}$. ∎

From Proposition 3 we see that convergence of the Newton iteration is a consequence of the boundedness of the function $F^{s}$, which again is ultimately determined by the structure of the multi-index set $\mathcal{A}$.

5.2 Analysis of convergence

Due to a lack of general mathematical theory for multi-parametric eigenvalue problems we rely on a slightly unconventional approach in analyzing our algorithm. First of all, we restrict ourselves to asymptotic analysis since the underlying problem is nonlinear and thus hard to analyze globally. Second, we will analyze the solutions pointwise in the parameter space and deduce convergence theorems from classical eigenvalue perturbation bounds.

5.2.1 Characterization of the dominant fixed point

The classical inverse iteration converges to the dominant eigenpair of the inverse matrix. In a somewhat similar fashion the spectral inverse iteration tends to converge to a certain fixed point, which we shall refer to as the dominant fixed point. Here we will establish a connection between this dominant fixed point of the spectral inverse iteration and the dominant eigenpair of the inverse of the parametric matrix under consideration. This connection is obtained by considering the fixed point as a pointwise perturbation of the eigenvalue problem of the parametric matrix.

If $\mathbf{u}_{\mathcal{A}}\in W_{\mathcal{A}}^{N}$ is a fixed point of the Algorithm 1, then there exists a pair $(s,\mathbf{v})\in W_{\mathcal{A}}\times W_{\mathcal{A}}^{N}$ such that $\mathbf{u}_{\mathcal{A}}=\mathbf{M}^{-1}P_{\mathcal{A}}(\mathbf{K}\mathbf{v})$ and

[TABLE]

We call $\mathbf{u}_{\mathcal{A}}$ the dominant fixed point if, whenever $(\tilde{s},\tilde{\mathbf{v}})\not=(s,\mathbf{v})$ also solves the system (50), then $s(y)>\tilde{s}(y)$ for all $y\in\Gamma$. For any fixed $y\in\Gamma$ we may write (50) as

[TABLE]

The following Lemma will be helpful in establishing a connection between the eigenpair of interest and the system (51).

Lemma 2**.**

Denote by $||\cdot||$ the standard Euclidean norm on $\mathbb{R}^{N}$. Assume that $S\in\mathbb{R}^{N\times N}$ can be diagonalized as

[TABLE]

where $\lambda_{1}\in\mathbb{R}$, $\Lambda=\operatorname*{diag}(\lambda_{2},\ldots,\lambda_{N})$ is real, and $(x\ X)$ is orthogonal. Assume also that $\lambda_{1}>\lambda_{2}\geq\ldots\geq\lambda_{N}$ and denote $\hat{\lambda}:=\lambda_{1}-\lambda_{2}$. Let $\rho\in\mathbb{R}$ and $r\in\mathbb{R}^{N}$ be such that $|\rho|\leq 1/2$ and $||r||\leq\hat{\lambda}/8$. Then there exist $\kappa\geq 1/2$ and $\pi\in\mathbb{R}^{N-1}$ such that

• (i)

The pair $(s,w)$ given by $s=\lambda_{1}-\kappa^{-1}x^{T}r$ and $w=\kappa x+X\pi$ solves the system

[TABLE]

• (ii)

If $(\tilde{s},\tilde{w})\not=(s,w)$ also solves the system (53), then $s>\tilde{s}$ or $x^{T}\tilde{w}<0$.

• (iii)

There exists $C>0$ such that $|\kappa-1|\leq C(|\rho|+\hat{\lambda}^{-2}||r||^{2})$ and $||\pi||\leq C\hat{\lambda}^{-1}||r||$.

Proof.

• (i)

Let $s(\kappa)=\lambda_{1}-\kappa^{-1}x^{T}r$. For any $\kappa\geq 1/2$ we have $|\kappa^{-1}x^{T}r|\leq\hat{\lambda}/4$ so that

[TABLE]

and

[TABLE]

The function

[TABLE]

is strictly increasing for $\kappa\geq 1/2$ since

[TABLE]

One may also verify that $f(1/2)<0$ and $f(2)>0$. Thus, we may choose $\kappa>1/2$ such that $f(\kappa)=0$. For $w=\kappa x+X\pi$ we obtain

[TABLE]

so the equation $Sw=sw+r$ is equivalent to

[TABLE]

Choosing $s=s(\kappa)$ and $\pi=(\Lambda-sI)^{-1}X^{T}r$ we see that both equations are satisfied. Moreover

[TABLE]

• (ii)

Suppose $(\tilde{s},\tilde{w})$ also solves the system (53) and write $\tilde{w}=\tilde{\kappa}x+X\tilde{\pi}$ for some $\tilde{\kappa}\in\mathbb{R}$ and $\tilde{\pi}\in\mathbb{R}^{N-1}$. In the nontrivial case we have $\tilde{\kappa}=x^{T}\tilde{w}>0$. Assume first that $0\leq\tilde{\kappa}\leq 1/2$. We have

[TABLE]

Since $s\geq\lambda_{1}-\kappa^{-1}||r||$, we deduce that

[TABLE]

Now let $\tilde{\kappa}\geq 1/2$. If $(\tilde{s},\tilde{w})$ is to solve (53) then, as in part (i), we should have

[TABLE]

From the first equation we obtain $\tilde{s}=\lambda_{1}-\tilde{\kappa}^{-1}x^{T}r$. Due to $|\tilde{\kappa}^{-1}x^{T}r|\leq\hat{\lambda}/4$ the matrix $(\Lambda-\tilde{s}I)$ is invertible so the second equation gives $\tilde{\pi}=(\Lambda-\tilde{s}I)^{-1}X^{T}r$. Here $\tilde{\kappa}\geq 1/2$ must be chosen so that $f(\tilde{\kappa})=0$ and therefore $(\tilde{s},\tilde{w})=(s,w)$.

• (iii)

From $f(\kappa)=0$ and $\kappa\geq 1/2$ we deduce that

[TABLE]

and

[TABLE]

∎

Applying Lemma 2 to the system (51) pointwise for $y\in\Gamma$ we obtain the following result.

Proposition 4**.**

Let $\mathbf{u}_{\mathcal{A}}\in W_{\mathcal{A}}^{N}$ be the dominant fixed point of Algorithm 1 and denote by $(s,\mathbf{v})$ the associated pair in $W_{\mathcal{A}}\times W_{\mathcal{A}}^{N}$ that solves (50). Let $\mu_{\mathcal{A}}\in W_{\mathcal{A}}$ be such that $P_{\mathcal{A}}(s\mu_{\mathcal{A}})=1$. For $y\in\Gamma$ denote by $\hat{\lambda}(y)$ the gap between the two largest eigenvalues of $\mathbf{K}^{-1}(y)\mathbf{M}$. Assume that $\inf_{y\in\Gamma}s(y)>0$ and $\inf_{y\in\Gamma}\hat{\lambda}(y)>0$. For $y\in\Gamma$ define

[TABLE]

and

[TABLE]

If

[TABLE]

and

[TABLE]

then there exists $C>0$ such that

[TABLE]

and

[TABLE]

where $\mu_{h}\!:\Gamma\to\mathbb{R}$ is the smallest eigenvalue of $\mathbf{M}^{-1}\mathbf{K}(y)$ and $\mathbf{u}_{h}\!:\Gamma\to\mathbb{R}^{N}$ is the corresponding eigenvector normalized in $||\cdot||_{\mathbb{R}^{N}_{\mathbf{M}}}$ (and with appropriate sign).

Proof.

It is easy to see that the system (51) is equivalent to

[TABLE]

where $\mathbf{w}(y)=s^{-1}(y)\mathbf{M}^{\frac{1}{2}}\mathbf{v}(y)$. By Lemma 2 the solution with the pointwise largest $s(y)$ can be written as

[TABLE]

where $\kappa\!:\Gamma\to[1/2,\infty)$ and $\boldsymbol{\pi}\!:\Gamma\to\mathbb{R}^{N-1}$ are such that

[TABLE]

$\lambda_{1}(y)=\mu_{h}^{-1}(y)$ is the pointwise largest eigenvalue of $\mathbf{S}(y)$ and $\mathbf{x}(y)=\mathbf{M}^{\frac{1}{2}}\mathbf{u}_{h}(y)$ is the corresponding eigenvector. The matrix $(\mathbf{x}(y)\ \mathbf{X}(y))$ is orthonormal for every $y\in\Gamma$. A Taylor expansion of $s^{-1}(y)$ yields

[TABLE]

where $\xi(y)$ is such that $0\leq\xi(y)\leq\kappa^{-1}(y)\mathbf{x}^{T}(y)\mathbf{M}^{\frac{1}{2}}\mathbf{r}(y)$. Combining this with the equation

[TABLE]

obtained from the condition $P_{\mathcal{A}}(s\mu_{\mathcal{A}})=1$, we have altogether that

[TABLE]

Furthermore,

[TABLE]

from which it follows that

[TABLE]

∎

Remark 5**.**

Note that we have not proven the existence of a dominant fixed point of the Algorithm 1. The residuals $\mathbf{r}$ and $\rho$ in Proposition 4 depend on the pair $(s,\mathbf{v})\in W_{\mathcal{A}}\times W_{\mathcal{A}}^{N}$ and hence Lemma 2 by itself is not sufficient to guarantee the existence of a dominant fixed point.

5.2.2 Convergence of the dominant fixed point to a parametric eigenpair

The next step in our analysis is to bound the error between the dominant fixed point of Algorithm 1 and the dominant eigenpair of the inverse of the parametric matrix. To this end we will use the pointwise estimate obtained previously.

From Proposition 4 we may easily deduce the following result.

Theorem 2**.**

Let $\mathbf{u}_{\mathcal{A}}\in W_{\mathcal{A}}^{N}$ be the dominant fixed point of Algorithm 1 and denote by $(s,\mathbf{v})$ the associated pair in $W_{\mathcal{A}}\times W_{\mathcal{A}}^{N}$ that solves (50). Let $\mu_{\mathcal{A}}\in W_{\mathcal{A}}$ be such that $P_{\mathcal{A}}(s\mu_{\mathcal{A}})=1$. For $y\in\Gamma$ denote by $\hat{\lambda}(y)$ the gap between the two largest eigenvalues of $\mathbf{K}^{-1}(y)\mathbf{M}$. Assume that $s_{*}:=\inf_{y\in\Gamma}s(y)>0$, $\hat{\lambda}_{*}:=\inf_{y\in\Gamma}\hat{\lambda}(y)>0$ and that the quantity

[TABLE]

is small enough. Then there exists $C>0$ such that

[TABLE]

and

[TABLE]

where $\mu_{h}\!:\Gamma\to\mathbb{R}$ is the smallest eigenvalue of $\mathbf{M}^{-1}\mathbf{K}(y)$ and $\mathbf{u}_{h}\!:\Gamma\to\mathbb{R}^{N}$ is the corresponding eigenvector normalized in $||\cdot||_{\mathbb{R}^{N}_{\mathbf{M}}}$ (and with appropriate sign). Here $C$ depends only on $s_{*}$, $\hat{\lambda}_{*}$, $\mu_{h}^{*}:=\sup_{y\in\Gamma}\mu_{h}(y)$, $K_{*}=\sup_{y\in\Gamma}||\mathbf{K}^{-1}(y)||_{\mathbb{R}^{N}_{\mathbf{M}}}$, and $M_{*}=||\mathbf{M}^{-1}||_{\mathbb{R}^{N}_{\mathbf{M}}}$.

Proof.

With $\mathbf{r}$ defined as in Proposition 4 we have

[TABLE]

and

[TABLE]

The bounds (2) and (2) now follow from Proposition 4. ∎

By Proposition 1 the exact eigenvalue and eigenvector of problem (9) are analytic functions of the parameter vector $y\in\Gamma$. This suggests that the residuals on the right hand side of equations (2) and (2) can be asymptotically estimated from Proposition 2.

5.2.3 Convergence of the spectral inverse iteration to the dominant fixed point

The classical inverse iteration converges to the dominant eigenpair of the inverse matrix at a speed characterized by the gap between the two largest eigenvalues. Here we will establish a similar asymptotic result for the convergence of the spectral inverse iteration towards the dominant fixed point.

Fixed points of the spectral inverse iteration may be characterized using the tensor notation of Algorithm 2. Let $\hat{\mathbf{u}}_{\mathcal{A}}\in\mathbb{R}^{PN}$ be a fixed point of the algorithm, i.e., $\hat{\mathbf{u}}_{\mathcal{A}}=\mathbf{S}\hat{\mathbf{v}}$ and $(\hat{s},\hat{\mathbf{v}})\in\mathbb{R}^{P}\times\mathbb{R}^{PN}$ are such that

[TABLE]

Define a linear operator $\mathbf{R}(\hat{s},\hat{\mathbf{v}})\!:\mathbb{R}^{PN}\to\mathbb{R}^{PN}$ by

[TABLE]

The convergence of the spectral inverse iteration to the fixed point $\hat{\mathbf{u}}_{\mathcal{A}}$ can now be related to the ratio of the norms of $\Delta^{-1}(\hat{s})$ and $\mathbf{R}(\hat{s},\hat{\mathbf{v}})\mathbf{S}^{-1}$.

Theorem 3**.**

Let $\hat{\mathbf{u}}_{\mathcal{A}}\in\mathbb{R}^{PN}$ be a fixed point of the Algorithm 2 and $(\hat{s},\hat{\mathbf{v}})\in\mathbb{R}^{P}\times\mathbb{R}^{PN}$ a corresponding solution to (79). Assume that $\Delta(\hat{s})$ is invertible. Let $\hat{\mu}_{\mathcal{A}}=\Delta^{-1}(\hat{s})\hat{e}_{1}$, where $\hat{e}_{1}=\{\delta_{\alpha 0}\}_{\alpha\in\mathcal{A}}\in\mathbb{R}^{P}$. Set $\phi_{\min}:=||\Delta^{-1}(\hat{s})||^{-1}_{\mathbb{R}^{P}}$ and $\psi_{\max}:=||\mathbf{R}(\hat{s},\hat{\mathbf{v}})\mathbf{S}^{-1}||_{\mathbb{R}^{P}\otimes\mathbb{R}^{N}_{\mathbf{M}}}$. Then the iterates of Algorithm 2 satisfy

[TABLE]

for $\hat{\mathbf{u}}^{(k)}$ sufficiently close to $\hat{\mathbf{u}}_{\mathcal{A}}$. Furthermore, there exists $C>0$ such that

[TABLE]

Proof.

The partial derivative (Jacobian) of the function $F(\hat{s}+\hat{t},\hat{\mathbf{v}}+\hat{\mathbf{w}})$ with respect to $\hat{t}$ at $\hat{t}=0$ is given by $2\Delta(\hat{s})$. The implicit function theorem now guarantees that there is a unique differentiable function $\hat{t}(\hat{\mathbf{w}})$ defined in a neighbourhood of $\hat{\mathbf{w}}=0$ such that $F(\hat{s}+\hat{t}(\hat{\mathbf{w}}),\hat{\mathbf{v}}+\hat{\mathbf{w}})=0$. Computing the first order approximation of this function we see that for $\hat{\mathbf{w}}$ small enough

[TABLE]

where h.o.t. stands for higher order terms. From (82) we obtain

[TABLE]

Set $\hat{\mathbf{v}}^{(k)}=\mathbf{S}^{-1}\hat{\mathbf{u}}^{(k)}$ and $\hat{\mathbf{w}}^{(k)}=\hat{\mathbf{v}}^{(k)}-\hat{\mathbf{v}}$. Now

[TABLE]

Since $\mathbf{S}\hat{\mathbf{w}}^{(k)}=\hat{\mathbf{u}}^{(k)}-\hat{\mathbf{u}}_{\mathcal{A}}$ we have that

[TABLE]

Equations (80) and (81) now follow from (85) and the fact that $\hat{\mu}^{(k)}$ is asymptotically given as a linear function of $\hat{\mathbf{u}}^{(k)}$. ∎

Adapting Theorem 3 to the context of Algorithm 1 we obtain the following Corollary.

Corollary 1**.**

Let $\mathbf{u}_{\mathcal{A}}\in W_{\mathcal{A}}^{N}$ be a fixed point of the Algorithm 1 and $(s,\mathbf{v})\in W_{\mathcal{A}}\times W_{\mathcal{A}}^{N}$ a corresponding solution to (50). Let $\mu_{\mathcal{A}}\in W_{\mathcal{A}}$ be such that $P_{\mathcal{A}}(s\mu_{\mathcal{A}})=1$. Assume that $s_{*}:=\inf_{y\in\Gamma}s(y)>0$ and let $\psi_{\max}$ be as in Theorem 3. Then the iterates of Algorithm 1 satisfy

[TABLE]

for $\mathbf{u}^{(k)}$ sufficiently close to $\mathbf{u}_{\mathcal{A}}$. Furthermore, there exists $C>0$ such that

[TABLE]

Proof.

Interpret Theorem 3 in the context of Algorithm 1. The bound $\phi_{\min}\geq s_{*}$ is a consequence of Lemma 1. ∎

Obviously the previous Corollary has practical value only if $\psi_{\max}<s_{*}$. Here we will briefly discuss the value of $\psi_{\max}$ in the case that $(\hat{s},\hat{\mathbf{v}})\in\mathbb{R}^{P}\times\mathbb{R}^{PN}$ is associated to the dominant fixed point of Algorithm 2. Observe that the equation $\hat{\mathbf{z}}=\mathbf{R}(\hat{s},\hat{\mathbf{v}})\hat{\mathbf{w}}$ is equivalent to the system

[TABLE]

for all $y\in\Gamma$. We see that, if $\mathbf{w}=\mathbf{v}$ then $\mathbf{z}=0$, whereas, if $\langle\mathbf{w}(y),\mathbf{v}(y)\rangle_{\mathbb{R}^{N}_{\mathbf{M}}}=0$ for all $y\in\Gamma$ then $\mathbf{z}=\mathbf{w}$. Thus, the matrix $\mathbf{R}(\hat{s},\hat{\mathbf{v}})$ acts as a deflation that shrinks vectors that are close to $\mathbf{v}(y)$ and preserves vectors that are almost orthogonal to $\mathbf{v}(y)$. From Proposition 4 we know that $s^{-1}(y)$ is an approximation of the smallest eigenvalue of $\mathbf{M}^{-1}\mathbf{K}(y)$ and $\mathbf{v}(y)$ is an approximation of the corresponding eigenvector. By Lemma 1 the operator norm of $\mathbf{S}^{-1}$ is bounded by $\sup_{y\in\Gamma}\lambda_{1}^{-1}(y)$, where $\lambda_{1}(y)$ is the smallest eigenvalue of $\mathbf{M}^{-1}\mathbf{K}(y)$. Analogously, since the eigenvector corresponding to this smallest eigenvalue is deflated by $\mathbf{R}(\hat{s},\hat{\mathbf{v}})$, we expect the norm of $\mathbf{R}(\hat{s},\hat{\mathbf{v}})\mathbf{S}^{-1}$ to be bounded by a value close to $\sup_{y\in\Gamma}\lambda_{2}^{-1}(y)$, where $\lambda_{2}(y)$ is the second smallest eigenvalue of $\mathbf{M}^{-1}\mathbf{K}(y)$. With this reasoning, if the deflation is sufficient, there is $\psi_{\max}^{*}\in\mathbb{R}$ such that

[TABLE]

The speed of convergence of the spectral inverse iteration is now characterized by what is essentially the largest ratio of the two smallest eigenvalues of the parametric matrix $\mathbf{M}^{-1}\mathbf{K}(y)$.

5.2.4 Combined error analysis

Let $(\mu,u)\in L^{2}_{\nu}(\Gamma)\times L^{2}_{\nu}(\Gamma)\otimes H_{0}^{1}(D)$ be the smallest eigenvalue and the associated eigenfunction of the continuous problem (9). Let $(\mu_{h},u_{h})\in L^{2}_{\nu}(\Gamma)\times L^{2}_{\nu}(\Gamma)\otimes V_{h}$ be the corresponding eigenpair of the semi-discrete problem (13). Assume that there exists a dominant fixed point $\mathbf{u}_{\mathcal{A}}\in W_{\mathcal{A}}^{N}$ of Algorithm 1 and an associated eigenvalue approximation $\mu_{h,\mathcal{A}}:=\mu_{\mathcal{A}}\in W_{\mathcal{A}}$ as in Proposition 4. Denote by $\mathbf{u}^{(k)}\in W_{\mathcal{A}}^{N}$ the $k$:th iterate of Algorithm 1 and by $\mu_{h,\mathcal{A},k}:=\mu^{(k)}\in W_{\mathcal{A}}$ the associated solution to (42). Let $u_{h,\mathcal{A}}$ and $u_{h,\mathcal{A},k}$ denote the functions in $W_{\mathcal{A}}\otimes V_{h}$, whose coordinates are defined by the vectors $\mathbf{u}_{\mathcal{A}}$ and $\mathbf{u}^{(k)}$ respectively. The spatial, stochastic, and iteration errors may now be separated in the following sense:

[TABLE]

and

[TABLE]

Under sufficient conditions we may now bound each term in the equations (89) and (5.2.4) separately using the theory developed earlier in this section. The first term may be approximated using Theorem 1, the second term may be approximated using Theorem 2 and Proposition 2, and the third term may be approximated using Corollary 1 of Theorem 3 and the hypothesis (88). We therefore expect that, with an optimal choise the multi-index sets $\mathcal{A}_{\epsilon}$ for $\epsilon>0$, the output of the spectral inverse iteration converges to the exact solution according to

[TABLE]

and similarly

[TABLE]

for certain rates $r>0$ and $l>0$.

5.3 Numerical examples

We present numerical evidence to verify the equations (91) and (92). In each of the following examples we compute the smallest eigenvalue and the corresponding eigenfunction of the model problem (4) in the unit square $D=[0,1]^{2}$ using the Algorithm 1. We use the smallest eigenvector at $y=0$ as an initial guess. For the diffusion coefficient we assume the form (5) with $a_{0}:=1$ and

[TABLE]

where we set $\varsigma=3.2$. Now $||a_{m}||_{L^{\infty}(D)}\leq Cm^{-\varsigma}$ and $||a_{m}||_{W^{2,\infty}(D)}\leq Cm^{-\varsigma+2}$ so that the assumptions (6) – (8) for $s=2$ are satisfied with $p_{0}>\varsigma^{-1}$ and $p_{2}>(\varsigma-2)^{-1}$. We therefore expect the regions of analyticity in Proposition 1 to increase according to $\tau_{m}\geq Cm^{\varsigma-1}$.

The deterministic mesh is a uniform grid of second order quadrilateral elements in all computations. The discretization in the parameter space is obtained by defining the multi-index sets

[TABLE]

where we set

[TABLE]

as in the proof of Proposition 2. Multi-index sets of this form have been introduced in [6] and in [15] an algorithm for generating them has been suggested.

We use a matrix free formulation of the conjugate gradient method for solving the linear systems (46) and (48). The preconditioner is constructed using the mean of the parametric matrix in question [17] and as an initial guess we set the solution of the system from the previous iteration. We wish to note that in this setting only a very few iterations of the conjugate gradient method are needed at each step of the spectral inverse iteration.

In the lack of an exact solution we compute an overkill solution $(\mu_{*},u_{*})$ for which the number of deterministic degrees of freedom is $N=36741$, the parameter $\epsilon$ is chosen such that $\#\mathcal{A}_{\epsilon}=264$, and the number of iterations is $k=16$. This results in roughly $10^{7}$ total degrees of freedom. The number of active dimensions in the overkill solution is $M(\mathcal{A}_{\epsilon})=113$. All the numerical examples in this section have been computed using this overkill solution as a reference. The expected value and variance of the eigenfunction are presented in Figure 1.

5.3.1 Convergence in space

Keeping the number of stochastic degrees of freedom $\#\mathcal{A}_{\epsilon}=264$ and the number of iterations $k=16$ fixed, we may investigate the convergence of the solution $(\mu_{*,h},u_{*,h})$ as a function of the spatial discretization parameter $h$. This convergence for piecewise quadratic basis functions is illustrated in Figure 2. We observe algebraic convergence rates of order 3 and 4 for the eigenfunction and eigenvalue respectively, exactly as predicted by Theorem 1. Thus, the error behaves like $N^{-3/2}$ and $N^{-2}$ with respect to the number of deterministic degrees of freedom.

5.3.2 Convergence in the parameter domain

Keeping the number of spatial degrees of freedom $N=36741$ and the number of iterations $k=16$ fixed, we may investigate the convergence of the solution $(\mu_{*,\mathcal{A}_{\epsilon}},u_{*,\mathcal{A}_{\epsilon}})$ as a function of $\#\mathcal{A}_{\epsilon}$ as $\epsilon\to 0$. This convergence is illustrated in Figure 3. We observe approximate algebraic convergence rates of order $-r=-1.9$ with respect to the number of stochastic degrees of freedom $\#\mathcal{A}_{\epsilon}$.

In Figure 4 we have presented the norms of the Legendre coefficients of the overkill solution. The ordering of the coefficients is the same as the order in which they would appear in the multi-index set $\#\mathcal{A}_{\epsilon}$ as $\epsilon\to 0$. We see that the norms converge at the rate $-r-1/2=-2.4$ exactly as we would expect from the proof of Proposition 2. In Figure 5 we have presented the norms of the same Legendre coefficients sorted by decreasing magnitude. From this Figure we estimate that, with an optimal selection of the multi-index sets we could in fact observe a rate of convergence $-r=-2.3$ for the error of the solution. This ideal rate of convergence is somewhat faster than the asymptotic theoretical bound of $-r=-\varsigma+3/2=-1.7$ predicted by Proposition 2.

Interestingly we observe two well separated clusters of values in Figure 4b. It seems that many of the multi-indices that correspond to relatively large Legendre coefficients of the eigenfunction, account only for a marginal contribution to the eigenvalue.

5.3.3 Convergence of the iteration error

Keeping the number of spatial basis functions $N=36741$ and the parameter $\epsilon$ fixed so that $\#\mathcal{A}_{\epsilon}=264$, we may investigate the convergence of the solution $(\mu_{*,k},u_{*,k})$ as a function of the number of iterations $k$. This convergence is illustrated in Figure 6. Assuming that the variation in the eigenvalues within the parameter space is small, the value $\lambda_{1/2}$ defined in (88) may be approximated by the ratio of the two smallest eigenvalues of the problem at $y=0$. Thus, Figure 6 suggests that the error behaves asymptotically like $\lambda_{1/2}^{k}$, just as predicted by Corollary 1.

It is worth noting that, from the analysis of the classical inverse iteration, one might expect the eigenvalue to converge faster than the eigenfunction. In fact, the eigenvalue exhibits a faster rate of convergence at first and the error behaves like $\lambda_{1/2}^{2k}$. Comparing to the results of the previous example, we see that $k\approx 9$ represents a turning point after which the stochastic error in the eigenfunction starts to dominate the iteration error. Hence, for $k\geq 9$ the polynomial approximation in the parameter domain is insufficient to guarantee the degree of accuracy that is required for the eigenvalue to exhibit the faster rate of convergence that is otherwise characteristic to it.

5.3.4 Concluding remarks and comparison to sparse collocation

Using the finest levels of discretization, i.e., $N=9296$ degrees of freedom for approximation in space and $\#\mathcal{A}_{\epsilon}=121$ degrees of freedom for approximation in the parameter domain, and computing $k=9$ steps of the inverse iteration we obtain a solution for which the $L^{2}_{\nu}(\Gamma)\otimes L^{2}(D)$ error of the eigenfunction is approximately $3\cdot 10^{-6}$. The number of total degrees of freedom in this case is more than $10^{6}$ and the number of active dimensions is $M(\mathcal{A}_{\epsilon})=60$. The total computational time on a standard desktop machine is approximately five minutes, most of which is spent in the conjugate gradient method for the linear systems (46) and (48).

When the solution computed via the spectral inverse iteration is compared to the results of the non-composite version of the sparse collocation method introduced in [4] and employed in e.g. [1] (see equations (5.12) – (5.13) and (5.16) – (5.17)), the statistics of the two solutions seem to almost coincide. Again using the finest levels of discretization ($N=9296$ and $\#\mathcal{A}_{\epsilon}=121$) for both methods, the $L^{2}(D)$ errors of mean and variance of the eigenfunction are both less than $3\cdot 10^{-8}$ and the errors in the eigenvalue are less than $3\cdot 10^{-11}$ and $3\cdot 10^{-9}$ for the mean and variance respectively.

6 Spectral subspace iteration

In this section we extend the spectral inverse iteration to a spectral subspace iteration, with which we can compute dominant subspaces of the inverse of the parametric matrix under consideration. The underlying assumption is that the subspace is sufficiently smooth with respect to the parameters. Convergence of the spectral subspace iteration is verified through numerical experiments.

6.1 Algorithm description

As with the classical subspace iteration, the idea in the spectral version is to perform inverse iteration for a set of vectors and orthogonalize these vectors at each step. Orthogonality should here be understood in a sense that the vectors are orthogonal for all points in the parameter space $\Gamma$. This can be approximately achieved by performing the Gram-Schmidt orthogonalization process for the vectors in the Galerkin sense, i.e., by projecting each elementary operation to the basis $W_{\mathcal{A}}$.

Fix a finite set of multi-indices $\mathcal{A}\subset(\mathbb{N}_{0}^{\infty})_{c}$ and let $P=\#\mathcal{A}$. The spectral subspace iteration for the system (16) is now defined in Algorithm 3. Observe that, if the projections were precise, then the Algorithm would correspond to performing the classical subspace iteration pointwise on $\Gamma$. Orhtogonalization of the basis vectors via the Gram-Schmidt process is achieved in step (2). We expect Algorithm 3 to converge to an approximate basis for the $Q$-dimensional invariant subspace associated to the smallest eigenvalues of the system.

Algorithm 3** (Spectral subspace iteration).**

Fix $tol>0$ and let $\{\mathbf{u}^{(0,q)}\}_{q=1}^{Q}\subset W_{\mathcal{A}}^{N}$ be an initial guess for the basis of the subspace. For $k=1,2,\ldots$ do

1. (1)

For each $q=1,\ldots,Q$ solve $\mathbf{v}^{(q)}\in W_{\mathcal{A}}^{N}$ from the linear equation

[TABLE] 2. (2)

For $q=1,\ldots,Q$ do

1. (2.1)

Set

[TABLE] 2. (2.2)

Solve $s^{(q)}\in W_{\mathcal{A}}$ from the nonlinear equation

[TABLE] 3. (2.3)

Solve $\mathbf{u}^{(k,q)}\in W_{\mathcal{A}}^{N}$ from the linear equation

[TABLE] 3. (3)

Stop if a suitable criterion is satisfied and return $\{\mathbf{u}^{(k,q)}\}_{q=1}^{Q}\subset W_{\mathcal{A}}^{N}$ as the approximate basis for the subspace.

In general we can not expect the output vectors $\{\mathbf{u}^{(k,q)}\}_{q=1}^{Q}\subset W_{\mathcal{A}}^{N}$ of Algorithm 3 to converge to any particular eigenvectors of the system (16). However, we still expect them to approximately span the subspace associated to the smallest eigenvalues of the system. In view of Section 3, if a cluster of eigenvalues is sufficiently well separated from the rest of the spectrum, then we assume the associated subspace to be analytic with respect to the parameter vector $y\in\Gamma$. In this case we may expect optimal convergence of the projections in the Algorithm.

Remark 6**.**

In order to measure convergence of the Algorithm 3 we should be able to estimate the angle between subspaces over the parameter space $\Gamma$. It is not entirely trivial to perform this kind of a computation in practise. The numerical examples in Section 6.2 will hopefully give some more insight on this.

Remark 7**.**

As noted in Section 3, the smallest eigenvalue of the problem (9) is always simple, hence analytic. For more general problems this might not be the case. For instance, in the event of an eigenvalue crossing, the eigenmode corresponding to the pointwise smallest eigenvalue is not (in general) even a continuous function of the parameter vector $y$. In this case we can modify the Algorithm 3 by adding the step

• (2.0)

Set $\mathbf{v}^{(1)}=\sum_{q=1}^{Q}\mathbf{v}^{(q)}$

before step (2.1). This should ensure optimal convergence, since even if the eigenmodes change places, we still expect their sum to be smooth with respect to $y$.

Using the tensors defined in Section 5 we may write Algorithm 3 in the following form.

Algorithm 4** (Spectral subspace iteration in tensor form).**

Fix $tol>0$ and let $\{\hat{\mathbf{u}}^{(0,q)}\}_{q=1}^{Q}\subset\mathbb{R}^{PN}$ be an initial guess for the basis of the subspace. For $k=1,2,\ldots$ do

1. (1)

For each $q=1,\ldots,Q$ solve $\hat{\mathbf{v}}^{(q)}\in\mathbb{R}^{PN}$ from the linear system

[TABLE] 2. (2)

For $q=1,\ldots,Q$ do

1. (2.1)

Set

[TABLE] 2. (2.2)

Solve $\hat{s}^{(q)}\in\mathbb{R}^{P}$ from the nonlinear system

[TABLE]

with the initial guess $s^{(q)}_{\alpha}=||\hat{\mathbf{w}}^{(q)}||_{\mathbb{R}^{P}\otimes\mathbb{R}^{N}_{\mathbf{M}}}\delta_{\alpha 0}$ for $\alpha\in\mathcal{A}$. 3. (2.3)

Solve $\hat{\mathbf{u}}^{(k,q)}\in\mathbb{R}^{PN}$ from the linear system

[TABLE] 3. (3)

Stop if a suitable criterion is satisfied and return $\{\hat{\mathbf{u}}^{(k,q)}\}_{q=1}^{Q}\subset\mathbb{R}^{PN}$ as the approximate basis for the subspace.

6.2 Numerical examples

We use Algorithm 3 to compute the 3-dimensional subspace associated with the smallest eigenvalues of the model problem considered in Section 5.3. We let the deterministic mesh be a uniform grid of second order quadrilateral elements with $N=2465$ degrees of freedom. As an initial guess we use the smallest eigenvectors of the problem at $y=0$. In Figure 7 we have presented the four smallest eigenvalues of the problem as a function of $y_{1}$, when $y_{2},y_{3},\ldots$ are held constant. We observe an eigenvalue crossing due to which the eigenvectors corresponding to the pointwise second and third smallest eigenvalues are discontinuous as functions of $y$.

In order to investigate the convergence of the spectral subspace iteration, we attempt to estimate the angle between the exact invariant subspace and the approximate one computed by Algorithm 3. For any fixed $y\in\Gamma$ we let $\mathbf{v}_{1}(y),\ldots,\mathbf{v}_{Q}(y)$ be a set of $\mathbb{R}^{N}_{\mathbf{M}}$ orthonormal exact eigenvectors corresponding to the $Q$-smallest eigenvalues of the problem. We define

[TABLE]

where $\Theta^{(k)}(y)\in R^{Q\times Q}$ is a matrix with elements $\Theta^{(k)}_{ij}(y)=\langle\mathbf{u}^{(k)}_{i}(y),\mathbf{v}_{j}(y)\rangle_{\mathbb{R}^{N}_{\mathbf{M}}}$. Now $\theta_{k}(y)$ can be viewed as the cosine of the angle between the two subspaces at $y\in\Gamma$ (see for instance [9] formula (2.2)). Thus, convergence of the algorithm can be measured in terms of the statistics of $\theta_{k}$. In the following examples we have estimated the mean and variance of $\theta_{k}$ using the non-composite version of the sparse collocation operator employed in [1]. For the definition of the collocation operator we have used the overkill multi-index set of section 5.3 ($\#\mathcal{A}_{\epsilon}=264$).

Convergence of the spectral subspace iteration for $Q=3$ is illustrated in Figure 8. We see that the values $\arccos(\mathbb{E}[\theta_{k}])$ behave like $\lambda_{3/4}^{k}$, where $\lambda_{3/4}$ is the ratio of the third and fourth smallest eigenvalues of the problem. Simultaneously the values $\mathrm{Var}[\theta_{k}]$ converge to zero. These results suggest that the angle between the exact subspace and the approximation computed by Algorithm 3 converges to zero on $\Gamma$. Furthermore, the rate of convergence is characterised by the rate $\lambda_{3/4}^{k}$ much like for the classical subspace iteration. Note however, that with a fixed basis for polynomial approximation, i.e. a fixed multi-index set $\mathcal{A}_{\epsilon}$, only a certain accuracy for the output may be reached. Increasing the number of basis polynomials makes more accurate solutions achievable.

7 Conclusions and future prospects

We have presented a comprehensive error analysis for the spectral inverse iteration, when applied to solving the ground state of a stochastic elliptic operator. We have also proposed a method of spectral subspace iteration and, using numerical examples, shown its potential in computing approximate subspaces associated to possibly clustered eigenvalues. Further analysis, both numerical and theoretical, of this algorithm is left for future research.

The numerical examples suggest that our algorithms are both accurate and efficient. However, theoretical estimates for the computational complexity are not entirely trivial to obtain as this would require information on the structure of the tensor of coeffiecients $c_{\alpha\beta\gamma}$. Moreover, when iterative solvers are used, the optimal strategy is to increase the associated tolerances in the course of the iteration. We note that sparse products of the spatial and stochastic approximation spaces, as in [15], may be applied to further reduce the computational effort, and that matrix free algorithms also allow for easy parallelization.