In numerical magnetohydrodynamics (MHD), a major challenge is maintaining del . B = O. Constrained transport (CT) schemes achieve this but have been restricted to specific methods. For more general (meshless, moving-mesh, ALE) methods, 'divergence-cleaning' schemes reduce the del . B errors; however they can still be significant and can lead to systematic errors which converge away slowly. We propose a new constrained gradient (CG) scheme which augments these with a projection step, and can be applied to any numerical scheme with a reconstruction. This iteratively approximates the least-squares minimizing, globally divergence-free reconstruction of the fluid. Unlike 'locally divergence free' methods, this actually minimizes the numerically unstable del . B terms, without affecting the convergence order of the method. We implement this in the mesh-free code GIZMO and compare various test problems. Compared to cleaning schemes, our CG method reduces the maximum del . B errors by similar to 1-3 orders of magnitude (similar to 2-5 dex below typical errors if no del . B cleaning is used). By preventing large del . B at discontinuities, this eliminates systematic errors at jumps. Our CG results are comparable to CT methods; for practical purposes, the del . B errors are eliminated. The cost is modest, similar to 30 per cent of the hydro algorithm, and the CG correction can be implemented in a range of numerical MHD methods. While for many problems, we find Dedner-type cleaning schemes are sufficient for good results, we identify a range of problems where using only Powell or '8-wave' cleaning can produce order-of-magnitude errors.