Identify non-interacting variables by querying a black-box function I have a black-box function of many variables, $f(x_1, ...,x_d)$, so that I can do with $f$ is to evaluate the function at chosen points. I would like to find the maximal set of subsets of the variables that do not have interaction so that I can express the black-box function as $f = g_1(z_1) + ... + g_m(z_m)$, where, 1) each $z_j = (x_{j_1}, ...,x_{j_t})$ corresponds to a subset of the variables $\{x_1, ...,x_d\}$,  2) $z_i \cap z_j = \emptyset$ for $i \not=j$, and 3) members ${x_{j_1}, ...,x_{j_t}}$ of $z_j$ are all interacting so that $g_j(z_j)$ cannot be expanded further.
I don't really care about recovering each individual functions $g_j$, but only to identify all the non-interacting groups of variables $\{z_j\}$.
Question 1:
If $f(x_1, ...,x_d)$ can be expanded like described above, then is there a procedure to asymptotically correctly identify, $\{z_j\}$, by querying the black box function $f$? 
Question 2:
How efficient can such algorithm be in terms of sample complexity (the number of queries required to find the correct $\{z_j\}$ on "average", assuming we have some probability measure of space of $f$). 
Question 3:
If Question 1 and/or 2 is too vague to answer, what is the weakest condition that can be placed on $f$ and/or the variables ${x_1, ..., x_d}$, so that these questions can be answered. 
I'm not even sure which branches of math are involved here, so any suggestion about the tags of this question would be helpful too. 
 A: What you're trying to do is an awful lot like Prinicipal Components Analysis (PCA). As a first go at this you should try the following using PCA:


*

*Get a Monte Carlo simulation program (or write a routine in your favorite numerical computation package) and use it to randomly sample points from the domain of your function. If your function has an unbounded domain, then sample each according to a standard normal. This is OK for your purposes, because you aren't trying to create a probability model, but merely trying to query the internal structure of $f(x)$.

*Generate a large sample (>1000). Try for at least 100 sample points per dimension to make sure you have enough points. 

*Examine cross sections of your data by looking at various two-variable combinations (you can crate a scatter plot matrix, which looks like this). This should give you a lot of info about the behavior of the function, especially since its deterministic.

*If your function is stongly non-linear in some (all) dimensions, you will either need to create several bi-variate models and/or transform varaibles so that you can get a an approximately linear function in the transformed or modeled explanatory varables.

*Finally, run your data through a Principal Components Analysis. This analysis will determine "orthogoal" vectors in your (possibly transformed) domain, composed of linear combinations of the actual variables. To use your notation: $z_i = c_{1i}x_1+c_{2i}x_2....$ So that each $z_i$ is indpendent of the other z's. You can also opt to apply a Varimax rotation to the results, which will tell the package to maximize the differences in how much each "z" explains. This will effectievly increase the correlation of some z's with f, while decreasing others, so that you have an even smaller set of z's to deal with.

*Once you get the function approximately linear and you have your equations for each z in terms of a linear combination of x's, you can perform a multiple linear regression on the z's on the transformed f(x) to get the final function and then back-transform to get the original functional realtionship.
Of course, there is no guarantee that every possible $f(x)$ can actually be represented, even partially, as a linear combination of independent functions as you hope, unless you are happy with a vary small region of validity.
There may be other ways to do this, but this is the most "algorithmic" approach I could think of...it still requires a bit from you though.
