# Improving Accuracy of Interpretability Measures in Hyperparameter
Optimization via Bayesian Algorithm Execution

Julia Moosbauer<EMAIL_ADDRESS>
Institute of Statistics, Munich Center for Machine Learning (MCML)
Ludwig-Maximilians-Universität München Giuseppe Casalicchio
<EMAIL_ADDRESS>
Institute of Statistics, Munich Center for Machine Learning (MCML)
Ludwig-Maximilians-Universität München Marius Lindauer<EMAIL_ADDRESS>hannover.de
Institute of Artificial Intelligence
Leibniz University Hannover Bernd Bischl<EMAIL_ADDRESS>
Institute of Statistics, Munich Center for Machine Learning (MCML)
Ludwig-Maximilians-Universität München

###### Abstract

Despite all the benefits of automated hyperparameter optimization (HPO), most
modern HPO algorithms are black-boxes themselves. This makes it difficult to
understand the decision process which leads to the selected configuration,
reduces trust in HPO, and thus hinders its broad adoption. Here, we study the
combination of HPO with interpretable machine learning (IML) methods such as
partial dependence plots. These techniques are more and more used to explain
the marginal effect of hyperparameters on the black-box cost function or to
quantify the importance of hyperparameters. However, if such methods are
naively applied to the experimental data of the HPO process in a post-hoc
manner, the underlying sampling bias of the optimizer can distort
interpretations. We propose a modified HPO method which efficiently balances
the search for the global optimum w.r.t. predictive performance _and_ the
reliable estimation of IML explanations of an underlying black-box function by
coupling Bayesian optimization and Bayesian Algorithm Execution. On benchmark
cases of both synthetic objectives and HPO of a neural network, we demonstrate
that our method returns more reliable explanations of the underlying black-box
without a loss of optimization performance.

##  1 Introduction

The performance of machine learning (ML) models usually depends on many
decisions, such as the choice of a learning algorithm and its hyperparameter
configurations. Manually reaching these decisions is usually a tedious trial-
and-error process. Automated machine learning (AutoML), e.g., hyperparameter
optimization (HPO), can support developers and researchers in this regard. By
framing these decisions as an optimization problem and solving them using
efficient black-box optimizers such as Bayesian Optimization (BO), HPO is
demonstrably more efficient than manual tuning, and grid or random search
(Bergstra et al., 2011; Snoek et al., 2012; Turner et al., 2020; Bischl et
al., 2021). However, there is still a lack of confidence in AutoML systems and
a reluctance to trust the returned best configuration (Drozdal et al., 2020).
One reason for why some practitioners still today prefer manual tuning over
automated HPO is that existing systems lack the ability to convey an
understanding of hyperparameter importance and how certain hyperparameters
affect model performance (Hasebrook et al., 2022), helping them to understand
why a final configuration was chosen.

Desirable insights into hyperparameter effects or importance could in
principle be generated by applying methods of interpretable machine learning
(IML) to experimental data from the HPO process, specifically the final
surrogate model generated by BO based on this HPO-data. However, these methods
– even though possible from a technical perspective and used before (Hutter et
al., 2014; Van Rijn & Hutter, 2018; Young et al., 2018; Head et al., 2022) –
should be used with caution in this context. The main reason is a sampling
bias caused by the desire for efficient optimization during HPO (Moosbauer et
al., 2021): Efficient optimizers typically sample more configurations in
promising regions with potentially well-performing hyperparameter
configurations, while other regions are underrepresented. This sampling bias
introduces a surrogate-model bias in under-explored regions as the surrogate
model is subject to high uncertainty in these regions. Consequently,
explanations of HPO runs, such as partial dependence plots (PDPs) (Friedman,
2001), can be misleading as they also rely on artificially created evaluations
in under-explored regions. Moosbauer et al. (2021) had been the first to
address this issue, and had proposed an approach to identify well-explored,
rather small subregions in which PDPs can be estimated accurately. While this
is valuable, it still does not allow to accurately estimate hyperparameter
effects globally.

To anticipate these unintended effects of this sampling bias as effectively as
possible already during the HPO process, we propose a modified BO algorithm
that efficiently searches for the global optimum _and_ accurate IML estimates
of the underlying black-box function at the same time. We build on the concept
of Bayesian Algorithm Execution (BAX) (Neiswanger et al., 2021) to estimate
the expected information gain (EIG) (Lindley, 1956) of configurations w.r.t.
the output of an interpretation method. We ultimately couple BO with BAX and
propose BOBAX as an efficient method that searches for accurate
interpretations without a relevant loss of optimization performance. Our
proposed method is generic as it is applicable to any BO variant (e.g.,
different acquisition functions or probabilistic surrogate models). As IML
technique we focus on PDPs (Friedman, 2001), which estimate the marginal
effect(s) of features (in our case: hyperparameters) on the output by
visualizing a marginal 1D or 2D function. PDPs constitute an established IML
technique (Lemmens & Croux, 2006; Cutler et al., 2007; Wenger & Olden, 2012;
Zhang et al., 2018), have been in use for more than 20 years to analyze ML
models, and have recently gained further interest in IML and XAI, and are also
increasingly used to analyze hyperparameter effects in HPO and AutoML (Young
et al., 2018; Zela et al., 2018; Head et al., 2022). We point out that our
technique is in principle not limited to PDPs, but can be combined with any
IML technique which can be quantitatively estimated from a surrogate model.

In a benchmark study, we demonstrate how BOBAX consistently yields more
reliable estimates for marginal effects estimated via the partial dependence
method while maintaining the same level of optimization efficiency as commonly
used methods. Finally, we demonstrate how BOBAX can give reliable insights
into hyperparameter effects of a neural network during tuning yielding state-
of-the-art performance. We believe that through our generic method, the
potential of IML methods can be unlocked in the context of HPO, thus paving
the way for more interpretability of and trust into human-centered HPO. Our
contributions include:

  1. 1.

The direct optimization for an accurate estimation of IML statistics, e.g.,
marginal effects for single or multiple hyperparameters, as part of BO for
HPO, making HPO interpretable and more trustworthy;

  2. 2.

The combination of BO and Bayesian Algorithm Execution (BAX), dubbed BOBAX,
where BAX is used to guide the search towards more accurate estimation of IML
statistics;

  3. 3.

Thorough study of different variants of BOBAX and baselines on synthetic
functions; and

  4. 4.

Empirical evidence that budget allocation regarding IML estimates does not
come at the expense of significantly reduced optimization performance on a
deep learning HPO benchmark.

##  2 Background

In this section, we formalize HPO and BO as the context of our work. We also
give an overview of Bayesian Algorithm Execution (BAX) as it serves as basis
for our work.

##### Hyperparameter Optimization

The aim of HPO is to efficiently find a well-performing configuration of a
learning algorithm. HPO is therefore commonly formalized as finding the
minimizer $\bm{\lambda}^{\ast}\in\mathop{\sf
arg\,min}\nolimits_{\bm{\lambda}\in\Lambda}c(\bm{\lambda})$ of a _black-box_
cost function $c:\Lambda\to\mathds{R}$ which maps a hyperparameter
configuration
$\mbox{$\bm{\lambda}=\left(\lambda_{1},...,\lambda_{d}\right)$}\in\Lambda$ to
the validation error of the model trained by a learning algorithm run using
$\bm{\lambda}$. The hyperparameter space
$\Lambda=\Lambda_{1}\times...\times\Lambda_{d}$ can be mixed, containing
categorical and continuous hyperparameters. Particularly in the context of
AutoML, where whole machine learning pipeline configurations are optimized
over, $\Lambda$ may even contain hierarchical dependencies between
hyperparameters (Thornton et al., 2013; Olson & Moore, 2016).

##### Bayesian Optimization

BO is a black-box optimization algorithm which has become increasingly popular
in the context of HPO (Jones et al., 1998; Snoek et al., 2012). BO
sequentially chooses configurations
$\bm{\lambda}^{(1)},...,\bm{\lambda}^{(T)}$ that are evaluated
$c_{\bm{\lambda}^{(1)}},...,c_{\bm{\lambda}^{(T)}}$ to obtain an archive
$A_{T}=\left\\{\left(\bm{\lambda}^{(i)},c_{\bm{\lambda}^{(i)}}\right)\right\\}_{i=1,...,T}$.
To choose the next configuration $\bm{\lambda}^{(T+1)}$ as efficiently as
possible, a surrogate model $\hat{c}$ is estimated on the archive $A_{T}$, and
a new point is proposed based on an acquisition function that leverages
information from the surrogate model $\hat{c}$. Typically, we chose a
probabilistic model and estimate a distribution over $c$, denoted by
$p(c~{}|~{}A_{T})$. A common choice are Gaussian processes
$c\sim\mathcal{GP}\left(\mu,k\right)$, characterized by a mean function
$\mu:\Lambda\to\mathds{R}$ and a covariance function
$k:\Lambda\times\Lambda\to\mathds{R}$. Acquisition functions usually trade off
exploration (i.e., sampling in regions with few data points and high posterior
uncertainty) and exploitation (i.e., sampling in regions with low mean).
Common examples are the expected improvement (EI) (Jones et al., 1998), the
lower confidence bound (LCB) (Jones, 2001; Srinivas et al., 2010), entropy
search (Hennig & Schuler, 2012; Hernández-Lobato et al., 2014) and knowledge
gradient (Wu et al., 2017).

##### Marginal Effects of Hyperparameters

Practitioners of HPO are often interested in whether and how individual
hyperparameters affect model performance. Not only is there a desire to gain
model comprehension Hasebrook et al. (2022), also such insights can influence
decisions, for example whether to tune a hyperparameter or not (Probst et al.,
2019), or modify hyperparameter ranges. One interpretation measure that the
community is looking at (Hutter et al., 2014; Zela et al., 2018; Young et al.,
2018; Van Rijn & Hutter, 2018; Zöller et al., 2022) is the marginal effect of
one or multiple hyperparameters $\bm{\lambda}_{S}$, $S\subset\\{1,2,...,d\\}$
on model performance, which is defined as111To keep notation simple, we denote
$c(\bm{\lambda})$ as a function of two arguments
$(\bm{\lambda}_{S},\bm{\lambda}_{R})$ to differentiate components in the index
set $S$ from those in the complement $R=\\{1,2,...,d\\}\setminus S$. The
integral shall be understood as a multiple integral of $c$ where
$\bm{\lambda}_{j}$, $j\in R$, are integrated out.

$\displaystyle
c_{S}(\bm{\lambda}_{S}):=\mathds{E}_{\bm{\lambda}_{R}}\left[c(\bm{\lambda})\right]=\int_{\Lambda_{R}}c(\bm{\lambda}_{S},\bm{\lambda}_{R})~{}\textrm{d}\mathbb{P}(\bm{\lambda}_{R}).$
(1)

In the context of HPO, $\mathds{P}$ is typically assumed to be the uniform
distribution over $\Lambda_{R}$ since we are interested in how hyperparameter
values $\bm{\lambda}_{S}$ impact model performance uniformly across the
hyperparameter space (Hutter et al., 2014; Moosbauer et al., 2021). Since
computing Eq. (1) analytically is usually possible, the PDP method (Friedman,
2001) approximates the integral, as in Eq. (1), by Monte Carlo approximation.

##### Information-based Bayesian Algorithm Execution

Information-based Bayesian Algorithm Execution (BAX) extends the idea of using
entropy search for estimating global optima to estimating other properties of
a function $f:\mathcal{X}\to\mathds{R}$ (Neiswanger et al., 2021). Similar to
BO, BAX tries to sequentially choose points $\mathbf{x}^{(i)}\in\mathcal{X}$
in order to estimate the quantity of interest accurately with as few
evaluations as possible. It is assumed that the quantity of interest can be
computed as the output
$\mathcal{O}_{\mathcal{A}}:=\mathcal{O}_{\mathcal{A}}(f)$ of running an
algorithm $\mathcal{A}$ on $f$, e.g. top-k estimation on a finite set,
computing level sets or finding shortest paths.

Similarly to BO, BAX sequentially builds a probabilistic model
$p(f~{}|~{}A_{T})$, e.g., a GP, over an archive of evaluated points $A_{T}$.
Based on $p(f~{}|~{}A_{T})$, they derive the posterior distribution over the
algorithm output $p(\mathcal{O}_{\mathcal{A}}~{}|~{}A_{T})$. To build the
archive $A_{T}$ as efficiently as possible, they choose to evaluate the point
$\mathbf{x}^{(T+1)}$ which maximizes the expected information gain about the
algorithm output $\mathcal{O}_{\mathcal{A}}$

$\displaystyle\text{EIG}_{T}(\mathbf{x})$
$\displaystyle:=\mathbb{H}\left[\mathcal{O}_{\mathcal{A}}|A_{T}\right]-\mathbb{E}_{f_{\mathbf{x}}|A_{T}}\left[\mathbb{H}\left[\mathcal{O}_{\mathcal{A}}|A_{T+1}\right]\right],$
(2)

where $\mathbb{H}$ denotes the entropy, and
$A_{T+1}:=A_{T}\cup\left\\{\left(\mathbf{x},f_{\mathbf{x}}\right)\right\\}$
with $f_{\mathbf{x}}$ the (unrevealed) value of $f$ at $\mathbf{x}$.

Neiswanger et al. (2021) propose an acquisition function to approximate the
EIG in Eq. (2). In its simplest form, the algorithm output
$\mathcal{O}_{\mathcal{A}}$ in the EIG is replaced by the algorithm’s
execution path $e_{\mathcal{A}}$, i.e., the sequence of all evaluations the
algorithm $\mathcal{A}$ traverses, which thus gives full information about the
output. The expected information gain estimated based on the execution path
$e_{\mathcal{A}}$ is given by

$\displaystyle\text{EIG}_{T}^{e}(\mathbf{x})$
$\displaystyle=\mathbb{H}\left[e_{\mathcal{A}}|A_{T}\right]-\mathbb{E}_{f_{\mathbf{x}}|A_{T}}\left[\mathbb{H}\left[e_{\mathcal{A}}|A_{T+1}\right]\right]$
(3)
$\displaystyle=\mathbb{H}\left[f_{\mathbf{x}}|A_{T}\right]-\mathbb{E}_{e_{\mathcal{A}}|A_{T}}\left[\mathbb{H}\left[f_{\mathbf{x}}|A_{T},e_{\mathcal{A}}\right]\right].$

where they used the symmetry of the mutual information to come up with the
latter expression. The first term
$\mathbb{H}\left[f_{\mathbf{x}}|A_{T}\right]$ is the entropy of the posterior
predictive distribution at an input $\mathbf{x}$ and can be computed in closed
form. The second term can be estimated as follows: A number of
$n_{\text{path}}$ samples $\tilde{f}\sim p(f~{}|~{}A_{T})$ is drawn from the
posterior process. The algorithm $\mathcal{A}$ is run on each of the samples
$\tilde{f}$ to produce sample execution paths $\tilde{e}_{\mathcal{A}}$,
yielding samples $\tilde{e}_{\mathcal{A}}\sim p(e_{\mathcal{A}}~{}|~{}A_{T})$,
used to estimate the second term as described by Neiswanger et al. (2021).

##  3 Related Work

Interpretability in AutoML refers either to (1) the interpretation of the
resulting model returned by an AutoML system (Xanthopoulos et al., 2020;
Binder et al., 2020; Carmichael et al., 2021; Coors et al., 2021), or (2) the
interpretation of hyperparameter effects and importance (Moosbauer et al.,
2021). We focus on the latter, specifically the construction of accurate and
unbiased estimators for, e.g., hyperparameter effects in HPO.

There are HPO and AutoML frameworks that provide visualisations and
interpretability statistics as additional outputs, e.g., _Google Vizier_
(Golovin et al., 2017) and _xAutoML_ (Zöller et al., 2022) provide an
interactive dashboard visualizing the progress of the optimization and
insights via parallel coordinate plots and multi-dimensional scaling on the
optimizer footprint. Similarly, the HPO frameworks _optuna_ (Akiba et al.,
2019) or _scikit-optimize_ (Head et al., 2022) allow for quick and simple
visualization of optimization progress and results. However, such relatively
simple visualizations do not give a deeper understanding of which
hyperparameter influence model performance in what way.

In the context of HPO, practitioners are commonly interested on the marginal
effects of hyperparameters on model performance Hutter et al. (2014); Young et
al. (2018); Zela et al. (2018) or the importance of hyperparameters on model
performance (Hutter et al., 2014; Biedenkapp et al., 2017; Van Rijn & Hutter,
2018; Probst et al., 2019). The latter is often directly derived from marginal
effects of hyperparameters (Hutter et al., 2014). Established HPO frameworks
(Head et al., 2022; Akiba et al., 2019) as well as visualization toolboxes
Zöller et al. (2022) already make implementations of these methods accessible
to users, however they neither discuss nor address a distortion of those
arising due to a sampling bias. While all of these approaches have their
merits, none of them address the imprecision in the estimates of these
interpretive measures caused by sample bias that is present in the archive
sampled by BO, since BO tends to exploit promising regions while leaving other
regions unexplored. So far, only Moosbauer et al. (2021) explicitly proposed a
post-hoc method that is able to identify subspaces of the configuration space
in which accurate and unbiased PDPs can be computed. However, the method does
not provide more accurate global IML estimates. To our knowledge, we are the
first to propose a method that improves the sampling process of HPO to provide
more accurate global estimates of such IML methods.

##  4 BOBAX: Enhanced Estimation of Interpretability Measures for HPO

We present our main contribution: BOBAX that efficiently searches for accurate
marginal effect estimates of hyperparameters while maintaining competitive HPO
performance.

###  4.1 Expected Information Gain for Partial Dependence

We first derive the information gained with regards to the estimate of a
marginal effect of a hyperparameter $\bm{\lambda}_{S}$ _if_ we observe
performance $c_{\bm{\lambda}^{(T+1)}}$ for a hyperparameter configuration
$\bm{\lambda}^{(T+1)}$. To this end, we quantify and analyze how a marginal
effect is estimated in the context of HPO. Two types of approximations are
performed: First, instead of estimating the marginal effect with regards to
the true, but unknown and expensive objective $c$, we estimate the marginal
effect of the surrogate model $\hat{c}$ 222Constructed by BO, usually this
will be the final surrogate model of the BO run, but this can also be applied
interactively to intermediate models, with $\hat{c}$ denoting the posterior
mean of a probabilistic model $p(c~{}|~{}A_{T})$. Secondly, we use the partial
dependence method (Friedman, 2001) for efficient estimation of marginal
effects of $\hat{c}:\Lambda\to\mathds{R}$, which estimates Eq. (1) by Monte-
Carlo sampling:

$\displaystyle\varphi_{\bm{\lambda}_{S}}=\frac{1}{n}\sum_{i=1}^{n}\hat{c}\left(\bm{\lambda}_{S},\bm{\lambda}_{R}^{(i)}\right),$
(4)

with $\bm{\lambda}_{S}$ fixed and
$\bm{\lambda}_{R}^{(i)}\overset{i.i.d.}{\sim}\mathds{P}(\bm{\lambda}_{R})$ a
Monte-Carlo sample drawn from a uniform distribution $\mathds{P}$. To bound
the computational effort to compute the PDP, Eq. (4) is evaluated for a
(typically equidistant) set of grid points
$\\{\bm{\lambda}_{S}^{(j)}\\}_{j=1,...,G}$. The PDP visualizes
$\varphi_{\bm{\lambda}_{S}}$ against $\bm{\lambda}_{S}$.

To define the expected information gain for partial dependence
$\textrm{EIG}_{\textrm{PDP}}$, we have the partial dependence method in terms
of a formal execution path (see also Algorithm 1): We iterate over all grid
points, and compute the mean prediction $\hat{c}^{(g,i)}$. The execution path
$e_{\mathcal{A}}$ thus corresponds to the Cartesian product
$\left(\bm{\lambda}_{S}^{(g)},\bm{\lambda}_{R}^{(i)}\right)$ for
$g\in\\{1,...,G\\}$ and $i\in\\{1,...,n\\}$ of all grid points
$\bm{\lambda}_{S}^{(g)}$ and the Monte-Carlo samples $\bm{\lambda}_{R}^{(i)}$.

As proposed by Neiswanger et al. (2021) as one variant, we estimate the
information gained with regards to the execution path of $e_{\mathcal{A}}$
instead of estimating the execution path with regards to the algorithm output
$O_{\mathcal{A}}$. Note that Neiswanger et al. (2021) argued that the
criterion in Eq. (3) is in general suboptimal, if for example large parts of
the execution path $e_{\mathcal{A}}$ do not have an influence on the algorithm
output. We argue, however, that it is not applicable to our use-case since
every element in the execution path of the PD method contributes with equal
weight to the computation of the partial dependence. Figure 1 illustrates the
computation of the PD based on the execution path, as well as the computation
of the $\textrm{EIG}_{\textrm{PDP}}$.

Input  $G$, $\hat{c}$,
$\left(\bm{\lambda}_{R}^{(i)}\right)\overset{i.i.d.}{\sim}\mathds{P}(\bm{\lambda}_{R})$

$\left(\bm{\lambda}_{S}^{(1)},...,\bm{\lambda}_{S}^{(G)}\right)\leftarrow$
equidist. grid on $\Lambda_{S}$

for $g\in\\{1,2,...,G\\}$ do

for $i\in\\{1,2,...,n\\}$ do

$\bm{\lambda}^{(g,i)}\leftarrow\left(\bm{\lambda}_{S}^{(g)},\bm{\lambda}_{R}^{(i)}\right)$

$\hat{c}^{(g,i)}\leftarrow\hat{c}\left(\bm{\lambda}^{(g,i)}\right)$

$e_{\mathcal{A}}\leftarrow
e_{\mathcal{A}}\cup\left(\bm{\lambda}^{(g,i)},\hat{c}^{(g,i)}\right)$

end for

$\varphi_{\bm{\lambda}_{S}^{(g)}}\leftarrow\frac{1}{n}\sum_{i=1}^{n}\hat{c}^{(g,i)}$

end for

Return
$\left(\bm{\lambda}_{S}^{(g)},\varphi_{\bm{\lambda}_{S}^{(g)}}\right)$,
$g=1,...,G$

Algorithm 1 PD algorithm with explicit execution path $e_{A}$

Figure 1: Shown are the elements of $e_{\mathcal{A}}$ (blue) PD method. The
grey points show the configurations in the archive which are used by BO to
construct the surrogate model. The green configuration is sampled by EI
(showing more exploitation) while the orange point is the point maximizing the
information gained about the PD estimate.

###  4.2 BOBAX: Efficient Optimization and Search for More Accurate
Interpretability Measures

Given the $\textrm{EIG}_{\textrm{PDP}}$ for PD, the optimization for
interpretability of hyperparameter effects as part of BO is possible by using
the $\textrm{EIG}_{\textrm{PDP}}$ as acquisition function. However,
interpretability alone is rarely of primary interest in practice; rather, the
goal is to identify well-performing configurations and obtaining reasonable
interpretations at the same time. We propose a method, dubbed BOBAX, that
allows to efficiently search for explanations without a relevant loss of
optimization efficiency.

BOBAX is an interleaving strategy which performs BO, and iterates between
using the EI (or any other suited acquisition function) and the
$\textrm{EIG}_{\textrm{PDP}}$ as acquisition function. Although we have
investigated also more complex variants (see Appendix B.2), interleaving
$\textrm{EIG}_{\textrm{PDP}}$ in every $k$-th iteration is simple yet
efficient. The smaller $k$ is, the higher is the weight of optimizing for
accurate interpretations in a BO run. We note that this strategy can replace
other interleaving exploration strategies, such as random samples (Hutter et
al., 2011), since optimizing for interpretability can be seen as another
strategy to cover the entire space in an efficient manner.333One might have
also considered addressing this as a multi-objective problem since we have two
objectives: (i) finding the optimum and (ii) obtaining good PDPs. However,
usually post-hoc multi-objective optimizers construct a Pareto front of a set
of multiple candidate solution, which we are not interested in here. Instead,
in each iteration of BO, the optimizer has to choose a concrete trade-off
between both objectives. For dynamically balancing out this trade-off, please
also refer to the next section.

From a practitioner’s point of view, it may be reasonable to consider accuracy
of interpretations rather as a constraint than an objective function to
optimize for. As soon as this constraint is fulfilled, a user may want to
invest all remaining budget into optimization only. Therefore, we also propose
an adaptive variant of BOBAX, dubbed a-BOBAX, which performs the interleaving
strategy in BOBAX as described above in a first phase, and transitions into
optimization only in a second phase as soon as the constraint is fulfilled. To
allow a user to input a meaningful constraint, the constraint must itself be
interpretable by a user. Therefore, we define this constraint by a desired
average width of confidence intervals around PD estimates, using the
definition444Confidence intervals are defined as
$\varphi_{\bm{\lambda}_{S}^{(g)}}\pm
q_{1-\alpha/2}\cdot\hat{s}_{\bm{\lambda}_{S}^{(g)}}$ around the PD estimate.
$\hat{s}_{\bm{\lambda}_{S}^{(g)}}$ denotes the uncertainty of a PD estimate
for a grid point $g$. As default, we look at $\alpha=0.05$.  of Moosbauer et
al. (2021). As an example, a user may want to specify a tolerance $\pm 1\%$ in
validation accuracy in estimation of PDs (see green tolerance bands in Figure
3 for illustration).

###  4.3 Theoretical and Practical Considerations

##### Runtime Complexity

Since BOBAX comes with additional overhead, we discuss this here in more
detail. The computation of the expectation requires posterior samples of the
execution path $e_{\mathcal{A}}~{}\sim p(e_{\mathcal{A}}~{}|~{}A_{T})$. This
is achieved by sampling from the posterior GP $\tilde{c}~{}\sim
p(c~{}|~{}A_{T})$ and execution of $O_{\mathcal{A}}$ on those samples, which
may produce a computational overhead depending on the costs of running
$O_{\mathcal{A}}$. We assume that executing $O_{\mathcal{A}}$ is neglectable
in terms of runtime. However, to compute the entropy
$\mathbb{H}\left[c_{\bm{\lambda}}|A_{T},e_{\mathcal{A}}\right]$, the posterior
process needs to be trained based on $A_{T}\cup e_{\mathcal{A}}$ (which has
size $T+n\cdot G$). Thus, the overall runtime complexity is dominated by
$\mathcal{O}\left(n_{\text{path}}\cdot(T+n)^{3}\right)$, as we compute the
entropy $n_{\text{path}}$ times to approximate the expectation and since
training a GP is cubic in the number of data points. Therefore, we recommend
to keep an eye on the runtime overhead of the calculation of
$\textrm{EIG}_{\textrm{PDP}}$ in relation to evaluating $c$ (e.g., training
and evaluating an ML algorithm). Especially in the context of deep learning,
the evaluation of a single configuration is usually by orders of magnitude
higher than that of computing the $\textrm{EIG}_{\textrm{PDP}}$555In our case,
the computation of the $\textrm{EIG}_{\textrm{PDP}}$ was ranging from the
order of a few seconds to a few minutes.. Also, we would like to emphasize
that the implementation of our method is based on GPflow (Matthews et al.,
2017), which allows fast execution of GPs on GPUs. Since GPUs are typically in
use for training in the context of DL anyway, they can easily be leveraged in
between iterations to speed up the computation of the
$\textrm{EIG}_{\textrm{PDP}}$.

##### Marginal Effects for Multiple Hyperparameters

Until now we have assumed that a user specifies a single hyperparameter of
interest $\bm{\lambda}_{S}$ for which we will compute the PD. However, it is
difficult to prioritize the hyperparameter of interest a-priori. Fortunately,
it is possible to extend the execution path to compute
$\textrm{EIG}_{\textrm{PDP}}$ by the respective execution paths of the PDs
with regards to all variables
$e_{\mathcal{A}}=e_{\mathcal{A},\bm{\lambda}_{1}}\cup
e_{\mathcal{A},\bm{\lambda}_{2}}\cup...\cup e_{\mathcal{A},\bm{\lambda}_{d}}$.
We investigate the differences between $\textrm{EIG}_{\textrm{PDP}}$ for a
single hyperparameter vs. for multiple hyperparameters in more detail in
Appendix C; in the practical use-case (see Section 2), we compute the
$\textrm{EIG}_{\textrm{PDP}}$ for multiple hyperparameters.

##  5 Benchmark

In this section, we present experiments to demonstrate the validity of our
method. In particular, we look at:

##### Hypothesis H1

Performing BO with $\text{EIG}_{\textrm{PDP}}$ as acquisition function is more
efficient than random search in optimizing for accurate interpretations

##### Hypothesis H2

Through BOBAX the accuracy of marginal effect estimates is clearly improved
without a significant loss of optimization performance.

###  5.1 Experimental Setup

Figure 2: The first three plots show the estimated PD with $95\%$ confidence
interval (blue) based on the surrogate model $\hat{c}$ after $T=30$ iterations
vs. the true marginal effect (black). BAX and BOBAX yield more accurate
estimates for the PD as compared the BO with EI. The right plot shows the
cumulative regret for the three methods. BAX, which is not performing
optimization at all, is also clearly outperformed in optimization performance.
BOBAX reaches the optimization result of BO with EI only after a few more
iterations.

##### Objective Functions

We apply our method to synthetic functions which are treated as black-box
function during optimization: Branin ($d=2$), Camelback ($d=2$), Stylinski-
Tang ($d=3$), Hartmann3 ($d=3$) and Hartmann6 ($d=6$).

##### Algorithms

To investigate H1, we consider BO with $\text{EIG}_{\text{PDP}}$ as
acquisition function (BAX). For H2, we consider BOBAX as described in
Algorithm 2, where we iterate evenly ($k=2$) between EI and
$\text{EIG}_{\text{PDP}}$ as acquisition function. Following Neiswanger et al.
(2021) we set the number of execution path samples to $20$ to approximate the
expectation in Eq. (3) in both variants. As strong baseline for accurate PDs
we consider random search (RS) and BO with posterior variance as acquistion
function (PVAR) as a pure exploration case of LCB. As strong baseline for
optimization we consider BO with EI (BO-EI). Further variants of our methods
(e.g., different frequencies of interleaving) and additional baselines (such
as BO with LCB with different exploration factors, or BO with EI and random
interleaving) are described in Appendix C.

##### Evaluation

We evaluate the accuracy of PD estimates by comparing the PD
$\varphi_{S}^{(g)}$ (estimated based on $\hat{c}$) against the PD
$\tilde{\varphi}_{S}^{(g)}$ computed on the ground-truth objective function
$c$, approximated with the same sample $\bm{\lambda}_{R}^{(i)}$ and the same
grid size $G$. As measure we use the $L_{1}$ distance
$\textrm{d}_{\text{L}_{1}}:=\frac{1}{G}\sum_{g=1}^{G}\left|\varphi_{S}^{(g)}-\tilde{\varphi}_{S}^{(g)}\right|$
averaged over all grid points. To assess optimization performance, we report
the simple regret $c(\hat{\bm{\lambda}})-c(\bm{\lambda}^{\ast})$, where
$\bm{\lambda}^{\ast}$ denotes the theoretical optimum of a function, and
$\hat{\bm{\lambda}}\in\textrm{argmin}\left\\{c_{\bm{\lambda}}~{}|~{}(\bm{\lambda},c_{\bm{\lambda}})\in
A_{T}\right\\}$ is the best found configuration during optimization.

##### Further Configurations

A Gaussian process with a squared exponential kernel is used as surrogate
model for all BO variants, and PDs are estimated on the respective surrogate
models. For RS, a GP (with same configuration) is fitted on $A_{T}$ and the PD
is computed thereon. Acquisition function optimization is performed by
randomly sampling $1500$ configurations, evaluating the respective acquisition
function and returning the best. Each (BO) run is given a maximum number of
$30\cdot d$ function evaluations.

##### Reproducibility and Open Science

The implementation of methods as well as reproducible scripts for all
experiments are publicly made available. Each experiment is replicated $20$
times based $20$ different seeds fixed across all variants. More details on
the code and on computational can be found in Appendix E.

###  5.2 H1: More accurate interpretations

Our experiments support hypothesis H1, i.e., we can achieve more accurate PD
estimates more efficiently through targeted sampling via the
$\textrm{EIG}_{\textrm{PDP}}$. An example run on the Branin function shown in
Figure 2 illustrates the behavior of the methods that is observable across all
experiments: BAX is yielding clearly more accurate PDPs than BO with EI
already after few iterations. Figure 4 in Appendix C.2 supports that PDs
estimated on data produced by BO with EI might provide not only
quantitatively, but also qualitatively wrong information in terms of ranking
the values $\varphi_{\bm{\lambda}_{S}^{(g)}}$ differently than the ground-
truth. As expected, increased accuracy of interpretations through BAX comes to
the cost of optimization efficiency. Results aggregated across all problems
and replications confirm this behavior on a broader scale, see Table 1666We
note that the different functions live on different scales s.t. we normalized
it by showing relative metrics wrt baselines, such RS for PDP estimates and EI
for optimization regret.. BAX is producing more accurate PDPs than RS (which
can be assumed to converge against the true marginal effect) already at early
stages, and is strongly significantly ($\alpha=1\%$) outperforming RS with
less iterations. We conclude that both BAX and PVAR can contribute to
approximating the true marginal effect well, but BAX is converging faster. In
addition, BO with EI is significantly outperformed in terms of accuracy of
PDPs, which supports our assumption of lowered quality caused through a heavy
sampling bias.

Table 1: Left: $L1$ error of the estimated PDP w.r.t. the ground truth PDP,
relative to RS as baseline. Negative values mean a relative reduction of the
$L1$ error compared to random search. Right: Optimization relative to BO-EI as
baseline. Results are averaged across all $20$ replications. Best values are
bold, and values are underlined if not significantly worse than the best based
on a Post-Hoc Friedman test ($\alpha=1\%$), see also Demsar (2006); García et
al. (2010) and Appendix C.1 for more details.

|  Relative $\textrm{d}_{L_{1}}(\textrm{PDP})$ after
---|---
| 25% | 50% | 75% | 100%
| Max. iterations spent
RS | 0.00 | 0.00 | 0.00 | 0.00
BO-EI | 0.18 | 0.39 | 0.47 | 0.67
PVAR | 0.13 | -0.08 | 0.08 | 0.14
BAX | -0.17 | -0.20 | -0.07 | 0.00
BOBAX | -0.14 | -0.16 | -0.04 | 0.03

| Relative optimization regret after
---|---
| 25% | 50% | 75% | 100%
| Max. iterations spent
RS | 2.42 | 160.99 | 530.70 | 951.47
BO-EI | 0.00 | 0.00 | 0.00 | 0.00
PVAR | 3.38 | 232.14 | 741.69 | 1887.22
BAX | 2.27 | 242.062 | 602.15 | 1408.62
BOBAX | 1.68 | 5.04 | 4.73 | 3.26

###  5.3 H2: More accurate interpretations at no relevant loss of optimization
efficiency

Our experiments also support hypothesis H2, i.e., with BOBAX we can achieve
clearly more accurate PD estimates while maintaining a competitive level of
optimization efficiency. Table 1 compares the accuracy of PD estimates
(measured via $\textrm{d}_{L_{1}}$) and optimization regret as compared to
baselines RS and BO-EI, respectively, aggregated over all five objective
functions. (BO)BAX allows for more accurate PDPs than the other methods, with
diminishing relative distance to RS, while BO with EI is clearly outperformed.
On the other hand, it can be observed that BOBAX is giving optimization
performance comparable to BO with EI throughout the course of optimization,
whereas RS is clearly outperformed. So, BOBAX combines the best of both
worlds: good interpretability (even better than RS) and efficient optimization
(on par with BO-EI). Figure 5 in Appendix C.2 shows that this effect is
visible for all objective function, but the strength of the effect depends on
the objective functions.

We conclude that our experiments support that BOBAX makes no (or only little)
compromises in optimization performance, but yields clearly better estimates
of marginal effects at the same time.

##  6 Practical HPO Application

Figure 3: Comparing PDP evolution for number of iterations for EI and BAX. BAX returns fairly certain PDPs early on, whereas BO with EI requires much more time. Table 2: Iterations needed to reach the desired precision of PD estimate of $\pm 1.5$ balanced accuracy points, accuracy of the final PD based on the L1 error to the ground truth, as well as the final model performance reached. Results are averaged across all 30 replications and all 15 datasets. Best values are bold, and values are underlined if not significantly worse than the best based on a Post-Hoc Friedman test ($\alpha=1\%$), see also Demsar (2006); García et al. (2010) and Appendix C.1 for more details.  | Iterations to desired precision | Rel. $\textrm{d}_{L_{1}}$ (PDP) | 1 - Balanced Accuracy
---|---|---|---
RS | 14.91 | 0.49 | 22.56
BO-EI | 22.59 | 0.57 | 19.38
BAX | 9.85 | 0.51 | 23.97
a-BOBAX | 11.56 | 0.52 | 19.94

We demonstrate a-BOBAX on a concrete HPO scenario, following the setup of
Moosbauer et al. (2021). We tune common hyperparameters of a neural network
with regards to balanced validation accuracy on 15 different datasets
respresenting different tasks from different domains (see Tables 5, 6 in
Appendix D) using the interface provided by YAHPO gym (Pfisterer et al.,
2021). We compare RS, EI, BAX, and adaptive BOBAX (a-BOBAX). In a-BOBAX, we
set the desired width of confidence intervals to $\pm 1.5\%$ balanced accuracy
points; we emphasize thought, that this value can be set by the user. For
a-BOBAX, we compute the $\textrm{EIG}_{\textrm{PDP}}$ jointly for the PDPs of
_learning rate_ , _dropout_ , _max. number of units_ , _weight decay_ , and
_momentum_. The respective methods ran under the same conditions as in Section
5, but were replicated $30$ times.

Figure 3 shows how accuracy of the PD estimate increases over time for BO with
EI vs. BAX. We observe that BAX is clearly more efficient in returning an
accurate estimate, which is in line with the results we observed in Section 5.
As motivated in Section 4.2, a practitioner might prefer to rather ensure a
minimum accuracy of IML measures, and therefore, handle this rather as a
constraint than as an objective. Table 2 is showing the time to reach the
desired precision of $\pm 1.5\%$ for the PDP, as well as final accuracy of PDs
and final optimization performance, aggregated over all experiments an
replications. We observe that a-BOBAX is (i) significantly faster in reaching
the desired precision threshold, allowing a user to interact earlier with
confidence, (ii) is comparable to RS in terms of final accurate representation
of PDs, and (iii) comparable to BO-EI in terms of optimization performance.
Note that the effect again depends strongly on the respective dataset (see
Figures 7, 8, 9 in Appendix D.2).

##  7 Discussion and Conclusion

##### Findings

We proposed (adaptive) BOBAX, modifying Bayesian Optimization (BO) for black-
box optimization and HPO to enhance interpretability of the optimization
problem at hand. We achieved this by adapting BAX to optimize for accurate
marginal effects and then interleaved BO and BAX. We further showed that BOBAX
can significantly enhance the accuracy of PD estimates during an optimization
procedure, while not losing optimization performance.

##### Usage

If a user has some desired precision of the IML estimates in mind, a-BOBAX
allows them to make use of BAX only until this level is not reached yet and
will focus on the optimization quality afterwards. This simple, yet efficient
strategy allows to get the most out of the overall budget.

##### Critical View and Limitations

Even though the usage of EIG is beneficial to the quality of a PD estimate,
there are also examples where no significant improvement is observed. We
assume that this particularly holds for hyperparameters that have a simple
(and therefore easy-to-learn) effect on performance. Consequently, the
marginal effect is easily learned for any of the methods. In addition to using
the adaptive version of BOBAX, we recommend dropping these simple-to-learn
hyperparameters from the joint computation of the EIG (4.3) as soon as the
PDPs are sufficiently certain. Furthermore, our method comes at a
computational overhead, being slightly larger than traditional BO since
computing EIG with BAX costs a bit more compute time. In terms of application
to HPO, we expect that the cost for training and validating hyperparameter
configurations or architectures of neural networks will be much larger than
BOBAX’s overhead in most relevant cases.

##### Outlook

We believe that BOBAX will contribute in particular towards more human-
centered HPO, where developers can start inspecting intermediate results as
soon as desired confidence was reached and then adapt the configuration space
if necessary. Although we focused on PDPs as an interpretability method,
extending our BOBAX idea to other IML approaches would be straightforward and
opens up new follow up directions. As one next step, we envision extending
BOBAX to the multi-fidelity setting (Li et al., 2017; Falkner et al., 2018)
which is required for more expensive HPO and AutoML problems. Last but not
least, we emphasize that we developed BOBAX primarily for HPO problems, but it
can also be applied to any black-box optimization problem, e.g., in
engineering or chemistry.

#### Acknowledgments

Use unnumbered third level headings for the acknowledgments. All
acknowledgments, including those to funding agencies, go at the end of the
paper. Only add this information once your submission is accepted and
deanonymized.

## References

  * Akiba et al. (2019) Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2019.
  * Bergstra et al. (2011) James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl.  Algorithms for hyper-parameter optimization.  _Advances in neural information processing systems_ , 24, 2011.
  * Biedenkapp et al. (2017) Andre Biedenkapp, Marius Lindauer, Katharina Eggensperger, Frank Hutter, Chris Fawcett, and Holger H. Hoos.  Efficient parameter importance analysis via ablation with surrogates.  In Satinder Singh and Shaul Markovitch (eds.), _Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, February 4-9, 2017, San Francisco, California, USA_ , pp. 773–779. AAAI Press, 2017.
  * Binder et al. (2020) Martin Binder, Julia Moosbauer, Janek Thomas, and Bernd Bischl.  Multi-objective hyperparameter tuning and feature selection using filter ensembles.  In Carlos Artemio Coello Coello (ed.), _GECCO ’20: Genetic and Evolutionary Computation Conference, Cancún Mexico, July 8-12, 2020_ , pp. 471–479. ACM, 2020.
  * Bischl et al. (2021) Bernd Bischl, Martin Binder, Michel Lang, Tobias Pielok, Jakob Richter, Stefan Coors, Janek Thomas, Theresa Ullmann, Marc Becker, Anne-Laure Boulesteix, Difan Deng, and Marius Lindauer.  Hyperparameter optimization: Foundations, algorithms, best practices and open challenges.  _CoRR_ , abs/2107.05847, 2021.
  * Carmichael et al. (2021) Zachariah Carmichael, Tim Moon, and Sam Ade Jacobs.  Learning interpretable models through multi-objective neural architecture search.  _CoRR_ , abs/2112.08645, 2021.
  * Coors et al. (2021) Stefan Coors, Daniel Schalk, Bernd Bischl, and David Rügamer.  Automatic componentwise boosting: An interpretable automl system.  _CoRR_ , abs/2109.05583, 2021.
  * Cutler et al. (2007) D Richard Cutler, Thomas C Edwards Jr, Karen H Beard, Adele Cutler, Kyle T Hess, Jacob Gibson, and Joshua J Lawler.  Random forests for classification in ecology.  _Ecology_ , 88(11):2783–2792, 2007.
  * Demsar (2006) Janez Demsar.  Statistical comparisons of classifiers over multiple data sets.  _J. Mach. Learn. Res._ , 7:1–30, 2006.
  * Drozdal et al. (2020) Jaimie Drozdal, Justin Weisz, Dakuo Wang, Gaurav Dass, Bingsheng Yao, Changruo Zhao, Michael Muller, Lin Ju, and Hui Su.  Trust in AutoML.  _Proceedings of the 25th International Conference on Intelligent User Interfaces_ , Mar 2020.
  * Falkner et al. (2018) Stefan Falkner, Aaron Klein, and Frank Hutter.  BOHB: robust and efficient hyperparameter optimization at scale.  In Jennifer G. Dy and Andreas Krause (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 1436–1445. PMLR, 2018\.
  * Fisher et al. (2019) Aaron Fisher, Cynthia Rudin, and Francesca Dominici.  All models are wrong, but many are useful: Learning a variable’s importance by studying an entire class of prediction models simultaneously.  _J. Mach. Learn. Res._ , 20:177:1–177:81, 2019.
  * Friedman (2001) Jerome H Friedman.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ , pp. 1189–1232, 2001.
  * García et al. (2010) Salvador García, Alberto Fernández, Julián Luengo, and Francisco Herrera.  Advanced nonparametric tests for multiple comparisons in the design of experiments in computational intelligence and data mining: Experimental analysis of power.  _Inf. Sci._ , 180(10):2044–2064, 2010.
  * Golovin et al. (2017) Daniel Golovin, Benjamin Solnik, Subhodeep Moitra, Greg Kochanski, John Karro, and D. Sculley.  Google vizier: A service for black-box optimization.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Halifax, NS, Canada, August 13 - 17, 2017_ , pp. 1487–1495. ACM, 2017.
  * Hasebrook et al. (2022) Niklas Hasebrook, Felix Morsbach, Niclas Kannengießer, Jörg K. H. Franke, Frank Hutter, and Ali Sunyaev.  Why do machine learning practitioners still use manual tuning? A qualitative study.  _CoRR_ , abs/2203.01717, 2022.
  * Head et al. (2022) Tim Head, Manoj Kumar, Holger Nahrstaedt, Gilles Louppe, and Iaroslav Shcherbatyi.  scikit-optimize/scikit-optimize, April 2022.  URL https://doi.org/10.5281/zenodo.6451894.
  * Hennig & Schuler (2012) Philipp Hennig and Christian J. Schuler.  Entropy search for information-efficient global optimization.  _J. Mach. Learn. Res._ , 13:1809–1837, 2012.
  * Hernández-Lobato et al. (2014) José Miguel Hernández-Lobato, Matthew W. Hoffman, and Zoubin Ghahramani.  Predictive entropy search for efficient global optimization of black-box functions.  In Zoubin Ghahramani, Max Welling, Corinna Cortes, Neil D. Lawrence, and Kilian Q. Weinberger (eds.), _Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada_ , pp. 918–926, 2014\.
  * Hutter et al. (2011) Frank Hutter, Holger H. Hoos, and Kevin Leyton-Brown.  Sequential model-based optimization for general algorithm configuration.  In Carlos A. Coello Coello (ed.), _Learning and Intelligent Optimization - 5th International Conference, LION 5, Rome, Italy, January 17-21, 2011. Selected Papers_ , volume 6683 of _Lecture Notes in Computer Science_ , pp. 507–523. Springer, 2011.
  * Hutter et al. (2014) Frank Hutter, Holger H. Hoos, and Kevin Leyton-Brown.  An efficient approach for assessing hyperparameter importance.  In _Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014_ , volume 32 of _JMLR Workshop and Conference Proceedings_ , pp. 754–762. JMLR.org, 2014\.
  * Hvarfner et al. (2022) Carl Hvarfner, Danny Stoll, Artur L. F. Souza, Marius Lindauer, Frank Hutter, and Luigi Nardi.  $\pi$bo: Augmenting acquisition functions with user beliefs for bayesian optimization.  _CoRR_ , abs/2204.11051, 2022.  doi: 10.48550/arXiv.2204.11051.  URL https://doi.org/10.48550/arXiv.2204.11051.
  * Jones (2001) Donald R. Jones.  A taxonomy of global optimization methods based on response surfaces.  _Journal of Global Optimization_ , 21(4):345–383, 2001.
  * Jones et al. (1998) Donald R. Jones, Matthias Schonlau, and William J. Welch.  Efficient global optimization of expensive black-box functions.  _J. Glob. Optim._ , 13(4):455–492, 1998.
  * Lemmens & Croux (2006) Aurélie Lemmens and Christophe Croux.  Bagging and boosting classification trees to predict churn.  _Journal of Marketing Research_ , 43(2):276–286, 2006.
  * Li et al. (2017) Lisha Li, Kevin G. Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar.  Hyperband: A novel bandit-based approach to hyperparameter optimization.  _J. Mach. Learn. Res._ , 18:185:1–185:52, 2017.
  * Lindley (1956) D. V. Lindley.  On a measure of the information provided by an experiment.  _The Annals of Mathematical Statistics_ , 27(4):986–1005, 1956.  ISSN 00034851.
  * Matthews et al. (2017) Alexander G. de G. Matthews, Mark van der Wilk, Tom Nickson, Keisuke. Fujii, Alexis Boukouvalas, Pablo León-Villagrá, Zoubin Ghahramani, and James Hensman.  GPflow: A Gaussian process library using TensorFlow.  _Journal of Machine Learning Research_ , 18(40):1–6, apr 2017.
  * Moosbauer et al. (2021) Julia Moosbauer, Julia Herbinger, Giuseppe Casalicchio, Marius Lindauer, and Bernd Bischl.  Explaining hyperparameter optimization via partial dependence plots.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), _Advances in Neural Information Processing Systems_ , volume 34, pp. 2280–2291. Curran Associates, Inc., 2021.
  * Neiswanger et al. (2021) Willie Neiswanger, Ke Alexander Wang, and Stefano Ermon.  Bayesian algorithm execution: Estimating computable properties of black-box functions using mutual information.  In Marina Meila and Tong Zhang (eds.), _Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 8005–8015. PMLR, 2021.
  * Olson & Moore (2016) Randal S. Olson and Jason H. Moore.  TPOT: A tree-based pipeline optimization tool for automating machine learning.  In _Proceedings of the 2016 Workshop on Automatic Machine Learning_ , volume 64 of _JMLR Workshop and Conference Proceedings_ , pp. 66–74. JMLR.org, 2016.
  * Pfisterer et al. (2021) Florian Pfisterer, Lennart Schneider, Julia Moosbauer, Martin Binder, and Bernd Bischl.  YAHPO gym - design criteria and a new multifidelity benchmark for hyperparameter optimization.  _CoRR_ , abs/2109.03670, 2021.
  * Probst et al. (2019) Philipp Probst, Anne-Laure Boulesteix, and Bernd Bischl.  Tunability: Importance of hyperparameters of machine learning algorithms.  _Journal of Machine Learning Research_ , 20:53:1–53:32, 2019.
  * Snoek et al. (2012) Jasper Snoek, Hugo Larochelle, and Ryan P. Adams.  Practical Bayesian optimization of machine learning algorithms.  In _Advances in Neural Information Processing Systems 25_ , pp. 2960–2968, 2012.
  * Srinivas et al. (2010) Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias W. Seeger.  Gaussian process optimization in the bandit setting: No regret and experimental design.  In Johannes Fürnkranz and Thorsten Joachims (eds.), _Proceedings of the 27th International Conference on Machine Learning_ , pp. 1015–1022. Omnipress, 2010.
  * Thornton et al. (2013) Chris Thornton, Frank Hutter, Holger H. Hoos, and Kevin Leyton-Brown.  Auto-weka: combined selection and hyperparameter optimization of classification algorithms.  In _The 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pp. 847–855. ACM, 2013.
  * Turner et al. (2020) Ryan Turner, David Eriksson, Michael McCourt, Juha Kiili, Eero Laaksonen, Zhen Xu, and Isabelle Guyon.  Bayesian optimization is superior to random search for machine learning hyperparameter tuning: Analysis of the black-box optimization challenge 2020.  In _NeurIPS 2020 Competition and Demonstration Track_ , volume 133 of _Proceedings of Machine Learning Research_ , pp. 3–26. PMLR, 2020\.
  * Van Rijn & Hutter (2018) Jan N Van Rijn and Frank Hutter.  Hyperparameter importance across datasets.  In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pp. 2367–2376, 2018.
  * Wenger & Olden (2012) Seth J Wenger and Julian D Olden.  Assessing transferability of ecological models: an underappreciated aspect of statistical validation.  _Methods in Ecology and Evolution_ , 3(2):260–267, 2012.
  * Wu et al. (2017) Jian Wu, Matthias Poloczek, Andrew Gordon Wilson, and Peter I. Frazier.  Bayesian optimization with gradients.  In _Advances in Neural Information Processing Systems 30_ , pp. 5267–5278, 2017.
  * Xanthopoulos et al. (2020) Iordanis Xanthopoulos, Ioannis Tsamardinos, Vassilis Christophides, Eric Simon, and Alejandro Salinger.  Putting the human back in the AutoML loop.  In _Proceedings of the Workshops of the EDBT/ICDT 2020 Joint Conference_ , volume 2578 of _CEUR Workshop Proceedings_. CEUR-WS.org, 2020\.
  * Young et al. (2018) M. Todd Young, Jacob D. Hinkle, Arvind Ramanathan, and Ramakrishnan Kannan.  Hyperspace: Distributed bayesian hyperparameter optimization.  In _30th International Symposium on Computer Architecture and High Performance Computing, SBAC-PAD 2018, Lyon, France, September 24-27, 2018_ , pp. 339–347. IEEE, 2018.  doi: 10.1109/CAHPC.2018.8645954.  URL https://doi.org/10.1109/CAHPC.2018.8645954.
  * Zela et al. (2018) Arber Zela, Aaron Klein, Stefan Falkner, and Frank Hutter.  Towards automated deep learning: Efficient joint neural architecture and hyperparameter search.  _CoRR_ , abs/1807.06906, 2018.  URL http://arxiv.org/abs/1807.06906.
  * Zhang et al. (2018) Zhongheng Zhang, Marcus W Beck, David A Winkler, Bin Huang, Wilbert Sibanda, Hemant Goyal, et al.  Opening the black box of neural networks: methods for interpreting neural network models in clinical applications.  _Annals of translational medicine_ , 6(11), 2018.
  * Zimmer et al. (2021) Lucas Zimmer, Marius Lindauer, and Frank Hutter.  Auto-PyTorch Tabular: Multi-fidelity metalearning for efficient and robust AutoDL.  _IEEE TPAMI_ , 2021.  Preprint via Early Access.
  * Zöller et al. (2022) Marc-André Zöller, Waldemar Titov, Thomas Schlegel, and Marco F. Huber.  Xautoml: A visual analytics tool for establishing trust in automated machine learning.  _CoRR_ , abs/2202.11954, 2022.

##  Appendix A Appendix

##  Appendix B Additional methodological aspects

###  B.1 Interpretability methods beyond the PDP

BOBAX is generic in the sense that it can be applied to other IML methods than
the PDP that are of interest to the user, as long as a the execution path of
the respective method is accessible to BOBAX.

While we considered the partial dependence method to estimate main effects
(i.e., the marginal effect of a single hyperparameter $\bm{\lambda}_{s}$ on
estimated performance) in our experiments, Algorithm 1 can be extended to
estimate interaction effects of two hyperparameters $S=\\{s,s^{\prime}\\}$.
This is done by simply replacing the grid points in Algorithm 1 by a two-
dimensional grid
$\left(\bm{\lambda}_{s}^{(g)},\bm{\lambda}_{s^{\prime}}^{(g^{\prime})}\right)$
for all pairs $g,g^{\prime}\in\\{1,2,...,G\\}$ with
$\left(\bm{\lambda}_{s}^{(1)},...,\bm{\lambda}_{s}^{(G)}\right)$ and
$\left(\bm{\lambda}_{s^{\prime}}^{(1)},...,\bm{\lambda}_{s^{\prime}}^{(G)}\right)$
representing equidistant grids. With this modified execution path our method
is be straightforwardly applied to estimate interaction effects.

Also, other methods within IML can be optimized for with BOBAX; for example,
the hyperparameter importance via permutation feature importance (PFI) (Fisher
et al., 2019). Importance of a single hyperparameter $\bm{\lambda}_{S}$ is
computed by shuffling the the values of this hyperparameter in the $A_{T}$,
resulting in a modified archive $\tilde{A}_{T,\bm{\lambda}_{S}}$ and the
difference in errors of the model $\hat{c}$ on $A_{T}$ and on
$\tilde{A}_{T,\bm{\lambda}_{S}}$ is compared. The respective execution path
$e_{\mathcal{A}}$ is the joint set of all shuffled versions of the archive
$\tilde{A}_{T,\bm{\lambda}_{1}}\cup\tilde{A}_{T,\bm{\lambda}_{2}}\cup...\cup\tilde{A}_{T,\bm{\lambda}_{d}}$.

###  B.2 BOBAX and Variants

Algorithm 2 BOBAX

Input  $k$, $n_{\text{init}}$, $O_{\mathcal{A}}$

$A_{T}\leftarrow$ Sample initial design of size $n_{\text{init}}$ over
$\Lambda$

while stopping criterion not met do

if $T\mod k=0$ then

$\bm{\lambda}^{(T+1)}\leftarrow\arg\max_{\bm{\lambda}\in\Lambda}\textrm{EIG}_{\textrm{PDP}}(\bm{\lambda})$

else

$\bm{\lambda}^{(T+1)}\leftarrow\arg\max_{\bm{\lambda}\in\Lambda}\text{EI}(\bm{\lambda})$

end if

$c_{\bm{\lambda}^{(T+1)}}\leftarrow c(\bm{\lambda}^{(T+1)})$

$A_{T+1}\leftarrow
A_{T}\cup\\{\left(\bm{\lambda}^{(T+1)},c_{\bm{\lambda}^{(T+1)}}\right)\\}$

$T\leftarrow T+1$

end while

Return  $A_{T},O_{\mathcal{A}}(\hat{c})$

Algorithm 2 shows the BOBAX algorithm as introduced and discussed in the main
paper. We have investigated two more alternative acquisition functions to
trade-off interpretability and optimization efficiency. One is a probabilistic
variant of interleaving $\textrm{EIG}_{\textrm{PDP}}$, where in every
iteration

$\displaystyle\bm{\lambda}^{(T+1)}=\textrm{arg
max}_{\bm{\lambda}\in\Lambda}\begin{cases}\textrm{EIG}_{\textrm{PDP}}(\bm{\lambda})&\textrm{if
}p\leq\pi\\\ \textrm{EI}(\bm{\lambda})&\textrm{if }p>\pi\end{cases}$

where $p\sim\textrm{Unif}(0,1)$ and $\pi$ is a threshold set by a user. If
$\pi$ is set to $0.5$ this corresponds to the probabilistic variant of
Algorithm 2 with $k=2$. We call this variant
$\textrm{BOBAX}_{\textrm{prob}}^{\pi}$. This method also opens up the
possibility to reduce the relative amount search for interpretability (as a
kind of exploration) over time by an annealing strategy where the probability
$\pi$ is lowered over time.

As a second variant, we investigated a multiplicative variant of
$\textrm{EIG}_{\textrm{PDP}}$ and EI inspired by Hvarfner et al. (2022):

$\displaystyle\textrm{EIBAX}^{\beta}(\bm{\lambda})=\textrm{EI}(\bm{\lambda})\cdot\textrm{EIG}_{\textrm{PDP}}(\bm{\lambda})^{\beta/T},$

where the values of a sampled batch of
$\textrm{EIG}_{\textrm{PDP}}(\bm{\lambda})$ are min-max-scaled to $[0,1]$.
Note that in comparison to the interleaving strategy, this method has a
computational disadvantage since it requires to compute the
$\textrm{EIG}_{\textrm{PDP}}$ in _every_ iteration.

Note that in any of the variants above, the EI can be replaced by any other
acquisition function.

##  Appendix C Benchmark

###  C.1 Additional Details

##### Details on evaluation

We performed a statistical test to allow for conclusions as to whether the
methods compared ($\textbf{RS},\textbf{EI},\textbf{BAX},\textbf{BOBAX}$) are
performing significantly differently in terms of (1) quality of the PD
estimate measured by $\textrm{d}_{\textrm{L}_{1}}$, (2) optimization
performance as measured by regret in Table 1. We applied a _Friedman aligned
ranks test_ as described in (García et al., 2010) on the respective
performance values on different objective functions and replications to
conclude whether there is a difference between methods. Note that the chosen
test is recommended over the Friedman test by García et al. (2010) in
particular if the number of algorithms is low (four to five) because of an
increased power. We applied a post hoc test with Hommel correction for
multiple testing, and report statistical significance based on corrected
p-values. We rely on the implementation
scmamp777https://github.com/b0rxa/scmamp.

##### Comparison with additional baselines

As additional baselines, we are running BO with LCB
$\hat{c}(\bm{\lambda})+\tau\cdot\hat{s}^{2}(\bm{\lambda})$ acquisition
function with different values of $\tau\in\\{1,2,5\\}$, denoted by
$\textbf{LCB}^{1}$, $\textbf{LCB}^{2}$, $\textbf{LCB}^{5}$. Also, we are
running BO with interleaved random configurations every $k\in\\{2,5,10\\}$
iterations, denoted by $\textbf{BO-RS}^{2},\textbf{BO-RS}^{5},\textbf{BO-
RS}^{10}$. We are in addition considering different variations of the BOBAX
method as described in Section B.2: We consider $\textbf{EIBAX}^{20}$,
$\textbf{EIBAX}^{50}$, $\textbf{EIBAX}^{10}$, as well as
$\textbf{BOBAX}^{0.5}_{\text{prob}}$. Also, we have run BOBAX for different
degrees of random interleaving $k\in\\{2,5,10\\}$, denoted by
$\textbf{BOBAX}^{2},\textbf{BOBAX}^{5},\textbf{BOBAX}^{10}$.

Note that all (BAX) variants optimize for a PD for one variable only; we have
chosen the first variable as default. To support our claims in Section 4.3
that our method can be easily applied to jointly compute the PDP for multiple
variables, we are also comparing to one variant which computes the PDP for
_all_ variables, denoted by $\textbf{BAX}_{\textrm{all}}$ and compare it to
BAX.

##### Technical details

All experiments only require CPUs (and no GPUs) and were computed on a Linux
cluster (see Table 3).

Table 3: Description of the infrastructure used for the experiments in this
paper.  Computing Infrastructure
---
Type | Linux CPU Cluster
Architecture | 28-way Haswell-EP nodes
Cores per Node | 1
Memory limit (per core) | 2.2 GB

##### Implementation details

Our implementation of BOBAX is based on the implementation provided by
(Neiswanger et al., 2021)888https://github.com/willieneis/bayesian-algorithm-
execution, which in turn is based on the GPflow (Matthews et al., 2017)
implementation for Gaussian processes.

Note that we are not optimizing the hyperparameters of the GP (lengthscale,
kernel variance, and nugget effect) during BOBAX to eliminate one source of
variance between methods. Instead, similarly to (Neiswanger et al., 2021), we
are setting those parameters to sensible default values. These are determined
by the following heuristic executed prior to all experiments: For every
objective function, we perform maximum likelihood optimization of these GP
hyperparameters based on $200$ randomly sampled points, and choose the
configuration with the highest likelihood. This configuration is fixed across
all replications and methods. While this heuristic does not impact the
expressiveness of our statements since all methods are based on the same
kernel hyperparameters, we emphasize that choosing appropriate hyperparameters
is crucial for the performance of our method; therefore, a stable
implementation (as done in established BO libraries) is regarded a necessary
requirement for practical usage.

###  C.2 Additional Results

First of all, to provide some evidence for our claim that that BO with EI can
return inaccurate PDPs not only in absolute terms but also when considering
ranks, we have computed Spearman’s rank correlation of the respective PD
estimate with the ground truth objective (see Figure 4).

To evaluate many different algorithms based on two criteria (1) error in PDP
estimate $\textrm{d}_{\textrm{dL}_{1}}$ and (2) optimization regret in a
compressed way, we are looking at the ranks of different methods with regards
to both metrics, resulting in two ranks
$\textrm{rank}_{\textrm{d}_{\textrm{dL}1}}$,
$\textrm{rank}_{\textrm{regret}}$. For the sake of evaluation we assume that
interpretability and optimization efficiency are of equal importance and
therefore assign each method a combined rank of
$\frac{1}{2}\cdot\textrm{rank}_{\textrm{d}_{\textrm{dL}1}}+\frac{1}{2}\cdot\textrm{rank}_{\textrm{regret}}$.
We average the combined ranks of every method across replications and problem
instances. Table 4 shows the combined ranks for our proposed methods BAX and
BOBAX (introduced in Section 2) as well as all baselines.

Figure 6 compares the $\textrm{EIG}_{\textrm{PDP}}$ computed w.r.t. the PD of
a single variable vs. jointly for the PDs of all variables. We observe that
there is no drop in performance; in particular, we observe that the joint
computation performs comparably to the computation for a single variable when
evaluated on a single variable; and the joint computation performs better, if
the accuracy of _all_ PDPs is considered.

Table 4: The table shows the combined ranks $\frac{1}{2}\cdot\textrm{rank}_{\textrm{d}_{\textrm{dL}1}}+\frac{1}{2}\cdot\textrm{rank}_{\textrm{regret}}$ of different methods introduced in Section 4 as well as additional baselines introduced in Appendix C.1. Results are averaged across $20$ replications and across all problems. We observe that $\text{BOBAX}^{2}$ is best in terms of the combined rank.  | Combined ranks after
---|---
| 25% | 50% | 75% | 100%
| Max. iterations spent
BOBAX2 | 6.30 | 5.24 | 4.88 | 4.88
BOBAX5 | 6.36 | 5.91 | 5.23 | 5.08
BOBAX10 | 6.51 | 6.10 | 5.46 | 5.14
BO-RS2 | 7.65 | 7.02 | 5.88 | 5.49
BOBAX${}_{\textrm{prob}}^{0.5}$ | 6.96 | 6.39 | 5.92 | 5.72
BO-RS5 | 7.24 | 6.60 | 5.92 | 5.73
BO-RS10 | 7.37 | 6.64 | 6.18 | 5.78
EIBAX100 | 6.71 | 6.40 | 6.00 | 5.94
BAX | 6.77 | 6.95 | 6.32 | 6.18
EIBAX20 | 6.91 | 6.79 | 6.10 | 6.20
LCB5 | 8.93 | 6.65 | 6.09 | 6.22
EI | 7.53 | 7.44 | 6.67 | 6.28
EIBAX50 | 6.83 | 6.45 | 6.21 | 6.49
RS | 8.82 | 9.19 | 7.74 | 6.92
PVAR | 9.80 | 8.01 | 7.18 | 7.00
LCB2 | 8.94 | 7.30 | 7.30 | 7.54
LCB1 | 8.36 | 7.99 | 8.20 | 8.10
Figure 4: The figure shows Spearman’s rank correlation of the estimated PDP
vs. the iterations performed for the Branin function. It demonstrates that
looking at PDPs computed on data from BO with expected improvement can be even
wrong in terms of correlation (as compared to BAX and RS), which matters a lot
in the context of optimization.

Figure 5: Error of PD estimates measured via measured by $\textrm{d}_{L_{1}}$
(left) and optimization regret (right) for the different synthetic objectives.
While RS is clearly outperformed in terms of optimization efficiency, BOBAX
and BO with EI perform comparable on this problem instance.  Figure 6: The
performance of BOBAX with $\textrm{EIG}_{\textrm{PDP}}$ computed with regards
to the first variable only (blue) vs. the performance of BOBAX when
$\textrm{EIG}_{\textrm{PDP}}$ is computed for the joint execution paths of PD
estimates with regards to _all_ variables (orange). Left: Error of the PD
estimate for the _first_ variable (measured via
$\textrm{d}_{\textrm{L}_{1}}$). Right: Error of the PD estimate for the _all_
variables (measured via $\textrm{d}_{\textrm{L}_{1}}$). We observe that the
joint computation delivers more accurate PDs over _all_ variables. However, we
also observe that the difference is not dramatically big.

##  Appendix D Practical HPO Application

###  D.1 Additional Details

Table 5: Hyperparameter space of the LCBench (Zimmer et al., 2021) benchmark suite within YAHPO gym (Pfisterer et al., 2021); _batch size_ and _maximum number of layers_ have been set to defaults $512$ and $5$, respectively.  Name | Range | log | type
---|---|---|---
Max. number of units | $[64,512]$ | yes | int
Learning rate (SGD) | $[1\textrm{e}^{-4},1\textrm{e}^{-1}]$ | yes | float
Weight decay | $[1\textrm{e}^{-5},1\textrm{e}^{-1}]$ | no | float
Momentum | $[0.1,0.99]$ | no | float
Max. dropout rate | $[0.0,1.0]$ | no | float
Table 6: Datasets accessed via the _lcbench_ suite of YAHPO gym (Pfisterer et al., 2021); the underlying data for the surrogate benchmark was made available by (Zimmer et al., 2021).  ID | Name | Usecase | $n$ | $d$
---|---|---|---|---
3945 | KDDCup09_appetency | Prediciton of customer behavior | 50000 | 231
34539 | drug-directory | Drug classification | 120215 | 21
7593 | covertype | Forest cover type | 581012 | 55
126025 | adult | Salary prediction | 48842 | 15
126026 | nomao | Active-learning in real-world | 34465 | 119
126029 | bank-marketing | Bank direct marketing | 4521 | 17
146212 | shuttle |  | 58000 | 10
167104 | Australian | Credit approval | 690 | 15
167149 | kr-vs-kp | Chess game | 3196 | 37
167152 | mfeat-factors | Handwritten numerals | 2000 | 217
167161 | credit-g | Credit risk prediciton | 1000 | 21
167168 | vehicle | Classification of vehicles | 846 | 22
167185 | cnae-9 | Classification of free text | 1080 | 857
167200 | higgs | Higgs boson detection | 98050 | 29
189908 | Fashion-MNIST | Classification of Zalando’s article images | 70000 | 785

As practical HPO application we have chosen the use case of tuning
hyperparameters of a neural network (as shown in Table 5) on the different
classification tasks (listed in Table 6) with regards to _Balanced accuracy_
as performance measures. In BAX / a-BOBAX, we are computing the
$\textrm{EIG}_{\textrm{PDP}}$ jointly for the PDPs of all hyperparameters
listed in Table 5. Each run is replicated $10$ times. Otherwise, all other
settings correspond to the settings in Sections 5 and Appendix C. Note that
the benchmark provided via Yahpo Gym (Pfisterer et al., 2021) is a surrogate
benchmark, which not only supports efficient execution of a benchmark, but
also gives access to a (reasonably cheap-to-evaluate) empirical performance
model as ground truth objective; allowing us to compute the ground-truth PDP
(and thus, any measure of error of the PDP) based on this empirical
performance model.

###  D.2 Additional Results

Figures 7, 8, 9 shows a more granular representation of results for the HPO
usecase.

Figure 7: The figure compares error of the PDP estimate after the full budget
spent (in terms of $\textrm{dL}_{1}$; shown in the first row), the percentage
of iterations needed to reach the desired level of confidence (middle row), as
well as the final regret (last row) for the different methods a-BOBAX, EI, and
RS on the different datasets (columns) that we tuned for. In most cases,
a-BOBAX has a final error in PDP comparable to RS, but clearly better than
with EI, and reaches the desired level of confidence faster then the two other
methods. In terms of optimization performance, a-BOBAX and EI perform
comparably, and both clearly outperform RS.  Figure 8: The figure compares
error of the PDP estimate after the full budget spent (in terms of
$\textrm{dL}_{1}$; shown in the first row), the percentage of iterations
needed to reach the desired level of confidence (middle row), as well as the
final regret (last row) for the different methods a-BOBAX, EI, and RS on the
different datasets (columns) that we tuned for. In most cases, a-BOBAX has a
final error in PDP comparable to RS, but clearly better than with EI, and
reaches the desired level of confidence faster then the two other methods. In
terms of optimization performance, a-BOBAX and EI perform comparably, and both
clearly outperform RS.  Figure 9: The figure compares error of the PDP
estimate after the full budget spent (in terms of $\textrm{dL}_{1}$; shown in
the first row), the percentage of iterations needed to reach the desired level
of confidence (middle row), as well as the final regret (last row) for the
different methods a-BOBAX, EI, and RS on the different datasets (columns) that
we tuned for. In most cases, a-BOBAX has a final error in PDP comparable to
RS, but clearly better than with EI, and reaches the desired level of
confidence faster then the two other methods. In terms of optimization
performance, a-BOBAX and EI perform comparably, and both clearly outperform
RS.

##  Appendix E Code and Implementation

All code and data needed to reproduce the benchmark will be made publicly
available via a Github repository after completion of the review process.
During review phase, all code is uploaded as a supplementary material, or can
alternatively be downloaded from https://figshare.com/s/d6ef1b8f4c9c1e844229.
Please refer to the README.md file for further information about how to use
the code to reproduce results.

Note that our implementation is based on the implementation provided by
Neiswanger et al. (2021)999https://github.com/willieneis/bayesian-algorithm-
execution.

Raw and processed results can be downloaded from
https://figshare.com/s/4573a2546f1d8a535c12.