An algorithm is demonstrated to yield accurate linear viscoelastic properties in a nonequilibrium stochastic simulation by variance reduction. For noninteracting polymer dumbbells, multiple orders of magnitude improvements in computational accuracy for nonequilibrium viscoelastic properties are easily realized. This simple and robust control variate scheme uses control variables defined as the viscoelastic properties obtained from a parallel, equilibrium simulation. Results for strongly interacting colloidal particles demonstrate some advantages as well as limitations of the method.