Abstract. 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.Key words. parameterized systems, spectral methods, stochastic Galerkin AMS subject classifications. 65C99, 65F10 DOI. 10.1137/1007990461. 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 a 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 ∈ S be a set of input parameters from a d-dimensional tensor product parameter space S = S 1 ⊗ · · · ⊗ S d ; the range of each S i may be bounded or unbounded. We equip the parameter space with a bounded, separable, positive weight function ω : S → R + , where