Exploring and interpolating a function using machine-learning?
I would look into the field of "optimal experimental design" in bayesian inverse problems, particularly the recent work of Alen Alexandrian.

http://arxiv.org/abs/1410.5899

http://www4.ncsu.edu/~aalexan3/research.html

Essentially, one has an inner inverse problem for approximating the function based on point measurements of derived quantities, hosted within an outer optimization problem for choosing the points based on minimizing a combination of the error and the variance.

Furthermore, you don't need to do a full inner-outer solve procedure. Rather, you can use the KKT conditions for the inner problem as the constraint for the outer problem, and formulate a "meta" KKT system for the combined problem.

It is formulated in the language of PDE-constrained inverse problems, but would also apply to simpler situations like your problem (the "PDE" becomes the identity matrix..)