Stochastic reduced basis methods for solving large-scale linear random algebraic systems of equations, such as those obtained by discretizing linear stochastic partial differential equations in space, time, and the random dimension, are introduced. The fundamental idea employed is to represent the system response using a linear combination of stochastic basis vectors with undetermined deterministic coefficients (or random functions). We present a theoretical justification for employing basis vectors spanning the preconditioned stochastic Krylov subspace to approximate the response process. Subsequently, variants of the Bubnov–Galerkin scheme are employed to compute the undetermined coefficients, which allow explicit expressions for the response quantities to be derived. We also examine some theoretical properties of the projection scheme and procedures for computing the response statistics. Numerical studies are presented for static and dynamic analysis of stochastic structural systems. We demonstrate that significant improvements over the Neumann expansion scheme, as well as other relevant techniques in the literature, can be achieved.