We introduce a randomized algorithm for overdetermined linear least-squares regression. Given an arbitrary full-rank m × n matrix A with m ≥ n , any m × 1 vector b , and any positive real number ε, the procedure computes an n × 1 vector x such that x minimizes the Euclidean norm ‖ Ax − b ‖ to relative precision ε. The algorithm typically requires 𝒪((log( n )+log(1/ε)) mn + n 3 ) floating-point operations. This cost is less than the 𝒪( mn 2 ) required by the classical schemes based on QR -decompositions or bidiagonalization. We present several numerical examples illustrating the performance of the algorithm.