## Appendix A Surrogate Models
### Svr  
The SVR approach is a generalization of the SVM classification problem by introducing an \(\epsilon\)-sensitive region around the approximated function, also called an \(\epsilon\)-tube. The optimization task in SVR contains two steps: first, finding a convex \(\epsilon\)-insensitive loss function that need to be minimized, and second, finding the smallest \(\epsilon\)-tube that contains the most training instances.  
The convex optimization has a unique solution and is solved using numerical optimization algorithms. One of the main advantages of SVR is that the computational complexity does not depend on the dimensionality of the input space [46]. To deal with otherwise intractable constraints of the optimization problem, we introduce slack variables \(\xi_{i}\) and \(\xi_{i}^{*}\)[47]. The positive constant \(C\) determines the trade-off between the flatness of the function and the magnitude up to which deviations greater than \(\epsilon\) are allowed. The primal quadratic optimization problem of SVR is defined as the following:  
\[\underset{\omega,b}{\text{minimize}} \frac{1}{2}||\omega||^{2}+C\sum_{i=1}^{n}(\xi_{i}+\xi_{i}^{*}) \tag{16}\] \[\text{subject to the following}:\begin{cases}y_{i}-\omega^{T}x_{i}- b&\leq&\epsilon+\xi_{i}\\ \omega^{T}x_{i}+b-y_{i}&\leq&\epsilon+\xi_{i}^{*}\\ \xi_{i},\xi_{i}^{*}&\geq&0\end{cases} \tag{17}\]  
Here, \(\omega\) is the weight and \(b\) the bias to be adjusted. The constrained quadratic optimization problem can be solved by minimizing the Lagrangian with non-negative Lagrange multipliers \(\lambda_{i},\lambda_{i}^{*},\alpha_{i},\alpha_{i}^{*}\), \(i\in\{1,\ldots,n\}\):  
\[\begin{split}\mathcal{L}(\omega,\xi^{*},\xi,\lambda,\lambda^{*}, \alpha,\alpha^{*})=\frac{1}{2}||\omega||^{2}+C\sum_{i=1}^{n}\xi_{i}+\xi_{i}^{* }+\sum_{i=1}^{n}\alpha_{i}^{*}(y_{i}-\omega^{T}x_{i}-\epsilon-\xi_{i}^{*})\\ +\sum_{i=1}^{n}\alpha_{i}(-y_{i}+\omega^{T}x_{i}-\epsilon-\xi_{i} )-\sum_{i=1}^{n}\lambda_{i}\xi_{i}+\lambda_{i}^{*}\xi_{i}^{*}\end{split} \tag{18}\]The minimum of \(\mathcal{L}\) can be found by taking the partial derivatives with respect to the variables and making them equal to zero (_Karush-Kuhn-Tucker_ (KKT) conditions). With the final KKT condition, we can state the following:  
\[\begin{split}\alpha_{i}(-y_{i}+\omega^{T}x_{i}-\varepsilon-\xi_{i} )&=0\\ \alpha_{i}^{*}(-y_{i}+\omega^{T}x_{i}-\varepsilon-\xi_{i}^{*})& =0\\ \lambda_{i}\xi_{i}&=0\\ \lambda_{i}^{*}\xi_{i}^{*}&=0\end{split} \tag{103}\]  
The Lagrange multipliers that are zero correspond to the inside of the \(\varepsilon\)-tube, while the support vectors have non-zero Lagrange multipliers. The function estimate depends only on the support vectors, hence this representation is sparse. More specifically, we can derive the following function approximation to predict \(\hat{y}_{j}(x_{j})\):  
\[\hat{y}_{j}(x_{j})=\sum_{i=1}^{n_{SV}}(\alpha_{i}^{*}-\alpha_{i})x_{i}^{T}x_{j} \tag{104}\]  
with \(\alpha_{i},\alpha_{i}^{*}\in[0,C]\) and the number of support vectors \(n_{SV}\). For nonlinear SVR we replace \(\omega^{T}x_{i}\) in (12)-(15) by \(\omega^{T}\phi(x_{i})\) and the inner product in (16) by the kernel \(K(x_{i},x_{j})\).