For a given smooth multivariate function \(\phi = \phi(Y^1,\dots,Y^d)\) depending on \(d\) random variables \(Y^1,\dots,Y^d,\) the authors analyze the problem of approximating \(\phi\) by discrete least-squares projection on a multivariate polynomial space, starting from noise-free observation of \(\phi\) on random evaluations of \(Y^1,\dots,Y^d.\) In Section 2 of the paper is introduced the approximation problem as an \(L^2\) projection on a space of polynomials in \(d\) underlying variables; some common choices of polynomial spaces are also described. Further, the optimality of the random \(L^2\) projection, in terms of a best approximation constant, is proved, and the asymptotic behavior of this best approximation constant is analyzed, as the number of random evaluation points goes to infinity. In Section 3 the previous study is restricted to polynomial spaces in one variable and a theorem is proved, which provides a rule to select the number of random points as a function of the maximal polynomial degree, which makes the discrete random \(L^2\) projection nearly optimal with any prescribed confidence level. Section 4 gives the algebraic formulation of the random projection problem, and Section 5 complements the analysis with numerical tests, both in the one-dimensional case and in higher dimensions, respectively.