I am looking for iterative algorithms that can approximate a multivariate convex real-valued function $f(\vec{x})=y,\, \Bbb{R}^n\rightarrow \Bbb{R}$. The function is not known beforehand, but it is possible to compute its value $y_i$ for any $\vec{x}_i$ (and depending on the case, $(\nabla f(\vec{x}))_{\vec{x}_i}$ might also be available).
By probing the function for a large set of $\vec{x}_i$, one can then define hyperplanes that forms upper/lower bounds, based on the $(\vec{x}_i,y_i)$ pairs and the $(\nabla f(\vec{x}))_{\vec{x}_i}$ (if available). The question is: Is there an algorithm to "optimally" select these points? It becomes relevant when computing $f(\vec{x}_i)$ is computationally intensive.
I am aware of the Sandwich Algorithm by Rote, and am looking for something similar for the multivariate case (ideally also with some sort of proof regarding the quality of the approximation). I am also aware of the work of Kamenev (for instance, Chapter 6 of this book), but he assumes that the support function of the convex body can be obtained for any direction, which is not the case in my problem since the function is unknown.
Thanks in advance for your help.