From a multivariate function $f$, depending on $n\geq 1$ variables, which can be computed numerically, but which does not admit simple analytic expression, I would like to approximate numerically the quantity: $$ \frac{\partial^n f}{\partial x_1...\partial x_n}(x_1,...,x_n)$$ using, e.g., finite differences.
Intuitively, using finite differences, I would proceed like this: Let ($e_1,...,e_n$) be the canonical base of $\mathbb{R}^n$, and let $h\in\mathbb{R}_+^*$ be a small number. If $n = 1$ (univariate function); I would compute: $$ \frac{\partial f}{\partial x_1}(x) \approx \frac{f(x + h e_1) - f(x - h e_1)}{2h}$$ Now, if $n = 2$ (multivariate function with 2 variables), I would compute: $$ \frac{\partial^2 f}{\partial x_1\partial x_2}(x) \approx \frac{(f(x + h e_1 + h e_2) - f(x + h e_1 - h e_2)) - (f(x - h e_1 + h e_2) - f(x - h e_1 - h e_2))}{(2h)^2}$$ and so on for larger $n$. My problem is that this approximation involves $2^n$ terms, which is cumbersome for large $n$. Is anyone aware of a procedure / reference to obtain a good approximation without computing as much as $2^n$ evaluations of $f$ , or is this hopeless ?
It really depends on the structure of your function. For instance, for an archimedean copula, one needs $O(n^2)$ operations, for nested archimedean copulas one can still compute this derivative with an effort of about $O(n^4)$ operations.
Keywords are "cross-derivative" and "anova decomposition", the latter with the aim to approximate the function by a sum of terms in few variables.