Non-parametric estimation of functions as well as their derivatives by means of local-polynomial regression is a subject that was studied in the literature since the late 1970's. Given a set of noisy samples of a $\mathcal{C}^k$ smooth function, we perform a local polynomial fit, and by taking its $m$-th derivative we obtain an estimate for the $m$-th function derivative. The known optimal rates of convergence for this problem for a $k$-times smooth function $f:\mathbb{R}^d \to \mathbb{R}$ are $n^{-\frac{k-m}{2k + d}}$. However in modern applications it is often the case that we have to estimate a function operating to $\mathbb{R}^D$, for $D \gg d$ extremely large. In this work, we prove that these same rates of convergence are also achievable by local-polynomial regression in case of a high dimensional target, given some assumptions on the noise distribution. This result is an extension to Stone's seminal work from 1980 to the regime of high-dimensional target domain. In addition, we unveil a connection between the failure probability $\varepsilon$ and the number of samples required to achieve the optimal rates.