We develop a least-squares method for the velocity-flux first-order form of the Navier-Stokes equations. This form involves the velocity gradient tensor as a new dependent variable, called velocity flux. The method is based on minimization of a functional involving a negative norm for the residual of the momentum equation. To develop a practical method this norm is replaced by a discrete equivalent, as suggested in a recent paper by Bramble and Pasciak. Among the main results are optimal discretization error estimates for conforming finite element approximations, and well-posedness of the discrete algebraic problem.