Active Learning for Derivative-Based Global Sensitivity Analysis with Gaussian Processes

Thank you for your review of our paper and for your questions, which we answer below.

Weaknesses:

1. It is correct that we do use two performance metrics in our paper, RMSEs and NDCGs. As discussed also in the response to reviewer oJrN above, RMSE is the most effective performance metric for this purpose. RMSE is typically more important than NDCG since it shows the ability of the method to converge to the true values of DGSMs. It is also the most commonly used metric in the literature. As far as we are aware, NDCG has not been used before as a performance metric. However, our own experience has been that in some applications, the ranking of the variables is important for practitioners. For instance, parameter selection for downstream tasks. We thus included a ranking metric to provide another perspective into the results. However, ranking metrics sometimes fail to provide insights when several variables are equally important and have the similar DGSMs values. We discussed this case in the appendix and provided results in figure 6. Hence, RMSE generally provides more important insights into the efficacy of various methods. We will clarify in the paper that RMSE (the typical performance metric in the literature) is the most effective.

2. Estimating the importance of the variables of black box expensive functions is used in several real-world applications (systems safety [2], biomedicine[3], environmental models [4,5], hydrology [6], and more [1]) where the function evaluations are expensive and the functions are high dimensional, so sample efficiency from active learning is crucial. Prior work has discussed the budget efficiency issue repeatedly [e.g., 7,8,9]. Within our paper, we did apply our methods to 3 real-world problems: the Car Side Impact problem, and two Vehicle Safety problems. These problems use simulations of real safety applications. While the simulations are not costly, we would not be able to rigorously evaluate a large suite of methods as we do in the paper on problems with actual costly evaluations.

Questions:

1. The two figures represent the acquisition function values for the information gain over the gradient and the square of the gradient. The figures show that both acquisition functions score most inputs in the search space similarly, reflecting the fact that in this problem, the best points for learning df/dx also learn (df/dx)**2, and vice versa. This will, of course, not always be the case.

2. This is a good question that was also asked by reviewer oJrN; see response there for more details. In short, as the method is exploring new regions, it might add points that are not immediately useful, thus temporarily increasing the RMSE by causing an over- or under-fitting of the model. This can temporarily lead to higher RMSE. However, with more data, the model self-corrects and RMSE decreases.It is also important to note that even when this happens, our proposed approaches still outperform the baselines.

3. To answer your question, we ran experiments to see if our methods that learn the derivatives of f also learn f, or if a separate procedure would be necessary to also estimate f. Conceptually, learning the derivatives of f will also learn f (via integration, by the Fundamental Theorem of Calculus). This is especially true in our setting where we are learning the derivatives from observations of f itself, with a model that is explicitly modeling f alongside its derivatives. However, it is not necessarily the most sample-efficient approach for learning f. The results of the experiments are in Fig. 2 in the rebuttal results, where we evaluate methods on their ability to learn f, by using RMSE of f as the performance metric. Those results show that the derivative information gain (DIG) is actually better at learning f than the f_{VAR} acquisition function that selects points of maximum uncertainty in f. So, the direct answer to your question is no, it is not necessary to use a separate method to accurately estimate f, as estimating the derivatives will in practice also estimate f. The fact that DIG outperforms f_{VAR} can be attributed to DIG being less myopic than f_{VAR}.

[1] Saltelli, Andrea, et al. Global sensitivity analysis: the primer. John Wiley & Sons, 2008.

[2] Qian, Gengjian, et al. "Sensitivity analysis of complex engineering systems: Approaches study and their application to vehicle restraint system crash simulation." Reliability Engineering & System Safety 187 (2019): 110-118.

[3] Wirthl, Barbara, et al. "Global sensitivity analysis based on Gaussian‐process metamodelling for complex biomechanical problems." International journal for numerical methods in biomedical engineering 39.3 (2023): e3675.

[4] Pianosi, Francesca, et al. "Sensitivity analysis of environmental models: A systematic review with practical workflow." Environmental Modelling & Software 79 (2016): 214-232.

Dear Reviewer Uip4,

Thank you for your feedback which will help improve the clarity and contribution of the work. We will provide additional context wrt existing acquisition functions, gold-standard evaluation metrics, and the importance of the problem in evaluating expensive-to-compute simulators in the sciences and engineering,

Please let us know if there is anything else you feel that we should touch upon in the CR if accepted.

Finally, as a gentle reminder, we would appreciate if you could raise your score as stated.

Thank you for your constructive comments and helpful suggestions. Here are responses to the weaknesses and questions:

Weaknesses :

1. While the variance and IG criteria are known in the literature, formulating them in a computationally efficient and tractable fashion is not straightforward for DGSMs, which is our main technical contribution.

2. It is correct that the assumption that the GP derivative approximates the gradient of the true function of interest is vital for the success of the methods we develop. This assumption is motivated both practically and theoretically. On the practical side, GPs are indeed the standard surrogate model for sensitivity analysis of blackbox expensive functions; see e.g. [5] from the main text, and citations therein. On the theoretical side, GPs are favored because they are universal approximators, meaning that, with a suitable kernel, they can indeed approximate any arbitrary continuous target function [1]. Thus, requiring the GP gradients to provide a good approximation of the gradients of f is not, in theory, a limiting assumption.

Questions:

1. Thank you for the suggestion, we implemented and evaluated the suggested idea where we evaluate pairs of points near each other. We generate a quasi-random sequence and then interleave it with points that are small perturbations of the point before it. Specifically, for all n odd, we took x_n from the quasi-random sequence, and then took x_{n+1} = x_n + eps, where eps ~ MVN(0, I*eta) and eta=0.01 (with box bounds scaled to [0, 1]). That is, each dimension is given a small Gaussian perturbation on the value from the previous iteration. Fig. 1 in the rebuttal results shows the performance of this method, which is there called QR-P. In one problem (Ishigami2) it performs very well, though generally the results are comparable to QR without the interleaved pairing. This approach does rely on the perturbation hyperparameter eta, and the “right” choice likely depends on the problem. We will add this ablation study to our paper.

2. In line 135, we use x to refer to any possible value, so it can be m_X if we are evaluating the mean function on X or can be m_{x*} if we are evaluating it on a new input.

3. Thank you for the suggestion. We derived and implemented the information gain for the absolute value of the gradient based on an approximation of the entropy of the folded distribution. We show the results in Fig. 3 in the rebuttal results, on the illustrative problem of Fig. 1 from the main text. Please see the response to reviewer oJrN for full details. In short, the approximation uses a truncated Taylor series with an exp() term that is numerically unstable for values of the posterior mean and variance that we run into in practice. Thus, this acquisition function is not a suitable optimization target, but we will add the result and discussion of it to the paper.

[1] CA Micchelli, Y Xu, H Zhang (2006) “Universal Kernels”, Journal of Machine Learning Research 7:2651-2667.

Thank you for your positive feedback and for recognizing the strengths of our paper. Below we provide answers to your concerns.

Weaknesses:

1. Thank you for the suggestion to fill out the matrix of acquisition functions. For the baseline f_{V_r}, since our experiments used noiseless function evaluations, the observation at x reduces the GP posterior variance at x to 0; thus f_{V_r} is equivalent to f_{VAR} (max variance of f). We will make sure this is clear in the paper. Information gain for the absolute value of the gradient was indeed missing. We derived and implemented information gain for the absolute derivative based on an approximation of the entropy of the folded normal distribution [1]. We evaluated this acquisition on the illustrative test problem in Figure 1 from the main text; the result is in Fig. 3 in the rebuttal results. While generally, the results look similar to other derivative information gain acquisitions in Fig. 1 of the main text, there is a noticeable discrepancy at x=0.6. This is because of numerical issues. The approximation from [1] uses a truncated Taylor expansion that includes an exponential term that can blow up and become a poor approximation depending on the posterior values. To be useful for acquisition optimization, a new, more accurate and stable approximation of the folded normal distribution entropy will be necessary. We will add this result and discussion of future work to the paper.

2. RMSE is typically more important than NDCG since it shows the ability of the method to converge to the true values of DGSMs. It is also the most commonly used metric in the literature. As far as we are aware, NDCG has not been used before as a performance metric. However, our own experience has been that in some applications, the ranking of the variables is important for practitioners. For instance, parameter selection for downstream tasks. We thus included a ranking metric to provide another perspective into the results. However, ranking metrics sometimes fail to provide insights when several variables are equally important and have the similar DGSMs values. We discussed this case in the appendix and provided results in figure 6. Hence, RMSE generally provides more important insights into the efficacy of various methods.

Questions:

1. We used an Automatic Relevance Determination (ARD) RBF kernel.

2. Active learning involves a trade-off between exploring new regions of the input space and exploiting known regions. While the model is still exploring, adding a data point in one part of the space may cause a global adjustment in the model predictions that can cause it to err in another part of the space. This is especially possible early on when the model is still learning the kernel lengthscales, and has the capacity to over- or under-fit the data. With more exploration and data, the model self-corrects and RMSE decreases.

3. Thank you for the suggestion. This certainly seems possible, and in fact methods that ensemble acquisition functions have recently shown great promise for general-purpose black-box optimization tasks [2]. We will draw connections between those results and this work in our discussion of future work.

[1] M. Tsagris, C. Beneki, H. Hassani (2014) "On the folded normal distribution." Mathematics 2:12-18

[2] R Turner et al. "Bayesian optimization is superior to random search for machine learning hyperparameter tuning: Analysis of the black-box optimization challenge 2020." NeurIPS 2020 Competition and Demonstration Track. PMLR, 2021.

Dear reviewer,

Thank you again for your review and time.

A major source of concern was that we were missing two baselines. (1) We have clarified that we have in fact already evaluated $f_{V_r}$ (it is equivalent to $f_\text{VAR}$). (2) We found that IG for the absolute derivative for GPs is not a "baseline" considered by any of the prior literature, and it is in fact a non-trivial acquisition function to derive in our setting. We have identified an approximation of this quantity based on a Taylor expansion of the entropy of the folded normal, ran experiments, and found that the acquisition function contained discontinuities, rendering it an ill-defined AF for this task. Thus there is no straightforward "baseline" for this in the literature.

We have also clarified a number of minor points, such as the fact that NDCG is not a common evaluation metric for sensitivity analysis (which is why we did not include it in the main text), and your query about why active learning may decrease RMSE in some situations.

Thank you for your positive feedback and for supporting our paper. Below we provide answers to your concerns.

Questions

1. We agree that adding a general takeaway would be beneficial for practitioners. We will add the following paragraph to the discussion: “Our general recommendation is to use the information gain of the gradient (DIG) or the information gain of the squared gradient (DSqIG) for low dimensional problems (up to d=10). For high dimensional problems (d>10), we recommend using the variance of the absolute value of the gradient (DAbV) or the variance of the squared gradient (DSqV)”. This will complement the existing discussion in Sec. 5.1 on why the exploratory nature of variance-based approaches are important in high dimensions.

2. See the table below for the acquisition optimization time for each method across different problems. Generally information gain methods are more expensive, with the squared derivative information gain particularly expensive due to the use of a hypergeometric function in its approximation. However, even the maximum time for that method (214 seconds) is fast relative to the time-consuming function evaluation setting that is the focus of this paper. It is interesting to also note that running time generally increases with dimension but not always. This is because the optimization time will depend on the shape of the acquisition function surface and how hard the optimization is, which is orthogonal to the dimensionality of the problem.