We consider the numerical solution of the large systems of linear equations obtained from the stochastic Galerkin formulation of stochastic partial differential equations. We propose an iterative algorithm that exploits Kronecker product structure of the linear systems. Solutions in such settings often have a special low-rank structure, and Krylov subspace methods have been developed to take advantage of this property by truncating various terms constructed during the course of the iteration. We show here that the costs of truncation can be reduced by identifying an important subspace in the stochastic domain and compresses tensors of high rank on the fly during the iterations. The proposed reduction scheme is a multilevel method in that the important subspace can be identified in a coarse spatial grid setting. In addition, we construct preconditioners for the linear operators in tensor product format that accelerates the convergence of the low-rank iterative algorithms. The efficiency of the proposed method is illustrated by numerical experiments on benchmark problems.