We introduced a random symmetric Gauss-Seidel (RSGS) method, which was designed to handle large scale linear least squares problems involving tall coefficient matrices. This RSGS method projected the approximate residual onto the subspace spanned by two symmetric columns at each iteration. These columns were sampled from the coefficient matrix based on an effective probability criterion. Our theoretical analysis indicated that RSGS converged when the coefficient matrix had full column rank. Furthermore, numerical experiments demonstrated that RSGS outperformed the baseline algorithms in terms of iteration steps and CPU time.