This paper is aimed at the efficient numerical simulation of optimization problems governed by either steady-state or unsteady partial differential equations involving random coefficients. This class of problems often leads to prohibitively high dimensional saddle point systems with tensor product structure, especially when discretized with the stochastic Galerkin finite element method. Here, we derive and analyze robust Schur complement-based block-diagonal preconditioners for solving the resulting stochastic optimality systems with all-at-once low-rank solvers. Moreover, we illustrate the effectiveness of our solvers with numerical experiments.