This study proposes a unified stochastic framework for approximating and computing the gradient of every smooth function evaluated at non-independent variables, using $\ell_p$-spherical distributions on $\R^d$ with $d, p\geq 1$. The upper-bounds of the bias of the gradient surrogates do not suffer from the curse of dimensionality for any $p\geq 1$. Also, the mean squared errors (MSEs) of the gradient estimators are bounded by $K_0 N^{-1} d$ for any $p \in [1, 2]$, and by $K_1 N^{-1} d^{2/p}$ when $2 \leq p \ll d$ with $N$ the sample size and $K_0, K_1$ some constants. Taking $\max\left\{2, \log(d) \right\} < p \ll d$ allows for achieving dimension-free upper-bounds of MSEs. In the case where $d\ll p< +\infty$, the upper-bound $K_2 N^{-1} d^{2-2/p}/ (d+2)^2$ is reached with $K_2$ a constant. Such results lead to dimension-free MSEs of the proposed estimators, which boil down to estimators of the traditional gradient when the variables are independent. Numerical comparisons show the efficiency of the proposed approach.