  We propose and analyze a weighted greedy scheme for computing deterministic
sample configurations in multidimensional space for performing least-squares
polynomial approximations on $L^2$ spaces weighted by a probability density
function. Our procedure is a particular weighted version of the approximate
Fekete points method, with the weight function chosen as the (inverse)
Christoffel function. Our procedure has theoretical advantages: when linear
systems with optimal condition number exist, the procedure finds them. In the
one-dimensional setting with any density function, our greedy procedure almost
always generates optimally-conditioned linear systems. Our method also has
practical advantages: our procedure is impartial to compactness of the domain
of approximation, and uses only pivoted linear algebraic routines. We show
through numerous examples that our sampling design outperforms competing
randomized and deterministic designs when the domain is both low and high
dimensional.
