#  LEO: Learning Efficient Orderings for Multiobjective Binary Decision
Diagrams

Rahul Patel, Elias B. Khalil

###### Abstract

Approaches based on Binary decision diagrams (BDDs) have recently achieved
state-of-the-art results for multiobjective integer programming problems. The
variable ordering used in constructing BDDs can have a significant impact on
their size and on the quality of bounds derived from relaxed or restricted
BDDs for single-objective optimization problems. We first showcase a similar
impact of variable ordering on the Pareto frontier (PF) enumeration time for
the multiobjective knapsack problem, suggesting the need for deriving variable
ordering methods that improve the scalability of the multiobjective BDD
approach. To that end, we derive a novel parameter configuration space based
on variable scoring functions which are linear in a small set of interpretable
and easy-to-compute variable features. We show how the configuration space can
be efficiently explored using black-box optimization, circumventing the curse
of dimensionality (in the number of variables and objectives), and finding
good orderings that reduce the PF enumeration time. However, black-box
optimization approaches incur a computational overhead that outweighs the
reduction in time due to good variable ordering. To alleviate this issue, we
propose LEO, a supervised learning approach for finding efficient variable
orderings that reduce the enumeration time. Experiments on benchmark sets from
the knapsack problem with 3-7 objectives and up to 80 variables show that LEO
is $\sim$30-300% and $\sim$10-200% faster at PF enumeration than common
ordering strategies and algorithm configuration. Our code and instances are
available at https://github.com/khalil-research/leo.

## Introduction

In many real-world scenarios, one must jointly optimize over a set of
conflicting objectives. For instance, solving a portfolio optimization problem
in finance requires simultaneously minimizing risk and maximizing return. The
field of multiobjective optimization deals with solving such problems. It has
been successfully applied in novel drug design (Lambrinidis and Tsantili-
Kakoulidou 2021), space exploration (Song et al. 2018; Tangpattanakul,
Jozefowiez, and Lopez 2012), administrating radiotherapy (Yu et al. 2000),
supply chain network design (Altiparmak et al. 2006), among others. In this
paper, we specifically focus on multiobjective integer programming (MOIP),
which deals with solving multiobjective problems with integer variables and
linear constraints.

The goal of solving a multiobjective problem is to find the Pareto frontier
(PF): the set of feasible solutions that are not dominated by any other
solution, i.e., ones for which improving the value of any objective
deteriorates at least one other objective. The PF solutions provide the
decision-maker with a set of trade-offs between the conflicting objectives.
Objective-space search methods iteratively solve multiple related single-
objective problems to enumerate the PF but suffer from redundant computations
in which previously found solutions are encountered again or a single-
objective problem turns out infeasible. On the other hand, decision-space
search methods leverage branch-and-bound. Unlike the single-objective case
where one compares a single scalar bound (e.g., in mixed-integer linear
programming (MIP)), one needs to compare bound sets to decide if a node can be
pruned; this in itself is quite challenging. Additionally, other crucial
components of branch-and-bound such as branching variable selection and
presolve are still underdeveloped, limiting the usability of this framework.

Binary decision diagrams (BDDs) have been a central tool in program
verification and analysis (Bryant 1992; Wegener 2000). More recently, however,
they have been used to solve discrete optimization problems (Bergman et al.
2016b, a) that admit a recursive formulation akin to that of dynamic
programming. BDDs leverage this structure to get an edge over MIP by
efficiently encoding the feasible set into a network model which enables fast
optimization. To the best of our knowledge, Bergman and Cire (2016) were the
first to use BDDs to solve multiobjective problems, achieving state-of-the-art
results for a number of problems. The VO used to construct a BDD has a
significant impact on its size and consequently any optimization of the
diagram. However, the VO problem within BDD-based MOIP has not been addressed
in the literature. We address this gap by designing a novel learning-based BDD
VO technique for faster enumeration of PF.

Figure 1: Searching for variable orderings. The robot represents any automated
search procedure that tries to find a good variable ordering (Approach #1) or
property weights (Approach #2) by iteratively proposing a variable ordering,
measuring its run time on one or more instances, and repeating until a
satisfactory ordering has been found. Approach #1 depicts a naive approach
that searches for an efficient variable ordering directly in the variable
ordering space. The proposed Approach #2 conducts an indirect search by
controlling a small set of variable property weights (PWs).

We begin with the following hypothesis: VO has an impact on the PF enumeration
time and an “efficient” VO can reduce it significantly. Following an empirical
validation of this hypothesis, we show that such orderings can be found using
black-box optimization, not directly in the (exponentially large) space of
variable orderings, but rather indirectly in the space of constant-size
variable scoring functions. The scoring function is a weighted linear
combination of a fixed set of variable properties (or attributes), and the
indirect search is in the space of possible weight combinations. Figure 1
illustrates how a search in the property-weight space can alleviate the curse
of dimensionality to make VO search scalable to problems with many variables.
However, solving the black-box optimization problem may be prohibitively time-
consuming for any one instance. For variable ordering to be useful in
practice, the time required to produce a good VO should be negligible relative
to the actual PF enumeration time. To alleviate this issue, we train a
supervised machine learning (ML) model on the orderings collected using black-
box optimization. A trained model can then be used on unseen (test) instances
to predict variable orderings. Should such a model generalize well, it would
lead to reduced PF enumeration times. We refer to our approach as LEO
(Learning Efficient Orderings). Our key contributions can be summarized as
follows:

  1. 1.

We show that variable ordering can have a dramatic impact on solving times
through a case study of the multiobjective knapsack, a canonical combinatorial
problem.

  2. 2.

We show how black-box optimization can be leveraged to find efficient variable
orderings at scale.

  3. 3.

We design a supervised learning framework for predicting variable orderings
which are obtained with black-box optimization on a set of training instances.
Our ML models are invariant to permutations of variables and independent of
the number of variables, enabling fast training and the use of one ML model
across instances of different sizes.

  4. 4.

We perform an extensive set of experiments on the knapsack problem and show
that LEO is $\sim$30-300% and $\sim$10-200% faster than the best non-ML
ordering strategy and the SMAC algorithm configuration tool, respectively.

  5. 5.

We perform a feature importance analysis of the best class of ML models we
have found, extreme gradient boosted ranking trees. The analysis reveals that:
(a) a single ML model can be trained across instances with varying numbers of
objectives and variables; and (b) a single knapsack-specific feature that we
had not initially contemplated performs reasonably well on its own, though far
worse than our ML models.

## Preliminaries

### Multiobjective Optimization

An MOIP problem takes the form
$\mathcal{M}:=\min_{x}\left\\{\bar{z}(x):x\in\mathcal{X},\mathcal{X}\subset\mathbb{Z}^{n}_{+}\right\\}$,
where $x$ is the decision vector, $\mathcal{X}$ is a polyhedral feasible set,
and $\bar{z}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}$ a vector-valued
objective function representing the $p$ objectives. In this work, we focus on
the knapsack problem with binary decision variables, hence
$\mathcal{X}\subset\\{0,1\\}^{n}$.

Definition: Pareto dominance. Let
$x^{1},x^{2}\in\mathcal{X},\bar{y}^{1}=\bar{z}(x^{1}),\bar{y}^{2}=\bar{z}(x^{2})$,
then $\bar{y}^{1}$ dominates $\bar{y}^{2}$ if
$\bar{y}^{1}_{j}\leq\bar{y}^{2}_{j},\forall j\in[p]$ and $\exists
j\in[p]:\bar{y}^{1}_{j}<\bar{y}^{2}_{j}$.

Definition: Efficient set. A solution $x^{1}\in\mathcal{X}$ is called an
efficient solution if $\nexists~{}x^{2}\in\mathcal{X}$ such that $x^{2}$
dominates $x^{1}$. The set of all efficient solutions to a multiobjective
problem is called an efficient set $\mathcal{X}_{E}$.

Definition: Pareto frontier (PF). The set of images of the efficient solutions
in the objective space, i.e.,
$\mathcal{Z}_{N}=\\{\bar{z}(x):x\in\mathcal{X}_{E}\\}$, is called the PF.

The exact solution approaches to solving multiobjective problems focus on
efficiently enumerating its PF.

### BDDs for Multiobjective Optimization

A BDD is a compact encoding of the feasible set of a combinatorial
optimization problem that exploits the recursive formulation of the problem.
Formally, a BDD is a layered acyclic graph
$G=(n,\mathcal{N},\mathcal{A},\ell,d)$, composed of nodes in $\mathcal{N}$,
arcs in $\mathcal{A}$, a node-level mapping $\ell:\mathcal{N}\rightarrow[n+1]$
that maps each node to a decision variable, and an arc-label mapping
$d:\mathcal{A}\rightarrow\\{0,1\\}$ that assigns a value from the variable
domain of an arc’s starting node to the arc. Here, $n$ is the number of
variables of the multiobjective problem $\mathcal{M}$.

The nodes are partitioned into $n+1$ layers $L_{1},\dots,L_{n+1}$, where
$L_{l}=\\{u:\ell(u)=l,u\in\mathcal{N}\\}$. The first and last layers have one
node each, the root $\mathbf{r}$ and terminal nodes $\mathbf{t}$,
respectively. The width of a layer $L_{l}$ is equal to the number of nodes in
that layer, $|L_{l}|$. The width of a BDD $G$ is equal to the maximum layer
width, $\max_{l\in[n+1]}|L_{l}|$. An arc $a:=(r(a),t(a))\in\mathcal{A}$ starts
from the node $r(a)\in L_{l}$ and ends in node $t(a)\in L_{l+1}$ for some
$l\in[n]$. It has an associated label $d(a)\in\\{0,1\\}$ and a vector of
values $\bar{v}(a)\in\mathbb{R}^{p}_{+}$ that represents the contribution of
that arc to the $p$ objective values.

Let $\mathcal{P}$ represent all the paths from the root to the terminal node.
A path $e(a_{1},a_{2},\dots,a_{n})\in\mathcal{P}$ is equal to the solution
$x=(d(a_{1}),d(a_{2}),\dots,d(a_{n}))$ and the corresponding image of this
solution in the objective space is given by
$\bar{v}(e)=\sum_{i=1}^{n}\bar{v}(a_{i})$, where the sum is taken elementwise
over the elements of each vector $\bar{v}(a_{i})$. The BDD representation of
$\mathcal{M}$ is valid if
$\mathcal{Z}_{N}=\texttt{ND}\left(\bigcup_{e\in\mathcal{P}}\bar{v}(e)\right)$,
where $\mathcal{Z}_{N}$ is the PF of $\mathcal{M}$ and ND an operator to
filter out dominated objective vectors. We refer the readers to Bergman et al.
(2021) and Bergman and Cire (2016) for the detailed description of the
multiobjective knapsack BDD construction. In what follows, we assume access to
a BDD construction and PF enumeration procedure and will focus our attention
on the variable ordering aspect.

## Related Work

### Exact Approaches for Solving MOIP Problems

Traditional approaches to exactly solve MOIP can be divided into objective-
space search and decision-space search methods. The objective-space search
techniques (Kirlik and Sayın 2014; Boland, Charkhgard, and Savelsbergh 2014,
2015; Ozlen, Burton, and MacRae 2014) enumerate the PF by searching in the
space of objective function values. They transform a multiobjective problem
into a single one either by weighted sum aggregation of objectives or
transforming all but one objective into constraints. The decision-space search
approaches (Przybylski and Gandibleux 2017; Sourd and Spanjaard 2008; Parragh
and Tricoire 2019; Adelgren and Gupte 2022; Belotti, Soylu, and Wiecek 2013)
instead explore the set of feasible decisions. Both these approaches have
their own set of challenges as described in the introduction. We point the
reader to (Ehrgott, Gandibleux, and Przybylski 2016; Ehrgott 2006) for a more
detailed background.

Bergman et al. (2021) showcased that BDD-based approaches can leverage problem
structure and can be orders of magnitude faster on certain problem classes;
the use of valid network pruning operations along with effective network
compilation techniques were the prime factors behind their success. However,
they do not study the effect of VO on the PF enumeration time.

### Variable Ordering for BDD Construction

Finding a VO that leads to a BDD with a minimal number of nodes is an NP-
Complete problem (Bollig, Löbbing, and Wegener 1996; Bollig and Wegener 1996).
This problem has received significant attention from the verification
community as smaller BDDs are more efficient verifiers. Along with the size,
VO also affects the objective bounds; specifically, smaller exact BDDs are
able to obtain better bounds on the corresponding limited width
relaxed/restricted BDDs (Bergman et al. 2016b).

The VO techniques for BDD construction can be broadly categorized as exact or
heuristic. The exact approaches to VO (Friedman and Supowit 1987; Ebendt and
Drechsler 2007; Ebendt, Gunther, and Drechsler 2005; Bergman et al. 2012),
though useful for smaller cases, are not able to scale with problem size. To
alleviate the issue of scalability for larger problems, heuristic VO
techniques are used; the proposed methodology falls into this category. These
heuristics can be general or problem-specific (Rice and Kulhari 2008; Bergman
et al. 2016b; Lu et al. 2000; Aloul, Markov, and Sakallah 2003; Chung, Hajj,
and Patel 1993; Butler et al. 1991; Rudell 1993) but the literature has not
tackled the multiobjective optimization setting. A VO method can also be
classified as either static or dynamic. Static orderings are specified in
advance of BDD construction whereas dynamic orderings are derived
incrementally as the BDD is constructed layer-by-layer. The latter may or may
not be applicable depending on the BDD construction algorithm; see (Karahalios
and van Hoeve 2022) for a discussion of this distinction for Graph Coloring.
We focus on static orderings and discuss extensions to the dynamic setting in
the Conclusion.

We focus on ML-based heuristics for VO as they can be learned from a relevant
dataset instead of handcrafted heuristics developed through a tedious trial-
and-error process.

### Machine Learning for Variable Ordering

Grumberg, Livne, and Markovitch (2003) proposed a learning-based algorithm to
construct smaller BDDs for model verification. In particular, they create
random orders for a given instance in the training set and tag variable pairs
based on their impact on the resulting BDD size. Using this data, they learn
multiple pair precedence classifiers. For a new instance at test time, they
query each trained pair precedence classifier to construct a precedence table.
These tables are merged into one to derive the ordering. The success of this
method hinges on the selective sampling of informative variable pairs and the
ability to generate orders with sufficient variability in BDD size to increase
the chance of observing high-quality labels. As they leverage problem-specific
heuristics for selective sampling and rely on random sampling for generating
labels, this method is not applicable to our setting. Specifically, we do not
have a notion of how informative a given variable pair is. Additionally,
random orders may not produce high variability in PF enumeration time as
evidenced by our experiments.

Carbin (2006) uses active learning to address the VO problem for BDDs used in
program analysis with the goal of minimizing the run time, similar to that of
LEO. However, certain differences make it less amenable to our setting.
Specifically, the technique to generate the candidate VOs in Carbin (2006) is
grounded in program analysis and cannot be applied to our problem. Instead,
LEO leverages bayesian optimization through SMAC in conjunction with the novel
property-weight search space to generate VO candidates.

Drechsler, Göckel, and Becker (1996) use an evolutionary search approach to
learn a single variable ordering heuristic for a set of instances. The learned
heuristic is a sequence of BDD operations (e.g., swapping two variables in an
ordering) applied to instances from circuit design and verification, where the
BDD represents a boolean function that must be verified efficiently. This
approach can be seen as one of algorithm configuration, which we will compare
to using SMAC and ultimately outperform.

The work of Cappart et al. (2019) is the first learning-based method to
address VO problem for BDDs used in solving discrete optimization problems.
Specifically, they learn a policy to order variables of relaxed/restricted
BDDs to obtain tighter bounds for the underlying optimization problem using
reinforcement learning (RL). A key component of training an RL policy is
devising the reward that the agent receives after taking an action in a given
state and moving to the next state. Note that we care about reducing the time
to enumerate the PF, which is only available at the end of a training episode.
The absence of intermediate rewards in our setting makes RL inapplicable.

Karahalios and van Hoeve (2022) developed an algorithm portfolio approach to
select the best VO strategy from a set of alternatives when constructing
relaxed decision diagrams for the (single-objective) graph coloring problem.
Fundamental to a portfolio-based approach is the existence of a set of
strategies, such that each one of them is better than the others at solving
some problem instances. However, such an algorithm portfolio does not exist in
our case and part of the challenge is to discover good ordering strategies.

There have been some recent contributions (Wu et al. 2022; Liu, Yan, and Jin
2023) relating to solving multiobjective problems using deep learning and
graph neural networks. However, these approaches are not exact and thus beyond
the scope of this paper. We discuss potential extensions of LEO to the inexact
setting in the Conclusion.

Figure 2: Schematic of the proposed method LEO, which comprises of three
phases. In Phase 1, we generate property-weight labels by iteratively
improving them using the interplay between “Black-box optimizer” and “BDD
Manager”. Phase 2 focuses on building datasets of tuples of features extracted
by “Featurizer“ and best property weights obtained from Phase 1 for training
learning-to-rank models. In Phase 3 we do “Model selection” to select the
model with best generalization and use it to predict efficient variable order.
Finally, we use the predicted variable order to construct the BDD and
enumerate the Pareto frontier.

## Methodology

The proposed methodology, as depicted in Figure 2, is divided into three
phases. We apply our technique to the multiobjective knapsack problem (MKP)
which can be described as:

$\max_{x\in\\{0,1\\}^{n}}\left\\{\left\\{\sum_{i\in[n]}a^{p}_{i}x_{i}\right\\}_{p=1}^{P}:\sum_{i\in[n]}w_{i}x_{i}\leq
W\right\\}.$

Here, $n$ is the number of items/variables, $w_{i}$ and $a^{p}_{i}$ are the
weight and profit corresponding to each item $i\in[n]$ and objective
$p\in[P]$. Finally, the capacity of the knapsack is $W\in\mathbb{Z}_{+}$.
Next, we define an efficient variable ordering, a concept which will be useful
in describing LEO.

Let $\mathcal{O}$ be the set of all possible variable orderings and
$\Gamma(o),o\in\mathcal{O}$ denote the PF enumeration time over a BDD
constructed using order $o$. Let
$o^{\star}\equiv\operatorname*{arg\,min}_{o\in\mathcal{O}}\Gamma(o)$ be the
optimal VO. Finding $o^{\star}$ among all $n!$ possible permutations of $n$
variables is intractable, so we will aim for an efficient variable ordering
(EVO) $o^{e}$ that is as close as possible to $o^{\star}$. Note that our
approach is heuristic and does not come with approximation guarantees on the
enumeration time of $o^{e}$ relative to $o^{\star}$.

The objective in the first phase is to find, for each training instance, an
EVO that acts as a label for supervising an ML model. In the second phase,
each training instance is mapped to a set of features and an ML model is
trained. Finally, we perform model selection and use the chosen model to
predict EVOs, referred to as $\hat{o}^{e}$, that are then used to construct a
BDD and compute the PF for any test instance.

Property |  | Definition
---|---|---
weight |  | $w_{i}$
avg-value |  | $\sum_{p=1}^{P}a^{p}_{i}/P$
max-value |  | $\max\\{a^{p}_{i}\\}_{p=1}^{P}$
min-value |  | $\min\\{a^{p}_{i}\\}_{p=1}^{P}$
avg-value-by-weight |  | $\left(\sum_{p=1}^{P}a^{p}_{i}/P\right)/w_{i}$
max-value-by-weight |  | $\max\\{a^{p}_{i}\\}_{p=1}^{P}/w_{i}$
min-value-by-weight |  | $\min\\{a^{p}_{i}\\}_{p=1}^{P}/w_{i}$
Table 1: Properties of a variable $i$ of an MKP.

### Phase 1: Finding an EVO

Since finding an optimal VO that minimizes the PF enumeration time is NP-
complete, we devise a heuristic approach for finding an EVO, $o^{e}$.

To find $o^{e}$ for a given instance $\mathcal{I}$, we use black-box
optimization as the run time $\Gamma^{\mathcal{I}}$ cannot be described
analytically and optimized over. A naive approach to find $o^{e}$ would be to
search directly in the variable ordering space, as suggested in Approach #1 of
Figure 1. While this might work well for tiny problems, it will not scale with
the increase in the problem size as there are $n!$ possible orderings.

To alleviate this issue, we define a surrogate search space that removes the
dependence on the problem size. Specifically, we introduce score-based
variable ordering where we order the variables based on the decreasing order
of their total score. For a given problem class, we define a set of properties
$\mathcal{K}$ for its variables that capture some problem structure. For
example, in an MKP, the weight of an item can act as property. Table 1 lists
all properties of a variable of an MKP. Let $g_{ik}$ be the property score of
variable $i$ for some property $k\in\mathcal{K}$,
$\bar{w}=(w_{1},\cdots,w_{k})$ be the property weights in
$[-1,1]^{|\mathcal{K}|}$. Then, the score of a variable $i$ is defined as

$s_{i}\equiv\sum_{k\in\mathcal{K}}w_{k}\cdot\frac{g_{ik}}{\sum_{i\in[n]}{g_{ik}}}.$
(1)

We recover the variable ordering by sorting them in decreasing order of their
score. Thus, as depicted in Approach #2 in Figure 1, the search for $o^{e}$ is
conducted in the surrogate search space $[-1,1]^{|\mathcal{K}|}$, which only
depends on $|\mathcal{K}|$ and not on the number of variables $n$. Note that
defining the search space in this manner gives an additional layer of dynamism
in the sense that two instances with the same property weights can have
different variable orders.

With a slight abuse of notation, let $\Gamma(\bar{w})$ represent the time
taken to enumerate the PF using a VO obtained by property weights $\bar{w}$.
Given a problem instance, the black-box optimizer controls the property
weights and the BDD manager computes the PF enumeration time based on the VO
derived from these property weights. The black-box optimizer iteratively tries
different property weight configurations, maintaining a list of incumbents, a
process that is depicted in Phase 1 of Figure 2. We propose to use the
variable ordering obtained from the best incumbent property weight as the
label for learning task.

Type | Feature |  Description |  Count
---|---|---|---
Variable | Properties |  The variable properties described in Table 1 |  7
Rank |  The ranks corresponding to the heuristic orderings |  10
|  |  described in Table 3 |
| Value SD |  The standard deviation of the values |  1
Context | # objectives |  The number of objectives in the MKP problem |  1
# items |  The number of items in the MKP problem |  1
Capacity |  The capacity of the MKP |  1
Weight stats. |  The mean, min., max. and std. of the MKP weights |  4
Aggregate value |  The mean, min., max. and std. of the values |  12
stats. |  of each objective |
Total count |  |  |  37
Table 2: Features for an MKP.

### Phase 2: Dataset Generation and Model Training

In this phase, we begin by generating the training dataset. We give special
attention to designing our features such that the resulting models are
permutation-invariant and independent of the size of the problem. Note that
instead of this feature engineering approach, one could use feature learning
through graph neural networks or similar deep learning techniques, see
(Cappart et al. 2023) for a survey. However, given that our case study is on
the knapsack problem, we opt for domain-specific feature engineering that can
directly exploit some problem structure and leads to somewhat interpretable ML
models.

Suppose we are given $J$ problem instances, each having $n_{j}$ variables,
$j\in[J]$. Let $\alpha_{ij}$ denote the features of variable $i\in[n_{j}]$ and
$\beta_{j}$ denote the instance-level context features. Using the features and
the EVO computed in Phase 1, we construct a dataset
$\mathcal{D}=\left\\{(\alpha_{ij},\beta_{j},r_{i}({o^{e}}_{j}):i\in[n_{j}]:j\in[J]\right\\}$.
Here, $o^{e}_{j}$ is the EVO of an instance $j$ and
$r_{i}:\mathbb{Z}^{n_{j}}_{+}\rightarrow\mathbb{Z}_{+}$ a mapping from the EVO
to the rank of a variable $i$. For example, if $n_{j}=4$ for some instance $j$
and $o^{e}_{j}=(2,1,4,3)$, then $r_{1}(o^{e}_{j})=3$, $r_{2}(o^{e}_{j})=4$,
$r_{3}(o^{e}_{j})=1$, and $r_{4}(o^{e}_{j})=2$. For a complete list of
variable and context features, refer to Table 2.

Learning-to-rank (LTR) is an ML approach specifically designed for ranking
tasks, where the goal is to sort a set of items based on their relevance given
a query or context. It is commonly used in applications such as information
retrieval. We formulate the task of predicting the EVO as an LTR task and use
the pointwise and pairwise ranking approaches to solve the problem.

In the pointwise approach, each item (variable) in the training data is
treated independently, and the goal is to learn a model that directly predicts
the relevance score or label for each variable. This is similar to solving a
regression problem with a mean-squared error loss. Specifically, we train a
model $f_{\theta}(\alpha_{ij},\beta_{j})$ to predict $r_{i}(o^{e}_{j})$. Once
the model is trained, it can be used to rank items based on their predicted
scores $\hat{o}^{e}$. However, this approach does not explicitly consider the
relationships between variables in its loss function.

The pairwise approach aims to resolve this issue (Joachims 2002). Let

$\mathcal{T}_{j}=\left\\{(i_{1},i_{2}):r_{i_{1}}(o^{e}_{j})>r_{i_{1}}(o^{e}_{j}))\right\\}$

be the set of all variable tuples $(i_{1},i_{2})$ such that $i_{1}$ is ranked
higher than $i_{2}$ for instance $j$. Then, the goal is to learn a model that
maximizes the number of respected pairwise-ordering constraints. Specifically,
we train a model $g_{\phi}(\cdot)$ such that the number of pairs
$(i_{1},i_{2})\in\mathcal{T}_{j}$ for which
$g_{\phi}(\alpha_{i_{1}j},\beta_{j})>g_{\phi}(\alpha_{i_{2}j},\beta_{j})$ is
maximized. This approach is better equipped to solve the ranking problem as
the structured loss takes into account the pairwise relationships.

### Phase 3: Model Selection and Testing

We follow a two-step approach to perform model selection as the ranking task
is a proxy to the downstream task of efficiently solving a multiobjective
problem. Firstly, for each model class (e.g., decision trees, linear models,
etc.), we select the best model based on Kendall’s Tau (Kendall 1938), a
ranking performance metric that measures the fraction of violated pairwise-
ordering constraints, on instances from a validation set different from the
training set. Subsequently, we pit the best models from each type against one
another and select the winner based on the minimum average PF enumeration time
on the validation set. Henceforth, for previously unseen instances from the
test set, we will use the model selected in Phase 3 to predict the EVO and
then compute the PF.

## Computational Setup

Our code and instances are available at https://github.com/khalil-
research/leo. All the experiments reported in this manuscript are conducted on
a computing cluster with an Intel Xeon CPU E5-2683 CPUs. We use SMAC (Lindauer
et al. 2021) – a black-box optimization and algorithm configuration library –
for generating the labels. The ML models are built using Python 3.8, Scikit-
learn (Pedregosa et al. 2011), and XGBoost (Chen and Guestrin 2016), and
SVMRank (Joachims 2006). The “BDD Manager” is based on the implementation of
Bergman and Cire (2016), which is available at
https://www.andrew.cmu.edu/user/vanhoeve/mdd/.

### Instance Generation

We use a dataset of randomly generated MKP instances as described in (Bergman
and Cire 2016). The values $w_{i}$ and $a^{p}_{i}$ are sampled randomly from a
discrete uniform distribution ranging from 1 to 100. The capacity $W$ is set
to $\left\lceil 0.5\sum_{i\in I}w_{i}\right\rceil$. We generate instances with
sizes

$\mathcal{S}=\\{(3,60),(3,70),(3,80),(4,50),(5,40),(6,40),(7,40)\\},$

where the first and second component of the tuple specify the number of
objectives and variables, respectively. For each size, we generate 1000
training instances, 100 validation instances, and 100 test instances.

### Instrumenting SMAC

As a labeling tool: To generate EVOs for the learning-based models, we use
SMAC (Lindauer et al. 2021). In what follows, SmacI refers to the use of SMAC
as a black-box optimizer that finds an EVO for a given training instance.
Specifically, SmacI solves $\bar{w}^{e}_{j}=\min_{\bar{w}}\Gamma_{j}(\bar{w})$
for each instance $j$ in the training set; we obtain $o^{e}_{j}$ by
calculating variable scores – a dot product between $\bar{w}^{e}_{j}$ and the
corresponding property value – and sorting variables in the decreasing order
of their score.

As a baseline: The other more standard use of SMAC, which we refer to as
SmacD, is as an algorithm configuration tool. To obtain a SmacD ordering, we
use SMAC to solve
$\bar{w}^{e}_{D}=\min_{\bar{w}}\mathbb{E}_{\bm{j\sim[J]}}\left[\Gamma_{j}(\bar{w})\right]$
and obtain an order $o^{e}_{D_{j}}$ for instance $j$ using single property
weight vector $\bar{w}^{e}_{D}$. The expectation of the run time here simply
represents its average over the $|J|$ training instances. Note that we get
only one property weight vector for the entire dataset in the SmacD case
rather than one per instance as in the SmacI case. However, we obtain an
instance-specific VO when using the SmacD configuration as the underlying
property values change across instances.

Initialization: For both uses of SMAC, we use the MinWt ordering as a warm-
start by assigning all the property weights to zero except the weight
property, which is set to -1. This reduces the need for extensive random
exploration of the configuration space by providing a reasonably good ordering
heuristic.

Random seeds: As SMAC is a randomized optimization algorithm, running it with
multiple random seeds increases the odds of finding a good solution. We
leverage this idea for SmacI as its outputs will be used as labels for
supervised learning and are thus expected to be of very high quality, i.e., we
seek instance-specific parameter configurations that yield variable orderings
with minimal solution times. We run SmacD with a single seed and use its
average run time on the training set as a target to beat for SmacI. Since the
latter optimizes at the instance level, one would hope it can do better than
the distribution-level configuration of SmacD.

As such, for SmacI, we run SMAC on each instance with 5 seeds for all sizes
except (5, 40) and (6, 40), for which we used one seed. We start with one
seed, then average the run time of the best-performing seed per instance to
that of the average enumeration time of SmacD on the training set, and
relaunch SMAC with a new seed until a desired performance gap between SmacI
and SmacD is achieved.

Computational budget: In the SmacD setting, we run SMAC with a 12-hour time
limit, whereas in the SmacI case the time limit is set to 20 minutes per
instance except for sizes (3, 60), (4, 50), (5, 40), for which it is set to 5
minutes per instance. In both settings, SMAC runs on a 4-core machine with a
20GB memory limit. It can be observed that generating the labels can be
computationally expensive. This dictates the choice of sizes in the set of
instance sizes, $\mathcal{S}$. Specifically, we select instance sets with an
average running time of at most 100 seconds (not too hard) and at least 5
seconds (nor too easy) using the top-down compilation method described in
Bergman et al. (2021).

### Learning Models

We use linear regression, ridge regression, lasso regression, decision trees,
and gradient-boosted trees (GBT) with mean-squared error loss to build size-
specific pointwise ranking models. Similarly, we train support vector machines
and GBT with pairwise-ranking loss to obtain pairwise ranking models for each
size. In the experimental results that follow, the best learning-based method
will be referred to as ML. This turns out to be GBT trained with pairwise-
ranking loss. The GBT models that were selected achieved a Kendall’s Tau
ranging between 0.67 to 0.81 across all problem sizes on the validation set.
The model selection follows the procedure mentioned in Phase 3 of the
Methodology section.

In terms of features, we omit the context features in Table 2 when training
linear size-specific models, as these features take on the same value for all
variables of the same instance and thus do not contribute to the prediction of
the rank of a variable. The context features are used with non-linear models
such as decision trees and GBT.

We also train two additional GBT models – ML+A and ML+AC – with pairwise-
ranking loss, on the union of the datasets of all sizes. In particular, ML+A
is trained with only variable features, similar to ML, whereas ML+AC adds the
instance context features to the variable features.

### Baselines

To evaluate the performance of learning-based orderings, we compare to four
baselines:

  * –

Lex uses the (arbitrary) default variable ordering in which the instance was
generated.

  * –

MinWt orders the variables in increasing order of their weight values,
$w_{i}$. This is a commonly used heuristic for solving the single-objective
knapsack problem.

  * –

MaxRatio orders the variables in decreasing order of the property min-value-
by-weight detailed in Table 1, which is defined as
$\min\\{a^{p}_{i}\\}_{p=1}^{P}/w_{i}$. This rule has an intuitive
interpretation: it prefers variables with larger worst-case (the minimum in
the numerator) value-to-weight ratio. It is not surprising that this heuristic
might perform well given that it is part of a $\frac{1}{2}$-approximation
algorithm for the single-objective knapsack problem (Chekuri and Fox 2009).

  * –

SmacD, as described in an earlier paragraph. It produces one weight setting
for the property scores per instance set. This baseline can be seen as a
representative for the algorithm configuration paradigm recently surveyed by
Schede et al. (2022).

Heuristic ordering | Size
---|---
⟨sort⟩_⟨property-name⟩ | (3, 20) | (3, 40) | (3, 60) | (3, 80) | (5, 20) | (5, 30) | (5, 40) | (7, 20) | (7, 30) | (7, 40)
max_weight | 0.778 | 1.147 | 1.378 | 1.317 | 0.815 | 1.865 | 2.760 | 1.165 | 2.769 | 1.830
min_weight | 0.812 | 0.654 | 0.648 | 0.685 | 0.769 | 0.700 | 0.533 | 0.694 | 0.588 | 0.581
max_avg-value | 0.861 | 0.891 | 0.843 | 0.972 | 0.839 | 0.992 | 0.840 | 0.845 | 1.125 | 0.823
min_avg-value | 0.893 | 1.619 | 2.144 | 1.692 | 0.936 | 1.721 | 2.132 | 1.103 | 1.918 | 1.524
max_max-value | 0.858 | 1.065 | 1.126 | 1.203 | 0.857 | 1.137 | 1.060 | 0.899 | 1.408 | 0.992
min_max-value | 0.855 | 1.018 | 1.010 | 0.995 | 0.854 | 1.077 | 1.009 | 0.909 | 0.981 | 0.975
max_min-value | 0.851 | 0.705 | 0.537 | 0.590 | 0.813 | 0.816 | 0.617 | 0.819 | 0.862 | 0.767
min_min-value | 0.883 | 1.662 | 2.266 | 1.729 | 0.930 | 1.700 | 2.020 | 1.094 | 1.980 | 1.430
max_avg-value-by-weight | 0.838 | 0.832 | 0.988 | 1.156 | 0.802 | 0.896 | 0.887 | 0.751 | 0.917 | 0.788
max_max-value-by-weight | 0.837 | 0.837 | 1.037 | 1.192 | 0.804 | 0.871 | 0.824 | 0.732 | 0.790 | 0.735

Table 3: Ratio of the average time taken by a heuristic ordering to that of 5
random orderings across 250 MKP instances across sizes specified. A ratio of
less (greater) than 1 indicates that heuristic ordering performs better
(worse) than random ordering, on average. The name of a heuristic ordering
follows the structure ⟨sort⟩_⟨property-name⟩, where ⟨sort⟩ can either be min
or max and ⟨property-name⟩ can be equal to one of the properties defined in
Table 1. A value of max (min) for ⟨sort⟩ would sort the variables in the
descending (ascending) order of the property ⟨property-name⟩.

## Experimental Results

We examine our experimental findings through a series of questions that span
the impact of VO on PF enumeration time (Q1), the performance of SmacI as a
black-box optimizer that provides labels for ML (Q2), the performance of LEO
on unseen test instances and comparison to the baselines (Q3, Q4), a feature
importance analysis of the best ML models used by LEO (Q5), and an exploration
of size-independent, unified ML models (Q6).

#### Q1. Does variable ordering impact the Pareto frontier enumeration time?

To test the hypothesis that VO has an impact on the PF enumeration time, we
compare the run time of ten heuristic orderings against the expected run time
of random orderings. By the run time of an ordering, we mean the time taken to
compute the PF on a BDD constructed using that particular ordering.

The heuristic orderings are constructed using the variable properties. For
example, the heuristic ordering called max_avg-value-by-weight sorts the
variables in descending order of the ratio of a variable’s average value by
its weight. We estimate the expected run time of a random ordering by sampling
5 heuristic variable orderings uniformly at random from all possible $n!$
orderings and averaging their run time. Table 3 summarizes the results of this
experiment. We work with 250 MKP instances of the problem sizes (3, 20), (3,
40), (3, 60), (3, 80), (5, 20), (5, 30), (5, 40), (7, 20), (7, 30), and (7,
40). The values in the table are the ratios of the average run time of a
heuristic ordering to that of the 5 random orderings; values smaller than one
indicate that a heuristic ordering is faster than a random one, on average.
The best heuristic ordering for each size is highlighted.

First, it is clear that some heuristic orderings consistently outperform the
random ones across all problem sizes (min_weight, max_min-value), and by a
significant margin. In contrast, some heuristic orderings are consistently
worse than random. For example, the heuristic ordering max_weight, min_avg-
value, and min_min-value consistently underperform when the number of
variables is more than 20.

Second, the choice of MinWt as a baseline was motivated by the results of this
experiment as min_weight wins on most sizes as highlighted in Table 3.

Altogether, this experiment validates the hypothesis that VO has an impact on
the PF enumeration time and that there exists EVO that can significantly
reduce this time.

Size | Method | Count | GMean | Min | Max
---|---|---|---|---|---
(5, 40) | ML | 100 | 1.43 | 1.26 | 40.91
ML+AC | 100 | 1.54 | 1.06 | 41.29
ML+A | 100 | 1.38 | 1.28 | 41.24
SmacD | 100 | 2.49 | 1.38 | 51.22
MaxRatio | 100 | 2.48 | 1.44 | 58.87
| MinWt | 100 | 2.42 | 1.44 | 52.61
| Lex | 100 | 6.23 | 1.86 | 109.18
(6, 40) | ML | 100 | 14.45 | 2.57 | 377.89
ML+AC | 100 | 13.62 | 2.01 | 427.80
ML+A | 100 | 13.67 | 2.29 | 355.64
SmacD | 100 | 18.62 | 2.95 | 434.34
MaxRatio | 100 | 19.84 | 3.05 | 664.88
| MinWt | 100 | 18.82 | 2.91 | 427.06
| Lex | 100 | 42.23 | 6.01 | 713.58
(7, 40) | ML | 100 | 41.97 | 2.26 | 1490.58
ML+AC | 99 | 45.39 | 2.27 | 1153.07
ML+A | 99 | 44.70 | 2.81 | 1130.82
SmacD | 99 | 64.84 | 2.95 | 1521.43
MaxRatio | 98 | 86.49 | 4.14 | 1700.44
| MinWt | 99 | 60.08 | 3.16 | 1561.36
| Lex | 95 | 102.23 | 7.66 | 1604.17
(4, 50) | ML | 100 | 5.75 | 2.54 | 79.42
ML+AC | 100 | 5.97 | 2.22 | 77.75
ML+A | 100 | 6.04 | 2.47 | 73.17
SmacD | 100 | 12.17 | 4.04 | 218.24
MaxRatio | 100 | 14.23 | 3.24 | 167.95
| MinWt | 100 | 18.37 | 4.67 | 553.15
| Lex | 100 | 32.25 | 6.67 | 395.67
(3, 60) | ML | 100 | 3.15 | 3.88 | 34.13
ML+AC | 100 | 3.07 | 3.44 | 31.99
ML+A | 100 | 3.23 | 3.56 | 31.14
SmacD | 100 | 4.07 | 3.53 | 43.58
MaxRatio | 100 | 9.24 | 5.76 | 64.21
| MinWt | 100 | 11.61 | 5.36 | 75.98
| Lex | 100 | 20.19 | 8.12 | 195.25
(3, 70) | ML | 100 | 15.05 | 7.92 | 64.55
ML+AC | 100 | 14.94 | 7.41 | 65.67
ML+A | 100 | 15.35 | 7.68 | 80.82
SmacD | 100 | 16.64 | 7.52 | 96.94
MaxRatio | 100 | 33.54 | 12.01 | 215.63
| MinWt | 100 | 44.43 | 12.17 | 194.07
| Lex | 100 | 71.55 | 19.87 | 308.21
(3, 80) | ML | 100 | 43.85 | 15.29 | 286.11
ML+AC | 100 | 42.89 | 15.43 | 262.57
ML+A | 100 | 44.64 | 15.35 | 262.77
SmacD | 100 | 53.33 | 15.66 | 632.92
MaxRatio | 100 | 100.31 | 29.64 | 741.82
| MinWt | 100 | 129.40 | 30.35 | 794.10
| Lex | 99 | 196.31 | 56.52 | 1332.19

Table 4: PF enumeration time evaluation results on the test set with a time
limit of 1800 seconds. Each row represents aggregate statistics for a given
instance “Size” and “Method”. “Count” stands for the number of instances on
which the algorithm ran successfully without hitting the time or memory
limits. “GMean”, “Min” and “Max” denotes the geometric mean, minimum and
maximum of the enumeration time computed across “count” many instances. We use
a shift of 5 to compute the “GMean”.

#### Q2. Can SmacI find good variable orderings for training instances?

Figure 3: SmacI performance w.r.t. wallclock time across different sizes. For
a given size, we first of all find the best seed, i.e., the seed which finds
an incumbent property-weight having minimum PF enumeration time, for each
instance. Then for a given wallclock time, we calculate the mean and standard
error of the PF enumeration time across instances for the incumbent property-
weight up to that wallclock time. The mean and standard error are plotted as a
blue line and blue error band around it, respectively. The horizontal black
line represents the average PF enumeration time of SmacD configuration on the
validation set.

Having established in Q1 that the search for high-quality variable orderings
is justified, we now turn to a key component of Phase 1: the use of SMAC as a
black-box optimizer that produces “label” variable orderings for subsequent ML
in Phase 2. Figure 3 shows the average run time, on the training instances, of
the incumbent property-weight configurations found by SMAC. These should be
interpreted as standard optimization convergence curves where we seek small
values (vertical axis) as quickly as possible (horizontal axis). Indeed, SMAC
performs well, substantially improving over its initial warm-start solution
(the MinWt ordering). The flat horizontal lines in the figure show the average
run time of the single configuration found by SmacD. The latter is ultimately
outperformed by the SmacI configurations on average, as desired.

#### Q3. How does LEO perform in comparison to baseline methods?

Table 4 presents the PF enumeration across different problem sizes and methods
on the test set. In this question, we will examine only the first of three ML
models, namely the first one referred to as ML in the table; the two other
models will be explored in Q7. We can observe that ML consistently outperforms
all baselines in terms of the geometric mean (GMean) of the PF enumeration
time, across all problem sizes.

SmacD acts as a strong baseline for the ML method, consistently being the
second-best in terms of average GMean, except for size (7, 40). The methods
MinWt and MaxRatio have a nice connection to the single-objective knapsack
problem (Chekuri and Fox 2009), as discussed earlier; this might explain why
these heuristics reduce the number of intermediate solutions being generated
in the BDDs during PF enumeration, helping the algorithm terminate quickly.
They closely follow SmacD in terms of GMean metric with sizes having more than
3 objectives; however, they are almost twice as worst as SmacD on instances
with 3 objectives.

Interestingly, when the number of objectives is larger than 3, MinWt
outperforms MaxRatio (and vice versa). This underscores the relationship
between the number of objectives and heuristics, i.e., one heuristic might be
preferred over another depending on the structure of the problem. This also
highlights why the ML method performs better than the heuristics as the
feature engineering helps create an ensemble of them rather than using only
one of them. Lastly, we observe that method Lex has the worst performance
across all sizes.

Figure 4: Box plots of time to enumerate the PF for different sizes. Figure 5:
Performance profile of different methods in terms of the fraction of instances
solved w.r.t. time for various sizes.

To complement Table 4, we present box plots and performance profiles in Figure
4 and Figure 5, respectively. Note how the PF enumeration time distribution
for the ML method is much more concentrated, with a smaller median and only a
few outliers compared to other methods. Note that the instances that had a
timeout are omitted from this analysis. This performance improvement leads to
a larger fraction of instances being solved in a smaller amount of time as
highlighted in Figure 5.

#### Q4. What explains LEO’s performance?

Size | Method | Nodes(%) | Width(%) | Checks(%) | GMean(%)
---|---|---|---|---|---
(5, 40) | ML | 81.46 | 88.13 | 36.28 | 22.91
SmacD | 69.49 | 79.11 | 56.15 | 39.95
MaxRatio | 88.80 | 94.10 | 52.49 | 39.75
MinWt | 69.49 | 79.11 | 56.15 | 38.76
(6, 40) | ML | 78.35 | 85.82 | 39.36 | 34.22
SmacD | 70.12 | 79.53 | 52.02 | 44.09
MaxRatio | 90.63 | 95.01 | 53.39 | 46.97
MinWt | 70.12 | 79.53 | 52.02 | 44.56
(7, 40) | ML | 81.21 | 88.70 | 56.99 | 41.06
SmacD | 71.18 | 80.18 | 79.12 | 63.43
MaxRatio | 72.11 | 81.37 | 100.93 | 84.61
MinWt | 70.11 | 79.52 | 69.13 | 58.77
(4, 50) | ML | 90.60 | 96.68 | 21.63 | 17.84
SmacD | 76.03 | 87.94 | 42.93 | 37.73
MaxRatio | 88.33 | 96.27 | 48.97 | 44.13
MinWt | 70.45 | 84.91 | 68.81 | 56.97
(3, 60) | ML | 93.11 | 98.36 | 20.18 | 15.58
SmacD | 91.70 | 98.81 | 26.14 | 20.14
MaxRatio | 86.32 | 96.42 | 48.45 | 45.76
MinWt | 68.68 | 86.68 | 62.06 | 57.51
(3, 70) | ML | 92.30 | 98.57 | 18.21 | 21.03
SmacD | 90.81 | 98.71 | 21.52 | 23.26
MaxRatio | 86.20 | 97.16 | 49.04 | 46.88
MinWt | 70.33 | 89.75 | 68.03 | 62.10
(3, 80) | ML | 91.09 | 98.69 | 19.94 | 22.34
SmacD | 92.34 | 99.27 | 27.77 | 27.17
MaxRatio | 84.45 | 97.34 | 57.17 | 51.10
MinWt | 68.79 | 90.58 | 74.72 | 65.91

Table 5: Analyzing the impact of BDD topology on the PF enumeration time on
the test set. For a given “Size” and “Method”, the metric value indicates the
percentage of mean performance by that “Method” compared to the Lex method.
The “Nodes” (“Width”) represents the percentage of the mean of the number of
nodes (mean width) of the BDDs generated by “Method” to that of Lex.
Similarly, “Checks” and “Time” indicates the percentage of the average Pareto-
dominance checks and the average geometric mean of the PF enumeration time,
respectively. For all metrics, values less than 100 are preferred. For
example, a metric with a value of $\delta$% indicates the “Method” is
$(100-\delta)$% better than Lex. Figure 6: Cumulative number of intermediate
solutions generated at a particular layer across all sizes.

Traditionally, smaller-sized exact BDDs are sought after for efficiently
performing tasks such as model checking or computing objective function
bounds. Hoping to find a similar connection in the multiobjective setting, we
analyze the relationship between the topology of BDDs generated by different
methods and the time to compute the PF in Table 5. Considering the fact that
Lex performs the worst among all methods and the impact of BDD size on the
downstream task, we expect that the size of the BDDs generated by Lex to be
bigger. This holds true across different sizes and methods as it can be
observed in Table 5 that Nodes and Width values are less than 100.

Extrapolating, one would expect BDDs generated using ML orderings to be the
smallest, as it achieves the best performance in terms of enumeration time.
Counter-intuitively, that is not the case: MinWt generates the smallest-sized
BDDs on average. For instance, the value of Nodes for size (3, 80) is 68.79%
for MinWt, which is lower than 91.09% for ML. However, the reduction in size
does not translate to improvements in the running time, as we already know
that ML performs best in terms of time.

We can decipher the performance gains of ML by studying the “Checks” metric.
This metric can be thought of as a proxy of the work done by the PF
enumeration algorithm; indeed, this metric is positively correlated with Time.
This phenomenon can be further studied in Figure 6, which shows the mean
cumulative solutions generated up to a given layer in the BDD. Clearly, ML has
the least number of intermediate solutions generated, which also translates to
fewer Pareto-dominance checks and smaller Checks.

To summarize, smaller-sized BDDs can be efficient in enumerating the PF.
However, reducing the BDD size beyond a certain threshold would inversely
affect its performance as it leads to more intermediate solutions being
generated, increasing the number of Pareto-dominance checks. This also
validates the need for task-specific methods like LEO that are specifically
geared towards reducing the run time rather than optimizing a topological
metric of the BDD.

#### Q5. How interpretable are the decisions made by LEO?

LEO uses GBT with the pairwise ranking loss for learning to order variables.
To obtain feature importance scores, we count the number of times a particular
feature is used to split a node in the decision trees. We then normalize these
scores by dividing them by the maximum score, resulting in values in $[0,1]$
for each size.

Figure 7: Heatmap of feature importance scores normalized across different
sizes. The features are described in Table 2. The property-based features
assume the same names as the properties detailed in Table 1. The rank-based
features of a variable have a prefix “rk_” before the heuristic ordering;
refer to Table 3 for more details. Finally, the “Value SD” feature from Table
2 is named std-value in this plot.

Figure 7 is a heatmap of feature importance scores for different sizes. We can
note that the min-value-by-weight feature is important across all sizes,
especially for cases with more than 3 objectives. In fact, the choice of VO
heuristic MaxRatio was driven by feature importance scores and is a case in
point of how learning-based methods can assist in designing good heuristics
for optimization. Furthermore, the real-valued features are more important
than the categorical rank features for problems with more than 3 objectives.

For problems with 3 objectives, the std-value feature deviation is extremely
crucial. Also, rank feature rk_max_avg-value-by-weight receives higher
importance than some of the real-valued features. It is also interesting to
observe that certain rank-based features are consistently ignored across all
sizes.

In a nutshell, the heatmap of feature importance scores helps in interpreting
the features that govern the decisions made by LEO, which also happens to be
in alignment with a widely used heuristic to solve the single-objective
knapsack problem.

#### Q6. Can we train a single size-independent ML model?

We answer in the affirmative: in fact, the methods ML+A and ML+AC in Table 4
refer to two models trained on the union of all training datasets of different
sizes. As described in the section “Learning Models”, ML+A uses the same
variable features as ML, whereas ML+AC adds the instance context features
described in Table 2. These two models perform remarkably well, sometimes even
outperforming the size-specific ML models that are being tested in-
distribution.

The main takeaway here is that our ML model architecture, which is independent
of the number of variables in an instance, enables the training of unified
models that are size-independent. In turn, these models perform very well.
Combined with the simple features and types of ML models we have used, this
finding points to a great potential for the use of ML in multiobjective BDD
problems. The appendix includes feature importance plots for both of these
models. ML+AC seems to make great use of the context features in its
predictions, as expected.

## Conclusion

LEO is the first machine learning framework for accelerating the BDD approach
for multiobjective integer linear programming. We contribute a number of
techniques and findings that may be of independent interest. For one, we have
shown that variable ordering does impact the solution time of this approach.
Our labeling method of using a black-box optimizer over a fixed-size parameter
configuration space to discover high-quality variable orderings successfully
bypasses the curse of dimensionality; the application of this approach to
other similar labeling tasks may be possible, e.g., for backdoor set discovery
in MIP (Khalil, Vaezipoor, and Dilkina 2022). An additional innovation is
using size-independent ML models, i.e., models that do not depend on a fixed
number of decision variables. This modeling choice enables the training of
unified ML models, which our experiments reveal to perform very well. Through
a comprehensive case study of the knapsack problem, we show that LEO can
produce variable orderings that significantly reduce PF enumeration time.

There are several exciting directions for future work that can build on LEO:

  1. –

The BDD approach to multiobjective integer programming has been applied to a
few problems other than knapsack (Bergman et al. 2021). Much of LEO can be
directly extended to such other problems, assuming that their BDD construction
is significantly influenced by the variable ordering. One barrier to such an
extension is the availability of open-source BDD construction code. As the use
of BDDs in optimization is a rather nascent area of research, it is not
uncommon to consider a case study of a single combinatorial problem, as was
done for example for Graph Coloring in (Karahalios and van Hoeve 2022) and
Maximum Independent Set in (Bergman et al. 2012).

  2. –

Our method produces a static variable ordering upfront of BDD construction.
While this was sufficient to improve on non-ML orderings for the knapsack
problem, it may be interesting to consider dynamic variable orderings that
observe the BDD construction process layer by layer and choose the next
variable accordingly, as was done in (Cappart et al. 2019).

  3. –

We have opted for rather interpretable ML model classes but the exploration of
more sophisticated deep learning approaches may enable closing some of the
remaining gap in training/validation loss, which may improve downstream
solving performance.

  4. –

Beyond the exact multiobjective setting, extending LEO to a heuristic that
operates on a restricted BDD may provide an approximate PF much faster than
full enumeration. We believe this to be an easy extension of our work.

## References

  * Adelgren and Gupte (2022) Adelgren, N.; and Gupte, A. 2022.  Branch-and-bound for biobjective mixed-integer linear programming.  _INFORMS Journal on Computing_ , 34(2): 909–933.
  * Aloul, Markov, and Sakallah (2003) Aloul, F. A.; Markov, I. L.; and Sakallah, K. A. 2003.  FORCE: a fast and easy-to-implement variable-ordering heuristic.  In _Proceedings of the 13th ACM Great Lakes symposium on VLSI_ , 116–119.
  * Altiparmak et al. (2006) Altiparmak, F.; Gen, M.; Lin, L.; and Paksoy, T. 2006.  A genetic algorithm approach for multi-objective optimization of supply chain networks.  _Computers & industrial engineering_, 51(1): 196–215.
  * Belotti, Soylu, and Wiecek (2013) Belotti, P.; Soylu, B.; and Wiecek, M. M. 2013.  A branch-and-bound algorithm for biobjective mixed-integer programs.  _Optimization Online_ , 1–29.
  * Bergman et al. (2021) Bergman, D.; Bodur, M.; Cardonha, C.; and Cire, A. A. 2021.  Network models for multiobjective discrete optimization.  _INFORMS Journal on Computing_.
  * Bergman and Cire (2016) Bergman, D.; and Cire, A. A. 2016.  Multiobjective optimization by decision diagrams.  In _International Conference on Principles and Practice of Constraint Programming_ , 86–95. Springer.
  * Bergman et al. (2016a) Bergman, D.; Ciré, A. A.; van Hoeve, W.; and Hooker, J. N. 2016a.  Discrete Optimization with Decision Diagrams.  _INFORMS J. Comput._ , 28(1): 47–66.
  * Bergman et al. (2016b) Bergman, D.; Cire, A. A.; Van Hoeve, W.-J.; and Hooker, J. 2016b.  _Decision diagrams for optimization_ , volume 1.  Springer.
  * Bergman et al. (2012) Bergman, D.; Cire, A. A.; Van Hoeve, W.-J.; and Hooker, J. N. 2012.  Variable ordering for the application of BDDs to the maximum independent set problem.  In _International Conference on Integration of Artificial Intelligence (AI) and Operations Research (OR) Techniques in Constraint Programming_ , 34–49. Springer.
  * Boland, Charkhgard, and Savelsbergh (2014) Boland, N.; Charkhgard, H.; and Savelsbergh, M. 2014.  The triangle splitting method for biobjective mixed integer programming.  In _Integer Programming and Combinatorial Optimization: 17th International Conference, IPCO 2014, Bonn, Germany, June 23-25, 2014. Proceedings 17_ , 162–173. Springer.
  * Boland, Charkhgard, and Savelsbergh (2015) Boland, N.; Charkhgard, H.; and Savelsbergh, M. 2015.  A criterion space search algorithm for biobjective integer programming: The balanced box method.  _INFORMS Journal on Computing_ , 27(4): 735–754.
  * Bollig, Löbbing, and Wegener (1996) Bollig, B.; Löbbing, M.; and Wegener, I. 1996.  On the effect of local changes in the variable ordering of ordered decision diagrams.  _Information Processing Letters_ , 59(5): 233–239.
  * Bollig and Wegener (1996) Bollig, B.; and Wegener, I. 1996.  Improving the variable ordering of OBDDs is NP-complete.  _IEEE Transactions on computers_ , 45(9): 993–1002.
  * Bryant (1992) Bryant, R. E. 1992.  Symbolic boolean manipulation with ordered binary-decision diagrams.  _ACM Computing Surveys (CSUR)_ , 24(3): 293–318.
  * Butler et al. (1991) Butler, K. M.; Ross, D. E.; Kapur, R.; and Mercer, M. R. 1991.  Heuristics to compute variable orderings for efficient manipulation of ordered binary decision diagrams.  In _Proceedings of the 28th ACM/IEEE Design Automation Conference_ , 417–420.
  * Cappart et al. (2023) Cappart, Q.; Chételat, D.; Khalil, E. B.; Lodi, A.; Morris, C.; and Velickovic, P. 2023.  Combinatorial Optimization and Reasoning with Graph Neural Networks.  _Journal of Machine Learning Research_ , 24(130): 1–61.
  * Cappart et al. (2019) Cappart, Q.; Goutierre, E.; Bergman, D.; and Rousseau, L. 2019.  Improving Optimization Bounds Using Machine Learning: Decision Diagrams Meet Deep Reinforcement Learning.  In _The Thirty-Third AAAI Conference on Artificial Intelligence, AAAI 2019, The Thirty-First Innovative Applications of Artificial Intelligence Conference, IAAI 2019, The Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019, Honolulu, Hawaii, USA, January 27 - February 1, 2019_ , 1443–1451. AAAI Press.
  * Carbin (2006) Carbin, M. 2006.  Learning effective BDD variable orders for BDD-based program analysis.  Technical report, Citeseer.
  * Chekuri and Fox (2009) Chekuri, C.; and Fox, K. 2009.  UIUC CS 598CSC: Approximation Algorithms, Lecture Notes on Knapsack.  URL: https://courses.engr.illinois.edu/cs598csc/sp2009/lectures/lecture˙4.pdf. Last visited on 2023/07/04.
  * Chen and Guestrin (2016) Chen, T.; and Guestrin, C. 2016.  XGBoost: A Scalable Tree Boosting System.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, 785–794. New York, NY, USA: ACM.  ISBN 978-1-4503-4232-2.
  * Chung, Hajj, and Patel (1993) Chung, P.-Y.; Hajj, I.; and Patel, J. H. 1993.  Efficient variable ordering heuristics for shared ROBDD.  In _1993 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 1690–1693. IEEE.
  * Drechsler, Göckel, and Becker (1996) Drechsler, R.; Göckel, N.; and Becker, B. 1996.  Learning heuristics for OBDD minimization by evolutionary algorithms.  In _Parallel Problem Solving from Nature—PPSN IV: International Conference on Evolutionary Computation—The 4th International Conference on Parallel Problem Solving from Nature Berlin, Germany, September 22–26, 1996 Proceedings 4_ , 730–739. Springer.
  * Ebendt and Drechsler (2007) Ebendt, R.; and Drechsler, R. 2007.  Exact BDD minimization for path-related objective functions.  In _Vlsi-Soc: From Systems To Silicon_ , 299–315. Springer.
  * Ebendt, Gunther, and Drechsler (2005) Ebendt, R.; Gunther, W.; and Drechsler, R. 2005.  Combining ordered best-first search with branch and bound for exact BDD minimization.  _IEEE Transactions on Computer-aided Design of Integrated Circuits and Systems_ , 24(10): 1515–1529.
  * Ehrgott (2006) Ehrgott, M. 2006.  A discussion of scalarization techniques for multiple objective integer programming.  _Annals of Operations Research_ , 147(1): 343–360.
  * Ehrgott, Gandibleux, and Przybylski (2016) Ehrgott, M.; Gandibleux, X.; and Przybylski, A. 2016.  Exact methods for multi-objective combinatorial optimisation.  _Multiple criteria decision analysis: State of the art surveys_ , 817–850.
  * Friedman and Supowit (1987) Friedman, S. J.; and Supowit, K. J. 1987.  Finding the optimal variable ordering for binary decision diagrams.  In _Proceedings of the 24th ACM/IEEE Design Automation Conference_ , 348–356.
  * Grumberg, Livne, and Markovitch (2003) Grumberg, O.; Livne, S.; and Markovitch, S. 2003.  Learning to order BDD variables in verification.  _Journal of Artificial Intelligence Research_ , 18: 83–116.
  * Joachims (2002) Joachims, T. 2002.  Optimizing search engines using clickthrough data.  In _Proceedings of the eighth ACM SIGKDD international conference on Knowledge discovery and data mining_ , 133–142.
  * Joachims (2006) Joachims, T. 2006.  Training linear SVMs in linear time.  In _Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining_ , 217–226.
  * Karahalios and van Hoeve (2022) Karahalios, A.; and van Hoeve, W.-J. 2022.  Variable ordering for decision diagrams: A portfolio approach.  _Constraints_ , 27(1-2): 116–133.
  * Kendall (1938) Kendall, M. G. 1938.  A new measure of rank correlation.  _Biometrika_ , 30(1/2): 81–93.
  * Khalil, Vaezipoor, and Dilkina (2022) Khalil, E. B.; Vaezipoor, P.; and Dilkina, B. 2022.  Finding Backdoors to Integer Programs: A Monte Carlo Tree Search Framework.  In _Thirty-Sixth AAAI Conference on Artificial Intelligence, AAAI 2022, Thirty-Fourth Conference on Innovative Applications of Artificial Intelligence, IAAI 2022, The Twelveth Symposium on Educational Advances in Artificial Intelligence, EAAI 2022 Virtual Event, February 22 - March 1, 2022_ , 3786–3795. AAAI Press.
  * Kirlik and Sayın (2014) Kirlik, G.; and Sayın, S. 2014.  A new algorithm for generating all nondominated solutions of multiobjective discrete optimization problems.  _European Journal of Operational Research_ , 232(3): 479–488.
  * Lambrinidis and Tsantili-Kakoulidou (2021) Lambrinidis, G.; and Tsantili-Kakoulidou, A. 2021.  Multi-objective optimization methods in novel drug design.  _Expert Opinion on Drug Discovery_ , 16(6): 647–658.
  * Lindauer et al. (2021) Lindauer, M.; Eggensperger, K.; Feurer, M.; Biedenkapp, A.; Deng, D.; Benjamins, C.; Sass, R.; and Hutter, F. 2021.  SMAC3: A Versatile Bayesian Optimization Package for Hyperparameter Optimization.  arXiv:2109.09831.
  * Liu, Yan, and Jin (2023) Liu, S.; Yan, X.; and Jin, Y. 2023.  End-to-end Pareto set prediction with graph neural networks for multi-objective facility location.  In _International Conference on Evolutionary Multi-Criterion Optimization_ , 147–161. Springer.
  * Lu et al. (2000) Lu, Y.; Jain, J.; Clarke, E.; and Fujita, M. 2000.  Efficient variable ordering using aBDD based sampling.  In _Proceedings of the 37th Annual Design Automation Conference_ , 687–692.
  * Ozlen, Burton, and MacRae (2014) Ozlen, M.; Burton, B. A.; and MacRae, C. A. 2014.  Multi-objective integer programming: An improved recursive algorithm.  _Journal of Optimization Theory and Applications_ , 160: 470–482.
  * Parragh and Tricoire (2019) Parragh, S. N.; and Tricoire, F. 2019.  Branch-and-bound for bi-objective integer programming.  _INFORMS Journal on Computing_ , 31(4): 805–822.
  * Pedregosa et al. (2011) Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; and Duchesnay, E. 2011.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12: 2825–2830.
  * Przybylski and Gandibleux (2017) Przybylski, A.; and Gandibleux, X. 2017.  Multi-objective branch and bound.  _European Journal of Operational Research_ , 260(3): 856–872.
  * Rice and Kulhari (2008) Rice, M.; and Kulhari, S. 2008.  A survey of static variable ordering heuristics for efficient BDD/MDD construction.  _University of California, Tech. Rep_ , 130.
  * Rudell (1993) Rudell, R. 1993.  Dynamic variable ordering for ordered binary decision diagrams.  In _Proceedings of 1993 International Conference on Computer Aided Design (ICCAD)_ , 42–47. IEEE.
  * Schede et al. (2022) Schede, E.; Brandt, J.; Tornede, A.; Wever, M.; Bengs, V.; Hüllermeier, E.; and Tierney, K. 2022.  A survey of methods for automated algorithm configuration.  _Journal of Artificial Intelligence Research_ , 75: 425–487.
  * Song et al. (2018) Song, Z.; Chen, X.; Luo, X.; Wang, M.; and Dai, G. 2018.  Multi-objective optimization of agile satellite orbit design.  _Advances in Space Research_ , 62(11): 3053–3064.
  * Sourd and Spanjaard (2008) Sourd, F.; and Spanjaard, O. 2008.  A multiobjective branch-and-bound framework: Application to the biobjective spanning tree problem.  _INFORMS Journal on Computing_ , 20(3): 472–484.
  * Tangpattanakul, Jozefowiez, and Lopez (2012) Tangpattanakul, P.; Jozefowiez, N.; and Lopez, P. 2012.  Multi-objective optimization for selecting and scheduling observations by agile earth observing satellites.  In _Parallel Problem Solving from Nature-PPSN XII: 12th International Conference, Taormina, Italy, September 1-5, 2012, Proceedings, Part II 12_ , 112–121. Springer.
  * Wegener (2000) Wegener, I. 2000.  _Branching Programs and Binary Decision Diagrams: Theory and Applications_ , volume 4.  SIAM.
  * Wu et al. (2022) Wu, Y.; Song, W.; Cao, Z.; Zhang, J.; Gupta, A.; and Lin, M. 2022.  Graph Learning Assisted Multi-Objective Integer Programming.  _Advances in Neural Information Processing Systems_ , 35: 17774–17787.
  * Yu et al. (2000) Yu, Y.; Zhang, J.; Cheng, G.; Schell, M.; and Okunieff, P. 2000.  Multi-objective optimization in radiotherapy: applications to stereotactic radiosurgery and prostate brachytherapy.  _Artificial Intelligence in Medicine_ , 19(1): 39–51.

##  Appendix A Appendix

### Feature importance plots for size-independent ML models

Figure 8: Feature importance plot for ML+A model. Refer to the caption of
Figure 7 for an explanation of the feature naming conventions. Figure 9:
Feature importance plot for ML+AC model. Features that start with “c_” are
context features; see Table 2.


a}(s-a)^{b}(\mathcal{M}\varphi)(-sL)Z_{\mathfrak{a}}(sL,f)$ is equal to

$\int_{\mathbf{T}_{\operatorname{NS}}(\mathbb{A}_{F})/T_{\operatorname{NS}}(F)}\mathbf{1}_{\operatorname{\overline{Eff}}_{(X;\mathbf{A})}}(\log|y|_{\mathbf{T}_{\operatorname{NS}}})|y|_{\mathbf{T}_{\operatorname{NS}}}^{-L}d\tau_{\mathbf{T}_{\operatorname{NS}}}(y)\int_{\mathbf{Y}(\mathbb{A}_{F})_{\mathbf{A}}}f(\pi_{\mathfrak{a}}(x))\varphi(||x||)dx.$

Make the change of variables $z=y\cdot x$ and $t=y\cdot||z||$. We then have
the equalities of measures $dx=|y|_{\Theta}^{-\mathbf{1}}dz$ and
$d\tau_{\mathbf{T}_{\operatorname{NS}}}(t)=d\tau_{\mathbf{T}_{\operatorname{NS}}}(y)$.

We have
$|y|_{\Theta}^{-\mathbf{1}}=|y|_{\mathbf{T}_{\operatorname{NS}}}^{K_{(X;\mathbf{A})}}=|y|_{\mathbf{T}_{\operatorname{NS}}}^{-aL}$
by (22). By definition we also have $||z||=y^{-1}\cdot||x||$ and
$\left|t\cdot||z||^{-1}\right|_{\mathbf{T}_{\operatorname{NS}}}=\mathbf{H}(t\cdot
z)$. The above integral is thus equal to

$\int_{\mathbf{T}_{\operatorname{NS}}(\mathbb{A}_{F})}\varphi(t)\int_{\mathbf{Y}(\mathbb{A}_{F})_{\mathbf{A}}/T_{\operatorname{NS}}(F)}\mathbf{1}_{\operatorname{\overline{Eff}}_{(X;\mathbf{A})}}(\log\mathbf{H}(t\cdot
z))f(\pi_{\mathfrak{a}}(z))\mathbf{H}(t\cdot
z)^{-(a+1)L}dzd\tau_{\mathbf{T}_{\operatorname{NS}}}(t).$

Now note that since $K_{(X;\mathbf{A})}=-aL$ and the definition
$\mathbf{H}(t\cdot z)^{-aL}dz$ is equal to the Tamagawa measure
$\tau_{\mathbf{Y},\mathbf{H}(t\cdot z)^{-aL}}$. Note now that
$\mathbf{H}(t\cdot z)=|t|_{\mathbf{T}_{\operatorname{NS}}}\mathbf{H}(z)$, if
one combines this with Lemma 6.8 we deduce that

$\begin{split}&\int_{\mathbf{Y}(\mathbb{A}_{F})_{\mathbf{A}}/T_{\operatorname{NS}}(F)}\mathbf{1}_{\operatorname{\overline{Eff}}_{(X;\mathbf{A})}}(\log\mathbf{H}(t\cdot
z))f(\pi_{\mathfrak{a}}(z))\mathbf{H}(t\cdot z)^{-(a+1)L}dz=\\\
&|t|_{\mathbf{T}_{\operatorname{NS}}}^{-aL}\int_{\mathbf{Y}(\mathbb{A}_{F})/T_{\operatorname{NS}}(F)}\mathbf{1}_{\operatorname{\overline{Eff}}_{(X;\mathbf{A})}}(\log\mathbf{H}(z))f(\pi_{\mathfrak{a}}(z))\mathbf{H}(z)^{-L}d\tau_{\mathbf{Y},\mathbf{H}^{-aL}}(z).\end{split}$

Combining the above we see that $(\mathcal{M}\varphi)(-aL)\lim_{s\to
a}(s-a)^{b}Z_{\mathfrak{a}}(sL,f)$ is equal to

$\int_{\mathbf{T}_{\operatorname{NS}}(\mathbb{A}_{F})}\varphi(t)|t|_{\mathbf{T}_{\operatorname{NS}}}^{-aL}d\tau_{\mathbf{T}_{\operatorname{NS}}}(t)\int_{\mathbf{Y}(\mathbb{A}_{F})}\mathbf{1}_{\operatorname{\overline{Eff}}_{(X;\mathbf{A})}}(\log\mathbf{H}(z))f(\pi_{\mathfrak{a}}(z))\mathbf{H}(z)^{-L}d\tau_{\mathbf{Y},\mathbf{H}^{-aL}}(z).$

By Remark 7.2 we may assume that
$(\mathcal{M}\varphi)(-aL)=\int_{\mathbf{T}_{\operatorname{NS}}(\mathbb{A}_{F})}\varphi(t)|t|_{\mathbf{T}_{\operatorname{NS}}}^{-aL}d\tau_{\mathbf{T}_{\operatorname{NS}}}(t)\neq
0$. The lemma follows after dividing both sides by
$(\mathcal{M}\varphi)(-aL)$. ∎

Combining this lemma, Lemma 7.11 and (24) we deduce the following theorem.

######  Theorem 8.6.

The series $Z(sL,f)$ converges absolutely for $\mathrm{Re}(s)$ sufficiently
large. There exists an $\varepsilon>0$ such that $Z(sL,f)$ has a meromorphic
continuation to the region $\mathrm{Re}(s)>a-\varepsilon$. Moreover, this
meromorphic continuation is holomorphic up to a single pole of order $b$ at
$s=a$ and $\lim_{s\to a}(s-a)^{b}Z(sL,f)$ is equal to

$\sum_{\mathbf{a}\in\operatorname{H}^{1}(F,T_{\operatorname{NS}})}\int_{\mathbf{Y}(\mathbb{A}_{F})}\mathbf{1}_{\operatorname{\overline{Eff}}_{(X;\mathbf{A})}}(\log\mathbf{H}(z))f(\pi_{\mathfrak{a}}(z))\mathbf{H}(z)^{-L}d\tau_{\mathbf{Y},\mathbf{H}^{-aL}}(z).$

This theorem implies that the integral Manin conjecture hold with our
prediction for the constant, see Conjecture 6.6.

######  Theorem 8.7.

Let $L\in\operatorname{Pic}(X;\mathbf{A})$ be a line bundle such that
$aL=-K_{(X;\mathbf{A})}$ for $a>0$. Let
$b=\emph{rk}\operatorname{Pic}(X;\mathbf{A})$. Let $H$ be a height function on
$L$, and $\mathbf{H}$ an extension to $\mathbf{Y}$. Let
$f:(X;\mathbf{A})(\mathbb{A}_{F})\to\mathbb{C}$ be any compactly supported
continuous function.

Assume that $c((X;\mathbf{A}),f,H)$ is non-zero, which we define to be the
finite sum

$\frac{1}{a(b-1)!}\sum_{\mathbf{a}\in\operatorname{H}^{1}(F,T_{\operatorname{NS}})}\int_{\mathbf{Y}(\mathbb{A}_{F})}\mathbf{1}_{\operatorname{\overline{Eff}}_{(X;\mathbf{A})}}(\log\mathbf{H}(z))f(\pi_{\mathfrak{a}}(z))H^{-1}(z)d\tau_{\mathbf{Y},H^{-a}}(z).$

We then have as $T\to\infty$ that

$\sum_{\begin{subarray}{c}P\in T(F)\\\ H(P)\leq T\end{subarray}}f(P)\sim
c((X;\mathbf{A}),f,H)T^{a}(\log T)^{b-1}.$

###### Proof.

By Proposition 6.14 it suffices to prove this for $\mathbf{H}$ and $f$ as in
§7.1 This case follows from Theorem 8.6 and a Tauberian theorem, e.g. [Ten15,
Thm. II.7.13]. ∎

Theorem 1.2 in the introduction follows from this theorem. Indeed, in that
case $a=1$. Also $f$ is non-negative and the measure
$\tau_{\mathbf{Y},H^{-a}}$ is supported on
$\mathbf{Y}(\mathbb{A}_{F})_{\mathbf{A}}$. So $c((X;\mathbf{A}),f,H)$ will be
non-zero as long as $f$ is non-zero on any element of the image of
$\mathbf{Y}(\mathbb{A}_{F})_{\mathbf{A}}$ along one of the
$\pi_{\mathfrak{a}}$. This image is equal to
$U(\mathbb{A}_{F})_{\mathbf{A}}^{\operatorname{Br}_{1}}$ by Theorem 5.10.

We can deduce the following corollary, implying Theorem 1.3.

######  Corollary 8.8.

The algebraic Brauer-Manin obstruction is the only obstruction to strong
approximation for $(X;\mathbf{A})$.

###### Proof.

Let $V\subset(X;\mathbf{A})(\mathbb{A}_{F})$ be an open such that $V\cap
U(\mathbb{A}_{F})_{\mathbf{A}}^{\operatorname{Br}_{1}}\neq\emptyset$. Let
$f:(X;\mathbf{A})(\mathbb{A}_{F})\to\mathbb{R}_{\geq 0}$ be continuous
function with compact support $\text{supp}(f)\subset V$ and such that
$\text{supp}(f)\cap
U(\mathbb{A}_{F})_{\mathbf{A}}^{\operatorname{Br}_{1}}\neq\emptyset$. We have
$c((X;\mathbf{A}),f,H)\neq 0$ by the same reasoning as above so Theorem 8.7
implies that $T(F)\cap V\neq\emptyset$.

We deduce that the closure $\overline{T(F)}$ in
$(X;\mathbf{A})(\mathbb{A}_{F})$ contains
$U(\mathbb{A}_{F})_{\mathbf{A}}^{\operatorname{Br}_{1}}$. ∎

###  8.4. Other heights

Let $L\in\operatorname{\overline{Eff}}_{(X;\mathbf{A})}^{\circ}$ now be any
big height. For each $\alpha\in\mathcal{R}$ let $H_{L,\alpha}$ be the
restriction of $H$ to $(U_{L};\mathbf{B})$ and $\mathbf{H}_{L,\alpha}$ the
height corresponding to an extension of the adelic metric defining $H$ to the
universal torsor $\mathbf{Y}_{L}\to(U_{L};\mathbf{B})$. This exists because of
Lemma 4.12.

######  Lemma 8.9.

The limit $\lim_{s\to a}(s-a)^{b}Z_{\mathfrak{a}}(sL,f)$ is equal to

$\sum_{\alpha\in\mathcal{R}}\int_{\mathbf{Y}_{L}(\mathbb{A}_{F})}\mathbf{1}_{\operatorname{\overline{Eff}}_{(U_{L};\mathbf{B})}}(\log\mathbf{H}_{L,\alpha}(z))f(\pi_{\mathfrak{a}}(z\alpha))H_{L,\alpha}^{-1}(z)d\tau_{\mathbf{Y},H_{L,\alpha}^{-a}}(z).$

###### Proof.

Let $W:=U\setminus\left(\cup_{Z_{i}\in D_{\operatorname{adj}}}Z_{i}\right)$
and for each place $v$ let $W_{v}:=U_{A_{v}}\setminus\left(\cup_{Z_{i}\in
D_{\operatorname{adj}}\cup B_{c}}Z_{i}\right)$ and let
$\mathcal{W}\subset\mathcal{U}$ be the complement of the closures of the
$Z_{i}$ for $Z_{i}\in D_{\operatorname{adj}}$. Let $S\subset S_{1}\subset
S_{2}\subset\cdots\subset\Omega_{F}$ be a sequence of subsets of places such
that $\cup_{n}S_{n}=\Omega_{F}$. For each $n\in\mathbb{N}$ and
$v\in\Omega_{F}$ let $g_{v,n}:W_{v}(F_{v})\to\mathbb{R}_{\geq 0}$ be a
compactly supported function which is smooth, resp. locally constant, if $v$
is archimedean, resp. non-archimedean, satisfying the following properties:
$g_{v,n-1}\leq g_{v,n}\leq f_{v}$, $g_{v,n}$ is the indicator function of
$\mathcal{W}(\mathcal{O}_{v})$ if $v\not\in S_{n}$ and
$\lim_{n\to\infty}g_{v,n}=f_{v}$.

Let $g_{n}=\prod_{v}g_{n,v}$ and $h_{n}=\left(\prod_{v\in
S_{n}}f_{v}-\prod_{v\in S_{n}}g_{n,v}\right)\prod_{v\not\in S_{n}}f_{v}$. The
function $g_{n}$, resp. $h_{n}$, is compactly supported continuous on
$\prod_{v}W_{v}(F_{v})$, resp. $(X;\mathbf{A})(\mathbb{A}_{F})$. Note that
$g_{n}\leq f\leq g_{n}+h_{n}$ and $\lim_{n\to\infty}g_{n}=f$ and
$\lim_{n\to\infty}h_{n}=0$.

Let $Z_{\mathfrak{a}}(sL,g_{n})=\sum_{P\in T(F)}g_{n}(P)H(P)^{-sL}$ and
$Z_{\mathfrak{a}}(sL,h_{n})=\sum_{P\in T(F)}h_{n}(P)H(P)^{-sL}$. Both
integrals converge absolutely for $\mathrm{Re}(s)>a$ by comparison with
$Z_{\mathfrak{a}}(sL,f)$ and Corollary 7.12.

We claim that $\liminf_{s\to
a,n\to\infty}|(s-a)^{b}Z_{\mathfrak{a}}(sL,h_{n})|=0$. Note that $h_{n}$ can
be written as a finite linear combination of functions as considered in
§7.1.2. By Lemma 8.4 we find that the limit of $Z_{\mathfrak{a}}(sL,h_{n})$ as
$s\to a$ is

$\sum_{\alpha\in\mathcal{R}}\int_{\mathbf{Y}_{L}(\mathbb{A}_{F})}\mathcal{H}_{\varphi}(\alpha
x;sL)g_{n}(\pi_{\mathfrak{a}}(\alpha x))dx.$

The above integral goes to $0$ as $n\to\infty$ by monotone convergence from
which the claim follows. We deduce from $g_{n}\leq f\leq g_{n}+h_{n}$ that
$Z_{\mathfrak{a}}(sL,f)=\lim_{n\to\infty}\lim_{s\to
a}Z_{\mathfrak{a}}(sL,g_{n})$.

Note that $U_{L,B_{v}}\subset W_{v}$ is a closed immersion. The function
$g_{n}$ is compactly supported on $\prod_{v}W_{v}(F_{v})$ so for each
$\alpha\in\mathcal{R}$ its restriction to the closed subset
$(U_{L};\mathbf{B})(\mathbb{A}_{F})\pi_{\mathfrak{a}}(\alpha)$ is also
compactly supported. Moreover, since it is compactly supported its restriction
to $(U_{L};\mathbf{B})(\mathbb{A}_{F})\pi_{\mathfrak{a}}(\alpha)$ is zero for
all but finitely many $\alpha\in\mathcal{R}$.

By Lemma 6.4 we can apply Lemma 8.5 to $(U_{L},\mathbf{B})$ and deduce that
for all $\alpha\in\mathcal{R}$ the limit $\lim_{s\to a}\sum_{P\in
T_{L}(F)\pi_{\mathfrak{a}}(\alpha)}g_{n}(P)H(P)^{-sL}$ is equal to

$\int_{\mathbf{Y}_{L}(\mathbb{A}_{F})}\mathbf{1}_{\operatorname{\overline{Eff}}_{(U_{L};\mathbf{B})}}(\log\mathbf{H}_{L,\alpha}(z))g_{n}(\pi_{\mathfrak{a}}(z\alpha))H_{L,\alpha}^{-1}(z)d\tau_{\mathbf{Y},H_{L,\alpha}^{-a}}(z).$

We conclude by summing over $\alpha\in\mathcal{R}$, noting that this is a
finite sum, and taking the increasing limit as $n\to\infty$. ∎

Combining this lemma, Lemma 7.11 and (24) we deduce the following theorem.

######  Theorem 8.10.

The series $Z(sL,f)$ converges absolutely for $\mathrm{Re}(s)$ sufficiently
large. There exists an $\varepsilon>0$ such that $Z(sL,f)$ has a meromorphic
continuation to the region $\mathrm{Re}(s)>a-\varepsilon$. Moreover, this
meromorphic continuation is holomorphic up to a single pole of order $b$ at
$s=a$ and $\lim_{s\to a}(s-a)^{b}Z(sL,f)$ is equal to

$\sum_{\mathbf{a}\in\operatorname{H}^{1}(F,T_{\operatorname{NS}})}\sum_{\alpha\in\mathcal{R}}\int_{\mathbf{Y}_{L}(\mathbb{A}_{F})}\mathbf{1}_{\operatorname{\overline{Eff}}_{(U_{L};\mathbf{B})}}(\log\mathbf{H}_{L,\alpha}(z))f(\pi_{\mathfrak{a}}(z\alpha))H_{L,\alpha}^{-1}(z)d\tau_{\mathbf{Y},H_{L,\alpha}^{-a}}(z).$

The same reasoning as for Theorem 8.7 shows that this implies the following.

######  Theorem 8.11.

Let $L\in\operatorname{\overline{Eff}}_{(X;\mathbf{A})}^{\circ}$ be big. Let
$H$ be a height function on $L$, and $\mathbf{H}_{L,\alpha}$ an extension of
$H$ to the universal torsor
$\pi_{\mathfrak{a},\alpha}:\mathbf{Y}_{L}\to(U_{L},\mathbf{B})\cdot\pi_{\mathfrak{a}}(\alpha)$.
Let $f:(X;\mathbf{A})(\mathbb{A}_{F})\to\mathbb{C}$ be any compactly supported
continuous function.

Assume that $c((X;\mathbf{A}),f,H)$ is non-zero, which we define to be the sum

$\frac{1}{a(b-1)!}\sum_{\mathbf{a}\in\operatorname{H}^{1}(F,T_{\operatorname{NS}})}\sum_{\alpha\in\mathcal{R}}\int_{\mathbf{Y}_{L}(\mathbb{A}_{F})}\mathbf{1}_{\operatorname{\overline{Eff}}_{(U_{L};\mathbf{B})}}(\log\mathbf{H}_{L,\alpha}(z))f(\pi_{\mathfrak{a}}(z\alpha))H_{L,\alpha}^{-1}(z)d\tau_{\mathbf{Y},H_{L,\alpha}^{-a}}(z).$

We then have as $T\to\infty$ that

$\sum_{\begin{subarray}{c}P\in T(F)\\\ H(P)\leq T\end{subarray}}f(P)\sim
c((X;\mathbf{A}),f,H)T^{a}(\log T)^{b-1}.$

We deduce Theorem 1.1 by using the partition of unity arguments at the
beginning of §6.

######  Remark 8.12.

If $L$ is adjoint rigid then $\pi_{\mathfrak{a}}(\alpha)\in
T(F)=T_{L}(F)\subset(U_{L},\mathbf{B})$. The maps
$\pi_{\mathfrak{a},\alpha}:\mathbf{Y}_{L}\to(U_{L},\mathbf{B})$ as
$\mathbf{a}$ ranges over all elements of
$\operatorname{H}^{1}(F,T_{\operatorname{NS}})$ and $\alpha$ ranges over all
elements of $\mathcal{R}$ then correspond to all universal torsors of
$(U_{L},\mathbf{B})$ by Lemma 5.9. Theorem 8.11 thus agrees with Conjecture
6.6 in this case.

## References

  * [BCHM10] Caucher Birkar, Paolo Cascini, Christopher D. Hacon, and James McKernan.  Existence of minimal models for varieties of log general type.  J. Amer. Math. Soc., 23(2):405–468, 2010.
  * [Bir62] B. J. Birch.  Forms in many variables.  Proc. Roy. Soc. London Ser. A, 265:245–263, 1961/62.
  * [BM90] Victor V. Batyrev and Yuri I. Manin.  Sur le nombre des points rationnels de hauteur borné des variétés algébriques.  Math. Ann., 286(1-3):27–43, 1990.
  * [Bou11] David Bourqui.  Fonction zêta des hauteurs des variétés toriques non déployées.  Mem. Amer. Math. Soc., 211(994):viii+151, 2011.
  * [BR95] Mikhail Borovoi and Zeév Rudnick.  Hardy-Littlewood varieties and semisimple groups.  Invent. Math., 119(1):37–66, 1995.
  * [Bro07] T. D. Browning.  An overview of Manin’s conjecture for del Pezzo surfaces.  In Analytic number theory, volume 7 of Clay Math. Proc., pages 39–55. Amer. Math. Soc., Providence, RI, 2007.
  * [BT95] Victor V. Batyrev and Yuri Tschinkel.  Rational points of bounded height on compactifications of anisotropic tori.  Internat. Math. Res. Notices, (12):591–635, 1995.
  * [BT96] Victor V. Batyrev and Yuri Tschinkel.  Height zeta functions of toric varieties.  J. Math. Sci., 82(1):3220–3239, 1996.  Algebraic geometry, 5.
  * [BT98a] Victor V. Batyrev and Yuri Tschinkel.  Manin’s conjecture for toric varieties.  J. Algebraic Geom., 7(1):15–53, 1998.
  * [BT98b] Victor V. Batyrev and Yuri Tschinkel.  Tamagawa numbers of polarized algebraic varieties.  Number 251, pages 299–340. 1998.  Nombre et répartition de points de hauteur bornée (Paris, 1996).
  * [Cho19] Dylon Chow.  The Distribution of Integral Points on the Wonderful Compactification by Height.  ProQuest LLC, Ann Arbor, MI, 2019.  Thesis (Ph.D.)–University of Illinois at Chicago.
  * [CLT10a] Antoine Chambert-Loir and Yuri Tschinkel.  Igusa integrals and volume asymptotics in analytic and adelic geometry.  Confluentes Math., 2(3):351–429, 2010.
  * [CLT10b] Antoine Chambert-Loir and Yuri Tschinkel.  Integral points of bounded height on toric varieties.  ArXiv preprint: arXiv:1006.3345, 2010.
  * [CLT12] Antoine Chambert-Loir and Yuri Tschinkel.  Integral points of bounded height on partial equivariant compactifications of vector groups.  Duke Math. J., 161(15):2799–2836, 2012.
  * [CTS87] Jean-Louis Colliot-Thélène and Jean-Jacques Sansuc.  La descente sur les variétés rationnelles. II.  Duke Math. J., 54(2):375–492, 1987.
  * [CTS21] Jean-Louis Colliot-Thélène and Alexei N. Skorobogatov.  The Brauer-Grothendieck group, volume 71 of Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge. A Series of Modern Surveys in Mathematics [Results in Mathematics and Related Areas. 3rd Series. A Series of Modern Surveys in Mathematics].  Springer, Cham, 2021.
  * [CX18] Yang Cao and Fei Xu.  Strong approximation with Brauer-Manin obstruction for toric varieties.  Ann. Inst. Fourier (Grenoble), 68(5):1879–1908, 2018.
  * [DW22] Ulrich Derenthal and Florian Wilsch.  Integral points on singular del pezzo surfaces.  Journal of the Institute of Mathematics of Jussieu, pages 1–36, 2022.
  * [DY22] Ratko Darda and Takehiko Yasuda.  Torsors for finite group schemes of bounded height.  ArXiv preprint: arXiv:2207.03642, 2022.
  * [FMT89] Jens Franke, Yuri I. Manin, and Yuri Tschinkel.  Rational points of bounded height on Fano varieties.  Invent. Math., 95(2):421–435, 1989.
  * [Hur11] Mathieu Huruguen.  Toric varieties and spherical embeddings over an arbitrary field.  J. Algebra, 342:212–234, 2011.
  * [Lou15] Daniel Loughran.  Rational points of bounded height and the Weil restriction.  Israel J. Math., 210(1):47–79, 2015.
  * [Lou18] Daniel Loughran.  The number of varieties in a family which contain a rational point.  J. Eur. Math. Soc. (JEMS), 20(10):2539–2588, 2018.
  * [Mil06] James S. Milne.  Arithmetic duality theorems.  BookSurge, LLC, Charleston, SC, second edition, 2006.
  * [NSW08] Jürgen Neukirch, Alexander Schmidt, and Kay Wingberg.  Cohomology of number fields, volume 323 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, second edition, 2008.
  * [Ono61] Takashi Ono.  Arithmetic of algebraic tori.  Ann. of Math. (2), 74:101–139, 1961.
  * [Ono63] Takashi Ono.  On the Tamagawa number of algebraic tori.  Ann. of Math. (2), 78:47–73, 1963.
  * [Pey95] Emmanuel Peyre.  Hauteurs et mesures de Tamagawa sur les variétés de Fano.  Duke Math. J., 79(1):101–218, 1995.
  * [Pey98] Emmanuel Peyre.  Terme principal de la fonction zêta des hauteurs et torseurs universels.  Number 251, pages 259–298. 1998.  Nombre et répartition de points de hauteur bornée (Paris, 1996).
  * [Poo17] Bjorn Poonen.  Rational points on varieties, volume 186 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2017.
  * [Sal98] Per Salberger.  Tamagawa measures on universal torsors and points of bounded height on Fano varieties.  Astérisque, (251):91–258, 1998.  Nombre et répartition de points de hauteur bornée (Paris, 1996).
  * [Ski97] C. M. Skinner.  Forms over number fields and weak approximation.  Compositio Math., 106(1):11–29, 1997.
  * [Sko01] Alexei Skorobogatov.  Torsors and rational points, volume 144 of Cambridge Tracts in Mathematics.  Cambridge University Press, Cambridge, 2001.
  * [Sta18] The Stacks Project Authors.  Stacks Project.  https://stacks.math.columbia.edu, 2018.
  * [TBT13] Ramin Takloo-Bighash and Yuri Tschinkel.  Integral points of bounded height on compactifications of semi-simple groups.  Amer. J. Math., 135(5):1433–1448, 2013.
  * [Ten15] Gérald Tenenbaum.  Introduction to analytic and probabilistic number theory, volume 163 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, third edition, 2015.  Translated from the 2008 French edition by Patrick D. F. Ion.
  * [Wil22] Florian Wilsch.  Integral points of bounded height on a certain toric variety.  ArXiv preprint: arXiv:2202.10909, 2022.
  * [Wil23] Florian Wilsch.  Integral points of bounded height on a log Fano threefold.  Int. Math. Res. Not. IMRN, (8):6780–6808, 2023.