# Interpretable Modeling of Single-cell perturbation Responses to Novel Drugs
Using Cycle Consistence Learning

Wei Huang, Aichun Zhu, Hui Liu*

###### Abstract

Phenotype-based screening has attracted much attention for identifying cell-
active compounds. Transcriptional and proteomic profiles of cell population or
single cells are informative phenotypic measures of cellular responses to
perturbations. In this paper, we proposed a deep learning framework based on
encoder-decoder architecture that maps the initial cellular states to a latent
space, in which we assume the effects of drug perturbation on cellular states
follow linear additivity. Next, we introduced the cycle consistency
constraints to enforce that initial cellular state subjected to drug
perturbations would produce the perturbed cellular responses, and, conversely,
removal of drug perturbation from the perturbed cellular states would restore
the initial cellular states. The cycle consistency constraints and linear
modeling in latent space enable to learn interpretable and transferable drug
perturbation representations, so that our model can predict cellular response
to unseen drugs. We validated our model on three different types of datasets,
including bulk transcriptional responses, bulk proteomic responses, and
single-cell transcriptional responses to drug perturbations. The experimental
results show that our model achieves better performance than existing state-
of-the-art methods.

## Introduction

Notwithstanding target-based drug discovery has made substantial advancements
in recent years, intervention of a specific target (protein or RNA) by
compounds has been insufficient in establishing the systematic correlation
with organism-level therapeutic effects or side effects. Consequently, the
failure rate of leading compounds generated from target-based screening to
approved drugs remains high. As such, there has been a renewed interest in
phenotypic drug discovery for the identification of cell-active compounds.

Transcriptional profiles serve as a robust and informative phenotypic measure
of cellular responses to perturbations. Large-scale compendia have been
established to examine the drug-induced phenotypic alterations across cancer
cell lines, including large-scale pharmacologic perturbation studies and cell
viability measurements upon different drug treatments. For example, the L1000
platform (Subramanian et al. 2017) has been developed for high-throughput
profiling of mRNA responses of cancer cell lines to diverse perturbations. In
parallel, proteomic responses of cancer cell lines to a diverse array of
clinically relevant drugs have also been generated using the reverse-phase
protein arrays (RPPAs). The cancer perturbed proteomics atlas (CPPA) profiled
large-scale drug responses of more than 200 clinically relevant proteins that
covered major targets for cancer therapy (Zhao et al. 2020). Transcriptional
and proteomic profiling reflect the multi-level regulatory state transition
upon external perturbations, providing sound measurements of cellular response
that greatly facilitate phenotype-based drug screening.

The single-cell RNA sequencing (scRNA-seq) can identify subtle changes in gene
expression and tumor heterogeneity at single-cell resolution, which is
important for distinguishing distinctive effects of certain perturbation on
cell subpopulations and identifying cellular subsets resistant to specific
drug. For example, sci-Plex (Srivatsan et al. 2020) uses nuclear hashing to
quantify global transcriptional responses to thousands of independent
perturbations at single-cell resolution. The single-cell profiling enables
interrogation of phenotypic heterogeneity at a level which has been hitherto
inaccessible. However, these high-throughput screening technologies is still
limited relative to the vast combinatorial landscape of all cell type (cell
number)-perturbation pairs. Therefore, computational model trained on the
observed experimental data to predict cellular responses within various
cellular contexts is of great importance.

In particular, the model capacity to predict cellular responses to unseen
perturbations (Yu and Welch 2022), particularly useful for drug repurposing,
may have significant medical implications. This requires the model to
effectively capture the intricate interactions between chemical components and
cellular molecules that may intrigue cascade biochemical reactions and
eventually drive the transition of molecular phenotype. The cellular response
to drug perturbation is controlled by underlying biological network, and the
cellular state transition is actually nonlinear. To build an interpretable
model, we assume the effects of drug perturbation on cellular states follow
linearly additive rule. Based on this assumption, we introduced the cycle
consistency constraints to enforce that initial cellular state subjected to
drug perturbations would produce the expected cellular responses, and,
conversely, removal of drug perturbation from the perturbed cellular states
would restore the initial cellular states. The cycle consistency constraints
and linear modeling enable to learn interpretable and transferable drug
perturbation representations, so that our model can predict cellular response
to unseen drugs. We validated our model on three different types of datasets,
including bulk transcriptional drug responses, bulk proteomic drug responses
and single-cell drug responses. The experimental results show that our model
achieves better performance than existing state-of-the-art methods.

We think this work has at least three contributions as below:

  * •

To our best knowledge, we are the first to introduce cycle consistency loss
into learning cellular responses to drug perturbations, which enables our
model to learn expressive and transferable drug representations.

  * •

We model the cellular response from two opposite perspectives, requiring the
model to simultaneously predict the cell state transitions from unperturbed
state (control) to perturbed state (treatment) and vice versa, which enforce
the encoder networks to capture the essential feature of drug perturbations to
cellular state.

  * •

We not only evaluate the proposed model on both bulk and single-cell
transcriptional responses, but also introduce a proteomic drug response
dataset to evaluate our model. The experimental results on these benchmark
datasets demonstrate the superior performance of our model. To our best
knowledge, this is the first to apply proteomic data to evaluate prediction
model of cellular response to novel drug perturbation.

## Related Works

### Predicting cellular response to perturbation

Some computational methods have been developed to predict cellular responses
to perturbations. Among them, mechanistic modeling has been leveraged to
predict cell viability or the abundance of specific proteins. While these
models are powerful at interpreting interactions, they typically require
longitudinal data, which is often unavailable in practice. Furthermore, most
mechanistic models do not scale well to genome-wide measurements or high-
dimensional scRNA-seq data, making them less suitable for predicting high-
dimensional responses. However, thanks to the development of deep learning
technology in recent years, this tool has been increasingly applied in the
analysis and interpretation of scRNA-seq data (Hetzel et al. 2021; Lopez,
Gayoso, and Yosef 2020).

Some machine learning methods that have been proposed for predicting cellular
responses to drug treatments. These methods include deep variational
autoencoder (Jia et al. 2021), kernelized Bayesian matrix factorization
(Madhukar et al. 2019), matrix factorization with similarity regularization
(Gao, Yang, and Wang 2021), convolutional neural network (Zhang, Lu, and Zang
2022). These methods leverage different techniques such as imputing drug
response through low embedding of multiple genes (Roohani, Huang, and Leskovec
2022), incorporating prior knowledge of pathway-drug associations (Chawla et
al. 2022), leveraging mutational signatures (Robichaux et al. 2021; Aissa et
al. 2021), and using gene expression data for prediction (Sharifi-Noghabi et
al. 2021; He et al. 2022).

### Linear model in latent space

The linear additive model in the latent space is widely used in deep learning
for interpretability. Specifically, these models use matrix factorization
techniques or deep generative models to predict drug response using linear
models in a low-dimensional latent space representation. CPA (Lotfollahi et
al. 2021) and chemCPA (Hetzel et al. 2022) are most related to our work, as
they combine the interpretability of linear models with the flexibility of
deep-learning approaches for single-cell response modeling. Although these
models generate easy-to-interpret embeddings for drugs and cells, their
accuracy is still insufficient to drive drug discovery.

### Cycle consistency

Cycle consistency is a concept used in many computer vision problems that
involve processing multiple entities. It is a way of seeking global agreement
by enforcing consistency between local relationships. One common way of
enforcing cycle consistency is through the use of a cycle consistency loss,
which is a type of loss function that encourages forward and backward
consistency between mappings. It was first proposed in CycleGAN (Zhu et al.
2017). In this essay, Cycle consistency refers to the property that the image
output by the first generator can be used as input to the second generator,
and the output of the second generator should match the original image. The
reverse is also true. This reduces the space of possible mapping functions and
helps to ensure that the mappings are consistent with each other. We are the
first to apply cycle consistency in predicting cellular responses to drug
perturbation, which enable our model to learn informative and transferable
drug perturbation features within various cellular contexts.

### Domain adaptation

Domain adaptation refers to the process of adapting a model from the source
domain to the target domain, where the data distributions of the source and
target domains are different. For general domain adaptation methods, they can
be classified into methods based on domain distribution differences (Tzeng et
al. 2019; Chen et al. 2019; Arjovsky, Chintala, and Bottou 2017), adversarial
learning-based methods (Wang, He, and Katabi 2020; Pei et al. 2018; Long et
al. 2018), reconstruction-based methods (Zheng et al. 2017), and sample
generation-based methods (Sankaranarayanan et al. 2018). In predicting drug
response, domain adaptation can be used to solve the problem of data
distribution mismatch. For example, knowledge learned from bulk-seq data can
be transferred to scRNA-seq data or patient data to predict drug response in
single cells or patients (Chen et al. 2022; Lotfollahi, Wolf, and Theis 2019;
Ma et al. 2021) which has important clinical significance.

## Cycle-consistent linear modeling

### Encoder-decoder architecture

We employ an encoder-decoder architecture to integrate the linear model and
cycle-consistent constraints to predict cellular drug response. For
simplicity, we refer to the proposed method as cycleCDR. Figure 1 shows the
illustrate diagram of our proposed learning framework. The molecular
signatures standing for cellular state serve as the input of an encoder, which
maps the cellular features of control samples into a latent space. The decoder
endeavors to yield expected cellular states, depending on the embedding to be
decoded is unperturbed or perturbed by drugs in the latent space.

Figure 1: The illustrative diagram of our learning framework. Our model
consists of two autoencoders, and they collaboratively learn and mutually
improve each other’s predictive performance. The effect of drugs on cell state
is linear model in the latent space so that our model is interpretable.

Formally, we define the encoder as $F_{\theta}$, which maps the cell state
$x\in R^{m}$ to an $l$-dimensional latent vector $z\in R^{l}$. The decoder is
defined as $G_{\varphi}$ that convert the latent vector into output space.
Denote by $x$ and $y$ the initial (untreated cells) and perturbed cellular
state (treated cells), and $z$ represents the mapped latent representation
whose dimension is equal to the size of autoencoder bottleneck layer. Our
encoder/decoder networks are fully or densely connected neural networks with
rectified linear unit (ReLU) activation function, $\theta$ and $\varphi$ are
the learnable parameters of the encoder and decoder. First, we require that
the encoder-decoder architecture functions as an autoencoder to map the input
molecular signatures to latent vector, and then recover the input signals. We
optimize their parameters to minimize the reconstruction loss as follows:

$\mathcal{L}_{reco}=\sum_{i=1}^{N}||x_{i}-G_{\varphi}(F_{\theta}(x_{i}))||_{2}^{2}$
(1)

The encoder converts the cellular states into low-dimensional but informative
representations in the latent space. More importantly, in the latent space we
can linearly model the cellular respose to drug perturbation (see Subsection
Linear modeling in latent space).

### Drug perturbation encoder

The graph attention network (GAT) was used to encode the drug perturbation
into representation in the latent space. We used the SMILES of the drug to
obtain the molecular graph $\mathcal{M}=<V,E>$, where $V$ is the set of nodes
(atoms) and $E$ is the set of edges (chemical bonds). Assuming that $h_{i}$ is
the embedding of node $i$ and $W$ is a learnable weight matrix, the attention
score $\alpha_{ij}$ between node $i$ and its first-order neighbor node $j$ can
be calculated using the following equation:

$\alpha_{ij}=\frac{\exp(elu(a^{T}(Wh_{i},Wh_{j})))}{\sum_{k\in\mathcal{N}(i)}\exp(elu(a^{T}(Wh_{i},Wh_{k})))}$
(2)

where $a$ is a learnable vector, $elu$ is the exponential linear unit
activation function, and $\mathcal{N}(i)$represents the first-order neighbor
of node $i$. The attention score $\alpha_{ij}$ was actually the softmax
normalized message between node $i$ and its neighbors. Once the attention
scores were computed, the output feature of node $i$ was computed by
aggregating its neighbor features weighted by corresponding attention scores:

$h_{i}=\sigma(a_{ii}Wh_{i}+\sum\limits_{j\in N(i)}\alpha_{ij}Wh_{j})$ (3)

where $\sigma(.)$ is the ReLU activation function.

### Linear modeling in latent space

Inspired by linear model in latent space, we assume that the effect of drug
perturbation on cellular states follow the linear additivity in the latent
space. Given that the cell state and drug disturbance are mapped to the same
latent space, we construct a linear and easy-to-interpret model of the
cellular drug response. Assuming that the $j$-th drug is mapped to $h_{j}$ via
the drug perturbation encoder, the drug-induced cellular response in the
latent space is defined as:

${z}_{ij}^{(+)}={z}_{i}+{h}_{j}$ (4)

in which ${z}_{ij}^{(+)}$ represents the perturbed cellular state of the
$i$-th cell by $j$-th drug in the latent space. Assuming that the actual cell
state induced by the drug is $y_{ij}$, the decoder $G_{\varphi}$ should yield
$\hat{y}_{ij}=G_{\varphi}({z}_{i}+{h}_{j})$ that approximates to $y_{ij}$ as
close as possible.

Based on the linear additivity, the cellular response to drug perturbation can
be modeled from the opposite direction. By subtracting the drug representation
in the latent space, the drug effect should be removed and the perturbed
cellular state could be restored to the perturbation-free state. Assume the
encoder $F_{\theta^{{}^{\prime}}}$ maps the actual cell response $y_{ij}$ to a
vector $z_{ij}^{\prime}=F_{\theta^{\prime}}(y_{ij})$ into the common latent
space, elimination of the drug interference can be directly subtracted by
$h_{j}$. The restored representation of cellular state in the latent space can
be defined as:

${z}_{ij}^{(-)}={z}_{ij}^{\prime}-{h}_{j}$ (5)

where $z_{ij}^{(-)}$ represents the latent representation of the $i$-th cell
type (single cell) eliminating the perturbation of the $j$-th drug.
Accordingly, another decoder $G_{\varphi^{\prime}}$ is used to map
$z_{ij}^{(-)}$ back to the cellular state that should be as close as possible
to $x_{i}$.

For scenarios with paired data (unperturbed vs. perturbed), we use the mean
squared error loss function as below:

$\begin{split}\mathcal{L}_{{MSE}}=\frac{1}{N*K}\sum_{i=1}^{N}\sum_{j=1}^{K}[(y_{ij}-G_{\varphi}(z_{i}+h_{j}))^{2}\\\
+(x_{i}-G_{\varphi^{\prime}}({z}_{ij}^{\prime}-{h}_{j}))^{2}]\end{split}$ (6)

The first term corresponds to loss in predicting the drug-perturbed cellular
state from the initial state, while the second term corresponds to the loss in
restoring the initial state from the drug-perturbed cell state.

However, for single-cell transcriptional response to drug perturbation, we
lose the one-to-one correspondence of individual cells before and after drug
treatment, because the cell body is destroyed in the single-cell RNA
sequencing assay. Therefore, we leverage adversarial learning to align the
data distribution between two domains (unperturbed state vs. perturbed state).
In this situation, the decoder $G_{\varphi}$ functions as a generator that
produces the perturbed cellular states. We introduce a discriminator
$D_{\psi}$ to distinguish between actual and generated cellular responses. The
adversarial loss is defined as follows:

$\begin{split}\mathcal{L}_{{GAN}}(G_{\varphi},{D}_{\psi},X,Y)={E}_{\psi}[\log{D}_{\psi}(y_{ij})]\\\
+{E}_{\varphi,\theta}[\log(1-{D}_{\psi}(G_{\varphi}({z}_{i}+{h}_{j})))]\end{split}$
(7)

where the generator $G_{\varphi}$ tries to produce cellular response
$G_{\varphi}({z}_{i}+{h}_{j})$ that look similar to actual response, while
$D_{\psi}$ aims to distinguish between generated samples and real samples.
$G_{\varphi}$ aims to minimize this objective against an adversary $D_{\psi}$
that tries to maximize it, i.e.,
$\min_{G_{\varphi}}\max_{D_{\psi}}{L}_{GAN}(G_{\varphi},{D}_{\psi},X,Y)$.
Correspondingly, another discriminator $D_{\psi^{\prime}}$ is introduced to
distinguish the actual initial state from generated cell states by
$G_{\varphi^{\prime}}$, and adversarial loss is defined as follows:

$\begin{split}\mathcal{L}_{GAN}(G_{\varphi^{\prime}},{D}_{\psi^{\prime}},{Y},{X})={E}_{\psi^{\prime}}[\log{D}_{\psi^{\prime}}(x_{i})]\\\
+{E}_{\varphi^{\prime},\theta^{\prime}}[\log(1-{D}_{\psi^{\prime}}(G_{\psi^{\prime}}({z}_{ij}^{\prime}-{h}_{j}))]\end{split}$
(8)

### Cycle-consistent loss

Although adversarial learning-based domain adaption could align the data
distribution of source domain and target domain, but a encoder network has
adequate capacity to map a specific set of cellular states to any random
permutation of cells within the target domain. Any of the learned mappings can
induce an output distribution that matches the target distribution. Thus,
adversarial losses alone cannot guarantee that the learned function can map an
individual input $x_{i}$ to a desired output $y_{i}$. To further reduce the
space of possible mappings, we require that the learned mapping functions
should be cycle-consistent.

Formally, the initial cellular state $x_{i}$ is mapped to the a drug-perturbed
state $\hat{y}_{i}j$ by drug $j$, which should be used to restore the initial
state. We thus require
$x_{i}\xrightarrow{{G}_{\varphi}(F_{\theta}(x_{i})+h_{j})}\widehat{y}_{ij}\xrightarrow{{G}_{\varphi^{\prime}}(F_{\theta^{\prime}}(\widehat{y}_{ij})-h_{j})}\widehat{x}_{i}\approx
x_{i}$. Correspondingly, the perturbed cellular state is mapped to an
unperturbed state that should be mapped to the corresponding perturbed state,
namely, we require
$y_{ij}\xrightarrow{G_{\varphi^{\prime}}(F_{\theta^{\prime}}(y_{ij})-h_{j})}\widehat{x}_{i}\xrightarrow{G_{\varphi}(F_{\theta}(\hat{x_{i}})+h_{j})}\widehat{y}_{ij}\approx
y_{ij}$. Thus, we define the cycle-consistency loss function as follows:

$\mathcal{L}_{{cyc}}=\sum_{i=1}^{N}(x_{i}-\hat{x}_{i})^{2}+\sum_{i=1}^{N}\sum_{j=1}^{K}(y_{ij}-\hat{y}_{ij})^{2}$
(9)

### Full objective

For the transcriptional or proteomic responses measured on cell population, we
have the paired data so that we define a relatively simple objective function
as below:

$\mathcal{L}=\mathcal{L}_{reco}+\mathcal{L}_{MSE}$ (10)

For the single-cell cellular response to drug perturbation, we have no paired
data and thus define the full objective function as:

$\mathcal{L}=\mathcal{L}_{reco}+\mathcal{L}_{GAN}+\lambda\mathcal{L}_{{cyc}}$
(11)

in which
$\mathcal{L}_{GAN}=\mathcal{L}_{{GAN}}(G_{\varphi},{D}_{\psi},X,Y)+\mathcal{L}_{GAN}(G_{\varphi^{\prime}},{D}_{\psi^{\prime}},{Y},{X})$,
$\lambda$ is the tradeoff parameter standing for the importance of cycle
consistence constraints. In addition, similar to Taigman et al. [49], we
regularize the generator produce nearly an identity mapping when real samples
of the target domain (perturbed) are provided as the input to the generator.
The identity loss is defined as
$\mathcal{L}_{identity}=\sum_{i}^{N}\sum_{j}^{K}(G_{\varphi}(F_{\theta}(y))-y)^{2}|$
and implicitly used in the full objective.

The two autoencoders included in our model function collaboratively and
mutually promote their predictive performance. In our ablation experiments, we
attempted to use a single autoencoder for the cycle consistence leaning,
namely $F_{\theta}=F_{\theta^{\prime}}$ and
$G_{\varphi}=G_{\varphi^{\prime}}$, and found that their performance is
comparative. We also evaluated the ablated model performance without
reconstruction loss, as well as the ablated model only adversarial loss alone
by removal of cycle consistency loss.

## Evaluation experiments

### Experimental settings

All the encoders and decoders for gene expression profiles were implemented
using a multi-layer MLP with Relu activation function and batch normalization.
The discriminator used in adversarial learning is implemented using a multi-
layer MLP to distinguish the domain of the gene expression profile. The
dimension of the bottleneck layer was set to 128. The drug molecular graph
encoder comprised of two layers of GAT with an additional fully connected
layer to adjust the embedding dimension. The multi-head attention mechanism
was applied to the first layer and the number of heads was set to 10. The
global max pooling over the node-level features is used to obtain the graph
embedding. The parameters of the GAT encoder are initialized using the
pretrained model based on attention-wise masked graph contrastive learning
(Liu et al. 2022). The hyperparameter $\lambda$ is set to 10 that used in
cycleGAN (Zhu et al. 2017).

To evaluate the performance of our model to predict cellular responses to drug
perturbations, we use the coefficient of determination $r^{2}$ as an
evaluation metric. For the evaluation of transcriptomic response prediction,
we compute the $r^{2}$ score and explained variance (EV) over all genes to
evaluate the performance in predicting transcriptome-wide response. Besides,
since the expression levels of most genes in the perturbed cells remains
similar to their control state, the performance metrics based on all genes are
high but does not reflect the true predictive capacity of a model. In
contrast, the differentially expressed genes (DEGs) can more faithfully
reflect the actual effect of drug interference on cellular state. Therefore,
we also calculated the $r^{2}$ and EV metrics based on 50 most significantly
differentially expressed genes.

### Performance evaluation on bulk transcriptional response

The L1000 dataset is a large collection of gene expression profiles that
measure the responses of human cell lines to various compounds. It consists of
approximately 1,400,000 gene expression profiles on the responses of about 50
human cell lines to one of about 20,000 compounds across a range of
concentrations. We obtained data from the L1000 website, and used only the
expression profiles treated by 10$\mu$M drug concentrations. The gene
expression values of technical replicates were averaged. As a result, we
obtained drug response expression profiles spanning 17,775 drugs and 42 cell
line types, including a total of 45,763 expression profiles over 978 landmark
genes.

Figure 2: Performance comparison of cycleCDR with the baseline model and
chemCPA method on L1000 bulk transcriptional response dataset.

To benchmark our model performance, we constructed a baseline model that
directly calculates the $r^{2}$ scores by discarding all perturbation
information. This baseline has been also adopted by previous study (Hetzel et
al. 2022). Also, we compared our model with another perturbation prediction
model chemCPA. Figure 2 shows the mean and median $r^{2}$ scores obtained by
each model on the test set. Compared to the baseline and chemCPA, our model
achieved the highest $r^{2}$ scores on all genes and DEG set. The result
strongly supports that our model achieves the current state-of-the-art level.
To visually demonstrate the predictive ability of our model, we used the UMAP
method to visualize the actual gene expression profiles and predicted gene
expression profiles, as shown in Figure 3 (a). It can be seen that our model’s
predicted values are very close to the actual values, indicating that our
model effectively captures the effect of drug perturbation on gene expression.
Moreover, Figure 4 shows the explained variance of the baseline and our model,
Figure 3 (b) visualize the boxplots of $r^{2}$ scores regrading 10 drugs with
the most perturbed cellular responses in the test set. The boxplots are
grouped by drugs, and each group contains the prediction results of the model
for different cell lines under the same drug perturbation.

(a)  (b)

Figure 3: Performance evaluation on L1000 bulk transctiptional response data.
(a) UMAP visualization of predicted and actual expression profiles of 2,289
test samples. (b) The boxplots of $r^{2}$ scores regrading 10 drugs having the
most number of perturbed cellular responses included in the test set. Figure
4: Explained variance comparison of cycleCDR with baseline model on L1000 bulk
transcriptional response dataset.

### Performance evaluation on proteomic response

Most molecular and targeted drugs achieve their pharmacological effects by
affecting the function of proteins and protein complexes. The CPPA portal
provides a set of large-scale proteomic expression levels measured by RPPA
assays. The CPPA dataset includes 549 clinically relevant protein levels of
126 human cell lines perturbed by 99 drugs, which enable us to evaluate our
model in predicting the proteomic response to drug perturbation.

We discard the data without dosage information and obtained 1,760 drug-cell
line combinations spanning 538 proteins. The processed dataset were randomly
divided into training ($n$=1,408), validation ($n$=246), and test ($n$=106)
sets. As the CPPA data contains the measured proteomic profiles upon different
drug dosage, we added an encoder for drug dosage to convert the drug dosage
into an embedding and perform element-wise multiplication with the drug
embedding. We show in Figure 5 the achieved mean and median $r^{2}$ and EV
metrics on the test set. Compared to the baseline model, our model showed
significant performance superiority.

Figure 5: Performance evaluation on bulk proteomic response dataset (CPPA).

### Performance evaluation on single-cell transcriptional response

The sci-Plex3 single-cell resolution transcriptome response dataset has
measured the transcriptional response of 3 human cancer cell lines to 188
compounds through high-throughput screening. For performance comparison, we
used the sci-Plex3 dataset processed by chemCPA to evaluate the predictive
ability of the model. Figure 6 shows the mean and median $r^{2}$ scores of our
model on the dataset. Compared to both chemCPA with and without pretraining on
L1000 bulk data, our model significantly achieved better performance. As shown
in Figure 7, the UMAP visualization verified that the predicted
transcriptional profiles are very close to the actual ones, indicating that
our model can effectively capture single-cell transcriptional drug response.
Figure 7 shows the obtained $r^{2}$ scores of three major cell lines upon drug
perturbations.

Figure 6: Performance comparison of on sci-Plex3 single-cell transcriptional
response dataset.

(a)  (b)

Figure 7: Performance evaluation on sci-plex3 single-cell transctiptional
response data. (a) UMAP visualization of predicted and actual expression
profiles. (b) The boxplots of $r^{2}$ scores regrading three cell lines.

### Model ablation

We first evaluate the effect of single autoencoder or dual autoencoder
architecture on model performance. Table 1 shows the results of model ablation
experiments performed on the sci-plex3 dataset. Although the dual autoencoder
structure can improve the performance, its improvement is relatively weak. In
addition, we also examined the effects of reconstruction loss and cycle
consistence loss on the model separately, and found that cycle consistence
loss contributes significantly to improve the model, especially for the
prediction of DEGS. Most interestingly, when using the single autoencoder the
combination usage of reconstruction loss and cycle consistence loss actually
reduces the performance of the model, which may be caused by the limited
expressive capacity of single autoencoder.

Table 1: Performance evaluation of ablated models Single AE | Dual AE | Reco Loss | Cycle Loss | mean $r^{2}$ (all) | mean $r^{2}$ (DEGs)
---|---|---|---|---|---
$\checkmark$ |  | $\checkmark$ |  | 0.53 | 0.61
$\checkmark$ |  |  | $\checkmark$ | 0.54 | 0.63
$\checkmark$ |  | $\checkmark$ | $\checkmark$ | 0.58 | 0.67
| $\checkmark$ | $\checkmark$ |  | 0.65 | 0.59
| $\checkmark$ |  | $\checkmark$ | 0.66 | 0.66
| $\checkmark$ | $\checkmark$ | $\checkmark$ | 0.72 | 0.68

## Conclusion

In this paper, we proposed a novel deep network model to learn cellular
response to drug perturbations. We integrated the cycle consistency loss and
linear model in latent space into a end-to-end learning framework, which
enables our model to learn expressive and transferable drug representations.
We evaluate the proposed model on both bulk and single-cell transcriptional
responses, as well as a proteomic drug response dataset. The experimental
results on these benchmark datasets demonstrate the superior performance of
our model.

## References

  * Aissa et al. [2021] Aissa, A. F.; Islam, A. B.; Ariss, M. M.; Go, C. C.; Rader, A. E.; Conrardy, R. D.; Gajda, A. M.; Rubio-Perez, C.; Valyi-Nagy, K.; Pasquinelli, M.; et al. 2021\.  Single-cell transcriptional changes associated with drug tolerance and response to combination therapies in cancer.  _Nature communications_ , 12(1): 1628.
  * Arjovsky, Chintala, and Bottou [2017] Arjovsky, M.; Chintala, S.; and Bottou, L. 2017.  Wasserstein generative adversarial networks.  In _International conference on machine learning_ , 214–223. PMLR.
  * Chawla et al. [2022] Chawla, S.; Rockstroh, A.; Lehman, M.; Ratther, E.; Jain, A.; Anand, A.; Gupta, A.; Bhattacharya, N.; Poonia, S.; Rai, P.; et al. 2022.  Gene expression based inference of cancer drug sensitivity.  _Nature communications_ , 13(1): 5680.
  * Chen et al. [2019] Chen, C.; Chen, Z.; Jiang, B.; and Jin, X. 2019.  Joint domain alignment and discriminative feature learning for unsupervised deep domain adaptation.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 33, 3296–3303.
  * Chen et al. [2022] Chen, J.; Wang, X.; Ma, A.; Wang, Q.-E.; Liu, B.; Li, L.; Xu, D.; and Ma, Q. 2022\.  Deep transfer learning of cancer drug responses by integrating bulk and single-cell RNA-seq data.  _Nature Communications_ , 13(1): 6494.
  * Gao, Yang, and Wang [2021] Gao, L.-G.; Yang, M.-Y.; and Wang, J.-X. 2021.  Collaborative matrix factorization with soft regularization for drug-target interaction prediction.  _Journal of Computer Science and Technology_ , 36: 310–322.
  * He et al. [2022] He, D.; Liu, Q.; Wu, Y.; and Xie, L. 2022.  A context-aware deconfounding autoencoder for robust prediction of personalized clinical drug response from cell-line compound screening.  _Nature Machine Intelligence_ , 4(10): 879–892.
  * Hetzel et al. [2022] Hetzel, L.; Boehm, S.; Kilbertus, N.; Günnemann, S.; Theis, F.; et al. 2022\.  Predicting cellular responses to novel drug perturbations at a single-cell resolution.  _Advances in Neural Information Processing Systems_ , 35: 26711–26722.
  * Hetzel et al. [2021] Hetzel, L.; Fischer, D. S.; Günnemann, S.; and Theis, F. J. 2021.  Graph representation learning for single-cell biology.  _Current Opinion in Systems Biology_ , 28: 100347.
  * Jia et al. [2021] Jia, P.; Hu, R.; Pei, G.; Dai, Y.; Wang, Y.-Y.; and Zhao, Z. 2021.  Deep generative neural network for accurate drug response imputation.  _Nature Communications_ , 12(1): 1740.
  * Liu et al. [2022] Liu, H.; Huang, Y.; Liu, X.; and Deng, L. 2022.  Attention-wise masked graph contrastive learning for predicting molecular property.  _Briefings in Bioinformatics_ , 23(5): bbac303.
  * Long et al. [2018] Long, M.; Cao, Z.; Wang, J.; and Jordan, M. I. 2018.  Conditional adversarial domain adaptation.  _Advances in neural information processing systems_ , 31.
  * Lopez, Gayoso, and Yosef [2020] Lopez, R.; Gayoso, A.; and Yosef, N. 2020.  Enhancing scientific discoveries in molecular biology with deep generative models.  _Molecular Systems Biology_ , 16(9): e9198.
  * Lotfollahi et al. [2021] Lotfollahi, M.; Susmelj, A. K.; De Donno, C.; Ji, Y.; Ibarra, I. L.; Wolf, F. A.; Yakubova, N.; Theis, F. J.; and Lopez-Paz, D. 2021.  Learning interpretable cellular responses to complex perturbations in high-throughput screens.  _BioRxiv_ , 2021–04.
  * Lotfollahi, Wolf, and Theis [2019] Lotfollahi, M.; Wolf, F. A.; and Theis, F. J. 2019.  scGen predicts single-cell perturbation responses.  _Nature methods_ , 16(8): 715–721.
  * Ma et al. [2021] Ma, J.; Fong, S. H.; Luo, Y.; Bakkenist, C. J.; Shen, J. P.; Mourragui, S.; Wessels, L. F.; Hafner, M.; Sharan, R.; Peng, J.; et al. 2021.  Few-shot learning creates predictive models of drug response that translate from high-throughput screens to individual patients.  _Nature Cancer_ , 2(2): 233–244.
  * Madhukar et al. [2019] Madhukar, N. S.; Khade, P. K.; Huang, L.; Gayvert, K.; Galletti, G.; Stogniew, M.; Allen, J. E.; Giannakakou, P.; and Elemento, O. 2019.  A Bayesian machine learning approach for drug target identification using diverse data types.  _Nature communications_ , 10(1): 5221.
  * Pei et al. [2018] Pei, Z.; Cao, Z.; Long, M.; and Wang, J. 2018.  Multi-adversarial domain adaptation.  In _Thirty-second AAAI conference on artificial intelligence_.
  * Robichaux et al. [2021] Robichaux, J. P.; Le, X.; Vijayan, R.; Hicks, J. K.; Heeke, S.; Elamin, Y. Y.; Lin, H. Y.; Udagawa, H.; Skoulidis, F.; Tran, H.; et al. 2021.  Structure-based classification predicts drug response in EGFR-mutant NSCLC.  _Nature_ , 597(7878): 732–737.
  * Roohani, Huang, and Leskovec [2022] Roohani, Y.; Huang, K.; and Leskovec, J. 2022.  GEARS: Predicting transcriptional outcomes of novel multi-gene perturbations.  _BioRxiv_ , 2022–07.
  * Sankaranarayanan et al. [2018] Sankaranarayanan, S.; Balaji, Y.; Castillo, C. D.; and Chellappa, R. 2018.  Generate to adapt: Aligning domains using generative adversarial networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 8503–8512.
  * Sharifi-Noghabi et al. [2021] Sharifi-Noghabi, H.; Harjandi, P. A.; Zolotareva, O.; Collins, C. C.; and Ester, M. 2021.  Out-of-distribution generalization from labelled and unlabelled gene expression data for drug response prediction.  _Nature Machine Intelligence_ , 3(11): 962–972.
  * Srivatsan et al. [2020] Srivatsan, S. R.; McFaline-Figueroa, J. L.; Ramani, V.; Saunders, L.; Cao, J.; Packer, J.; Pliner, H. A.; Jackson, D. L.; Daza, R. M.; Christiansen, L.; et al. 2020.  Massively multiplex chemical transcriptomics at single-cell resolution.  _Science_ , 367(6473): 45–51.
  * Subramanian et al. [2017] Subramanian, A.; Narayan, R.; Corsello, S. M.; Peck, D. D.; Natoli, T. E.; Lu, X.; Gould, J.; Davis, J. F.; Tubelli, A. A.; Asiedu, J. K.; et al. 2017.  A next generation connectivity map: L1000 platform and the first 1,000,000 profiles.  _Cell_ , 171(6): 1437–1452.
  * Tzeng et al. [2019] Tzeng, E.; Hoffman, J.; Zhang, N.; Saenko, K.; and Darrell, T. 2019.  Deep domain confusion: Maximizing for domain invariance. arXiv 2014.  _arXiv preprint arXiv:1412.3474_.
  * Wang, He, and Katabi [2020] Wang, H.; He, H.; and Katabi, D. 2020.  Continuously indexed domain adaptation.  _arXiv preprint arXiv:2007.01807_.
  * Yu and Welch [2022] Yu, H.; and Welch, J. D. 2022.  PerturbNet predicts single-cell responses to unseen chemical and genetic perturbations.  _bioRxiv_ , 2022–07.
  * Zhang, Lu, and Zang [2022] Zhang, C.; Lu, Y.; and Zang, T. 2022.  CNN-DDI: a learning-based method for predicting drug–drug interactions using convolution neural networks.  _BMC bioinformatics_ , 23(1): 1–12.
  * Zhao et al. [2020] Zhao, W.; Li, J.; Chen, M.-J. M.; Luo, Y.; Ju, Z.; Nesser, N. K.; Johnson-Camacho, K.; Boniface, C. T.; Lawrence, Y.; Pande, N. T.; et al. 2020\.  Large-scale characterization of drug responses of clinically relevant proteins in cancer cell lines.  _Cancer cell_ , 38(6): 829–843.
  * Zheng et al. [2017] Zheng, H.; Fu, J.; Mei, T.; and Luo, J. 2017.  Learning multi-attention convolutional neural network for fine-grained image recognition.  In _Proceedings of the IEEE international conference on computer vision_ , 5209–5217.
  * Zhu et al. [2017] Zhu, J.-Y.; Park, T.; Isola, P.; and Efros, A. A. 2017.  Unpaired image-to-image translation using cycle-consistent adversarial networks.  In _Proceedings of the IEEE international conference on computer vision_ , 2223–2232.


# BoGraph: Structured Bayesian Optimization From Logs for Expensive Systems
with Many Parameters

Sami Alabed  0000-0001-8716-526X University of Cambridge, UK
The Alan Turing Institute<EMAIL_ADDRESS>and  Eiko Yoneki  University of
Cambridge, UK<EMAIL_ADDRESS>
###### Abstract.

Current auto-tuners struggle with computer systems due to their large complex
parameter space and high evaluation cost. We propose BoGraph, an auto-tuning
framework that builds a graph of the system components before optimizing it
using causal structure learning. The graph contextualizes the system via
decomposition of the parameter space for faster convergence and handling of
many parameters. Furthermore, BoGraph exposes an API to encode experts’
knowledge of the system via performance models and a known dependency
structure of the components. We evaluated BoGraph via a hardware design case
study achieving $5x-7x$ improvement in energy and latency over the default in
a variety of tasks.

##  1\. Introduction

Computer systems facilitate the execution of a wide range of diverse workloads
by exposing tunable configurations to meet users’ demands. Examples of these
configurations are hardware designs (Gem5Simulatorbinkert2011, ) and database
knobs (InquiryMachineLearningbasedvanaken2021, ). Tuning systems configuration
is tedious due to the many parameters involved and the long evaluation time.
For example, in our case study on hardware design optimization, there are
$2^{64}$ unique designs. This complexity necessitates for auto-tuners to aid
practitioners. However, current research into auto-tuners focuses on
optimizing neural network hyperparameters (MicrosoftNniopensource2020, ;
RandomSearchHyperparameterbergstra2012, ;
PracticalBayesianOptimizationsnoek2012, ) and encounters significant
challenges when applied to tuning computer systems specifically.

###  1.1. BoGraph

BoGraph builds on Structured Bayesian Optimization (SBO) (Dalibard2017, )
principles to contextualize the system and finds optimal system configurations
in a few trials. SBO requires experts to encode their system knowledge using
hand-crafted probabilistic models. BoGraph differentiates itself from previous
works by simplifying the process of using SBO. It discovers probabilistic
models from logs using its novel structure discovery pipeline. BoGraph’s API
allows experts to encode contextual knowledge of the system, such as models or
dependencies between components, for a more robust and interpretable auto-
tuner. Next, we summarize the challenges BoGraph targets in optimizing
computer systems.

###  1.2. Challenges in systems optimization

####  1.2.1. Evaluation is slow and expensive

Evaluating a single configuration on a system takes several minutes (see 6(a))
to several hours (Floratou2017, ; InquiryMachineLearningbasedvanaken2021, ;
rlDevicePlacement, ). This is costly especially when optimizing cloud
instances (Alipourfard2017, ) or databases
(InquiryMachineLearningbasedvanaken2021, ). Hence, the need to find optimized
configurations with the fewest evaluations. As a result, tuners that rely on
many iterations are too costly to use, e.g., reinforcement learning
(ReinforcementLearningIntroductionsutton1998, ), random search
(RandomSearchHyperparameterbergstra2012, ), evolutionary search (Ansel2010, ),
hill-climbing (Ansel2014, ), or population-search
(PopulationBasedTrainingjaderberg2017, ). BoGraph leverages all available
information from the system logs and captures it using Gaussian Process (GP)
(GaussianProcessesMachinerasmussen2008, ) a sample-efficient model to
contextualize the system. Additionally, it allows incorporating existing
experts’ knowledge to speed up the process further.

####  1.2.2. Large number of parameters

Using GPs comes at the cost of increased computational complexity and
difficulty to scale to $D>10$ parameters without sacrificing their efficiency
(WhenGaussianProcessliu2019, ). Standard Bayesian Optimization
(PracticalBayesianOptimizationsnoek2012, ) methods often rely on GPs, which do
not scale to many parameters. BoGraph enables the GPs to scale by using the
system’s dependency graph to naturally reduce each model dimension as
illustrated in Figure 1.

(a) The model has to consider all the parameters jointly.

(b) Decomposition reduces the effective dimensions naturally.

Figure 1. The structure reduces the effective dimensions.

The Figure shows the system decomposability. On the left is a standard
optimization approach that considers all parameters jointly, with a max of
five dimensions. On the right, a decomposed method breaks down optimization
along the available metrics line, leading to fewer dimensions considered for
each optimization target.

####  1.2.3. Complex parameters dependency

Computer systems parameters often have a complex relationship; a parameter can
depend on another parameter, e.g. the memory’s structure and size. This
relationship is non-linear and difficult to learn using traditional methods.
BoGraph exposes an API for experts to encode a parametric probabilistic model
that captures parameters’ complex relationship.

####  1.2.4. Designing bespoke models is complicated

Hand-crafting models that capture the system’s complexity requires system and
machine learning expertise. Furthermore, validating and evaluating the models
is a trial-and-error process. Making tuners that rely on expert models such as
BOAT (Dalibard2017, ) and Causal Bayesian Optimization
(CausalBayesianOptimizationaglietti2020, ) difficult to use. BoGraph allows
encoding models if they exist, and it provides a pipeline that learns these
models if necessary. It captures statistical dependency between parameters and
system metrics then uses GPs to model their interactions.

###  1.3. Contributions

We utilized BoGraph to optimize the design space ($\approx 2^{64}$) of an
accelerator simulator, gem5-Aladdin (CodesigningAcceleratorsSocshao2016, ;
Gem5AladdinSoCSimulatorharvard-acc2016, ). BoGraph learned the system’s
components dependencies from logs, then modeled them using a graph of
probabilistic models. BoGraph found configurations in twenty iterations that
took the next best auto-tuner a hundred iterations. Moreover, while improving
the energy-latency utilization by $5-7x$ over the default, no other method
came close to its performance.

We summarize our contributions as follows:

  * •

An auto-tuner for computer systems that handle many parameters and find
optimal configurations quickly.

  * •

A SBO framework automatically learns a graph of probabilistic models from the
system’s logs.

  * •

An API for experts to express the system’s dependencies or encode
probabilistic models of its components.

  * •

Optimized all the accelerator design choices in gem5-Aladdin and improved
energy-latency by $5-7$x.

##  2\. Background

###  2.1. Bayesian optimization

Bayesian Optimization (BO) (TakingHumanOutshahriari2015, ) is an iterative
sample-efficient method to optimize configurations. Formulated as
$x^{*}=argmin(\alpha(\mathbb{M}(f))$, where $x^{*}$ is the optimal
configuration, $\alpha$ is the acquisition function, $\mathbb{M}$ is a
probabilistic model of the system, and $f$ is the system. Its efficiency comes
from generating samples using $\mathbb{M}$ that $\alpha$ optimizes to find
configurations that improve the system without direct interaction with it. The
choice of $\mathbb{M}$ and $\alpha$ directly impact the required evaluations
to find optimal configurations.

###  2.2. Gaussian Process

Gaussian Process (GP) (GaussianProcessesMachinerasmussen2008, ) is a flexible
sample-efficient model commonly used in BO. GP captures the impact of the
configurations against the objective in the mean function and co-variance in
its kernel. In practice, GP struggles with problems that have $D>10$
(BayesianOptimizationBillionwang2016, ) as it scales cubically to the data and
dimensions. Its runtime requires an expensive matrix inversion operation for
training and sampling, hence its scalability issues. The computational
complexity of the GP renders it ineffective for high-dimensional problems.

###  2.3. Probabilistic DAG

BoGraph models the system using a directed-acyclic-graph (DAG), where the
nodes are probabilistic models that capture a component of the system, and the
edges are the conditional dependency between the models. This representation,
a DAG of probabilistic models, is known as Bayesian Network
(ProbabilisticGraphicalModelskoller2009, ). Bayesian Network has the advantage
of factorizing the joint distribution of a node into a local distribution that
depends only on the parent: $P(X;G)=\prod_{j=i}^{p}P(X_{i}|X_{parent(i)})$
where $X_{i}$ is the ith node in the graph and parent(i) is the direct parent
of the node. Nodes factorization circumvents the curse of dimensionality and
allows each node to scale independently.

##  3\.  BoGraph Framework

Figure 2.  BoGraph’s optimization loop. The expert provides specifications of
the problem. Then, BoGraph learns a system structure from the logs and
combines expert knowledge to build a probabilistic model. Finally, it
generates samples from the model and optimizes them to propose optimized
configurations to the system.

BoGraph is an auto-tuning framework for systems with many configurations. It
builds on Structured Bayesian Optimization (SBO) (Dalibard2017, ) by utilizing
contextual knowledge to converge quickly. BoGraph novelty is to allow a hybrid
method for injecting expert knowledge into the system as well as learning
system structure from logs. Using its pipeline, it identifies statistical
links between the system components and uses them to build a system’s
probabilistic DAG.

Figure 2 shows BoGraph’s main components. The expert encodes any knowledge of
the system and the target objectives. Then, BoGraph scans the system logs to
learn the statistical dependency of the system’s components. Next, it combines
the expert’s knowledge with the learned dependency producing a probabilistic
DAG. Finally, it suggests configurations to evaluate by generating many
samples from the model and optimizing them using an acquisition function. Its
modular design facilitates the use of any sampler and acquisition functions,
with BoTorch (BoTorchFrameworkEfficientbalandat2020, ) being the default.

###  3.1. Structure discovery from logs

Figure 3.  A) A snippet of gem5’s (the environment) logs. B) The logs after
summarizing and grouping. C) The structure resulted by calculating a
statistical correlation between components (respecting parameters-component
causality) and applying the expert-defined edges (the entire graph is too
large for the paper).

BoGraph’s first component is the structure discovery process shown in Figure 3
that learns a dependency DAG from the system logs. The pipeline first
summarizes the logs, then discovers the statistical dependency using structure
learning (ProbabilisticGraphicalModelskoller2009, ; DagsNoTearszheng2018, ),
and finally connects any expert-defined edges to produce a dependency DAG of
the system components.

#### Summarizing the logs

Computer systems generate logs used for monitoring and debugging by experts
(SiteReliabilityEngineeringbeyer2016, ), BoGraph uses the logs to learn
components dependencies. However, the logs need to be summarized first as
there are often many log entries causing the curse of dimensionality. For
example, in the case study described in section 4, gem5-Aladdin reports more
than $1000$ log entry (Figure 3A), BoGraph pipeline summarizes them to $7$
entries (Figure 3B).

#### Preprocessing

First, BoGraph removes uninformative metrics by pruning ones with low variance
between executions. Then, they are standardized to ensure that different
measuring units do not skew the learning algorithm. Then, BoGraph exploits the
structured property of logs to summarize them. The logging protocol encodes
where the log was generated in the system
(CharacterizingLoggingPracticesyuan2012, ; WhereDevelopersLogfu2014, ). For
example, in the case study on gem5, the structure is
Component.SubComponent.Measurement, as shown in Figure 3A, BoGraph groups all
measurements under SubComponent in one group as they influence the same
component in the system, then they are compressed into a single statistic
using Factor Analysis (FA) (EasyGuideFactorkline2014, ). This compression
works because logs generated at the SubComponent level are dependent on each
other. The user can decide on the granularity of the grouping based on their
need. A fine granularity is helpful to understand the system, but it is more
expensive to evaluate. We used the highest level of granularity which mapped
gem5 to six main categories: L2, Datapath, Mem Ctrls, Membus, Tol2Bus, and
CPU. This process applies to all system logs, including previous system
evaluations. At the end of this stage, BoGraph produces a summary of the logs
as seen in Figure 3B for the structure learning algorithm.

###  3.2. Structure learning

Defining systems’ structure is an error-prone process; hence the need to learn
it from data, an active research area in the Bayesian community
(SurveyBayesianNetworkscanagatta2019, ). These methods propose several DAGs
and score them based on a trade-off between data fit and the complexity of the
graph. BoGraph uses a structure learning algorithm, NoTears
(DagsNoTearszheng2018, ). NoTears approximates the DAG score using a
differentiable function, then minimizes its error, which is easily
parallelizable on GPUs (Causalnex2021, ). It is simple to extend BoGraph with
other structure learning algorithms by implementing the BoGraph’s
StructureLearning interface. We encode a set of restrictions in these
algorithms that ensures the system’s parameters are not allowed to have
parents, and the objectives nodes are not allowed to have children to ensure
the causal flow of the optimization task. Then apply any expert-defined edges
resulting in Figure 3C. Optionally, BoGraph takes an update scheduler that
decides when to perform structure learning. By default, it performs structure
learning every quarter of the optimization budget.

###  3.3. Graph of a Probabilistic Models

Figure 4.  BoGraph maps the node in the structure to probabilistic models such
as GPs or expert’s models, creating probabilistic DAG.

BoGraph builds the probabilistic DAG by placing a GP node on intermediate
nodes, and then it connects any expert-defined models as illustrated in Figure
4. The probabilistic DAG nodes’ are compartmentalized, where each node
predicts a single value conditioned on its parents directly and can be trained
independently. As a result, any node can be an optimization target, providing
free multi-objective modeling and enabling a parallelizable training. We will
investigate both in a follow-up project.

Algorithm 1 Sampling the probabilistic DAG

Input: sampler $s$, DAG model $M$, system $SY$

initialize a samples cache: $Q:=Dict()$

get nodes along the path to the objectives: $P=M.path(SY.obj)$

for node $n$ in $P.nodes$ do

if $n\in SY.params$ then

sample parameters: $p_{samples}=s.sample(n.space())$

cache generated samples: $Q\cup(n.key,p_{samples})$

else

condition on parents: $C=n.predict(Q[n.parents])$

cache node results conditioned on parents: $Q\cup(n.key,C)$

end if

end for

Return: $Q[SY.objectives]$

Sampling process. Allowing general probabilistic models comes at the cost that
closed-form acquisition function optimization is not possible; instead,
BoGraph samples the model and then optimizes over the samples. This design
gives it the freedom to use any probabilistic model inside the graph, such as
Bayesian Neural Networks (BayesianNeuralNetworksmackay1995, ). Algorithm 1
shows the process to generate a prediction and samples from a probabilistic
DAG. First, BoGraph finds all the nodes with a direct path from the parameters
to the objective. Then, it creates a cache to hold all the samples and avoid
recomputing them for nodes with multiple children. Next, each node condition
on the parents’ prediction of the samples, and finally output the posterior of
the objectives.

###  3.4. Structured Bayesian Optimization

Algorithm 2 BoGraph Bayesian Optimization loop

Input: acquisition function $\alpha$, system $E$, budget

repeat

build DAG model from logs: $DAG_{M}=BoGraph{}{}(E.logs())$

select next to point to evaluate: $x_{n+1}=\alpha(s.sample(DAG_{M}))$

submit $x_{n+1}$ to be evaluated: $E.update(x_{n+1})$

until $budget=0$

Optimization loop. Algorithm 2 details the optimization loop. BoGraph’s
probabilistic DAG combines several models’ posterior, producing a complicated
posterior that cannot be optimized directly. Instead, we use quasi-acquisition
functions that optimize the estimate from samples rather than the model
directly. BoTorch (BoTorchFrameworkEfficientbalandat2020, ) provides a suite
of utility functions for BO frameworks, including quasi-acquisition functions
and samplers. We use BoTorch’s Sobol sequence sampler to generate diverse
samples of the posterior, then use the quasi-Expected Improvement acquisition
function to optimize the estimates and produce optimization candidates.

###  3.5.  BoGraph API

BoGraph API streamlines the process of expressing system designers’ knowledge.
Listing 5 shows several of the functions BoGraph that enable experts to encode
their knowledge in two ways: The first is a DAG API based on networkX
(ExploringNetworkStructurehagberg2008, ) to express the dependencies between
parameters and metrics supporting many methods to store and load DAGs. The
second is an advanced API that expresses experts’ knowledge using Pyro
(PyroDeepUniversalbingham2018, ), a probabilistic programming language. By
default, BoGraph uses a GP from GPyTorch
(GPyTorchBlackboxMatrixMatrixgardner2019, ) a performant GP implementation as
the probabilistic model if no model is defined. The user can override the
kernel of the GP to inject any prior information about the objective function
distribution; this is only recommended for advanced users.

⬇

\parclass ExpertModel(|\textcolor{red}{PyroNode}|): # User extends PyroNode to
build models

def model(x1: Tensor, obs: Tensor = None):

i = pyro.sample(”intercept”, dist.Uniform(0., 5.))

linear_mean = i + pyro.sample(”bias_x1”, dist.Normal(0.,1.))

… # a standard Bayesian model using Pyro

prior_dag = | \BoGraph{}{}|Prior() # | \BoGraph{}{}| uses ‘networkx‘ API

prior_dag.add_edges_from([’x1’, ’f1(x1)’]) # f1(x1) depends on x1

prior_dag.add_model(”f1(x1)”, ExpertModel()) # add expert’s model

prior_dag.add_model(”y”, GPyTorchModel()) # standard GP model

prior_dag.add_tabu_edge(”x1”, ”y”) # Edges known to not be correlated

|\BoGraph{}{}|.optimize(objective=Forrester2D, prior=prior_dag)

Figure 5.  BoGraph API showing expert knowledge injection

###  3.6. Summary

BoGraph simplifies the process of integrating SBO in systems for faster
optimizer convergence. The main loop performs structure learning of the
system’s components by processing the logs and transforming it into a
probabilistic model. BoGraph exposes an API that allows experts to encode
dependency structure and probabilistic models.

##  4\. Evaluation

System-on-a-chip (SoC) designers often experiment on a simulator
(Gem5Simulatorbinkert2011, ; AladdinPrertlPowerperformanceshao2014, ;
EyerissV2Flexiblechen2019, ) to reduce the energy and latency of integrated
chips (ExploringExplorationTutorialpimentel2016, ;
DeterminingOptimalCoherencybhardwaj2019, ). We optimized Energy-Delay Product
(EDP)(EnergyperformanceTradeoffsProcessorazizi2010, ) of gem5-Aladdin
(CodesigningAcceleratorsSocshao2016, ) where
$EDP=energy*(\frac{1}{sim\\_seconds})^{2}$ by tuning all the design parameters
in Table 1.

Table 1.  The gem5-aladdin parameters from their sweeper (Gem5AladdinSoCSimulatorharvard-acc2021, ). The $\approx 2^{64}$ unique design is enormous hence the need for BoGraph. The parameters spaces in red are pruned using a power function.  Parameter | Min | Max
---|---|---
Cycle Time | 1 | 10
Pipelining | 0 | 1
Cache Size (power of two) | 14 | 17
Cache Assoc (power of two) | 0 | 5
Cache Line Sz(power of two) | 5 | 7
Cache Hit Latency | 1 | 5
Cache Queue Size | 16 | 256
Cache Bandwidth | 2 | 18
Tlb hit latency | 1 | 32
Tlb miss latency | 1 | 32
Tlb entries (power of two) | 0 | 5
Tlb assoc (power of two) | 0 | 5
Tlb page size | 1024 | 32768
Tlb max outstanding walks | 2 | 16
Tlb bandwidth | 2 | 16
L2cache Size (power of two) | 17 | 21
Enable L2 | 0 | 1
Pipelined Dma | 0 | 1
Ready Mode | 0 | 1
Ignore Cache Flush | 0 | 1

#### Evaluation goals

We evaluated the following claims:

  * •

Executing BoGraph is faster than a full system evaluation, 4.1.

  * •

BoGraph converges towards optimal configurations faster than other state-of-
the-art methods, 4.2.

  * •

BoGraph produces interpretable system structures that reduces the maximum
dimension, 4.3.

#### Setup

The experiment ran on RTX3060TI, Intel i7 8700k, and 32GB RAM. BoGraph used
GPs from GPyTorch v1.2 (GPyTorchBlackboxMatrixMatrixgardner2019, ), samplers
and acquistion functions from BoTorch v0.5
(BoTorchFrameworkEfficientbalandat2020, ), and CausalNex v0.11 (Causalnex2021,
) for structure learning. The gem5-Aladdin simulator was compiled using MESI
Two Level Aladdin protocol with MathExprPowerModel as power module.

#### Benchmark

Used a mix of data and compute workloads from MachSuite
(MachsuiteBenchmarksAcceleratorreagen2014, ) as is common in SoC literature
(CodesigningAcceleratorsSocshao2016, ;
ExploringExplorationTutorialpimentel2016, ).

#### Expert knowledge.

The EDP formula.

#### Baselines

We used these baselines in the experiments:

  * •

Default (Gem5AladdinSoCSimulatorharvard-acc2016, ): the settings provided by
the simulator.

  * •

Random (RandomSearchHyperparameterbergstra2012, ; MicrosoftNniopensource2020,
): useful when small subset of the dimensions significantly impact the
objective.

  * •

PBTTuner (PopulationBasedTrainingjaderberg2017a, ; MicrosoftNniopensource2020,
): a bag of model approach to finding a collection of optimal configurations.

  * •

BoTorch (BoTorchFrameworkEfficientbalandat2020, ): shows the effectiveness of
BoGraph as they use the same samplers and acquisition functions.

  * •

DeepGP (DeepGaussianProcessesdamianou2013, ): Provides unsupervised structure
discovery in the latent space. We configured a two-layer GPs with 64 inducing
points (Our GPU’s max VRAM).

###  4.1. Execution time

We compared the execution time of the simulator for each task and compared it
to the overhead of using an auto-tuner in Figure 6. The results show that the
optimizers are cheaper than a full system evaluation. Hence the need for auto-
tuners that perform better than simple parameter sweeps.

(a) Tasks’ average execution time, with variance from $1800$ runs.

(b) Auto-tuners’ average time to train and propose a configuration, with
variances from $300$ runs.

Figure 6.  Optimizers’ execution time is dominated by the gem5-Aladdin task
execution time.

###  4.2. EDP Optimization

Figure 7.  The best found EDP in $100$ steps. The y-axis shows X-factors
improvement over default settings and the variance of three runs.

Figure 7 examines the improvement in the EDP from the default configurations
across all tasks. BoGraph outperformed every other method and improved over
the default by $5x-7x$ factors. Examining that further, Table 2 shows that
BoGraph found configurations with the least power consumption for the ELLPACK
task. Other tuners are stuck in local minima that do not improve the energy as
they are unaware the EDP is impacted by power. BoGraph is aware that the
latency and power impact EDP thus can avoid the local minima. The same story
applies to the other tasks in the benchmark.

Table 2.  ELLPACK’s performance for the median best configuration proposed by the models, lower values are better.  model | log EDP | Latency | Power (mWatts)
---|---|---|---
BoGraph | ${\mathbf{8.12\pm 0.28}}$ | $\mathbf{11.06\pm 0.19}$ | $\mathbf{27.37\pm 11.59}$
DeepGP | $10.14\pm 0.27$ | $11.06\pm 3.92$ | $203.15\pm 46.34$
PBTTuner | $10.34\pm 0.60$ | $10.81\pm 4.46$ | $264.36\pm 147.79$
BoTorch | $10.44\pm 0.62$ | $11.06\pm 3.90$ | $108.46\pm 102.97$
Random | $10.74\pm 0.19$ | $18.32\pm 4.81$ | $143.94\pm 194.83$
Default | $14.50\pm 0$ | $94.80\pm 0$ | $220.34\pm 0$

#### Convergence

Figure 8.  Auto-tuners’ improvement over default, with a median best
configuration and the error bars showing a minimum of 3 runs.

Figure 8 examines the convergence rate of each optimizer as minimizing the
evaluations reduces the cost significantly (Challenge 1.2.1). BoGraph reached
the best configuration the quickest and continued to improve while other
methods plateaued. Confirming that tuning computer systems is challenging for
most auto-tuners. BoGraph improvement jumps at every 20th step in Figure 8 are
caused by the structure learning algorithm, as its update frequency defaults
to a quarter of the evaluation budget.

###  4.3. Structure discovery

Figure 9.  BoGraph DAG max-dimension (lower is better) of three repeated runs.
The max dimension is defined as $MD=\forall_{n\in G}max(indegree(n))$, each
node becomes a probabilistic model.

There are over $1000$ metrics reported in a gem5 statistics file. BoGraph
shields the end-user from the difficult task of defining a probabilistic model
while providing an insightful structure to the designer as shown in Figure 3C.
Furthermore, Figure 9 shows that learning structure does reduce the maximum
dimension for the models in the graph.

##  5\. Oppurtinities and limitation.

###  5.1. Oppurtinities

BoGraph is still evolving, with aims to apply to different problem spaces,
evaluate it on a multi-objective problem, and improve its computational
engine.

#### Multi-objective optimization.

Multi-objective modeling is an interesting opportunity for computer system
objectives often compete (latency vs. throughput, energy vs. latency).
Modeling multiple objectives is trivial in BoGraph as a result of using DAGs
as an internal data structure as any node can be sampled. Then using off-the-
shelf multi-objective utility functions such as expected hypervolume
improvement function (DifferentiableExpectedHypervolumedaulton2020, ) to
optimize the samples from multiple nodes.

#### Computational overhead.

The overhead of using structured DAG is high compared to the standard
approaches. During training and sampling, BoGraph performs sequential
operations without exploiting the independence between nodes to batch and
distribute the workload. Leveraging parallelization will reduce the
computational time of BoGraph.

#### Application domains.

We have shown BoGraph ability to find meaningful structure from logs and use
that to optimize gem5-Aladdin. Other computer systems will benefit from
BoGraph ability to make sense of logs in aiding the optimizer. We are working
on a case study on optimizing the latency of PostgreSQL
(PostgreSQLIntroductionConceptsmomjian2001, ), BoGraph is using the reported
370 metrics to find links between the metrics and the latency objective and
tuning of thirty parameters.

###  5.2. Limitations

#### BoGraph’s pipeline overhead

BoGraph maintains the system trace (previous evaluation, configurations, and
logs) on disk, reading and writing these traces have both a memory and
execution time overhead. This is a design choice we took for the following
reasons: First, the cost of the real system dominates the execution of the
pipeline (subsection 4.1). Secondly, this enables BoGraph in the future to
learn from distributed instances. Finally, it provides a fault-tolerance
mechanism for when the instance has to restart.

Mitigation. Every stage and component of the BoGraph pipeline is configurable.
For example, the user can decide to reduce the frequency of a full causal
structure discovery or include preprocessing steps that reduce the computation
further by extending the preprocessing interface. This flexibility offers some
mitigation to the pipeline overhead. Furthermore, to mitigate this overhead,
BoGraph can run a helper instance that runs the pipeline and communicate only
the summarized information to the BoGraph main loop.

#### Log annotation

A key contribution of BoGraph is providing a pipeline to make use of system
logs and metrics to provide additional context that aids the system model.
Relying on the user to annotate their logs is one limitation of this work. The
annotation process involves providing a parser of the metrics or system logs;
in gem5 and PostgresSQL experiments, it is a one-line regex expression. Logs
in computer systems follow a standard practice
(CharacterizingLoggingPracticesyuan2012, ) as they are often produced to be
ingested into metric monitoring applications such as Prometheus
(PrometheusRunningInfrastructurebrazil2018, ). However, suppose the system is
very different where the logs do not correlate with the system. In that case,
it is possible to disable the structure learning pipeline in BoGraph and rely
on a manual structure, resulting in a similar performance to the standard
structure Bayesian Optimization method (Dalibard2017, ).

##  6\. Related work

Search-based auto-tuners. Reinforcement learning
(ReinforcementLearningIntroductionsutton1998, ), random search
(RandomSearchHyperparameterbergstra2012, ), evolutionary search (in
PetaBricks) (Ansel2010, ), hill-climbing (in OpenTuner) (Ansel2014, ), or
population-search (PopulationBasedTrainingjaderberg2017, ) are simple auto-
tuner to use and scale to many parameters. However, they require many system
evaluations, waste valuable expert knowledge, and are highly dependent on
initializing seed (unstable) (DeepReinforcementLearninghenderson2017, ). These
drawbacks render them unsuitable for our problem due to Challenge 1.2.1.

Dimensionality reduction. OtterTune (VanAken2017, ) used Factor Analysis
(EasyGuideFactorkline2014, ) to reduce the system dimensions. However, as the
OtterTune experiment expanded to a large system
(InquiryMachineLearningbasedvanaken2021, ), it made apparent that computer
system parameters have complex interactions that standard dimensionality
reduction methods fail to capture (Challenge 1.2.3). Principled Component
Analysis (PrincipalComponentAnalysiswold1987, ) suffer from a similar problem.

BO surrogate models. RandomForest
(ClassificationRegressionRandomForestliaw2002, ) used in SMAC
(EvaluationSequentialModelbasedhutter2013, ) and HyperMapper
(PracticalDesignSpacenardi2019a, ) scales to many parameters. Unfortunately,
non-Bayesian methods underestimate the models’ variance
(TakingHumanOutshahriari2015, ) meaning they explore less frequently and miss
on finding the optimal configurations. DeepGPs
(DeepGaussianProcessesdamianou2013, ) use several layers of GPs that reduce
the objective dimensionality by finding a latent structure of it. However, it
requires many experiments to find a meaningful latent structure as it ignores
existing expert knowledge.

Parallel based auto-tuners. AutoTVM (LearningOptimizeTensorchen2018, ),
Hyperband (HyperbandNovelBanditbasedli2017, ) and BOHB
(BOHBRobustEfficientfalkner2018, ) are techniques that depend on the
environment being cheap to terminate and easily parallelizable. Computer
systems are expensive and slow to evaluate; engineering them to launch many
parallel instances is non-trivial without impacting the measurement
sensitivity.

Experts-driven methods. Ernest (ErnestEfficientPerformancevenkataraman2016, ),
and RubberBand (RubberBandCloudbasedHyperparametermisra2021, ) use hand-
crafted models to measure the proposed configurations’ effectiveness quickly
leading to quick optimization. While Causal Bayesian Optimization
(CausalBayesianOptimizationaglietti2020, ) assumes that a full causal graph of
the system exists to speed up finding optimal system configurations.

##  7\. Conclusion

BoGraph is a framework for optimizing systems with many parameters and slow
execution. BoGraph contextualizes the system through learning association
between its metrics and parameters. Furthermore, it provides an API to encode
experts’ knowledge as probabilistic models. The contextualization leads to
tuning more parameters in fewer experiments. Using BoGraph we optimized the
design of an accelerator to improve its energy-latency objective by $2x$ in
${1}/{4}$ of the evaluations of the next best method. For future work, we are
applying BoGraph to PostgreSQL, showing applicability to a range of computer
systems.

###### Acknowledgements.

We would like to thank the reviewers for their valuable feedback. A special
thanks to Thomas Vanderstichele, Mihai Bujanca, Yomna El-Serafy, Ibrahim
Abdou, and Rana ElRashidy for their comments that improved the readability of
the paper. This research was supported by the Alan Turing Institute.

## References

  * (1) Virginia Aglietti, Xiaoyu Lu, Andrei Paleyes, and Javier González.  Causal Bayesian Optimization.  In International Conference on Artificial Intelligence and Statistics, pages 3155–3164. PMLR, 2020.
  * (2) Omid Alipourfard, Hongqiang Harry Liu, Jianshu Chen, Shivaram Venkataraman, Minlan Yu, and Ming Zhang.  CherryPick: Adaptively Unearthing the Best Cloud Configurations for Big Data Analytics, 2017.
  * (3) Jason Ansel and Cy Chan.  PetaBricks.  XRDS: Crossroads, The ACM Magazine for Students, 17(1):32, 2010\.
  * (4) Jason Ansel, Shoaib Kamil, Kalyan Veeramachaneni, Jonathan Ragan-Kelley, Jeffrey Bosboom, Una-May O’Reilly, and Saman Amarasinghe.  OpenTuner.  In Proceedings of the 23rd International Conference on Parallel Architectures and Compilation - PACT ’14, pages 303–316, 2014\.
  * (5) Omid Azizi, Aqeel Mahesri, Benjamin C Lee, Sanjay J Patel, and Mark Horowitz.  Energy-performance tradeoffs in processor architecture and circuit design: A marginal cost analysis.  ACM SIGARCH Computer Architecture News, 38(3):26–36, 2010.
  * (6) Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy.  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  arXiv:1910.06403 [cs, math, stat], December 2020.
  * (7) Paul Beaumont, Ben Horsburgh, Philip Pilgerstorfer, Angel Droth, Richard Oentaryo, Steven Ler, Hiep Nguyen, Gabriel Azevedo Ferreira, Zain Patel, and Wesley Leong.  CausalNex, 2021.
  * (8) James Bergstra and Yoshua Bengio.  Random search for hyper-parameter optimization.  Journal of machine learning research, 13(Feb):281–305, 2012.
  * (9) Betsy Beyer, Chris Jones, Jennifer Petoff, and Niall Richard Murphy.  Site Reliability Engineering: How Google Runs Production Systems.  ” O’Reilly Media, Inc.”, 2016.
  * (10) Kshitij Bhardwaj, Marton Havasi, Yuan Yao, David M. Brooks, José Miguel Hernández Lobato, and Gu-Yeon Wei.  Determining optimal coherency interface for many-accelerator socs using bayesian optimization.  IEEE Computer Architecture Letters, 18(2):119–123, 2019.
  * (11) Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman.  Pyro: Deep Universal Probabilistic Programming.  October 2018.
  * (12) Nathan Binkert, Bradford Beckmann, Gabriel Black, Steven K Reinhardt, Ali Saidi, Arkaprava Basu, Joel Hestness, Derek R Hower, Tushar Krishna, and Somayeh Sardashti.  The gem5 simulator.  ACM SIGARCH computer architecture news, 39(2):1–7, 2011.
  * (13) Brian Brazil.  Prometheus: Up & Running: Infrastructure and Application Performance Monitoring.  ” O’Reilly Media, Inc.”, 2018.
  * (14) Tianqi Chen, Lianmin Zheng, Eddie Yan, Ziheng Jiang, Thierry Moreau, Luis Ceze, Carlos Guestrin, and Arvind Krishnamurthy.  Learning to optimize tensor programs.  arXiv preprint arXiv:1805.08166, 2018.
  * (15) Yu-Hsin Chen, Tien-Ju Yang, Joel Emer, and Vivienne Sze.  Eyeriss v2: A Flexible Accelerator for Emerging Deep Neural Networks on Mobile Devices.  arXiv:1807.07928 [cs], May 2019.
  * (16) Valentin Dalibard, Michael Schaarschmidt, and Eiko Yoneki.  BOAT: Building auto-tuners with structured Bayesian optimization.  In Proceedings of the 26th International Conference on World Wide Web - WWW ’17, pages 479–488, New York, New York, USA, 2017\. ACM Press.
  * (17) Andreas Damianou and Neil D Lawrence.  Deep gaussian processes.  In Artificial Intelligence and Statistics, pages 207–215. PMLR, 2013.
  * (18) Samuel Daulton, Maximilian Balandat, and Eytan Bakshy.  Differentiable expected hypervolume improvement for parallel multi-objective Bayesian optimization.  arXiv preprint arXiv:2006.05078, 2020.
  * (19) Stefan Falkner, Aaron Klein, and Frank Hutter.  BOHB: Robust and efficient hyperparameter optimization at scale.  arXiv preprint arXiv:1807.01774, 2018.
  * (20) Avrilia Floratou, Ashvin Agrawal, Bill Graham, Sriram Rao, and Karthik Ramasamy.  Dhalion: Self-regulating stream processing in heron.  Proceedings of the VLDB Endowment, 10(12):1825–1836, 2017.
  * (21) Qiang Fu, Jieming Zhu, Wenlu Hu, Jian-Guang Lou, Rui Ding, Qingwei Lin, Dongmei Zhang, and Tao Xie.  Where do developers log? an empirical study on logging practices in industry.  In Companion Proceedings of the 36th International Conference on Software Engineering, pages 24–33, 2014.
  * (22) Jacob R. Gardner, Geoff Pleiss, David Bindel, Kilian Q. Weinberger, and Andrew Gordon Wilson.  GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration.  arXiv:1809.11165 [cs, stat], January 2019.
  * (23) Aric Hagberg, Pieter Swart, and Daniel S Chult.  Exploring network structure, dynamics, and function using NetworkX.  Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008.
  * (24) gem5-aladdin harvard-acc.  Gem5-Aladdin SoC Simulator.  Harvard Architecture, Circuits, and Compilers, October 2016.
  * (25) gem5-aladdin-param harvard-acc.  Gem5-Aladdin SoC Simulator.  Harvard Architecture, Circuits, and Compilers, October 2021.
  * (26) Peter Henderson, Riashat Islam, Philip Bachman, Joelle Pineau, Doina Precup, and David Meger.  Deep reinforcement learning that matters.  arXiv preprint arXiv:1709.06560, 2017.
  * (27) Frank Hutter, Holger Hoos, and Kevin Leyton-Brown.  An evaluation of sequential model-based optimization for expensive blackbox functions.  In Proceedings of the 15th Annual Conference Companion on Genetic and Evolutionary Computation, pages 1209–1216, 2013.
  * (28) Max Jaderberg, Valentin Dalibard, Simon Osindero, Wojciech M. Czarnecki, Jeff Donahue, Ali Razavi, Oriol Vinyals, Tim Green, Iain Dunning, and Karen Simonyan.  Population based training of neural networks.  arXiv preprint arXiv:1711.09846, 2017.
  * (29) Max Jaderberg, Valentin Dalibard, Simon Osindero, Wojciech M. Czarnecki, Jeff Donahue, Ali Razavi, Oriol Vinyals, Tim Green, Iain Dunning, Karen Simonyan, Chrisantha Fernando, and Koray Kavukcuoglu.  Population Based Training of Neural Networks.  November 2017.
  * (30) Paul Kline.  An Easy Guide to Factor Analysis.  Routledge, 2014.
  * (31) Daphne Koller and Nir Friedman.  Probabilistic Graphical Models: Principles and Techniques.  MIT press, 2009.
  * (32) Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar.  Hyperband: A novel bandit-based approach to hyperparameter optimization.  The Journal of Machine Learning Research, 18(1):6765–6816, 2017\.
  * (33) Andy Liaw and Matthew Wiener.  Classification and regression by randomForest.  R news, 2(3):18–22, 2002.
  * (34) Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfei Cai.  When Gaussian Process Meets Big Data: A Review of Scalable GPs.  arXiv:1807.01065 [cs, stat], April 2019.
  * (35) David JC MacKay.  Bayesian neural networks and density networks.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 354(1):73–80, 1995.
  * (36) Azalia Mirhoseini, Hieu Pham, Quoc V. Le, Benoit Steiner, Rasmus Larsen, Yuefeng Zhou, Naveen Kumar, Mohammad Norouzi, Samy Bengio, and Jeff Dean.  Device Placement Optimization with Reinforcement Learning.  CoRR, abs/1706.0, June 2017.
  * (37) Ujval Misra, Richard Liaw, Lisa Dunlap, Romil Bhardwaj, Kirthevasan Kandasamy, Joseph E. Gonzalez, Ion Stoica, and Alexey Tumanov.  RubberBand: Cloud-based hyperparameter tuning.  In Proceedings of the Sixteenth European Conference on Computer Systems, pages 327–342, Online Event United Kingdom, April 2021\. ACM.
  * (38) Bruce Momjian.  PostgreSQL: Introduction and Concepts, volume 192.  Addison-Wesley New York, 2001.
  * (39) Luigi Nardi, David Koeplinger, and Kunle Olukotun.  Practical design space exploration.  In 2019 IEEE 27th International Symposium on Modeling, Analysis, and Simulation of Computer and Telecommunication Systems (MASCOTS), pages 347–358. IEEE, 2019.
  * (40) OpenSource.  Microsoft/nni.  Microsoft, December 2020.
  * (41) Andy D Pimentel.  Exploring exploration: A tutorial introduction to embedded systems design space exploration.  IEEE Design & Test, 34(1):77–90, 2016.
  * (42) Carl Edward Rasmussen and Christopher K. I. Williams.  Gaussian Processes for Machine Learning.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, Mass., 3. print edition, 2008.
  * (43) Brandon Reagen, Robert Adolf, Yakun Sophia Shao, Gu-Yeon Wei, and David Brooks.  Machsuite: Benchmarks for accelerator design and customized architectures.  In 2014 IEEE International Symposium on Workload Characterization (IISWC), pages 110–119. IEEE, 2014.
  * (44) Mauro Scanagatta, Antonio Salmerón, and Fabio Stella.  A survey on Bayesian network structure learning from data.  Progress in Artificial Intelligence, pages 1–15, 2019.
  * (45) Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas.  Taking the human out of the loop: A review of Bayesian optimization.  Proceedings of the IEEE, 104(1):148–175, 2015.
  * (46) Yakun Sophia Shao, Brandon Reagen, Gu-Yeon Wei, and David Brooks.  Aladdin: A pre-rtl, power-performance accelerator simulator enabling large design space exploration of customized architectures.  In 2014 ACM/IEEE 41st International Symposium on Computer Architecture (ISCA), pages 97–108. IEEE, 2014.
  * (47) Yakun Sophia Shao, Sam Likun Xi, Vijayalakshmi Srinivasan, Gu-Yeon Wei, and David Brooks.  Co-designing accelerators and soc interfaces using gem5-aladdin.  In 2016 49th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO), pages 1–12. IEEE, 2016.
  * (48) Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical bayesian optimization of machine learning algorithms.  In Advances in Neural Information Processing Systems, pages 2951–2959, 2012.
  * (49) Richard S. Sutton and Andrew G. Barto.  Reinforcement Learning : An Introduction.  MIT Press, Cambridge Mass., 1998.
  * (50) Dana Van Aken, Andrew Pavlo, Geoffrey J. Gordon, and Bohan Zhang.  Automatic Database Management System Tuning Through Large-scale Machine Learning.  In Proceedings of the 2017 ACM International Conference on Management of Data - SIGMOD ’17, pages 1009–1024, New York, New York, USA, 2017. ACM Press.
  * (51) Dana Van Aken, Dongsheng Yang, Sebastien Brillard, Ari Fiorino, Bohan Zhang, Christian Bilien, and Andrew Pavlo.  An inquiry into machine learning-based automatic configuration tuning services on real-world database management systems.  Proceedings of the VLDB Endowment, 14(7):1241–1253, 2021.
  * (52) Shivaram Venkataraman, Zongheng Yang, Michael Franklin, Benjamin Recht, and Ion Stoica.  Ernest: Efficient performance prediction for large-scale advanced analytics.  In 13th {}USENIX{} Symposium on Networked Systems Design and Implementation ({}NSDI{} 16), pages 363–378, 2016.
  * (53) Ziyu Wang, Frank Hutter, Masrour Zoghi, David Matheson, and Nando de Feitas.  Bayesian optimization in a billion dimensions via random embeddings.  Journal of Artificial Intelligence Research, 55:361–387, 2016.
  * (54) Svante Wold, Kim Esbensen, and Paul Geladi.  Principal component analysis.  Chemometrics and intelligent laboratory systems, 2(1-3):37–52, 1987\.
  * (55) Ding Yuan, Soyeon Park, and Yuanyuan Zhou.  Characterizing logging practices in open-source software.  In 2012 34th International Conference on Software Engineering (ICSE), pages 102–112. IEEE, 2012.
  * (56) Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing.  Dags with no tears: Continuous optimization for structure learning.  arXiv preprint arXiv:1803.01422, 2018.