We present a randomized primal-dual algorithm that solves the problem $min_{x} max_{y} y^top A x$ to additive error $epsilon$ in time $mathrm{nnz}(A) + sqrt{mathrm{nnz}(A)n}/epsilon$, for matrix $A$ with larger dimension $n$ and $mathrm{nnz}(A)$ nonzero entries. This improves the best known exact gradient methods by a factor of $sqrt{mathrm{nnz}(A)/n}$ and is faster than fully stochastic gradient methods in the accurate and/or sparse regime $epsilon le sqrt{n/mathrm{nnz}(A)}$. Our results hold for $x,y$ in the simplex (matrix games, linear programming) and for $x$ in an $ell_2$ ball and $y$ in the simplex (perceptron / SVM, minimum enclosing ball). Our algorithm combines Nemirovskis conceptual prox-method and a novel reduced-variance gradient estimator based on sampling from the difference between the current iterate and a reference point.