  In this work, we study a random orthogonal projection based least squares
estimator for the stable solution of a multivariate nonparametric regression
(MNPR) problem. More precisely, given an integer $d\geq 1$ corresponding to the
dimension of the MNPR problem, a positive integer $N\geq 1$ and a real
parameter $\alpha\geq -\frac{1}{2},$ we show that a fairly large class of
$d-$variate regression functions are well and stably approximated by its random
projection over the orthonormal set of tensor product $d-$variate Jacobi
polynomials with parameters $(\alpha,\alpha).$ The associated uni-variate
Jacobi polynomials have degree at most $N$ and their tensor products are
orthonormal over $\mathcal U=[0,1]^d,$ with respect to the associated
multivariate Jacobi weights. In particular, if we consider $n$ random sampling
points $\mathbf X_i$ following the $d-$variate Beta distribution, with
parameters $(\alpha+1,\alpha+1),$ then we give a relation involving $n, N,
\alpha$ to ensure that the resulting $(N+1)^d\times (N+1)^d$ random projection
matrix is well conditioned. Moreover, we provide squared integrated as well as
$L^2-$risk errors of this estimator. Precise estimates of these errors are
given in the case where the regression function belongs to an isotropic Sobolev
space $H^s(I^d),$ with $s> \frac{d}{2}.$ Also, to handle the general and
practical case of an unknown distribution of the $\mathbf X_i,$ we use
Shepard's scattered interpolation scheme in order to generate fairly precise
approximations of the observed data at $n$ i.i.d. sampling points $\mathbf X_i$
following a $d-$variate Beta distribution. Finally, we illustrate the
performance of our proposed multivariate nonparametric estimator by some
numerical simulations with synthetic as well as real data.
