This paper develops a multilevel least-squares approach for the numerical solution of the complex scalar exterior Helmholtz equation. This second-order equation is first recast into an equivalent first-order system by introducing several "field" variables. A combination of scaled L 2 and H −1 norms is then applied to the residual of this system to create a least-squares functional. It is shown that, in an appropriate Hilbert space, the homogeneous part of this functional is equivalent to a squared graph norm, that is, a product norm on the space of individual variables. This equivalence to a norm that decouples the variables means that standard finite element discretization techniques and standard multigrid solvers can be applied to obtain optimal performance. However, this equivalence is not uniform in the wavenumber k, which can signal degrading performance of the numerical solution process as k increases. To counter this difficulty, we obtain a result that characterizes the error components causing performance degradation. We do this by defining a finite-dimensional subspace of these components on whose orthogonal complement k-uniform equivalence is proved for this functional and an analogous functional that is based only on L 2 norms. This subspace equivalence motivates a nonstandard multigrid method that attempts to achieve optimal convergence uniformly in k. We report on numerical experiments that empirically confirm k-uniform optimal performance of this multigrid solver. We also report on tests of the error in our discretization that seem to confirm optimal accuracy that is free of the so-called pollution effect.