Many applications require multi-dimensional numerical integration, often in the form of a cubature formula. These cubature formulas are desired to be positive and exact for certain finite-dimensional function spaces (and weight functions). Although there are several efficient procedures to construct positive and exact cubature formulas for many standard cases, it remains a challenge to do so in a more general setting. Here, we show how the method of least squares can be used to derive provable positive and exact formulas in a general multi-dimensional setting. Thereby, the procedure only makes use of basic linear algebra operations, such as solving a least squares problem. In particular, it is proved that the resulting least squares cubature formulas are ensured to be positive and exact if a sufficiently large number of equidistributed data points is used. We also discuss the application of provable positive and exact least squares cubature formulas to construct nested stable high-order rules and positive interpolatory formulas. Finally, our findings shed new light on some existing methods for multivariate numerical integration and under which restrictions these are ensured to be successful.