# Bayesian Optimization over Permutation Spaces

Aryan Deshwal, Syrine Belakaria, Janardhan Rao Doppa, Dae Hyun Kim

###### Abstract

Optimizing expensive to evaluate black-box functions over an input space
consisting of all permutations of $d$ objects is an important problem with
many real-world applications. For example, placement of functional blocks in
hardware design to optimize performance via simulations. The overall goal is
to minimize the number of function evaluations to find high-performing
permutations. The key challenge in solving this problem using the Bayesian
optimization (BO) framework is to trade-off the complexity of statistical
model and tractability of acquisition function optimization. In this paper, we
propose and evaluate two algorithms for BO over Permutation Spaces (BOPS).
First, BOPS-T employs Gaussian process (GP) surrogate model with Kendall
kernels and a Tractable acquisition function optimization approach based on
Thompson sampling to select the sequence of permutations for evaluation.
Second, BOPS-H employs GP surrogate model with Mallow kernels and a Heuristic
search approach to optimize expected improvement acquisition function. We
theoretically analyze the performance of BOPS-T to show that their regret
grows sub-linearly. Our experiments on multiple synthetic and real-world
benchmarks show that both BOPS-T and BOPS-H perform better than the state-of-
the-art BO algorithm for combinatorial spaces. To drive future research on
this important problem, we make new resources and real-world benchmarks
available to the community.

##  1 Introduction

Optimizing expensive black-box functions is a common problem in many science
and engineering applications (Frazier and Wang 2016; Yang, Wu, and Arnold
2019; Zhou et al. 2020; Deshwal, Simon, and Doppa 2021). An important class of
such problems involve optimizing over the space of all permutations of a set
of objects. For example, in the design of integrated circuits (ICs), we need
to find placement of functional blocks to optimize performance via expensive
simulations. Another example is to find a schedule for the jobs in additive
manufacturing (parts can be built layer-by-layer) to optimize throughput via
simulations. Bayesian optimization (BO) (Shahriari et al. 2016) has proven to
be a successful approach for optimizing black-box functions in a sample-
efficient manner. The key idea is to learn a surrogate statistical model,
e.g., Gaussian process, and use it to select the sequence of inputs for
function evaluation to uncover high-performing inputs.

A large body of the BO literature focuses on continuous input spaces
(Greenhill et al. 2020) with few recent works (Oh et al. 2019; Baptista and
Poloczek 2018; Dadkhahi et al. 2020; Doppa 2021) tackling the challenging
setting of input spaces over discrete structures, e.g., sets, sequences, and
graphs. Unlike continuous spaces, discrete spaces come with many unique
challenges such as difficulty of defining a general representation, non-
smoothness, etc. which require specialized treatment of different types of
combinatorial structures (BakIr et al. 2007). In this paper, we consider the
understudied space of all permutations of $d$ objects, which is equivalently
characterized by the non-Abelian Symmetric group $\mathcal{S}_{d}$ of
cardinality $d!$, consisting of all bijections from a set to itself. In
contrast, all of the existing work on combinatorial BO is focused on the input
space with categorical/binary variables which corresponds to the direct sums
of abelian cyclic groups $\mathbb{Z}/c\mathbb{Z}$ of cardinality $c^{d}$,
where each of the $d$ variables can take one of the $c$ categories.
Combinatorial BO methods such as COMBO (Oh et al. 2019) cannot exploit the
special characteristics of permutation spaces. For example, they account for
$d^{d}$ space $\gg$ $d!$ permutation space.

The key challenge to devise effective solutions for BO over Permutation Spaces
(BOPS) is to trade-off the complexity of statistical model and tractability of
search to select permutations with high utility/acquisition values (e.g.,
expected improvement). We propose and evaluate two algorithms for BOPS by
addressing this challenge. First, BOPS-T uses a simpler statistical model in
the form of GP with Kendall kernels (Jiao and Vert 2015). By considering the
weight space view of the GP model and employing Thompson sampling as the
acquisition strategy, we show that the acquisition function optimization is a
Quadratic Assignment problem which can be solved using Tractable Semi-definite
programming relaxation based solvers. Second, BOPS-H uses a complex
statistical model in the form of GP with Mallow kernels (Jiao and Vert 2015;
Mania et al. 2018). We employ expected improvement as the acquisition strategy
and perform Heuristic search (local search with multiple restarts) to select
the permutations for function evaluation. We analyze the theoretical
properties of BOPS-T and in terms of regret bounds and show that it achieves a
sublinear (time) regret. Our comprehensive experiments on both synthetic and
real-world benchmarks show that both BOPS-T and BOPS-H perform better than
COMBO (Oh et al. 2019).

Contributions. The main contributions of this paper are:

  * •

Two BO algorithms for permutation spaces, namely, BOPS-T and BOPS-H, that make
varying trade-offs between the complexity of statistical model and
tractability of search for selecting permutations for evaluation.

  * •

We theoretically analyze BOPS-T and show that it achieves sublinear (in time)
regret bounds of $\mathcal{O}^{*}(d^{3/2}\sqrt{T})$, where $d$ and $T$ refers
to dimensionality and time (no. of iterations) respectively and
$\mathcal{O}^{*}$ denotes upto logarithmic factors.

  * •

We evaluate the efficacy of our algorithms and compare with the state-of-the-
art COMBO algorithm on multiple synthetic and real-world benchmarks.

  * •

We make new resources and benchmarks based on important real-world problems
available to the BO community to drive future research on this under-studied
problem. The resources and source code of BOPS algorithms are available at
https://github.com/aryandeshwal/BOPS/.

##  2 Problem Formulation

In this paper, we consider optimization problems with the input space
consisting of all permutations over $d$ objects. Given
$[1,d]:=\\{1,2,\cdots,d\\}$, indexing the $d$ objects, a permutation is
defined as a bijective mapping $\pi:[1,d]\mapsto[1,d]$. The set of all
permutations along with the composition binary operation
$\left((\pi_{1}\circ\pi_{2})(x)=\pi_{1}(\pi_{2}(x)\right)$ is known as the
Symmetric group $\mathcal{S}_{d}$ which has a cardinality $|\mathcal{S}_{d}|$=
$d!$.

Let $f:\mathcal{S}_{d}\mapsto\mathbb{R}$ be a black-box objective function
that is expensive to evaluate. Our goal is to optimize $f$ while minimizing
the number of function evaluations:

$\displaystyle\pi^{*}=arg\min_{\pi\in\mathcal{S}_{d}}f(\pi)$ (2.1)

For a concrete example problem, consider the domain of design and optimization
of integrated circuits (ICs). There are many applications in IC design, where
we need to optimize over permutations of functional blocks of different
granularity (small cells to processing cores). Some example objectives include
performance and manufacturing cost. We need to perform expensive computational
simulations to evaluate each candidate permutation.

We solve this problem using the Bayesian optimization (BO) framework. A BO
algorithm is composed of two main components: 1) a surrogate statistical model
built using past evaluations of the expensive black-box objective function;
and 2) an acquisition function to capture the expected utility of evaluating a
new input, which is optimized at each BO iteration to find the next best input
for evaluation. Algorithm 1 shows the generic BOPS approach. We provides
concrete instantiations for these two components by handling the unique
challenges of permutation spaces. In each BO iteration, we select the most
promising permutation for evaluation (line 3) and update the statistical
models using all the training examples (line 6). The best found permutation is
returned at the end of maximum BO iterations (line 8).

Algorithm 1 Bayesian Optimization over Permutation Spaces (BOPS)

requires: black-box objective $f(\pi)$, Gaussian process $GP(0,k)$ with Kernel
over permutations $k$, Acquisition function $AF$

1: ${D}_{0}\leftarrow$ small initial data; and $GP_{0}\leftarrow
GP(\mu_{D_{0}},k_{D_{0}})$

2: for $j$=$1,2,\dots$ do

3: Acquisition function optimization to select the next permutation: $\pi_{j}$
= $arg\min_{\pi\in S_{d}}AF(\pi)$

4: Evaluate the selected permutation $\pi_{j}$ to get $f(\pi_{j})$

5: Aggregate training data: ${D_{j}}\leftarrow
D_{j-1}\cup\\{\pi_{j},f(\pi_{j})\\}$

6: Update GP posterior: $GP_{j}\leftarrow GP(\mu_{D_{j}},k_{D_{j}})$

7: end for

8: return $\pi_{best}$ = $arg\min\\{f(\pi_{1}),f(\pi_{2})\cdots\\}$

##  3 Related Work

Bayesian optimization (BO) methods for input spaces over discrete structures
can be broadly classified into two categories based on the amount of available
training data: large and small data settings. Our paper focuses on the small-
data setting due to the unavailability of large training sets and the
objective function evaluations being expensive.

BO methods for large-data setting. Most of these methods employ the large set
of training examples to learn embeddings of discrete structures in a latent
space and perform continuous BO in the latent space (Gómez-Bombarelli et al.
2018; Lu et al. 2018; Eissman et al. 2018; Kajino 2019). This setting is
especially relevant for domains such as molecular design (Griffiths and
Hernández-Lobato 2017; Hernández-Lobato et al. 2017) and biological sequence
design (Angermueller et al. 2020, 2019).

BO methods for small-data setting. Prior work in the small-data setting
consider input spaces over binary/categorical structures: $x_{1}\times
x_{2}\cdots\times x_{d}$ space, where each $x_{i}$ can take all possible
values from a set of categories $\\{0,1,\cdots,c-1\\}$, where $c$=2 and $c>2$
for binary and categorical variables respectively. This results in a search
space of size $c^{d}$, which has a different structure than our considered
setting of the space of all permutations of a set of $d$ objects. Prior
methods in this setting can be classified based on their choice of surrogate
models (linear models (Baptista and Poloczek 2018; Dadkhahi et al. 2020),
random forests (Hutter, Hoos, and Leyton-Brown 2011), Gaussian processes (Oh
et al. 2019; Deshwal and Doppa 2021)) or acquisition function optimization
approaches (heuristic search (Hutter, Hoos, and Leyton-Brown 2011; Oh et al.
2019; Dadkhahi et al. 2020; Deshwal and Doppa 2021), principled mathematical
optimization (Baptista and Poloczek 2018; Deshwal, Belakaria, and Doppa 2020,
2021)), and combination thereof (Garrido-Merchán and Hernández-Lobato 2020;
Luong et al. 2019; Deshwal et al. 2020; Swersky et al. 2020).

Gaussian process models for discrete structures. COMBO (Oh et al. 2019) is the
state-of-the-art method for small-data setting. It employs Gaussian process
(GP) with diffusion kernels (Kondor and Lafferty 2002) defined over a
combinatorial graph representation of the input space. However, COMBO’s graph
representation cannot exploit the special characteristics of permutation
spaces. For example, for a space of $d$ objects, it requires $d$ subgraphs
each with at least $d$ nodes to account for $d^{d}$ combinatorial space which
is unnecessarily large for our usage. Another line of work considers GP models
with kernels for specific structures including sets (double sum/embedding
kernels (Buathong, Ginsbourger, and Krityakierne 2020; Kim et al. 2021)),
molecular graphs (optimal transport based kernel (Korovina et al. 2020)), and
strings (string kernels (Moss et al. 2020; Lodhi et al. 2002)). We employ
Kendall kernels and Mallow kernels (Jiao and Vert 2015; Mania et al. 2018) for
permutation spaces which were shown to have superior performance on
permutation based classification tasks.

##  4 BO Algorithms for Permutation Spaces

In this section, we provide two algorithms for BO over permutation spaces that
make varying trade-offs between the complexity of statistical model and
tractability of acquisition function optimization. First, BOPS-T employs a
simple statistical model with an efficient Semi-definite programming (SDP)
relaxation based optimization method. Second, BOPS-H employs a complex
statistical model and performs heuristic search for optimizing the acquisition
function. We employ Gaussian processes (GPs) (Rasmussen and Williams 2006) as
the surrogate model in both algorithms. GPs are effective statistical models
commonly used for BO as they provide a principled framework for uncertainty
estimation. They are fully characterized by a kernel $k$ (Kanagawa et al.
2018) which intuitively captures the similarity between two candidate inputs
from the same input space.

###  4.1 BOPS-T Algorithm

Surrogate model. The similarity between any two permutations
($\pi,\pi^{\prime}$) can be naturally defined by considering the number of
pairs of objects ordered in the same way or in opposite ways. This is captured
by the notion of the number of discordant pairs $n_{d}(\pi,\pi^{\prime})$,
also known as Kendall-tau distance (Kendall 1938). $n_{d}(\pi,\pi^{\prime})$
counts the number of pairs of objects ordered oppositely by $\pi$ and
$\pi^{\prime}$ as defined below:

$\displaystyle\begin{split}n_{d}(\pi,\pi^{\prime})=\sum_{i<j}[1_{\pi(i)>\pi(j)}&1_{\pi^{\prime}(i)<\pi^{\prime}(j)}\\\
&+1_{\pi(i)<\pi(j)}1_{\pi^{\prime}(i)>\pi^{\prime}(j)}]\end{split}$ (4.1)

A related notion of concordant pairs $n_{c}(\pi,\pi^{\prime})$ counts the
number of object pairs ordered similarly by $\pi$ and $\pi^{\prime}$:

$\displaystyle n_{c}(\pi,\pi^{\prime})=\binom{d}{2}-n_{d}(\pi,\pi^{\prime})$
(4.2)

Kendall kernels (Jiao and Vert 2015) are positive-definite kernels defined
over permutations using the notion of discordant and concordant pairs as
follows:

$\displaystyle
k(\pi,\pi^{\prime})=\frac{n_{c}(\pi,\pi^{\prime})-n_{d}(\pi,\pi^{\prime})}{\binom{d}{2}}$
(4.3)

Because of their proven effectiveness over permutations (Jiao and Vert 2015,
2018), we propose using Kendall kernels with GPs as surrogate model in our
BOPS-T algorithm.

For our surrogate model, we consider the weight-space formulation of the GP.
This weight-space formulation is essential for the SDP based acquisition
function optimization approach described in the next section. In the weight-
space view, we can reason about GPs as a weighted sum of basis functions
$\phi$ = $\\{{\phi_{i}(\cdot)}\\}$, i.e.,

$\displaystyle w^{T}\phi(\cdot);\quad w\sim N(0,I)$ (4.4)

where $N(\cdot)$ represents multi-variate Gaussian distribution and $I$ is the
identity matrix. Every kernel has a canonical feature map (as per the Moore-
Aronszajn theorem (Aronszajn 1950)) $\phi:\mathcal{S}_{d}\mapsto H_{k}$,
$H_{k}$ being its associated Reproducing Kernel Hilbert Space (RKHS), that is
employed as the basis function in 4.4. The feature map expression for Kendall
kernel (constructed by (Jiao and Vert 2015)) is given below:

$\displaystyle\phi(\pi)=\\{\sqrt{\binom{d}{2}^{-1}}\left(1_{\pi(i)>\pi(j)}-1_{\pi(i)<\pi(j)}\right)\\}_{(1\leq
i<j\leq d)}$ (4.5)

Acquisition function and optimizer. In order to sequentially select the next
permutation for evaluation guided by the learned surrogate model, we employ
Thompson sampling as our acquisition function. Thompson sampling is a
powerful, practitioner-friendly, and parameter-free approach for appropriately
balancing the exploration vs. exploitation dilemma (Russo and Van Roy 2014) in
sequential bandit optimization. The key idea is to sample a function from the
surrogate model’s posterior and select its optimizer as the next permutation
for evaluation. In the weight-space view of GPs, this corresponds to sampling
a weight vector $\hat{w}$ from its posterior and solving the following
optimization problem:

$\displaystyle\pi_{next}=arg\min_{\pi\in\mathcal{S}_{d}}\hat{w}^{T}\phi(\pi)$
(4.6)

It should be noted that the sampled weight vector $\hat{w}$ is an exact
function defined by GP (with Kendall kernel) over permutation spaces and has
no approximation error when compared to the function space approach. This is
in contrast to the common practice of using Thompson sampling over continuous
spaces, where random Fourier features based weight-space representation of GPs
is used which inevitably results in approximation error because of sampling a
finite number of features from an infinite feature space.

We now show that the above acquisition function optimization problem (4.6) is
a Quadratic Assignment Problem (QAP) (Burkard et al. 1998). To observe that,
the objective in 4.6 is written in an equivalent form in terms of $P_{d}$, the
set of all possible permutation matrices $P$ of size $d\times d$, as follows:

$\displaystyle\min_{P\in P_{d}}Tr(WPAP^{T})$ (4.7)

where $Tr$ is the matrix trace operation and $A$ is a $d\times d$ matrix
defined as follows:

$\displaystyle A=\left\\{\begin{array}[]{ll}1&\mbox{if }i<j\\\ -1&\mbox{if
}i>j\\\ 0&\mbox{if }i==j\end{array}\right.\forall i,j\in[1,d]$

and $W$ is another $d\times d$ matrix given as follows:

$\displaystyle
W=\left\\{\begin{array}[]{ll}w_{\frac{(i-1)}{2}(2d-i)+(j-i)}&\mbox{if }i<j\\\
0&\mbox{if }i\geq j\end{array}\right.\forall i,j\in[1,d]$

Concretely, the equivalence of objectives in 4.6 and 4.7 can be seen as
follows:

$\displaystyle Tr(WPAP^{T})=\sum_{i=1}^{d}(WPAP^{T})_{ii}\quad$ (4.8)

Equation 4.8 is the definition of the trace of a matrix. Now, considering each
entry $(WPAP^{T})_{ii}$ in 4.8:

$\displaystyle(WPAP^{T})_{ii}$
$\displaystyle=\sum_{j=1}^{d}W_{ij}\cdot(PAP^{T})_{ji}$ (4.9)
$\displaystyle=\sum_{j>i}^{d}w_{\frac{(i-1)}{2}(2d-i)+(j-i)}\cdot(PAP^{T})_{ji}$
(4.10) $\displaystyle=\sum_{j>i}^{d}w_{\frac{(i-1)}{2}(2d-i)+(j-i)}\cdot
A_{\pi(j)\pi(i)}$ (4.11)

where 4.10 follows from the definition of $W$ and 4.11 follows from the fact
that pre-multiplying (post-multiplying) by a permutation matrix permutes the
rows (columns) of $A$. Using 4.11 in 4.8:

$\displaystyle
Tr(WPAP^{T})=\sum_{i=1}^{d}\sum_{j>i}^{d}w_{\frac{(i-1)}{2}(2d-i)+(j-i)}\cdot
A_{\pi(j)\pi(i)}$ (4.12)

By noting that $A_{\pi(j)\pi(i)}$ is exactly the feature map in 4.5 (upto
multiplication by a constant $\sqrt{\binom{d}{2}^{-1}}$ which doesn’t change
the optimal solution), the equivalence between 4.6 and 4.7 is established.

Although, in general, Quadratic assignment probem is NP-hard (Sahni and
Gonzalez 1976), we leverage existing Semi-definite programming (SDP) based
strong relaxations (Zhao et al. 1998) to obtain good approximate solutions to
the acquisition function optimization problem. Using the invariance of the
trace under cyclic permutations and vectorization identity
$(vec(APW)=(W^{T}\otimes A)vec(P))$, objective in 4.7 is standardized as:

$\displaystyle\min_{P,Q}$ $\displaystyle((W^{T}\otimes A)Q)$ (4.13)
$\displaystyle P$ $\displaystyle\in P_{n}$ $\displaystyle Q$
$\displaystyle=vec(P)vec(P)^{T}$

where $vec(P)$ is the column-wise vectorization of $P$. We leverage the
clique-based SDP relaxation approach of (Ferreira, Khoo, and Singer 2018)
which can exploit matrix sparsity (e.g., zeros in the upper-triangular matrix
$W$) for solving 4.13. The key idea is to enforce semi-definiteness only over
groups of $Q$’s entries (i.e., cliques) to get a relaxation that can be solved
using fast and accurate algorithms.

###  4.2 BOPS-H Algorithm

Surrogate model. We propose to employ Mallows kernel which plays a role on the
symmetric group $\mathcal{S}_{d}$ similar to the Gaussian (RBF) kernel on the
Euclidean space. Given a pair of permutations $\pi$ and $\pi^{\prime}$, the
Mallows kernel is defined as the exponentiated negative of the number of
discordant pairs $n_{d}(\pi,\pi^{\prime})$ between $\pi$ and $\pi^{\prime}$
i.e.

$\displaystyle k_{m}{\pi,\pi^{\prime}}=\exp(-ln_{d}(\pi,\pi^{\prime}))$ (4.14)

where $l\geq 0$ is a hyper-parameter of the kernel similar to the length-scale
hyper-parameter of the Gaussian kernels on Euclidean space. A key measure of
the expressivity of a kernel is based on a property called universality which
captures the notion of whether the RKHS of the kernel is rich enough to
approximate any function on a given input space arbitrary well. It was
recently shown (Mania et al. 2018) that Mallows kernel is universal over the
space of permutations in contrast to the Kendall kernel discussed in the
previous section. Therefore, Mallows kernels are more powerful than Kendall
kernels and allows us to capture richer structure in permutations when used to
learn GP based surrogate models. Indeed, our experiments also demonstrate
empirically the superior modeling capability of Mallows Kernel.

Acquisition function and optimizer. Unlike Kendall kernel, the feature space
of Mallows Kernel is exponentially large (Mania et al. 2018) making it
practically inefficient to sample functions from the GP posterior (in the
weight-space style as described earlier).Therefore, we propose to employ
expected improvement (EI) as our acquisition function. The additional
complexity of GP based statistical model with Mallows kernel makes the
acquisition function optimization problem $\pi_{next}$ =
$arg\min_{\pi\in\mathcal{S}_{d}}AF(\pi)$ is intractable for EI. Therfore, we
propose to perform Heuristic search in the form of local search with multiple
restarts that has been shown to be very effective in practice for solving
combinatorial optimization problems. To search over only valid permutations
$\pi\in\mathcal{S}_{d}$, at each local search step, we consider only those
neighbors which are permutations of the current state. Otherwise, we will be
searching over a huge combinatorial space with both valid (permutations) and
invalid structures (non-permutations) as done by COMBO: may not result in
producing a permutation from its acquisition function optimization procedure.
Indeed, we observed this behavior in our experiments with COMBO. We use the
modified local search procedure over permutations for both COMBO and our
BOPS-H algorithm in experiments.

##  5 Theoretical Analysis for BOPS-T

In this section, we analyze the theoretical properties of our BOPS-T algorithm
in terms of regret metric (Srinivas et al. 2009), which is a commonly used
measure for analyzing BO algorithms. Note that there is no prior regret bound
analysis for BO algorithms for EI even in continuous spaces. Hence, we leave
the analysis of BOPS-H algorithm for future work. Let simple regret $R$ be
defined as follows:

$\displaystyle R=\sum_{t=1}^{T}(f(\pi_{t})-f(\pi^{*}))$ (5.1)

where $\pi_{t}$ is the permutation picked by the BO algorithm at time
(iteration) $t$. In our case of using Thompson sampling as an acquisition
function, it is natural to consider the expected form of this regret (Russo
and Van Roy 2014) where the expectation is taken over the distribution of
functions as given by the GP prior with Kendall kernel. We analyze this
expected form of regret, also known as Bayesian regret:

$\displaystyle\mathcal{BR}=\sum_{t=1}^{T}\mathbb{E}(f(\pi_{t})-f(\pi^{*}))$
(5.2)

where the expectation is over the distribution of functions $f\sim GP(0,k)$.
The below theorem bounds the Bayesian regret of our BOPS-T algorithm:

###### Theorem 5.1

Let $f\sim GP(0,k)$ with Kendall kernel $k$ (4.3), the Bayesian regret of the
BOPS-T algorithm after $T$ observations $y_{i}$ =
$f(\pi_{i})+\epsilon_{i},i\in\\{1,2,\cdots T\\}$ with $\epsilon_{i}$ being
Gaussian distributed i.i.d. noise $\epsilon_{i}\sim N(0,\sigma^{2})$ is :
$\mathcal{BR}$ = $\mathcal{O}^{*}(d^{3/2}\sqrt{T})$, where $\mathcal{O}^{*}$
denotes upto log factors.

Proof. The key quantity in bounding the regret of Bayesian optimization with
Gaussian processes (also known as GP bandits) is an information-theoretic
quantity called as maximum information gain $\gamma_{T}$ (Srinivas et al.
2009) that depends on the kernel $k$ and intuitively captures the maximum
information that can be gained about $f$ after $T$ observations, i.e.,

$\displaystyle\gamma_{T}=\max_{A\subset\mathcal{S}_{d},|A|=T}I(\mathbf{y_{A}};f)$
(5.3)

where $I$ is the mutual information and $A$ is a subset of permutations with
corresponding function evaluations $\mathbf{y_{A}}$.

(Russo and Van Roy 2014) proved the Bayesian regret for Thompson sampling by
characterizing it in terms of upper confidence bound based results from
(Srinivas et al. 2009):

###### Proposition 1

(Proposition 5 (Russo and Van Roy 2014)). If $|X|<\infty$, $\\{f(x):x\in X\\}$
follows a multivariate Gaussian distribution with marginal variances bounded
by 1, the Bayesian regret for Thompson sampling based bandit policy is given
as:

$\displaystyle\mathcal{BR}=1+2\sqrt{T\gamma_{T}\ln{(1+\sigma^{2})^{-1}}\ln{\left(\frac{(T^{2}+1)|X|}{\sqrt{2\pi}}\right)}}$
(5.4)

where $X$ is the action space.

This proposition is directly applicable in our setting because the action
space, being the cardinality of the symmetric group $\mathcal{S}_{d}$, is
finite (i.e., $|\mathcal{S}_{d}|$ = $d!$) and the function
$\\{f(\pi):\pi\in\mathcal{S}_{d}\\}$ follows a multivariate Gaussian
distribution (by definition of Gaussian process with Kendall kernel). We
compute the specific terms in the right-hand side of 5.4 that are applicable
in our setting to prove the regret bound.

The maximum information gain for kernels with finite feature maps can be
computed in the weight-space form (Sec 4.1) as a special case of linear kernel
(Srinivas et al. 2009).

$\displaystyle\gamma_{T}\leq C\log|I+\sigma^{-2}K|$ (5.5)

where $C$ = $1/2\cdot(1-1/e)^{-1}$ is a constant, $K$ is a $T\times T$ matrix
with each entry $K_{ij}$ = $k(\pi_{i},\pi_{j})$. As per kernel trick,

$\displaystyle K=\Phi^{T}\Phi$ (5.6)

where $\Phi$ is a matrix with
$\Sigma^{1/2}\phi(\pi_{i}),i\in\\{1,2,\cdots,T\\}$ as the columns (4.5).
Therefore,

$\displaystyle\gamma_{T}\leq C\ln|I+\sigma^{-2}\Phi^{T}\Phi|$ (5.7)

By Schur’s complement:

$\displaystyle\gamma_{T}\leq C\ln|I+\sigma^{-2}\Phi^{T}\Phi|\leq
C\ln|I+\sigma^{-2}\Phi\Phi^{T}|$ (5.8)

By Hadamard’s inequality:

$\displaystyle\gamma_{T}$ $\displaystyle\leq C\ln|I+\sigma^{-2}\Phi\Phi^{T}|$
(5.9) $\displaystyle\leq
C\sum_{i=1}^{\binom{d}{2}}\ln(1+\sigma^{-2}\lambda_{i})$ (5.10)

where $\\{\lambda_{1},\lambda_{2},\cdots\\}$ is the eigenvalue set of the
matrix $\Phi\Phi^{T}$.

By Gershgorin circle theorem (Varga 2010), all the eigenvalues of a matrix is
upper bounded by the maximum absolute sum of rows, i.e. $\lambda_{i}\leq
d^{2}T$ with the assumption that $\|\Sigma^{1/2}\phi(\pi)\|\leq 1$.

$\displaystyle\gamma_{T}=O(d^{2}\ln(d^{2}T))$ (5.11)

Now, using Stirling’s approximation, we can bound the $\ln(|X|)$ term in 5.4,
where $|X|=|\mathcal{S}_{d}|$ in our case:

$\displaystyle\ln(|\mathcal{S}_{d}|)=O(d\ln d)$ (5.12)

Plugging 5.11 and 5.12 in 5.4 and ignoring constants, we get the following
expression:

$\displaystyle\mathcal{BR}$ $\displaystyle=O(\sqrt{Td^{2}\ln d^{2}T(\ln
T^{2}+d\ln d)})$ (5.13) $\displaystyle\mathcal{BR}$
$\displaystyle=O(\sqrt{(Td^{2}\ln d^{2}T\ln T^{2}+Td^{3}\ln d^{2}T\ln d)})$
(5.14) $\displaystyle\mathcal{BR}$
$\displaystyle=\mathcal{O}^{*}(d^{3/2}\sqrt{T})$ (5.15)

Hence, ignoring log factors, our proposed BOPS-T algorithm achieves sublinear
(time) regret.

##  6 Experiments and Results

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: Results comparing BOPS-T, BOPS-H, and COMBO (best objective function
value vs. number of BO iterations) on both synthetic and real-world
benchmarks: (Top row) QAP, TSP, CP; and (Bottom row) FP1, FP2, and HMD.

In this section, we describe the benchmarks and experimental setup followed by
results and discussion.

###  6.1 Benchmarks

We employ diverse and challenging benchmarks for black-box optimization over
permutations for our experiments. We have the following two synthetic
benchmarks.

1) Quadratic assignment problem (QAP). QAPLIB (Burkard, Karisch, and Rendl
1997) is a popular library that contains multiple QAP instances. Each QAP
instance contains a cost matrix ($A$) and distance matrix ($B$) sized $n\times
n$, where $n$ is the number of input dimensions. The goal is to find the best
permutation that minimizes the quadratic assignment objective $Tr(APBP^{T})$,
where $P$ is an $n\times n$ permutation matrix. We use input space with $n=15$
dimensions in our experiments.

2) Traveling salesman problem (TSP). TSP problems are derived from low-
dimensional variants of the printed circuit board (PCB) problems from the
TSPLIB library (Reinelt 1995). The overall goal is to find the route of
drilling holes in the PCB that minimizes the time taken to complete the job.
We use input space with $d=10$ dimensions from the data provided in the
library.

We perform experiments on three important real-world applications from the
domain of computer-aided design of integrated circuits (ICs). These
applications are characterized by permutations over functional blocks at
different levels of granularity that arise in different stages of design and
optimization of ICs. Importantly, even tiny improvements in solution has huge
impact (e.g., improved performance over the lifespan of the IC or reduced cost
for manufacturing large samples of the same IC). A big challenge in the
combinatorial BO literature is the availability of challenging real-world
problems to evaluate new approaches. Hence, we provide our three real-world
benchmarks as a new resource to allow rapid development of the field.

3) Floor planning (FP). We are given $k$ rectangular blocks with varying width
and height, where each block represents a functional module performing certain
task. Each placement of the given blocks is called a floor-plan. Our goal is
to find the floor plan that minimizes the manufacturing cost per chip. We use
two variants of this benchmark with 10 blocks (FP1 and FP2) that differ in the
functionality of the blocks.

4) Cell placement (CP). We are given 10 rectangular cells with same height and
a netlist that contains the connection information among the cells. The goal
is to place the 10 rectangular cells for optimizing the performance of the
circuit. Intuitively, shorter nets have shorter delays, so placements with
shorter wire-length will result in higher performance.

5) Heterogeneous manycore design (HMD). This is a manycore architecture
optimization problem from the rodinia benchmark (Che et al. 2009). We are
given 16 cores of three types: 2 CPUs, 10 GPUs, and 4 memory units. They are
connected by a mesh network (each core is connected to its four neighboring
cores) to facilitate data transfer. The goal is to place the given 16 cores to
optimize the energy delay product (EDP) objective that captures both latency
and energy, two key attributes of a manycore chip.

(a)

(b)

(c)

(d)

(e)

Figure 2: Results comparing the three surrogate models of BOPS-T, BOPS-H and
COMBO on negative log-likelihood (NLL) metric computed on a test set on five
benchmarks: (Top row) QAP, TSP, CP; and (Bottom row) FP1, FP2.

###  6.2 Experimental Setup

Configuration of algorithms. We compare our proposed BOPS-T and BOPS-H
algorithms with the state-of-the-art combinatorial BO algorithm COMBO (Oh et
al. 2019). COMBO employs a diffusion kernel based GP surrogate model and
optimizes expected improvement acquisition function using local search with
restarts to select inputs for evaluation. Each local search step considers all
neighbors of the current structure in the combinatorial graph (i.e.,
structures with Hamming distance one). We modify COMBO’s local search
procedure (https://github.com/QUVA-Lab/COMBO) to consider only those neighbors
which are permutations of the current state thereby helping COMBO to avoid
searching a large combinatorial space with huge number of invalid structures
(non-permutations).

We used the SDP relaxation based QAP solver code from
(https://github.com/fsbravo/csdp) for implementing BOPS-T. BOPS-H is built
using popular GPyTorch (Gardner et al. 2018) and BoTorch (Balandat et al.
2020) libraries. We used 10 restarts for local search based EI optimization
for BOPS-H. BOPS-T, BOPS-H, and COMBO are initialized with the same 20 random
permutations in each experiment.

Evaluation metric. We plot the objective function value of the best
permutation over different BO iterations. Each experiment is repeated 20 times
and we plot the mean of the best objective value plus and minus the standard
error.

###  6.3 Results and Discussion

In this section, we present and discuss our experimental results along
different dimensions.

Figure 1 shows the results for BO performance (best objective value vs. number
of function evaluations / BO iterations) of BOPS-T, BOPS-H, and COMBO on all
six benchmarks. Below we discuss these results in detail.

BOPS-T vs. BOPS-H. Recall that BOPS-T and BOPS-H makes varying trade-offs
between the complexity of statistical model and tractability of acquisition
function optimization: BOPS-T uses simple model and tractable search; and
BOPS-H employs complex model and heuristic search. From Figure 1, we can
observe that BOPS-H performs significantly better than BOPS-T on all six
benchmarks.

BOPS vs. COMBO. From the results shown in Figure 1, we make the following
observations: 1) BOPS-H performs significantly better than both BOPS-T and
COMBO on all six benchmarks; and 2) BOPS-T is comparable or slightly better
than COMBO on all benchmarks except TSP and CP.

We hypothesize that the performance of different BO algorithms, namely,
BOPS-H, BOPS-T, and COMBO is proportional to the quality of their surrogate
models in terms of making predictions on unknown permutations and their
uncertainty quantification ability. To verify this hypothesis, we compare the
three surrogate models quantitatively in terms of their performance on the
log-likelihood metric.

Comparison of surrogate models. We compare the three surrogate models on the
log-likelihood metric (Murphy 2012) because it captures both the prediction
and uncertainty quantification of a model which are essential for the
effectiveness of BO. We plot the negative log-likelihood (NLL) of the three
surrogate models on a testing set of 50 instances as a function of the
increasing size of training data. Each experiment is replicated with 10
different training sets and each method is evaluated using the median of the
NLL metric on 10 different test sets of 50 permutations each. Figure 2 shows
the results on all benchmarks except HMD. We do not show results on HMD since
each function evaluation is much more expensive when compared to all other
benchmarks, and we are generating multiple replications of the training and
testing sets (10 $\times$ 10 = 100 runs). We make the following observations
from Figure 2: 1) BOPS-H shows the best performance among the three surrogate
models; 2) BOPS-T does better than COMBO on all benchmarks other than cell-
placement. Since both COMBO and BOPS-H employ the same acquisition function
(EI) and optimizer (local search), it is evident that the gains in the BO
performance comes from the superior surrogate model of BOPS-H.

##  7 Conclusions

We proposed and evaluated two effective Bayesian optimization algorithms with
varying trade-offs for optimizing expensive black-box functions over the
challenging input space of permutations. The results point to a key conclusion
that it is important to use an appropriate model that exploits the specific
structure of permutation spaces, which is different than the generic
combinatorial space over categorical variables. We characterized the
importance of this problem setting by describing three important real-world
applications from the domain of computer-aided design of integrated circuits.
Furthermore, we make all these benchmarks available to drive future research
in this problem space. Future work includes studying extensions to handle
high-dimensional spaces (Oh, Gavves, and Welling 2018) and multiple objectives
(Belakaria, Deshwal, and Doppa 2019).

Acknowledgments. This research is supported in part by NSF grants IIS-1845922,
OAC-1910213,and SII-2030159.

## References

  * Angermueller et al. (2020) Angermueller, C.; Belanger, D.; Gane, A.; Mariet, Z.; Dohan, D.; Murphy, K.; Colwell, L.; and Sculley, D. 2020.  Population-Based Black-Box Optimization for Biological Sequence Design.  In _International Conference on Machine Learning_ , 324–334. PMLR.
  * Angermueller et al. (2019) Angermueller, C.; Dohan, D.; Belanger, D.; Deshpande, R.; Murphy, K.; and Colwell, L. 2019.  Model-based reinforcement learning for biological sequence design.  In _International Conference on Learning Representations_.
  * Aronszajn (1950) Aronszajn, N. 1950.  Theory of reproducing kernels.  _Transactions of the American mathematical society_ , 68(3): 337–404.
  * BakIr et al. (2007) BakIr, G.; Hofmann, T.; Schölkopf, B.; Smola, A. J.; and Taskar, B. 2007.  _Predicting structured data_.  MIT press.
  * Balandat et al. (2020) Balandat, M.; Karrer, B.; Jiang, D. R.; Daulton, S.; Letham, B.; Wilson, A. G.; and Bakshy, E. 2020.  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  In _Advances in Neural Information Processing Systems 33_.
  * Baptista and Poloczek (2018) Baptista, R.; and Poloczek, M. 2018.  Bayesian Optimization of Combinatorial Structures.  In _Proceedings of the 35th International Conference on Machine Learning_ , 462–471.
  * Belakaria, Deshwal, and Doppa (2019) Belakaria, S.; Deshwal, A.; and Doppa, J. R. 2019.  Max-value entropy search for multi-objective bayesian optimization.  In _Advances in Neural Information Processing Systems (NeurIPS)_.
  * Buathong, Ginsbourger, and Krityakierne (2020) Buathong, P.; Ginsbourger, D.; and Krityakierne, T. 2020.  Kernels over sets of finite sets using rkhs embeddings, with application to bayesian (combinatorial) optimization.  In _International Conference on Artificial Intelligence and Statistics_ , 2731–2741. PMLR.
  * Burkard et al. (1998) Burkard, R. E.; Cela, E.; Pardalos, P. M.; and Pitsoulis, L. S. 1998.  The quadratic assignment problem.  In _Handbook of combinatorial optimization_ , 1713–1809. Springer.
  * Burkard, Karisch, and Rendl (1997) Burkard, R. E.; Karisch, S. E.; and Rendl, F. 1997.  QAPLIB–a quadratic assignment problem library.  _Journal of Global optimization_ , 10(4): 391–403.
  * Che et al. (2009) Che, S.; Boyer, M.; Meng, J.; Tarjan, D.; Sheaffer, J. W.; Lee, S.-H.; and Skadron, K. 2009.  Rodinia: A Benchmark suite for Heterogeneous Computing.  In _IEEE International Symposium on Workload Characterization (IISWC)_ , 44–54. IEEE.
  * Dadkhahi et al. (2020) Dadkhahi, H.; Shanmugam, K.; Rios, J.; Das, P.; Hoffman, S. C.; Loeffler, T. D.; and Sankaranarayanan, S. 2020.  Combinatorial Black-Box Optimization with Expert Advice.  In _ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD)_ , 1918–1927.
  * Davis (1991) Davis, L. 1991.  _Handbook of genetic algorithms_.  CumInCAD.
  * Deshwal, Belakaria, and Doppa (2020) Deshwal, A.; Belakaria, S.; and Doppa, J. R. 2020.  Scalable Combinatorial Bayesian Optimization with Tractable Statistical models.  _CoRR_ , abs/2008.08177.
  * Deshwal, Belakaria, and Doppa (2021) Deshwal, A.; Belakaria, S.; and Doppa, J. R. 2021.  Mercer Features for Efficient Combinatorial Bayesian Optimization.  In _AAAI_ , 7210–7218.
  * Deshwal et al. (2020) Deshwal, A.; Belakaria, S.; Doppa, J. R.; and Fern, A. 2020.  Optimizing discrete spaces via expensive evaluations: A learning to search framework.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 3773–3780.
  * Deshwal and Doppa (2021) Deshwal, A.; and Doppa, J. R. 2021.  Combining Latent Space and Structured Kernels for Bayesian Optimization over Combinatorial Spaces.  In _Thirty-fifth Conference on Neural Information Processing Systems (NeurIPS)_.
  * Deshwal, Simon, and Doppa (2021) Deshwal, A.; Simon, C.; and Doppa, J. R. 2021.  Bayesian optimization of nanoporous materials.  _Molecular Systems Design and Engineering, Royal Society of Chemistry_.
  * Doppa (2021) Doppa, J. R. 2021.  Adaptive Experimental Design for Optimizing Combinatorial Structures.  In _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence (IJCAI)_ , 4940–4945.
  * Eissman et al. (2018) Eissman, S.; Levy, D.; Shu, R.; Bartzsch, S.; and Ermon, S. 2018.  Bayesian optimization and attribute adjustment.  In _Proc. 34th Conference on Uncertainty in Artificial Intelligence_.
  * Ferreira, Khoo, and Singer (2018) Ferreira, J. F. B.; Khoo, Y.; and Singer, A. 2018.  Semidefinite programming approach for the quadratic assignment problem with a sparse graph.  _Computational Optimization and Applications_ , 69(3): 677–712.
  * Frazier and Wang (2016) Frazier, P. I.; and Wang, J. 2016.  Bayesian optimization for materials design.  In _Information Science for Materials Discovery and Design_ , 45–75. Springer.
  * Gardner et al. (2018) Gardner, J. R.; Pleiss, G.; Bindel, D.; Weinberger, K. Q.; and Wilson, A. G. 2018\.  GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration.  In _Advances in Neural Information Processing Systems_.
  * Garrido-Merchán and Hernández-Lobato (2020) Garrido-Merchán, E. C.; and Hernández-Lobato, D. 2020.  Dealing with categorical and integer-valued variables in bayesian optimization with gaussian processes.  _Neurocomputing_ , 380: 20–35.
  * Gómez-Bombarelli et al. (2018) Gómez-Bombarelli, R.; Wei, J. N.; Duvenaud, D.; Hernández-Lobato, J. M.; Sánchez-Lengeling, B.; Sheberla, D.; Aguilera-Iparraguirre, J.; Hirzel, T. D.; Adams, R. P.; and Aspuru-Guzik, A. 2018.  Automatic Chemical Design using a Data-driven Continuous Representation of Molecules.  _ACS central science_ , 4(2): 268–276.
  * Greenhill et al. (2020) Greenhill, S.; Rana, S.; Gupta, S.; Vellanki, P.; and Venkatesh, S. 2020.  Bayesian optimization for adaptive experimental design: A review.  _IEEE Access_ , 8: 13937–13948.
  * Griffiths and Hernández-Lobato (2017) Griffiths, R.-R.; and Hernández-Lobato, J. M. 2017.  Constrained bayesian optimization for automatic chemical design.  _arXiv preprint arXiv:1709.05501_.
  * Hernández-Lobato et al. (2017) Hernández-Lobato, J. M.; Requeima, J.; Pyzer-Knapp, E. O.; and Aspuru-Guzik, A. 2017.  Parallel and distributed Thompson sampling for large-scale accelerated exploration of chemical space.  In _International conference on machine learning_ , 1470–1479. PMLR.
  * Hutter, Hoos, and Leyton-Brown (2011) Hutter, F.; Hoos, H. H.; and Leyton-Brown, K. 2011.  Sequential Model-based Optimization for General Algorithm Configuration.  In _International conference on Learning and Intelligent Optimization_ , 507–523.
  * Jiao and Vert (2015) Jiao, Y.; and Vert, J.-P. 2015.  The Kendall and Mallows kernels for permutations.  In _International Conference on Machine Learning_ , 1935–1944. PMLR.
  * Jiao and Vert (2018) Jiao, Y.; and Vert, J.-P. 2018.  The weighted kendall and high-order kernels for permutations.  _arXiv preprint arXiv:1802.08526_.
  * Kajino (2019) Kajino, H. 2019.  Molecular hypergraph grammar with its application to molecular optimization.  In _International Conference on Machine Learning_ , 3183–3191. PMLR.
  * Kanagawa et al. (2018) Kanagawa, M.; Hennig, P.; Sejdinovic, D.; and Sriperumbudur, B. K. 2018.  Gaussian processes and kernel methods: A review on connections and equivalences.  _arXiv preprint arXiv:1807.02582_.
  * Kendall (1938) Kendall, M. G. 1938.  A new measure of rank correlation.  _Biometrika_ , 30(1/2): 81–93.
  * Kim et al. (2021) Kim, J.; McCourt, M.; You, T.; Kim, S.; and Choi, S. 2021.  Bayesian optimization with approximate set kernels.  _Machine Learning_ , 110(5): 857–879.
  * Kondor and Lafferty (2002) Kondor, R.; and Lafferty, J. 2002.  Diffusion Kernels on Graphs and other Discrete Structures.  In _Proceedings of the 19th International Conference on Machine Learning_ , volume 2002, 315–322.
  * Korovina et al. (2020) Korovina, K.; Xu, S.; Kandasamy, K.; Neiswanger, W.; Poczos, B.; Schneider, J.; and Xing, E. 2020.  Chembo: Bayesian optimization of small organic molecules with synthesizable recommendations.  In _International Conference on Artificial Intelligence and Statistics_ , 3393–3403. PMLR.
  * Lodhi et al. (2002) Lodhi, H.; Saunders, C.; Shawe-Taylor, J.; Cristianini, N.; and Watkins, C. 2002\.  Text classification using string kernels.  _Journal of Machine Learning Research_ , 2(Feb): 419–444.
  * Lu et al. (2018) Lu, X.; Gonzalez, J.; Dai, Z.; and Lawrence, N. D. 2018.  Structured variationally auto-encoded optimization.  In _International conference on machine learning_ , 3267–3275. PMLR.
  * Luong et al. (2019) Luong, P.; Gupta, S.; Nguyen, D.; Rana, S.; and Venkatesh, S. 2019.  Bayesian optimization with discrete variables.  In _Australasian Joint Conference on Artificial Intelligence_ , 473–484. Springer.
  * Mania et al. (2018) Mania, H.; Ramdas, A.; Wainwright, M. J.; Jordan, M. I.; and Recht, B. 2018.  On kernel methods for covariates that are rankings.  _Electronic Journal of Statistics_ , 12(2): 2537–2577.
  * Moss et al. (2020) Moss, H. B.; Beck, D.; González, J.; Leslie, D. S.; and Rayson, P. 2020.  BOSS: Bayesian Optimization over String Spaces.  _arXiv preprint arXiv:2010.00979_.
  * Murphy (2012) Murphy, K. P. 2012.  _Machine learning: a probabilistic perspective_.  MIT press.
  * Oh, Gavves, and Welling (2018) Oh, C.; Gavves, E.; and Welling, M. 2018.  BOCK : Bayesian Optimization with Cylindrical Kernels.  In Dy, J. G.; and Krause, A., eds., _Proceedings of the 35th International Conference on Machine Learning (ICML)_ , volume 80 of _Proceedings of Machine Learning Research_ , 3865–3874. PMLR.
  * Oh et al. (2019) Oh, C.; Tomczak, J.; Gavves, E.; and Welling, M. 2019.  Combinatorial Bayesian Optimization using the Graph Cartesian Product.  In _Advances in Neural Information Processing Systems_ , 2910–2920.
  * Rasmussen and Williams (2006) Rasmussen, C. E.; and Williams, C. K. I. 2006.  _Gaussian processes for machine learning_.  Adaptive computation and machine learning. MIT Press.
  * Reinelt (1995) Reinelt, G. 1995.  Tsplib95.  _Interdisziplinäres Zentrum für Wissenschaftliches Rechnen (IWR), Heidelberg_ , 338: 1–16.
  * Russo and Van Roy (2014) Russo, D.; and Van Roy, B. 2014.  Learning to optimize via posterior sampling.  _Mathematics of Operations Research_ , 39(4): 1221–1243.
  * Sahni and Gonzalez (1976) Sahni, S.; and Gonzalez, T. 1976.  P-complete approximation problems.  _Journal of the ACM (JACM)_ , 23(3): 555–565.
  * Shahriari et al. (2016) Shahriari, B.; Swersky, K.; Wang, Z.; Adams, R. P.; and De Freitas, N. 2016.  Taking the human out of the loop: A review of bayesian optimization.  _Proceedings of the IEEE_ , 104(1): 148–175.
  * Srinivas et al. (2009) Srinivas, N.; Krause, A.; Kakade, S. M.; and Seeger, M. 2009.  Gaussian process optimization in the bandit setting: No regret and experimental design.  _arXiv preprint arXiv:0912.3995_.
  * Swersky et al. (2020) Swersky, K.; Rubanova, Y.; Dohan, D.; and Murphy, K. 2020.  Amortized bayesian optimization over discrete spaces.  In _Conference on Uncertainty in Artificial Intelligence_ , 769–778. PMLR.
  * Varga (2010) Varga, R. S. 2010.  _Geršgorin and his circles_ , volume 36.  Springer Science & Business Media.
  * Yang, Wu, and Arnold (2019) Yang, K. K.; Wu, Z.; and Arnold, F. H. 2019.  Machine-learning-guided directed evolution for protein engineering.  _Nature methods_ , 16(8): 687–694.
  * Zhao et al. (1998) Zhao, Q.; Karisch, S. E.; Rendl, F.; and Wolkowicz, H. 1998.  Semidefinite programming relaxations for the quadratic assignment problem.  _Journal of Combinatorial Optimization_ , 2(1): 71–109.
  * Zhou et al. (2020) Zhou, Z.; Belakaria, S.; Deshwal, A.; Hong, W.; Doppa, J. R.; Pande, P. P.; and Heo, D. 2020.  Design of Multi-Output Switched-Capacitor Voltage Regulator via Machine Learning.  In _DATE_.

##  8 Additional Experiments

(a) QAP (15 dimensions)

(b) TSP (10 dimensions)

(c) Cell Placement

(d) Floorplanning 1

(e) Floorplanning 2

Figure 3: Results comparing BOPS-T, BOPS-H, and COMBO and Genetic Algorithm
(GA) (best objective function value vs. number of BO iterations) on both
synthetic and real-world benchmarks:

Genetic algorithms (Davis 1991) are a commonly employed technique for
optimizing discrete structures. They have also been used for optimizing
acquisition functions (Moss et al. 2020) over discrete spaces in Bayesian
optimization. We present additional results comparing our proposed algorithms
with genetic algorithms in Figure 3. We used a genetic algorithm with
permutation based operators as implemented in Platypus
https://github.com/Project-Platypus/Platypus library. We tried few different
choices for the two hyper-parameters (population size and offspring size) and
picked population size to be 20 and offspring size to 10. Since we are in a
small data setting, large values for these hyper-parameters will not be
practically useful. As seen in the figure, BOPS-H still performs the best
outperforming even genetic algorithms demonstrating the importance
incorporating a model-guided search like Bayesian optimization for finding the
best structure.


1Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di
Milano - Milan, Italy
2ARIANET srl, Via Benigno Crespi 52, Milan, Italy
††thanks: This work was supported by the Italian Ministry of University and
Research (MUR) and the European Union (EU) under the PON/REACT project.

# Back to the Future: GNN-based NO2 Forecasting Via Future Covariates

###### Abstract

Due to the latest environmental concerns in keeping at bay contaminants
emissions in urban areas, air pollution forecasting has been rising the
forefront of all researchers around the world. When predicting pollutant
concentrations, it is common to include the effects of environmental factors
that influence these concentrations within an extended period, like traffic,
meteorological conditions and geographical information. Most of the existing
approaches exploit this information as past covariates, i.e., past exogenous
variables that affected the pollutant but were not affected by it. In this
paper, we present a novel forecasting methodology to predict NO2 concentration
via both past and future covariates. Future covariates are represented by
weather forecasts and future calendar events, which are already known at
prediction time. In particular, we deal with air quality observations in a
city-wide network of ground monitoring stations, modeling the data structure
and estimating the predictions with a Spatiotemporal Graph Neural Network
(STGNN). We propose a conditioning block that embeds past and future
covariates into the current observations. After extracting meaningful
spatiotemporal representations, these are fused together and projected into
the forecasting horizon to generate the final prediction. To the best of our
knowledge, it is the first time that future covariates are included in time
series predictions in a structured way. Remarkably, we find that conditioning
on future weather information has a greater impact than considering past
traffic conditions. We release our code implementation at
https://github.com/polimi-ispl/MAGCRN.

Index Terms—  Air Quality Forecasting, Spatiotemporal Data, Graph Neural
Network, Graph-based Forecasting, Urban Computing

##  1 Introduction

In the landscape of 21st-century environmental concerns, air pollution
forecasting has risen to the forefront, capturing widespread attention at a
global scale. Addressing this challenge, numerous researchers have delved into
the exploration of effective solutions, endeavoring to precisely forecast air
pollutant concentrations through a variety of methods. Among these approaches,
Deep Learning (DL) methods currently hold the predominant interest, marking
significant advancements in the efforts to forecasting air quality parameters
and mitigating the impact of air pollution [1, 2, 3].

Pollutant concentrations exhibit intricate correlations in both temporal and
spatial domains, and these correlations dynamically evolve over time [4, 5,
6]. Recently, DL-based approaches have demonstrated considerable success in
pollutant forecasting by capturing nonlinear temporal and spatial patterns
very effectively [2, 7].

In the last few years, there has been a surge in the popularity of applying
Spatiotemporal Graph Neural Networks (STGNNs) to pollutants forecasting [8, 9,
10, 11, 12, 13, 14] since they can process the data structure by modeling it
as a graph. Indeed, the underlying assumption is to incorporate the data
structure as an inductive bias [15, 16, 17]. Time-evolving observations coming
from ground monitoring stations inherently contains rich spatiotemporal
structure and spatiotemporal dynamics. The dependency structure of the
observations can be captured through pairwise relationships among the
stations. These representations form a graph where each station corresponds to
a node and the functional dependencies among the stations can be seen as
edges.

Most of the existing approaches capture spatial dependencies on a fixed graph
structure, assuming that the underlying relationship between entities is fixed
and pre-determined [18, 19, 20, 12]. However, the explicit graph structure may
not necessarily reflect the true dependency between measurements and existing
relationships may be missing due to incomplete data connections [21, 22].
Recently, graph learning methods have been proposed [23] that allow to learn
the graph structure directly from data. These techniques are promising,
pushing forward the forecasting abilities [21, 24, 25].

In the context of air pollutant forecasting like NO2, PM10 and PM2.5, it is
reasonable to assume that including auxiliary information helps the prediction
performance. This is motivated by the fact that air quality is often
influenced by exogenous variables like traffic and weather conditions. For
example, traffic-related emissions have been one of the top contributors to
air pollution in many cities around the world. It has been proved that these
emissions can deteriorate ambient air quality on a large spatial scale,
especially during the morning and evening rush hours in urban regions [4, 8,
18]. In addition, studies have shown that air pollutants vary under different
meteorological conditions [5, 26]. Indeed, the temperature affects the
atmospheric and ventilation conditions; humidity and precipitation can change
the deposition characteristics of particulate matter; wind speed promotes the
diffusion and spread of pollutants [27].

All this information could be included in the air quality forecasting as
covariates, i.e., exogenous variables that affect the pollutant to predict,
but are not affected by it. For instance, past covariates are represented by
traffic conditions, meteorological factors and geographical information.
Future covariates are all the future-related information that are known in
advance, like weather forecasts and calendar events (seasons, days of a week,
time of a day) [4].

Past covariates have been successfully exploited in [11, 13, 18, 12, 14, 19].
To the best of our knowledge, very few STGNN-based methods that exploit future
covariates in performing pollutant prediction tasks have been proposed [8, 10,
20], but they lack a structural way to consider future information into the
forecasting process.

In this paper, we explore the potential of both past and future covariates to
perform forecasting of NO2 concentration. To do so, we focus on STGNN-based
architectures that automatically learn the graph structure from data. We build
our method upon the STGNN module in [24], and we propose a conditioning block
to incorporate additional past and future information into the forecasting
process. Then, we condense the covariate information through a straightforward
fusion process for predicting NO2 concentration.

When tested on a recently released real-world air quality dataset, the
proposed approach outperforms existing state-of-art forecasting methods,
showing superior performances in predicting air quality parameters up to a
forecasting horizon of three days. Interestingly, the role of future
covariates reveals paramount for enhancing the prediction performances.

##  2 Problem Formulation

We consider a collection of $N$, hourly-sampled, 1-dimensional time series,
acquired from a network of $N$ ground monitoring stations. We denote with the
matrix $\mathbf{X}^{p}_{\\{t_{0},W\\}}\in\mathbb{R}^{N\times W}$ the stacked
$N$ observations acquired for $W$ consecutive past time instants [hours] up to
a final reference time instant $t_{0}$. For instance,
$\mathbf{X}^{p}_{\\{t_{0},W\\}}$ can represent the hourly evolution of NO2
concentration acquired at $N$ different stations in the last $W$ hours.

The aim of this work is to perform a multi-step forecasting, that is,
predicting for each time step in the forecasting horizon $H$ [hours] the
future observation $\mathbf{X}^{f}_{\\{t_{1},H\\}}\in\mathbb{R}^{N\times H}$
starting from a generic time instant $t_{1}\geq t_{0}$ for a time window of
$H$ hours. In general, $H$ can be different from $W$. To face this problem, we
exploit the knowledge of past observations $\mathbf{X}^{p}_{\\{t_{0},W\\}}$
combined with their covariate time series, i.e., the time-dependent exogenous
variables that influence the time series in $\mathbf{X}^{p}_{\\{t_{0},W\\}}$
but are not affected by them.

We refer to the term past covariates as the covariate time series known only
into the past in the same time window of the acquired observations, i.e., in a
temporal window $W$ up to the time instant $t_{0}$. An example of past
covariate of NO2 concentration is represented by the traffic data like the
number of vehicles and their speed, which are strictly linked to air quality
conditions. Also past cyclical events derived from calendar like the daily
time evolution and weekday information can be used as additional covariates
[28, 29]. For simplicity, we denote with $\mathcal{U}^{p}_{\\{t_{0},W\\}}$ the
set of time series associated with $P$ past covariates, being
$\mathcal{U}^{p}_{\\{t_{0},W\\}}=\\{\mathbf{U}^{p}_{i}\in\mathbb{R}^{N\times
W}\\},i\in[1,P]$.

With the term future covariates we refer to the covariate time series known
into the future in the forecasting horizon, i.e., from the time instant
$t_{1}$ for a time window $H$. This is the case for instance of future weather
forecasts, which are known in advance and can provide useful information for
predicting air quality conditions [5, 27]; future calendar events and
daily/weekly time series can be used as well to predict air quality. We define
the set of time series related to $F$ future covariates as
$\mathcal{U}^{f}_{\\{t_{1},H\\}}=\\{\mathbf{U}^{f}_{i}\in\mathbb{R}^{N\times
H}\\},i\in[1,F]$.

We formulate our forecasting problem as finding a function $F(\cdot)$ to
forecast the next $H$ time steps observations based on the knowledge of the
past temporal observations and past and future covariates. Therefore,

$\hat{\mathbf{X}}^{f}_{\\{t_{1},H\\}}=F_{\theta}(\mathbf{X}^{p}_{\\{t_{0},W\\}},\mathcal{U}^{p}_{\\{t_{0},W\\}},\mathcal{U}^{f}_{\\{t_{1},H\\}})$

where $\theta$ denotes all the learnable parameters of the forecasting model.

##  3 Proposed Methodology

As recently done in the literature of time series forecasting [30, 25, 24], we
model the tackled problem as a graph-based forecasting problem. Our network of
ground monitoring stations can be represented with a graph structure, i.e.,
ground stations (graph nodes) that exhibit a certain degree of correlation
(encoded in the graph edges).

In particular, we propose a modified version of the Adaptive Graph
Convolutional Recurrent Network (AGCRN) presented in [24] and successfully
used for predicting the future status of a traffic network from past data.
This network is part of a specific family of Graph Neural Networks (GNNs)
known as STGNNs, i.e., it is characterized by a graph topology that remains
constant, while node and edge features change over time. We decided to get
advantages of this architecture for two main reasons. First, differently from
others STGNNs, AGCRN captures specific spatial and temporal characteristics of
each node. This information improves the representation of station-level
peculiarities and it is crucial since a station’s measurement can greatly
depend on its physical location (urban or rural) [31]. Second, the hidden
graph topology is directly inferred from training data, which usually benefits
the forecasting performance [15].

Fig. 1: The training pipeline of the proposed MAGCRN. The AGCRN∗ used here
represents AGCRN without the last FC layer.

To face our air quality forecasting task, we propose to modify the standard
prediction paradigm of AGCRN by exploiting the knowledge of covariates.
Differently from the original version, our proposed Modified AGCRN (MAGCRN)
leverages the ability of conditioning the forecast not only on past but also
on future information. An example of future information that can influence our
predictions is represented by the future weather forecasts [5, 26]. By
conditioning NO2 forecasting on future weather forecasts, we can build a more
comprehensive and accurate model that takes into account the intricate
relationships between atmospheric conditions and air quality. It is worth
noticing that, in a real scenario, future weather forecasts are available to
the algorithm at prediction time.

The complete training pipeline of our proposed MAGCRN is shown in Fig. 1. We
can recognize three main steps:

  1. 1.

Conditioning, where the input observations $\mathbf{X}^{p}_{\\{t_{0},W\\}}$
are conditioned on past ($\mathcal{U}^{p}_{\\{t_{0},W\\}}$) and future
($\mathcal{U}^{f}_{\\{t_{1},H\\}}$) covariates separately, returning
$\mathbf{C}^{p}$ and $\mathbf{C}^{f}$ signals.

  2. 2.

Embedding estimation, where we extract two embedding signals $\mathbf{E}^{p}$
and $\mathbf{E}^{f}$ from the conditioned signals $\mathbf{C}^{p}$ and
$\mathbf{C}^{f}$.

  3. 3.

Past & Future fusion, where we combine the $\mathbf{E}^{p}$ and
$\mathbf{E}^{f}$ embeddings to construct the final prediction
$\hat{\mathbf{X}}^{f}_{\\{t_{1},H\\}}$.

Conditioning.  The most important part of our system is the fully-learnable
input conditioning module named $\operatorname{Cond}$, which conditions the
current observed data $\mathbf{X}^{p}_{\\{t_{0},W\\}}$ on both past and future
covariates. Given two generic data matrices $\mathbf{K}_{1}$ and
$\mathbf{K}_{2}$, the proposed conditioning module takes inspiration from
[32], and can be computed as

$\operatorname{Cond}(\mathbf{K}_{1},\mathbf{K}_{2})=\operatorname{MLP}(\phi(\operatorname{MLP}(\mathbf{K}_{1})))+\operatorname{MLP}(\phi(\operatorname{MLP}(\mathbf{K}_{2})))$

where the $\operatorname{MLP}$ module is based on a Multilayer Perceptron
(MLP) performing a linear transformation on the input data and $\phi$ is a
non-linearity carried out by a Rectified Linear Unit (ReLU). On the one hand,
the input data are conditioned on past covariates, obtaining
$\mathbf{C}^{p}=\operatorname{Cond}(\mathbf{X}^{p}_{\\{t_{0},W\\}},\mathcal{U}^{p}_{\\{t_{0},W\\}})$;
on the other hand, the same input is conditioned on future covariates,
obtaining
$\mathbf{C}^{f}=\operatorname{Cond}(\mathbf{X}^{p}_{\\{t_{0},W\\}},\mathcal{U}^{f}_{\\{t_{1},H\\}})$.
In particular, $\mathbf{C}^{p}\in\mathbb{R}^{N\times W\times Z}$ and
$\mathbf{C}^{f}\in\mathbb{R}^{N\times H\times Z}$, where $Z$ is the output
dimension introduced by the MLP-based conditioning modules. The structure of
the two $\operatorname{Cond}$ blocks is the same; however, the four involved
$\operatorname{MLP}$ blocks are different in the sense that their parameters
are learned independently from each other.

Fig. 2: The proposed $\operatorname{Cond}$ module.

In Fig. 2, we report an example of the $\operatorname{Cond}$ module in
conditioning the input window $\mathbf{X}^{p}_{\\{t_{0},W\\}}$ with past
information $\mathcal{U}^{p}_{\\{t_{0},W\\}}$.

Embedding estimation.  After the conditioning, $\mathbf{C}^{p}$ and
$\mathbf{C}^{f}$ are given as input to two identical AGCRN cells to
automatically learn the latent graph structure and, at the same time, capture
node-specific patterns. To perform a late fusion between past and future
conditioning contributions, these cells are modified by removing the last
Fully Connected (FC) layer. We define the outputs of the AGCRN cells as past
and future embeddings, $\mathbf{E}^{p}$ and $\mathbf{E}^{f}$, respectively,
with $\mathbf{E}^{p},\mathbf{E}^{f}\in\mathbb{R}^{N\times Z}$. Notice that the
temporal dimension has been embedded into the last features’ dimension.

Past & Future fusion.  The embeddings are then fused together to obtain the
final feature tensor $\mathbf{O}\in\mathbb{R}^{N\times Z}$, defined as

$\mathbf{O}=(1-\alpha)\mathbf{E}^{p}+\alpha\mathbf{E}^{f}$

where $\alpha$ is a mixing coefficient determining the balance between the
past and future information. If $\alpha=0$, the output is entirely determined
by past conditioning; if $\alpha=1$ the output relies only on future
conditioning. The resulting $\mathbf{O}$ vector is mapped into the forecasting
horizon $H$ by a linear FC layer, obtaining an estimation of the future NO2
concentration $\hat{\mathbf{X}}^{f}_{\\{t_{1},H\\}}\in\mathbb{R}^{N\times H}$.

The network is trained minimizing the Mean Absolute Error (MAE) between the
groundtruth $\mathbf{X}^{f}_{\\{t_{1},H\\}}$ and its forecast
$\hat{\mathbf{X}}^{f}_{\\{t_{1},H\\}}$.

In testing phase, given $\mathbf{X}^{p}_{\\{t_{0},W\\}}$ along with its
related past ($\mathcal{U}^{p}_{\\{t_{0},W\\}}$) and future
($\mathcal{U}^{f}_{\\{t_{1},H\\}}$) covariates, we exploit the trained MAGCRN
to perform the forecast, obtaining $\hat{\mathbf{X}}^{f}_{\\{t_{1},H\\}}$.

##  4 Experimental Setup

Fig. 3: The distribution of the 24 ground monitoring stations placed in the
city of Madrid. Each point represents one observation station.

Dataset. We conduct experiments on a recently released real-world dataset
presented in [18]. This dataset is composed of hourly-sampled air quality,
meteorological and traffic data from 24 different ground monitoring stations
in the city of Madrid (Spain), from January to June 2019. The station spatial
distribution is reported in Fig. 3. The released data include:

  * •

Air quality data — NO2 concentration $[\mu$g/m${}^{3}]$;

  * •

Meteorological data — wind speed [m/s], wind direction [rad], temperature
[°C], relative humidity [%], barometric pressure [mb], solar irradiance
[W/m2];

  * •

Traffic data — intensity [vehicles/h], occupancy time [%], load (degree of
congestion) [%], average traffic speed [km/h].

Our target data consists in the NO2 concentration. We consider as a past
covariates the measurements related to the traffic conditions. Indeed, this
information has been proven to be strongly correlated with air quality
variations [4, 33]. As for future covariates, we consider all the
meteorological data related to the forecasting horizon. In addition, calendar-
derived covariates are included in both past a future covariates set since
they are known in advance [4].

Network baselines.  In this work, we focus on STGNNs that directly learn the
adjacency matrix of the graph nodes. To compare our proposed methodology with
the state-of-the-art, we select three recently proposed STGNNs for time series
forecasting:

  * •

Graph WaveNet [21], which integrates diffusion graph convolutions with dilated
convolutions to capture spatial and temporal dependencies;

  * •

Gated GNN (GGNN) [34], that integrates gated mechanisms and residual
connections for effective information propagation;

  * •

AGCRN [24], the reference architecture which our method is built upon.

Notice that none of these STGNNs directly supports future covariates,
contrarily to our proposed methodology.

Moreover, we do not compare ourselves with state-of-the-art models exploiting
future covariates since they address the problem in a different manner. The
work in [8] adopts a hierarchical GNN by constructing a three-level hierarchy,
from cities to regions. Our data do not allow for this granularity as they
regard a more limited area (city-wide). Similarly, in [10] the authors focus
on regional-level predictions.

Training and testing specifications.  To train and test the selected STGNNs,
we split the dataset in chronological order by the ratio of 7:2:1 to generate
training, validation, and test data.

To train our model, we construct consecutive input/output windows pairs with
$W=24$ past information and $H=24$ future information, resulting in 3820
window pairs. We train our model in predicting 1-day ahead NO2 concentration,
thus we consider $t_{1}=t_{0}$. The employed MLP in the $\operatorname{Cond}$
module consists in a FC linear layer, with $Z=64$. Moreover, we equally weight
past and future covariates’ contributions, setting $\alpha=0.5$. We combine
the ADAM optimizer with the cosine annealing learning rate scheduler, setting
the maximum and minimum learning rate to 10-2 and 10-7 respectively, with a
maximum period of 20. We stop the learning phase if the validation error does
not decrease after 100 epochs.

At testing stage we always evaluate a forecasting horizon $H=24$, considering
429 window pairs. However, in addition to the 1-day ($t_{1}=t_{0}$)
forecasting horizon, we also evaluate the model in performing 2-day
($t_{1}=t_{0}+24$) and 3-day ($t_{1}=t_{0}+48$) hourly NO2 predictions without
retraining the network. We do this for testing the robustness and the
generalization capabilities of the proposed method to predict farther events
in the future. Then, we evaluate the model performance in terms of MAE, Root
Mean Squared Error (RMSE) and Mean Relative Error (MRE), i.e., the MAE
normalized by the $\ell_{1}$-norm of the target window.

##  5 Experimental Results

Table 1: Forecasting error by means of MAE, RMSE, MRE, considering the
baseline networks and our proposed MAGCRN. We can exploit two different groups
of covariates (i.e., past $\mathcal{U}^{p}$ and future $\mathcal{U}^{f}$) for
three different forecasting horizons, i.e., 1-24h (1-day), 25-48h (2-day),
49-72h (3-day). Best results in bold.

Model |  |  | Forecasting Error (1-day / 2-day / 3-day)
---|---|---|---
$\mathcal{U}^{p}_{\\{t_{0},W\\}}$ | $\mathcal{U}^{f}_{\\{t_{1},H\\}}$ | MAE | RMSE | MRE
Graph WaveNet | ✓ |  | 12.33 / 13.51 / 14.16 | 19.41 / 20.59 / 21.37 | 0.48 / 0.52 / 0.55
GGNN | ✓ |  | 11.53 / 13.93 / 14.01 | 18.15 / 20.75 / 20.74 | 0.44 / 0.54 / 0.54
AGCRN | ✓ |  | 11.56 / 13.17 / 13.69 | 18.31 / 20.14 / 20.70 | 0.45 / 0.51 / 0.53
MAGCRN (ours) | ✓ | ✓ | 9.43 / 10.30 / 10.85 | 14.64 / 15.56 / 16.19 | 0.36 / 0.40 / 0.42

General results.  Tab. 1 reports the overall results of MAGCRN and the
compared baselines in terms of MAE, RMSE and MRE for three different
forecasting horizons. As expected, the performances decrease as the
forecasting horizon moves far from the 1-day predictions. Nonetheless, our
approach, which jointly considers past and future covariates in learning the
graph structure and performing the desired forecast, consistently outperforms
all the baselines.

Since NO2 and air pollutants in general strongly depends on weather conditions
[5, 6], we believe that considering the future weather forecast as an
additional conditioning information benefits the prediction. The combination
of two different conditioning (past and future) brings more accurate
predictions, showing the inherent dependence that exists between the pollutant
and weather conditions.

Parameter Sensitivity.  In this section, we investigate the sensitivity of the
most important parameters of the proposed MAGCRN for our forecasting task.
Specifically, we explore the role of: (i) the embedding size $Z$; (ii) the
$\alpha$ factor used to balance the action between past and future
conditioning. For the sake of brevity, we only report the achieved MAE as a
function of the analyzed parameters, for different forecasting horizons. Each
time we vary a parameter, we set the other to its default value, i.e., $Z=64$
and $\alpha=0.5$. The results are reported in Fig. 4.

Fig. 4: MAE results vs network parameters, according to different time
forecasting horizons. In particular, $\blacktriangledown$ indicates 1-day time
horizon, $\blacktriangle$ the 2-day time horizon, $\blacktriangleleft$ the
3-day time horizon.

Fig. 4 (a) depicts the performances by varying $Z$, that directly influences
the generation of node-specific parameters in the STGNN. This parameter proves
crucial for the near future, i.e., in the 1-day forecasting horizon, with an
optimal value of $64$. On the contrary, $Z$ does not seem to impact more
distant forecasting horizons (especially the 3-day case), since the trend is
approximately stable for all the tested values.

In Fig. 4 (b) we vary the $\alpha$ factor used to balancing the action between
the past and future conditioning. A lower $\alpha$ promotes the conditioning
of past covariates in generating the forecast, up to the extreme case where
the model considers only past covariates $\mathcal{U}^{p}_{\\{t_{0},W\\}}$ ,
i.e., $\alpha=0$. On the other hand, a higher $\alpha$ favors the action of
the future conditioning, up to the extreme case in which $\alpha=1$, where we
discard $\mathcal{U}^{p}_{\\{t_{0},W\\}}$ and consider only
$\mathcal{U}^{f}_{\\{t_{1},H\\}}$ in generating the forecast. The value
$\alpha=0.5$ gives equal importance to past and future covariates in the
prediction process.

It is possible to notice a particular behavior, consistent among all the
analyzed forecasting horizons: assigning equal importance to past and future
information does not lead to the best results. The weather forecasts seem to
play a paramount role in the NO2 prediction. Indeed, the best performance
corresponds to giving more importance ($\alpha=0.8$) to the future weather
conditions rather than to the past information, which is provided by the
traffic conditions. The extreme case in which $\alpha=1$, i.e., considering
only the weather forecasts as covariates, does not bring the best results.
Nonetheless, for the forecasting horizons that are distant to the one adopted
in training (2-day and 3-day scenarios), the future conditioning becomes even
more crucial to achieve acceptable results.

##  6 Conclusions

This paper introduces an innovative forecasting methodology for predicting NO2
time series concentrations by effectively leveraging past and future
covariates. Differently from the solutions proposed in the literature, this is
the first instance of future covariates being systematically included in time
series predictions.

The proposed study addresses air quality observations in a city-wide network
of ground monitoring stations, utilizing a Spatiotemporal Graph Neural Network
(STGNN) to model the data structure and make predictions. We incorporate past
covariates (e.g., traffic conditions) and future covariates (e.g., weather
forecasts and upcoming calendar events) through an input conditioning module.
Then, we extract spatiotemporal representations of NO2 observations
conditioned to past and future covariates. We integrate these representations
and project them into the forecasting horizon for the final forecast. Our
findings emphasize the importance of including future information, especially
for longer forecasting horizons, surpassing the impact of past traffic
conditions.

Future work will delve more deeply into the role of the coefficient balancing
past and future covariates contributions, enabling it to be learned in the
training phase for enhanced performances. Moreover, we will investigate
potential improvements to the network architecture, integrating future
covariates into the network’s message-passing process and exploring dynamic
graph structures.

##  7 References

## References

  * [1] Konstantinos Benidis et al. “Deep Learning for Time Series Forecasting: Tutorial and Literature Survey” In _ACM Comput. Surv._ 55.6 ACM, 2022
  * [2] Guangyin Jin et al. “Spatio-Temporal Graph Neural Networks for Predictive Learning in Urban Computing: A Survey” In _arXiv_ , 2023 DOI: 10.48550/arxiv.2303.14483
  * [3] Jesus Rodrigo Cedeno Jimenez, Angelly de Jesus Pugliese Viloria and Maria Antonia Brovelli “Estimating Daily NO2 Ground Level Concentrations Using Sentinel-5P and Ground Sensor Meteorological Measurements” In _ISPRS International Journal of Geo-Information_ 12.3, 2023 DOI: 10.3390/ijgi12030107
  * [4] Bo Zhang et al. “Deep learning for air pollutant concentration prediction: A review” In _Atmospheric Environment_ 290, 2022, pp. 119347 DOI: 10.1016/j.atmosenv.2022.119347
  * [5] Daniel J. Jacob and Darrell A. Winner “Effect of climate change on air quality” Atmospheric Environment In _Atmospheric Environment_ 43.1, 2009, pp. 51–63 DOI: 10.1016/j.atmosenv.2008.09.051
  * [6] Yansui Liu, Yang Zhou and Jiaxin Lu “Exploring the relationship between air pollution and meteorological conditions in China under environmental governance” In _Scientific Reports_ 10, 2020 DOI: 10.1038/s41598-020-71338-7
  * [7] Rodrigo de Medrano, Víctor de Buen Remiro and José L. Aznarte “SOCAIRE: Forecasting and monitoring urban air quality in Madrid” In _Environmental Modelling & Software_ 143, 2021, pp. 105084 DOI: https://doi.org/10.1016/j.envsoft.2021.105084
  * [8] Jindong Han, Hao Liu, Haoyi Xiong and Jing Yang “Semi-Supervised Air Quality Forecasting via Self-Supervised Hierarchical Graph Neural Network” In _IEEE Transactions on Knowledge and Data Engineering_ 35.5, 2023, pp. 5230–5243 DOI: 10.1109/TKDE.2022.3149815
  * [9] Xue-Bo Jin et al. “Deep Spatio-Temporal Graph Network with Self-Optimization for Air Quality Prediction” In _Entropy_ 25.2, 2023 DOI: 10.3390/e25020247
  * [10] Jing Xu et al. “Dynamic graph neural network with adaptive edge attributes for air quality prediction: A case study in China” In _Heliyon_ 9.7, 2023
  * [11] Jindong Han, Hao Liu, Hengshu Zhu and Hui Xiong “Kill Two Birds With One Stone: A Multi-View Multi-Adversarial Learning Approach for Joint Air Quality and Weather Prediction” In _IEEE Transactions on Knowledge and Data Engineering_ 35.11, 2023, pp. 11515–11528 DOI: 10.1109/TKDE.2023.3236423
  * [12] Thu Hang Phung et al. “Unsupervised Air Quality Interpolation with Attentive Graph Neural Network” In _Proceedings of the 11th International Symposium on Information and Communication Technology_ , SoICT’22 ACM, 2022, pp. 103–110
  * [13] Chunyang Wang et al. “Modeling Inter-Station Relationships with Attentive Temporal Graph Convolutional Network for Air Quality Prediction” In _Proceedings of the 14th ACM International Conference on Web Search and Data Mining_ , WSDM ’21 ACM, 2021, pp. 616–634
  * [14] Yanlin Qi, Qi Li, Hamed Karimian and Di Liu “A hybrid model for spatiotemporal forecasting of PM2.5 based on graph convolutional neural network and long short-term memory” In _Science of The Total Environment_ 664, 2019, pp. 1–10 DOI: 10.1016/j.scitotenv.2019.01.333
  * [15] Andrea Cini, Ivan Marisca, Daniele Zambon and Cesare Alippi “Graph Deep Learning for Time Series Forecasting” In _arXiv_ , 2023 DOI: 10.48550/arxiv.2310.15978
  * [16] Antonio Longa et al. “Graph Neural Networks for temporal graphs: State of the art, open challenges, and opportunities” In _arXiv_ , 2023 DOI: 10.48550/arxiv.2302.01018
  * [17] Geert Leus et al. “Graph Signal Processing: History, development, impact, and outlook” In _IEEE Signal Processing Magazine_ 40.4, 2023, pp. 49–60 DOI: 10.1109/MSP.2023.3262906
  * [18] Ditsuhi Iskandaryan, Francisco Ramos and Sergio Trilles “Graph Neural Network for Air Quality Prediction: A Case Study in Madrid” In _IEEE Access_ 11, 2023, pp. 2729–2742 DOI: 10.1109/ACCESS.2023.3234214
  * [19] Peijiang Zhao and Koji Zettsu “MASTGN: Multi-Attention Spatio-Temporal Graph Networks for Air Pollution Prediction” In _2020 IEEE International Conference on Big Data_ , 2020, pp. 1442–1448 DOI: 10.1109/BigData50022.2020.9378156
  * [20] Shuo Wang et al. “PM2.5-GNN: A Domain Knowledge Enhanced Graph Neural Network For PM2.5 Forecasting” In _Proceedings of the 28th International Conference on Advances in Geographic Information Systems_ , SIGSPATIAL ’20 ACM, 2020, pp. 163–166
  * [21] Zonghan Wu et al. “Graph Wavenet for Deep Spatial-Temporal Graph Modeling” In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_ , IJCAI’19, 2019
  * [22] Ming Jin et al. “A Survey on Graph Neural Networks for Time Series: Forecasting, Classification, Imputation, and Anomaly Detection” In _arXiv_ , 2023 DOI: 10.48550/arxiv.2307.03759
  * [23] Andrea Cini, Daniele Zambon and Cesare Alippi “Sparse Graph Learning from Spatiotemporal Time Series” In _Journal of Machine Learning Research_ 24.242, 2023, pp. 1–36
  * [24] Lei Bai et al. “Adaptive Graph Convolutional Recurrent Network for Traffic Forecasting” In _Proceedings of the 34th International Conference on Neural Information Processing Systems_ , NeurIPS’20, 2020
  * [25] Zonghan Wu et al. “Connecting the Dots: Multivariate Time Series Forecasting with Graph Neural Networks” In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, KDD ’20 ACM, 2020, pp. 753–763
  * [26] J.-M. Bertrand et al. “Technical note: Improving the European air quality forecast of the Copernicus Atmosphere Monitoring Service using machine learning techniques” In _Atmospheric Chemistry and Physics_ 23.9, 2023, pp. 5317–5333 DOI: 10.5194/acp-23-5317-2023
  * [27] Daniel J. Jacob “Introduction to Atmospheric Chemistry” Princeton University Press, 1999
  * [28] M.W Gardner and S.R Dorling “Meteorologically adjusted trends in UK daily maximum surface ozone concentrations” In _Atmospheric Environment_ 34.2, 2000, pp. 171–176 DOI: 10.1016/S1352-2310(99)00315-5
  * [29] Massimo Stafoggia et al. “Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013–2015, using a spatiotemporal land-use random-forest model” In _Environment International_ 124, 2019, pp. 170–179
  * [30] Binqing Wu and Ling Chen “DSTCGCN: Learning Dynamic Spatial-Temporal Cross Dependencies for Traffic Forecasting” In _arXiv_ , 2023 DOI: 10.48550/arxiv.2307.00518
  * [31] Andrea Cini, Ivan Marisca, Daniele Zambon and Cesare Alippi “Taming Local Effects in Graph-based Spatiotemporal Forecasting” In _Proceedings of the 37th International Conference on Neural Information Processing Systems_ , NeurIPS’23, 2023
  * [32] Andrea Cini and Ivan Marisca “Torch Spatiotemporal”, 2022
  * [33] Chao Zhang, Yu Zheng, Xiuli Ma and Jiawei Han “Assembler: Efficient Discovery of Spatial Co-Evolving Patterns in Massive Geo-Sensory Data” In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2015
  * [34] Victor Garcia Satorras, Syama Sundar Rangapuram and Tim Januschowski “Multivariate Time Series Forecasting with Latent Graph Inference” In _arXiv_ arXiv, 2022 DOI: 10.48550/arXiv.2203.03423