Stochastic Galerkin finite element approximation of PDEs with random inputs leads to linear systems of equations with coefficient matrices that have a characteristic Kronecker product structure. By reformulating the systems as multi-term linear matrix equations, we develop an efficient solution algorithm which generalizes ideas from rational Krylov subspace approximation. Our working assumptions are that the number of random variables
characterizing the random inputs is modest, in the order of a few tens, and that the dependence on these variables is linear, so that it is sufficient to seek only a reduction in the complexity associated with the spatial component of the
approximation space. The new approach determines a low-rank approximation to the solution matrix by performing a projection onto a low-dimensional space and provides an efficient solution strategy whose convergence rate is
independent of the spatial approximation. Moreover, it requires far less memory than the standard preconditioned conjugate gradient method applied to the Kronecker formulation of the linear systems.