Based on a novel point of view on 1-dimensional Gaussian quadrature, we present a new approach to the computation of d-dimensional cubature formulae. It is well known that the nodes of 1-dimensional Gaussian quadrature can be computed as eigenvalues of the so-called Jacobi matrix. The d-dimensional analog is that cubature nodes can be obtained from the eigenvalues of certain mutually commuting matrices. These are obtained by extending (adding rows and columns to) certain noncommuting matrices A_1,...,A_d, related to the coordinate operators x_1,...,x_d, in R^d. We prove a correspondence between cubature formulae and "commuting extensions" of A_1,...,A_d, satisfying a compatibility condition which, in appropriate coordinates, constrains certain blocks in the extended matrices to be zero. Thus, the problem of finding cubature formulae can be transformed to the problem of computing (and then simultaneously diagonalizing) commuting extensions. We give a general discussion of existence and of the expected size of commuting extensions and describe our attempts at computing them, as well as examples of cubature formulae obtained using the new approach.