{'fraction_non_alphanumeric': 0.037366612407306925, 'fraction_numerical': 0.00991137637909206, 'mean_word_length': 4.543613394826549, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Factorizing the Stochastic Galerkin System

Recent work has explored solver strategies for the linear system of equations arising from a spectral Galerkin approximation of the solution of PDEs with parameterized (or stochastic) inputs. We consider the related problem of a matrix equation whose matrix and right hand side depend on a set of parameters (e.g. a PDE with stochastic inputs semidiscretized in space) and examine the linear system arising from a similar Galerkin approximation of the solution. We derive a useful factorization of this system of equations, which yields bounds on the eigenvalues, clues to preconditioning, and a flexible implementation method for a wide array of problems. We complement this analysis with (i) a numerical study of preconditioners on a standard elliptic PDE test problem and (ii) a fluids application using existing CFD codes; the MATLAB codes used in the numerical studies are available online.

1. Introduction. Complex engineering models are often described by systems of equations, where the model outputs depend on a set of input parameters. Given values for the input parameters, the model output may be computed by solving a discretized differential equation, which typically involves the solution (or series of solutions) of a system of linear equations for the unknowns. Each of these computations may be very expensive depending on factors such as grid resolution or physics components in the model. As the role of simulation gains prominence in areas such as decision support, design optimization, and predictive science, understanding the effects of variability in the input parameters on variability in the model output becomes more important. Exhaustive parameter studies (e.g. uncertainty quantification, sensitivity analysis, model calibration) can be prohibitively expensive -particularly for a large number of input parameters -and therefore accurate interpolation and robust surrogate models are essential. We consider the model problem of a matrix equation whose matrix and right hand side depend on a set of parameters. Let s ∈ [−1, 1] d be a set of input parameters from an idealized parameter space. We equip the parameter space with a bounded, separable, positive weight function ω : [−1, 1] d → R + , where ω(s) = ω 1 (s 1 ) · · · ω d (s d ).
(In a probabilistic context, this weight function represents a probability measure on the input parameter space.) Let A(s) be an N × N matrix-valued function where we assume that A(s) is invertible for all s ∈ [−1, 1] d , and let b(s) be the vector-valued function with N square-integrable components. We seek the vector valued function x(s) that satisfies Such parameterized matrix problems often arise as an intermediate step when computing an approximate solution of a complex model with multiple input parameters. They appear in such diverse fields as differential equations with random (or parameterized) inputs [1,31], electronic circuit design [23], image deblurring models [5], and ranking methods for nodes in a graph [25,7]. Given a parameterized matrix equation, one may wish to approximate the vector valued function that satisfies the equation or estimate its statistics. Once an approximation is constructed that is cheaper to evaluate than the true solution, statistics -such as mean and variance -of the approximation represent estimates of statistics of the true solution.
The approximation model is a vector of multivariate polynomials represented as a series of orthonormal polynomial basis functions where each basis function is a product of univariate orthonormal polynomials. We employ the standard multi-index notation; let α = (α 1 , . . . , α d ) ∈ N d be a multi-index, and define the basis polynomial π α (s) = π α1 (s 1 ) · · · π α d (s d ). (1. 2) The polynomial π αi (s i ) is the orthonormal polynomial of degree α i , where the orthogonality is defined with respect to the weight function ω i (s i ). Then for α, β ∈ N d , where equality between multi-indicies means component-wise equality. For a given index set I ⊂ N d with size |I| < ∞, the polynomial approximation can be written The N -vector x α is the coefficient of the series corresponding to π α (s). The N × |I| matrix X has columns x α , and the parameterized vector π(s) contains the basis polynomials. The goal of the approximation method is to compute the unknown coefficients X. Such polynomial models have become popular for approximating the solution of PDEs with random inputs; they appear under names such as polynomial chaos methods [36], stochastic finite element methods [18], stochastic Galerkin methods [2], and stochastic collocation methods [35,3]. The Galerkin methods compute the series coefficients such that the equation residual is orthogonal to the approximation space defined by the index set I; they typically employ a full polynomial basis of order n given by although such a basis set is not strictly necessary. The number of terms in this basis set is |I n | = n+d n , which grows exponentially with the number of input parameters d. The process of computing the coefficients involves solving a linear system of size N |I|× N |I|, which can be prohibitively large for even a moderate number of input parameters (6 to 10) and low order polynomials (degree > 5). For this reason, there has been a flurry of recent work on solver strategies for the matrix equations arising from the Galerkin methods [27,32,21,10,12,13], including papers on preconditioning [30,11,29,33]. Such work has relied on knowing the matrix valued coefficients A α of a series expansion of the parameterized matrix A(s), which typically come from a specific form of the coefficients of the elliptic operator in the PDE models. Another drawback of the Galerkin method is its limited ability to take advantage of existing solvers for the problem A(λ)x(λ) = b(λ) given a parameter point λ ∈ [−1, 1] d . In contrast, the pseudospectral and collocation methods can compute coefficients of the model (1.4) using only evaluations of the solution vector x(λ). The chosen points typically correspond to a multivariate quadrature rule, and computing the coefficients of the polynomial model (1.4) becomes equivalent to approximating its Fourier coefficients with the quadrature rule [3,8]. Thus, from the point of view of code reuse and rapid implementation, the pseudospectral/collocation methods have a distinct advantage.
We propose a variant of the Galerkin method that alleviates the drawbacks of limited code reuse and memory limitations. By formally replacing the integration in the Galerkin method by a multivariate quadrature rule -a step that is more often than not performed in practice -we derive a decomposition of the linear system of equations used to compute the Galerkin coefficients of (1.4); such a method has been called Galerkin with numerical integration in the context of numerical methods for PDEs [4]. This decomposition allows us to compute the Galerkin coefficients using only evaluations of the parameterized matrix A(λ) and parameterized vector b(λ) for points λ ∈ [−1, 1] d corresponding to a quadrature rule. In fact, if one requires only matrix-vector multiplies as in Krylov-based iterative methods (e.g. cg [20] or minres [26] solvers) for the Galerkin system of equations, this restriction can be relaxed to matrix-vector products A(λ)v for a given N -vector v; there is no need for the coefficients of an expansion of A(s) as in (1.6). Therefore the method takes full advantage of sparsity of the parameterized system resulting in reduced memory requirements. The decomposition also yields insights for preconditioning the Galerkin system that generalize existing work. Additionally, the decomposition immediately reveals bounds on the eigenvalues of the Galerkin system for the case of symmetric A(s).
The remainder of the paper is structured as follows. In Section 2, we derive the Galerkin method for a given problem and basis set. We then derive the decomposition of the matrix used to compute the Galerkin coefficients and examine its consequences including bounds on the eigenvalues of the matrix, strategies for simple implementation and code reuse, and insights into preconditioning. We then provide a numerical study of various preconditioners suggested by the decomposition using the common test case of an elliptic PDE with parameterized coefficients in Section 3; the codes for the numerical study are available in a Matlab suite of tools accessible online. To emphasize the advantages of code reuse, we apply the method to an engineering test problem using existing codes in Section 4. Finally, we conclude in Section 5 with summarizing remarks.

Notation.
For the remainder of the paper, we will use the bracket notation · to denote a discrete approximation to the integral with respect to the weight function ω. In other words, for functions f : where the points λ β = (λ β1 , . . . , λ β d ) ∈ [−1, 1] d and weights ν β ∈ R define a multivariate quadrature rule for multi-indicies in the set J ⊂ N d . We will abuse this notation by putting matrix and vector valued functions inside the brackets, as well. For example, A denotes the mean of A(s) computed with a quadrature rule.

A Spectral Galerkin Method.
We first review the Galerkin method [8] for computing the coefficients of the polynomial model (1.4). Define the residual and let x g (s) be the Galerkin approximation. We require that each component of the residual r i (x g , s) be orthogonal to the approximation space defined by the span of π α with α ∈ I, We can combine the equations in (2.2) using the matrix notation or equivalently, upon substituting the model x g (s) = Xπ(s), Using the vec notation [19, Section 4.5], we can rewrite (2.4) as where x = vec(X) is an N |I| × 1 constant vector equal to the columns of X stacked on top of each other. The constant matrix ππ T ⊗ A has size N |I| × N |I| and a distinct block structure; the α, β block of size N × N is equal to π α π β A for multiindicies α, β ∈ I. Similarly, the α block of the N |I| × 1 vector π ⊗ b is equal to bπ α , i.e. the Fourier coefficient of b associated with π α (s) approximated with the quadrature rule.
Much of the literature on PDEs with random inputs [27,30] points out an interesting block sparsity pattern that arises in the matrix ππ T ⊗ A when the parameterized matrix A(s) depends at most linearly on any componenets of s; this is a result the orthogonality of the bases π(s). For general analytic dependence on the parameters, such sparsity patterns do not appear [14]. However, one can always mimick the sparsity pattern of A(s) with a simple reordering of the variables. By taking the transpose of (2.4) and using the same vec operations, if we definex = vec(X T ) (i.e. the same unknowns reordered), then Notice that the matrix in (2.6) retains the sparsity of A(s) in its blocks, since each For the remainder of the paper, however, we will work with the form (2.5). By writing out the numerical integration rule for the integrals used to form ππ T ⊗ A , we uncover an interesting decomposition, which we state as a theorem.
Theorem 2.1. Let {(λ β , ν β )} with β ∈ J be a multivariate quadrature rule. The matrix ππ T ⊗ A can be decomposed as where I is the N × N identity matrix, and Q is a matrix of size |I| × |J | -one row for each polynomial basis and one column for each point in the quadrature rule. The Proof. Writing out the quadrature rules, Notice that the elements of the vector π(λ β ) are the polynomials π α (s) evaluated at the quadrature point λ β . If we define the vectors then we have Let Q be the matrix whose columns are q β , and define the block diagonal matrix A(λ) with diagonal blocks equal to A(λ β ). Then for an N × N identity matrix I, we can rewrite (2.10) as (2.7), as required.
As an aside, we note that if A(s) depends polynomially on the parameters s, then each integrand in the matrix ππ T ⊗ A is a polynomial in s. Therefore by the polynomial exactness, there is a Gaussian quadrature rule such that the numerical integration approach exactly recovers the true Galerkin matrix.
The elements of Q are the orthogonal polynomials evaluated at the quadrature points and multiplied by the square root of the quadrature weights. Thus, for two polynomials π α (s) and π β (s) with α, β ∈ I, define r α and r β to be the corresponding rows of Q; then If the quadrature rule is a tensor product Gaussian quadrature rule of sufficiently high order to exactly compute the polynomial integrand, then this implies QQ T = I, where I is the |I| × |I| identity matrix; we will assume from here on that the chosen quadrature rule yields this property. We will also assume that the number of basis polynomials is less than the number of points used in the quadrature rule, i.e. |I| ≤ |J |; otherwise the Galerkin matrix ππ T ⊗ A will be rank deficient.

12)
where θ(X) denotes the eigenvalues of a matrix X, and θ min (X) and θ max (X) denote the smallest and largest eigenvalues of X, respectively. Proof. Let {(λ β , ν β )} be a tensor product Gaussian quadrature rule of sufficiently high order so that QQ T = I. Letπ(s) be a vector containing the polynomial bases such that [π(s) T ,π(s) T ] T contains the tensor product polynomial basis corresponding to the tensor product quadrature grid {λ β }. LetQ be the matrix with columns 2.1. Iterative Solvers. Due to their size and sparsity, the preferred way of solving (2.5) is with Krylov based iterative solvers [27] that rely on matrix-vector products with the matrix ππ T ⊗ A . By employing the decomposition in Theorem 2.1, we can compute these using only multiplication of the parameterized matrix evaluated at the quadrature point by a given vector. More precisely, given a vector u = vec(U), suppose we wish to compute (2.14) We accomplish this in three steps using the properties of the Kronecker product: Step 1 can be thought of as pre-processing, and step 3 as post-processing. In practice, each row of Q may have a Kronecker structure corresponding to the tensor product quadrature rule. In this case, steps 1 and 3 can be computed accurately and efficiently using multiplication methods such as [15]. The second step requires only constant matrix-vector products where the matrix is A(s) evaluated at the quadrature points. Therefore we can take advantage of a memory-efficient, reusable interface for the matrix-vector multiplies that will exploit any sparsity in the matrix. We reiterate that this can be accomplished without any knowledge of the specific type of parameter dependence in A(s).

Preconditioning
Strategies. The number of iterations required to achieve a given convergence criterion (e.g. a sufficiently small residual) can be greatly reduced for Krylov-based iterative methods with a proper preconditioner. In general, preconditioning a system is highly problem dependent and begs for the artful intuition of the scientist. However, the structure revealed by the decomposition from Theorem 2.1 offers a number of useful clues.
Suppose we have an N × N matrix P that is easily invertible. We can construct a block-diagonal preconditioner I ⊗ P −1 , where I is the identity matrix of size |I| × |I|. If we premultiply the preconditioner against the factored form of ππ T ⊗ A , we get Notice that by the mixed product property and commutativity of the identity matrix, the block-diagonal preconditioner slips past Q⊗I to act directly on the parameterized matrix evaluated at the quadrature points. The blocks on the diagonal of the inner matrix product are P −1 A(λ β ) for β ∈ J . In other words, we can choose one constant matrix P to affect the parameterized system at all quadrature points. A reasonable and popular choice is the mean P = A ; see [29,27] for a detailed analysis of this preconditioner for stochastic finite element systems. Notice that this is also the first N × N block of the Galerkin matrix. However, if A(s) is very large or has some complicated parametric dependence, then forming the mean system and inverting it (or computing partial factors) for the preconditioner may be prohibitively expensive. If the dependence of A(s) on s is close to linear, then P = A ( s ) may be much easier to evaluate and just as effective.
One goal of the preconditioner is reduce the condition number of the matrix, and one way of achieving this is to reduce the spread of the eigenvalues. If we knew a priori which region of the parameter space produced the extrema of the parameterized eigenvalues of A(s) (e.g. the boundaries), then we could choose an appropriate parameter value to construct an effective preconditioner. Unfortunately, we only get to use one such evaluation. Therefore, if the largest possible value of the parameterized eigenvalues is very large, we may choose this parameter value. Alternatively, if the smallest eigenvalue is close to zero (for positive definite systems), then this may be a better option to reduce the condition number of the Galerkin system. In the next section, we explore a few choices for the preconditioner P on a standard test problem.
3. Preconditioning Study. Consider the following parameterized elliptic partial differential equation with homogeneous Dirichlet boundary conditions; variations of this problem can be found in [1,2,16], amongst others. We seek a solution u(x, s) that satisfies where u = 0 on the boundary of the square domain [0, 1] 2 , and the parameter space is the hypercube s ∈ [−1, 1] 4 equipped with a uniform measure. The logarithm of the elliptic coefficient is given by a truncated Karhunen-Loeve like expansion [24] of a zero-mean random field with covariance We discretize (3.1) using the finite element method implemented in Matlab's pde toolbox on an irregular mesh of N = 1,921 triangles. We solve the eigenvalue problem with the discrete covariance matrix to compute the values of the eigenfunctions ψ k (x) on the given mesh. The square roots of first ten eigenvalues σ k are plotted in Figure 3.1; we truncate the expansion at d = 4 to retain roughly 90% of the energy of the field. Given a point in the parameter space, the PDE Toolbox allows us to access the stiffness matrix. Therefore we can apply Matlab's minres solver using only matrix-vector multiplies against the stiffness matrix evaluated at the quadrature points. For the polynomial basis, we use the normalized multivariate product Legendre polynomials of order 5, which includes |I| = 126 basis functions. Therefore the number of unknowns is N |I| = 242,046. We use a tensor product Gauss-Legendre quadrature rule of order 12 in each parameter (20,736 points) to implicitly form the The numbers indicate the percentage of field energy captured when using that many terms in the expansion. Galerkin matrix; this is more than sufficient to maintain the orthogonality in the basis functions.
To test the preconditioning ideas suggested by the decomposition from Theorem 2.1 and equation (2.15), we try five different choices of P. For each P, we precompute the Cholesky factors to apply (2.15) efficiently.  Table 3.1 summarize the results, which we now explain in detail. In the first type of preconditioning, we set P = A(s r ) for a random point s r in [−1, 1] 4 . To sample the likely effectiveness of any random point, we select 25 random points and evaluate each choice of s r . Second, we set P = A(s max ) where A(s max ) is the matrix with the largest eigenvalue in [−1, 1] 4 . To compute s max , we use up to 100 iterations of the power method to estimate the largest eigenvalue at each point in the quadrature rule, that is for each λ β for β ∈ J ; we also evaluate A(s) for the tensor product of the endpoints as well. Third, we set P = A(s min ) where A(s min ) is the matrix with the smallest eigenvalue in [−1, 1] 4 . For this computation, we use Matlab's optimization routine fmincon and use eigs/arpack to evaluate the smallest eigenvalue [22]. Fourth, we set P = A(s mid ) where s mid is the midpoint of the domain, which is the origin for our experiments. Fifth, and finally, we set P = A , the mean preconditioner. To evaluate the mean, we use a 2nd order quadrature rule for a fast approximation and a 5th order quadrature rule for a more exact approximation.
In Figure 3.2, we show the convergence of the 2-norm of the residual for each of the preconditioning strategies, as well as no preconditioning. We only show the result from the 5th order mean based preconditioner as there was no appreciable difference in convergence. This plot clearly shows that the mean and midpoint preconditioners are excellent choices. Further note that any random point is more than two orders of magnitude better than no preconditioning at all. In fact, using a random point is better than using the point with largest eigenvalue and may compare with using the point with the smallest eigenvalue.
Next, Table 3.1 shows a few timing results from these experiments. We present the time taken to compute/setup the preconditioner, the average time between iterations,  with 24GB of RAM. Matlab used four cores for its own multithreading and we used the Parallel Computing toolbox's parfor construction to parallelize the matrix-vector product over 4 cores as well. By watching a processor performance meter, we observed high utilization of all four cores. We note, however, that during some of the random point evaluations, another process intruded and caused the system to swap memory to disk heavily. We believe this caused the inconsistent average iteration time. From these timing results, we conclude that preconditioning changes the iteration time only slightly -if at all. Furthermore, the setup times are often small when compared with the savings in runtime.
The Matlab codes for this study include a utility for generating realizations of random fields [9] and a suite of tools called PMPack (Parameterized Matrix Pack- age) for using spectral methods to approximate the solution of parameterized matrix equations. The codes for generating the results of the study can be found at: http://www.stanford.edu/~dgleich/publications/2010/codes/pmpack-sisc/ 4. Application -Heat Transfer with Uncertain Material Properties. As a proof of concept, we examine an application from computational fluid dynamics with uncertain model inputs. The flow solver used to compute the deterministic version of this problem -i.e. for a single realization of the model inputs -was developed at Stanford's Center for Turbulence Research as part of the Department of Energy's Predictive Science Academic Alliance Program; the numerical method used is described in [28] and is based on an implicit, second order spatial discretization. For this example, we slightly modified the codes to extract of the non-zero elements of the matrix and right hand side used in the computation of the temperature distribution. With access to the linear system, we were able to apply the Galerkin method to the stochastic version of the problem to approximate the statistics of solution.
4.1. Problem Set-up. The governing equation is the integral version of a twodimensional steady advection-diffusion equation. We seek a scalar field φ = φ(x, y) representing the temperature defined on the domain Ω that satisfies ∂Ω ρφ v(s) · dS = ∂Ω (Γ(s) + Γ t ) ∇φ · dS . (4.1) The density ρ is assumed constant. The velocity v is precomputed by solving the incompressible Reynolds averaged Navier-Stokes equations and randomly perturbed by three spatially varying oscillatory functions with different frequencies such that the divergence free constraint is satisfied; the magnitudes of the perturbations parameterized by s 1 , s 2 , and s 3 , respectively. We interpret the magnitudes as uniform random perturbations of the velocity field, which is simply an input to this model. The diffusion coefficient Γ = Γ(s 4 , s 5 , s 6 ) is similarly altered by a strictly positive, parameterized, spatially varying function that models random perturbation. The turbulent diffusion coefficient Γ t is evaluated according to the Spalart-Allmarass model [28]. The domain Ω is a channel with a series of cylinders; the computational mesh on the domain Ω contains roughly 10,000 nodes and is shown in Figure 4.1. Inflow and outflow boundary conditions are prescribed in the streamwise direction, and periodic boundary conditions are set along the y coordinate. Specified heat flux boundary conditions are applied on the boundaries of the cylinders to model a cooling system. The goal is to compute the expectation and variance of the scalar field φ = φ(x, y, s) over the domain Ω with respect to the variability introduced by the parameters. We use the Galerkin method to construct a polynomial approximation of φ s 1 s 2 s 3 s 4 s 5 s 6 3 1 1 8 5 5  along the coordinates induced by the parameters s using product Legendre polynomials. To show that this method applies for an arbitrary basis set, we choose the basis polynomials to reflect the solution's anisotropic dependence on the parameters [6]; it is neither the standard full polynomial or tensor product polynomial basis. The order of univariate polynomial associated with each parameter is given in Table 4.1. The basis contains 142 multivariate polynomials. To solve the Galerkin system, we use Matlab's BiCGstab [34] method (since the matrix is not symmetric) preconditioned by the diagonal elements of the parameterized matrix evaluated at the midpoint of the hypercube parameter space.
We plot the expectation and variance of φ over the domain Ω in Figure 4.2, which are computed in the standard way as functions of the Galerkin coefficients. As expected, the variance in φ occurs in the downstream portion of the domain as a result of the variability in the diffusion coefficient.

Summary.
We have examined the system of equations arising from a spectral Galerkin approximation of the vector valued solution x(s) to the parameterized matrix equation A(s)x(s) = b(s). Such problems appear when PDE models with parameterized (or random) inputs are discretized in space, and a Galerkin projection with an orthogonal polynomial basis is used for approximation in the parameter space. We showed how the system of equations used to compute the coefficients of the Galerkin approximation admits a factorization once the integration is formally replaced by a numerical quadrature rule -a common step in practice. The factorization involves (i) a matrix with orthogonal rows related to the chosen polynomial basis and numerical quadrature rule and (ii) a block diagonal matrix with nonzero blocks equal to A(s) evaluated at the quadrature points. Then matrix-vector products with the Galerkin system can be computed with only the action of A(s) on a vector at a point in the parameter space; this yields a reusable interface for implementing the Galerkin method. The factorization also reveals bounds on the spectrum of the Galerkin matrix, and the Kronecker structure of the factorization gives clues to successful preconditioners. We tested some ideas from these preconditioner clues on the standard test problem of an elliptic PDE with parameterized coefficients, and we saw that using A(s) evaluated at the midpoint of the parameter space was as effective as the popular mean-based preconditioner; the midpoint preconditioner is much easier to compute, in general. As a proof of concept, we applied the method to an engineering flow problem by slightly modifying an existing CFD code to retreive the matrix elements.
6. Acknowledgments. This material is based upon work supported by the Department of Energy [National Nuclear Security Administration] under Award Number NA28614.