  Computing cross-partial derivatives using fewer model runs is relevant in
modeling, such as stochastic approximation, derivative-based ANOVA, exploring
complex models, and active subspaces. This paper introduces surrogates of all
the cross-partial derivatives of functions by evaluating such functions at $N$
randomized points and using a set of $L$ constraints. Randomized points rely on
independent, central, and symmetric variables. The associated estimators, based
on $NL$ model runs, reach the optimal rates of convergence (i.e.,
$\mathcal{O}(N^{-1})$), and the biases of our approximations do not suffer from
the curse of dimensionality for a wide class of functions. Such results are
used for i) computing the main and upper-bounds of sensitivity indices, and ii)
deriving emulators of simulators or surrogates of functions thanks to the
derivative-based ANOVA. Simulations are presented to show the accuracy of our
emulators and estimators of sensitivity indices. The plug-in estimates of
indices using the U-statistics of one sample are numerically much stable.
