11institutetext: Heinz Nixdorf Institute and Department of Computer Science,
Paderborn University, Warburger Str. 100, 33100 Paderborn, Germany
11email: {alexander.tornede, marcel.wever<EMAIL_ADDRESS>
# Extreme Algorithm Selection with
Dyadic Feature Representation††thanks: This paper is published at the
Discovery Science 2020 conference. The final authenticated version is
available online at https://doi.org/10.1007/978-3-030-61527-7_21.

Alexander Tornede 11 0000-0002-2415-2186 Marcel Wever 11 0000-0001-9782-6818
Eyke Hüllermeier 11 0000-0002-9944-4108

###### Abstract

Algorithm selection (AS) deals with the automatic selection of an algorithm
from a fixed set of candidate algorithms most suitable for a specific instance
of an algorithmic problem class, e.g., choosing solvers for SAT problems.
Benchmark suites for AS usually comprise candidate sets consisting of at most
tens of algorithms, whereas in algorithm configuration (AC) and combined
algorithm selection and hyperparameter optimization (CASH) problems the number
of candidates becomes intractable, impeding to learn effective meta-models and
thus requiring costly online performance evaluations. In this paper, we
propose the setting of extreme algorithm selection (XAS), which, despite
assuming limited time resources and hence excluding online evaluations at
prediction time, allows for considering thousands of candidate algorithms and
thereby facilitates meta learning. We assess the applicability of state-of-
the-art AS techniques to the XAS setting and propose approaches leveraging a
dyadic representation, in which both problem instances and algorithms are
described in terms of feature vectors. We find this approach to significantly
improve over the current state of the art in various metrics.

###### Keywords:

extreme algorithm selection dyadic ranking surrogate model

##  1 Introduction

Algorithm selection (AS) refers to a specific recommendation task, in which
the choice alternatives are algorithms: Given a set of candidate algorithms to
choose from, and a specific instance of a problem class, such as SAT or
integer optimization, the task is to select or recommend an algorithm that
appears to be most suitable for that instance, in the sense of performing best
in terms of criteria such as runtime, solution quality, etc. Hitherto
practical applications of AS, as selecting a SAT solver for a logical formula,
typically comprise candidate sets consisting of at most tens of algorithms,
and this is also the order of magnitude that is found in standard AS benchmark
suites such as ASlib [2].

This is in contrast with the problem of combined algorithm selection and
hyperparameter optimization (CASH) [24] as considered in automated machine
learning (AutoML), where the number of potential candidates is very large and
potentially infinite [24, 6, 16]. Corresponding methods heavily rely on
computationally extensive search procedures combined with costly online
evaluations of the performance measure to optimize for, since learning
effective meta models for an instantaneous recommendation becomes infeasible.

In this paper, we propose _extreme algorithm selection_ (XAS) as a novel
setting in-between traditional AS and AC/CASH, which is motivated by
application scenarios characterized by

  * •

the demand for prompt recommendations in quasi real time,

  * •

an extremely large (though still finite) set of candidate algorithms.

An example is the scenario of “On-the-fly computing” [10], including “On-the-
fly machine learning” [17] as one of its instantiations, where users can
request online (machine learning) software services customized towards their
needs. Here, users are unwilling to wait for several hours until their service
is ready, but rather claim a result quickly. Hence, for providing a first
version of an appropriate service, costly search and online evaluations are
not affordable. As will be seen, XAS offers a good compromise solution:
Although it allows for the consideration of extremely many candidate
solutions, and even offers the ability to recommend configurations that have
never been encountered so far, it is still amenable to AS techniques and
avoids costly online evaluations.

In a sense, XAS relates to standard AS as the emerging topic of extreme
classification (XC) [1] relates to standard multi-class classification.
Similar to XC, the problem of learning from sparse data is a major challenge
for XAS: For a single algorithm, there are typically only observations for a
few instances. In this paper, we propose a benchmark dataset for XAS and
investigate the ability of state-of-the-art AS approaches to deal with this
sparsity and to scale with the size of candidate sets. Furthermore, to support
more effective learning from sparse data, we propose methods based on “dyadic”
feature representations, in which both problem instances and algorithms are
represented in terms of feature vectors. In an extensive experimental study,
we find these methods to yield significant improvements.

##  2 From Standard to Extreme Algorithm Selection

In the standard (per-instance) algorithm selection setting, first introduced
in [20], we are interested in finding a mapping
$s:\mathcal{I}\longrightarrow\mathcal{A}$, called algorithm selector. Given an
instance $i$ from the instance space $\mathcal{I}$, the latter selects the
algorithm $a^{*}$ from a set of candidate algorithms $\mathcal{A}$, optimizing
a performance measure
$m:\mathcal{I}\times\mathcal{A}\longrightarrow\mathbb{R}$. Furthermore, $m$ is
usually costly to evaluate. The optimal selector is called _oracle_ and is
defined as

$s^{*}(i)\mathrel{\vbox{\hbox{\footnotesize.}\hbox{\footnotesize.}}}=\arg\max_{a\in\mathcal{A}}\mathbb{E}\big{[}\,m(i,a)\,\big{]}$
(1)

for all $i\in\mathcal{I}$. The expectation operator $\mathbb{E}$ accounts for
any randomness in the application of the algorithm — in the non-deterministic
case, the result of applying $a$ to $i$, and hence the values of the
performance measure, are random variables.

Most AS approaches leverage machine learning techniques, in one way or another
learning a surrogate (regression) model
$\widehat{m}:\mathcal{I}\times\mathcal{A}\longrightarrow\mathbb{R}$, which is
fast to evaluate and thus allows one to compute a selector
$\widehat{s}:\mathcal{I}\longrightarrow\mathcal{A}$ by

$\widehat{s}(i)\mathrel{\vbox{\hbox{\footnotesize.}\hbox{\footnotesize.}}}=\arg\max_{a\in\mathcal{A}}\widehat{m}(i,a)\,.$
(2)

In order to infer such a model, we usually assume the existence of a set of
training instances $\mathcal{I}_{D}\subset\mathcal{I}$ for which we have
instantaneous access to the associated performances of some or often all
algorithms in $\mathcal{A}$ according to $m$.

The XAS setting distinguishes itself from the standard AS setting by two
important properties. Firstly, we assume that the set of candidate algorithms
$\mathcal{A}$ is extremely large. Thus, approaches need to be able to scale
well with the size of $\mathcal{A}$. Secondly, due to the size of
$\mathcal{A}$, we can no longer reasonably assume to have evaluations for each
algorithm on each training instance. Instead, we assume that the training
matrix spanned by the training instances and algorithms is only sparsely
filled. In fact, we might even have algorithms without any evaluations at all.
Hence, suitable approaches need to be able to learn from very few data and to
tackle the problem of “zero-shot learning” [29].

Similarly, the XAS setting differs from the AC and CASH settings in two main
points. Firstly, dealing with real-valued hyperparameters, the set of
(configured) algorithms $\mathcal{A}$ is generally assumed to be infinite in
both AC and CASH, whereas $\mathcal{A}$ is still finite (even if extremely
large) in XAS. More importantly, in both AC and CASH, one usually assumes
having time to perform online evaluations of solution candidates at
recommendation time. However, as previously mentioned, this is not the case in
XAS, where instantaneous recommendations are required. Hence, the XAS setting
significantly differs from the AS, AC, and CASH settings. A summary of the
main characteristics of these settings is provided in Table 1.

Table 1: Overview of the characteristics of the problem settings we
distinguish.

Characteristics/Setting | AS | XAS | AC | CASH
---|---|---|---|---
Size of $\mathcal{A}$ |  at most tens |  extremely many |  potentially infinite |  potentially infinite
Training data | complete | sparse | mostly not present | mostly not present
Online evaluations | no | no | yes | yes

##  3 Exploiting Instance Features

Instance-specific AS is based on the assumption that instances can be
represented in terms of feature information. For this purpose,
$f_{I}:\mathcal{I}\longrightarrow\mathbb{R}^{k}$ denotes a function
representing instances as $k$-dimensional, real-valued feature vectors, which
can be used to learn a surrogate model (2). This can be done based on
different types of data and using different loss functions.

###  3.1 Regression

The most common approach is to tackle AS as a regression problem, i.e., to
construct a regression dataset for each algorithm, where entries consist of an
instance representation and the associated performance of the algorithm at
question. Accordingly, the dataset associated with algorithm $a\in\mathcal{A}$
consists of tuples of the form $\big{(}f_{I}(i),m(i,a)\big{)}$, created for
those instances $i\in\mathcal{I}_{D}$ to which $a$ has been applied, so that a
performance evaluation $m(i,a)\in\mathbb{R}$ is available. Using this dataset,
a standard regression model $\widehat{m}_{a}$ can be learned per algorithm
$a$, and then used as a surrogate. The model can be realized as a neural
network or a random forest, and trained on loss functions such as root mean
squared or absolute error. For an overview of methods of this kind, we refer
to Section 6.

This approach has two main disadvantages. Firstly, it is not well suited for
the XAS setting, as it requires learning a huge number of surrogate models,
one per algorithm. Although these models can usually be trained very quickly,
the assumption of sparse training data in the XAS setting requires them to be
learned from only a handful of training examples — it is not even uncommon to
have algorithms without any performance value at all. Accordingly, the sparser
the data, the more drastically this approach drops in performance, as will be
seen in the evaluation in Section 5. Secondly, it requires precise real-valued
evaluations of the measure $m$ as training information, which might be costly
to obtain. In this regard, one may also wonder, whether regression is not
solving an unnecessarily difficult problem: Eventually, AS is only interested
in finding the best algorithm for a given problem instance, or, more
generally, in ranking the candidate algorithms in decreasing order of their
expected performance. An accurate prediction of absolute performances is a
_sufficient_ but not a _necessary_ condition for doing so.

###  3.2 Ranking

As an alternative to regression, one may therefore think of tackling AS as a
_ranking_ problem. More specifically, the counterpart of the regression
approach outlined above is called _label ranking_ (LR) in the literature [28].
Label ranking deals with learning to rank choice alternatives (referred to as
“labels”) based on given contexts represented by feature information. In the
setting of AS, contexts and labels correspond to instances and algorithms,
respectively. The type of training data assumed in LR consists of rankings
$\pi_{i}$ associated with training instances $i\in\mathcal{I}_{D}$, that is,
order relations of the form
$(f_{I}(i),a_{i,1})\succ\ldots\succ(f_{I}(i),a_{i,l_{i}})$, in which $\succ$
denotes an underlying preference relation; thus,
$(f_{I}(i),a)\succ(f_{I}(i),a^{\prime})$ means that, for instance $i$
represented by features $f_{I}(i)$, algorithm $a$ is preferred to (better
than) algorithm $a^{\prime}$. If $i$ is clear from the context, we also
represent the ranking by $a_{1}\succ\ldots\succ a_{l_{i}}$. Compared to the
case of regression, a ranking dataset of this form can be constructed more
easily, as it only requires qualitative comparisons between algorithms instead
of real-valued performance estimates.

A common approach to label ranking is based on the so-called Plackett-Luce
(PL) model [4], which specifies a parameterized probability distribution on
rankings over labels (i.e., algorithms in our case). The underlying idea is to
associate each algorithm $a$ with a latent utility function
$\widehat{m}_{a}:\mathcal{I}\longrightarrow\mathbb{R}_{+}$ of a context (i.e.,
an instance), which estimates how well an algorithm is suited for a given
instance. The functions $\widehat{m}_{a}$ are usually modeled as log-linear
functions

$\widehat{m}_{a}(i)=\exp\big{(}\boldsymbol{\theta}_{a}^{\top}\,f_{I}(i)\big{)}\,,$
(3)

where $\boldsymbol{\theta}_{a}\in\mathbb{R}^{k}$ is a real-valued,
$k$-dimensional vector, which has to be fit for each algorithm $a$. The PL
model specifies a probability distribution on rankings: given an instance
$i\in\mathcal{I}$, the probability of a ranking $a_{1}\succ\ldots\succ a_{z}$
over any subset $\\{a_{1},\ldots,a_{z}\\}\subseteq\mathcal{A}$ is

$\mathbb{P}(a_{1}\succ\ldots\succ
a_{z}\,|\,\boldsymbol{\Theta})=\prod_{n=1}^{z}\frac{\widehat{m}_{a_{n}}(i)}{\widehat{m}_{a_{n}}(i)+\ldots+\widehat{m}_{a_{z}}(i)}\,.$
(4)

A probabilistic model of that kind suggests learning the parameter matrix
$\boldsymbol{\Theta}=\\{\boldsymbol{\theta}_{a}\,|\,a\in\mathcal{A}\\}$ via
maximum likelihood estimation, i.e., by maximizing

$L(\boldsymbol{\Theta})=\prod_{i\in\mathcal{I}_{D}}\mathbb{P}(\pi_{i}\,|\,\boldsymbol{\Theta})$

associated with (4); this approach is explained in detail in [4]. Hence, the
associated loss function under which we learn is now of a probabilistic nature
(the logarithm of the PL-probability). It no longer focuses on the difference
between the approximated performance $\widehat{m}_{a}(i)$ and the true
performance $m(i,a)$, but on the ranking of the algorithms with respect to $m$
— putting it in the jargon of preference learning, the former is a “pointwise”
while the latter is a “listwise” method for learning to rank [3].

This approach potentially overcomes the second problem explained for the case
of regression, but not the first one: It still fits a single model per
algorithm $a$ (the parameter vector $\boldsymbol{\theta}_{a}$), which
essentially disqualifies it for the XAS setting.

###  3.3 Collaborative Filtering

This may suggest yet another approach, namely the use of collaborative
filtering (CF) [8], in the setting of AS originally proposed by [23]. In CF
for AS, we assume a (usually sparse) performance matrix
$R^{|\mathcal{I}_{D}|\times|\mathcal{A}|}$, where an entry $R_{i,a}=m(i,a)$
corresponds to the performance of algorithm $a$ on instance $i$ according to
$m$ if known, and $R_{i,a}=\,?$ otherwise. CF methods were originally designed
for large-scale settings, where products (e.g. movies) are recommended to
users, and data to learn from is sparse. Hence, they appear to fit well for
our XAS setting.

Similar to regression and ranking, model-based CF methods also learn a latent
utility function. They do so by applying matrix factorization techniques to
the performance matrix $R$, trying to decompose it into matrices
$U\in\mathbb{R}^{|\mathcal{I}_{D}|\times t}$ and
$V\in\mathbb{R}^{t\times|\mathcal{A}|}$ w.r.t. some loss function $L(R,U,V)$,
such that

$R\approx\widehat{R}=UV^{\top}\,,$ (5)

where $U$ ($V$) can be interpreted as latent features of the instances
(algorithms), and $t$ is the number of latent features. Accordingly, the
latent utility of a known algorithm $a$ for a known instance $i$ can be
computed as

$\widehat{m}_{a}(i)=U_{i,\bullet}V_{\bullet,a}^{\top}\,,$ (6)

even if the associated value $R_{i,a}$ is unknown in the performance matrix
used for training. The loss function $L(R,U,V)$ depends on the exact approach
used — examples include the root mean squared error and the absolute error
restricted by some regularization term to avoid overfitting. In [15], the
authors suggest a CF approach called Alors, which we will use in our
experiments later on. It can deal with unknown instances by learning a feature
map from the original instance to the latent instance feature space. Alors
leverages the CF approach CoFi${}^{\text{RANK}}$ [31] using the normalized
discounted cumulative gain (NDCG) [30] as loss function $L(R,U,V)$. Since the
NDCG is a ranking loss, it focuses on decomposing the matrix $R$ so as to
produce an accurate ranking of the algorithms. More precisely, it uses an
exponentially decaying weight function for ranks, such that more emphasis is
put on the top and less on the bottom ranks. Hence, it seems particularly well
suited for our use case.

##  4 Dyadic Feature Representation

As discussed earlier, by leveraging instance features, or learning such a
representation as in the case of Alors, the approaches presented in the
previous section can generalize over instances. Yet, none of them scales well
to the XAS setting, as they do not generalize over algorithms; instead, the
models are algorithm-specific and trained independently of each other. For the
approaches presented earlier (except for Alors), this does not only result in
a large number of models but also requires these models to be trained on very
few data. Furthermore, it is not uncommon to have algorithms without any
observation. A natural idea, therefore, is to leverage feature information on
algorithms as well.

More specifically, we use a feature function
$f_{A}:\mathcal{A}\longrightarrow\mathbb{R}^{d}$ representing algorithms as
$d$-dimensional, real-valued feature vectors. Then, instead of learning one
latent utility model per algorithm, the joint feature representation of a
“dyad” consisting of an instance and an algorithm, allows us to learn a single
joint model

$\widehat{m}:f_{I}(\mathcal{I})\times
f_{A}(\mathcal{A})\longrightarrow\mathbb{R}\,,$ (7)

and hence to estimate the performance of a given algorithm $a$ on a given
instance $i$ in terms of $\widehat{m}(f_{I}(i),f_{A}(a))$.

###  4.1 Regression

With the additional feature information at hand, instead of constructing one
dataset per algorithm, we resolve to a single joint dataset comprised of
examples $\Big{(}\psi\big{(}f_{I}(i),f_{A}(a)\big{)},m(i,a)\Big{)}$ with
dyadic feature information for all instances $i\in\mathcal{I}_{D}$ and
algorithms $a\in\mathcal{A}$ for which a performance value $m(i,a)$ is known.
Here,

$\psi:\mathbb{R}^{k}\times\mathbb{R}^{d}\longrightarrow\mathbb{R}^{q}$ (8)

is a joint feature map that defines how the instance and algorithm features
are combined into a single feature representation of a dyad. What is sought,
then, is a (parametrized) latent utility function
$\widehat{m}_{\boldsymbol{\theta}}:\mathbb{R}^{q}\longrightarrow\mathbb{R}$,
such that

$\widehat{m}_{\boldsymbol{\theta}}\Big{(}\psi\big{(}f_{I}(i),f_{A}(a)\big{)}\Big{)}\,$
(9)

is an estimation of the performance of algorithm $a$ on instance $i$.
Obviously, the choice of $\psi$ will have an important influence on the
difficulty of the regression problem and the quality of the model (9) induced
from the data $\mathcal{D}^{\mathit{REG}}$. The regression task itself comes
down to learning the parameter vector $\boldsymbol{\theta}$. In principle,
this can be done exactly as in Section 3.1, also using the same loss function.
Note that this is a generalization of the approach used by SMAC [11] for
predicting performances across instances in algorithm configuration. We allow
for a generic joint feature map $\psi$ and an arbitrary model for
$\widehat{m}_{\boldsymbol{\theta}}$, whereas SMAC limits itself to a
concatenation of features and trains a random forest for modeling
$\widehat{m}_{\boldsymbol{\theta}}$. Once again, it is noteworthy that SMAC by
itself is not applicable in the XAS setting, as it relies on costly online
evaluations.

###  4.2 Ranking

A similar adaptation can be made for the (label) ranking approach presented in
Section 3.2 [25]. Formally, this corresponds to a transition from the setting
of label ranking to the setting of _dyad ranking_ (DR) as recently proposed in
[21]. The first major change in comparison to the setting of label ranking
concerns the training data, where the rankings $\pi_{i}$ over subsets of
algorithms $\\{a_{i,1},\ldots,a_{i,l_{i}}\\}\subseteq\mathcal{A}$ for instance
$i$ are now of the form

$\psi\big{(}f_{I}(i),f_{A}(a_{i,1})\big{)}\succ\ldots\succ\psi\big{(}f_{I}(i),f_{A}(a_{i,l_{i}})\big{)}\,.$
(10)

Thus, we no longer represent an algorithm $a$ simply by its label ($a$) but by
features $f_{A}(a)$. Furthermore, like in the case of regression, we no longer
learn one latent utility function per algorithm, but a single model of the
form (9) based on a dyadic feature representation. In particular, we model
$\widehat{m}_{\boldsymbol{\theta}}$ as a feed-forward neural network, where
$\boldsymbol{\theta}$ represents its weights, which, as shown in [21], can be
learned via maximum likelihood estimation on the likelihood function implied
by the underlying PL model. Note that the use of a neural network is of
particular interest here, since it allows one to learn the underlying joint
feature map $\phi$ implicitly. Although both instance and algorithm features
are simply fed as a concatenated vector into the network, it can recombine
these features due to its structure and thus implicitly learn such a joint
feature representation.

In contrast to the methods presented in the previous section, the methods
based on dyadic feature information are capable of assigning a utility to
unknown algorithms. Thus, they are well suited for the XAS setting and in
principle even applicable when $\mathcal{A}$ is infinite, as long as a
suitable feature representation $f_{A}$ is available. Furthermore, as
demonstrated empirically in Section 5, the dyadic feature approaches are very
well suited for dealing with sparse performance matrices that are typical of
the XAS setting.

##  5 Experimental Evaluation

In our experiments, we evaluate well established state-of-the-art approaches
to algorithm selection as well as the proposed dyadic approaches in the XAS
setting. More specifically, we consider the problem of selecting a machine
learning classifier (algorithm) for a new classification dataset (instance) as
a case study related to the “on-the-fly machine learning” scenario [17].
Please note that this is just one amongst many conceivable instantiations of
the XAS setting, which is supposed to demonstrate the performance of the
presented methods. To this end, we first generate a benchmark and then use
this benchmark for comparison. The generated benchmark dataset as well as the
implementation of the approaches including detailed documentation is provided
on GitHub111https://github.com/alexandertornede/extreme_algorithm_selection.

###  5.1 Benchmark Dataset

In order to benchmark the generalization performance of the approaches
presented above in the XAS setting, we consider the domain of machine
learning. More precisely, the task is to select a classification algorithm for
an (unseen) dataset. Therefore, a finite set of algorithms $\mathcal{A}$ for
classification and a set of instances $\mathcal{I}$ corresponding to
classification datasets need to be specified. Furthermore, a performance
measure $m$ is needed to score the algorithms’ performance.

The set of candidate algorithms $\mathcal{A}$ is defined by sampling up to 100
different parameterizations of 18 classification algorithms from the machine
learning library WEKA [7], ensuring these parameterizations not being too
similar. An overview of the algorithms, their parameters and the number of
instantiations contained in $\mathcal{A}$ is given in Table 2. This yields
$|\mathcal{A}|=1,270$ algorithms in total. The last row of the table sums up
the items of the respective column, providing insights into the dimensionality
of the space of potential candidate algorithms.

The set of instances $\mathcal{I}$ is taken from the OpenML CC-18 benchmarking
suite222https://docs.openml.org/benchmark/#openml-cc18 (Excluding datasets
554, 40923, 40927, 40996 due to technical issues.) [27], which is a curated
collection of various classification datasets that are considered interesting
from a model selection resp. hyperparameter optimization point of view. This
property makes the datasets particularly appealing for the XAS benchmark
dataset, as it ensures more diversity across the algorithms.

Accordingly, the total performance matrix spanned by the algorithms and
classification datasets in principle features $1,270\cdot 71=88,900$ entries
for which the benchmark contains $55,919$ actual values and the rest are
unknown.

In the domain of machine learning, one is usually more interested in the
generalization performance of an algorithm than in the runtime. Therefore, $m$
is chosen to assess the solution performance of an algorithm. To this end, we
carry out a 5-fold cross validation and measure the mean accuracy across the
folds333The standard deviation of the performance values per dataset is on
average $0.101$, minimum $0.0064$ and maximum $0.33$.. As the measure of
interest, accuracy is a reasonable though to some extent arbitrary choice.
Note that in principle any other measure could have been used for generating
the benchmark as well.

Table 2: The table shows the types of classifiers used to derive the set
$\mathcal{A}$. Additionally, the number of numerical parameters (#num.P),
categorical parameters (#cat.P), and instantiations (n) is shown.

Learner |  0R |  1R |  BN |  DS |  DT |  IBk |  J48 |  JR |  KS |  L |  LMT |  MP |  NB |  PART |  REPT |  RF |  RT |  SMO
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
#num.P | 0 | 1 | 0 | 0 | 1 | 1 | 2 | 2 | 1 | 1 | 2 | 2 | 0 | 2 | 3 | 3 | 4 | 1
#cat.P | 0 | 0 | 2 | 0 | 3 | 3 | 6 | 2 | 2 | 0 | 5 | 6 | 2 | 2 | 2 | 2 | 4 | 2
n | 1 | 30 | 12 | 1 | 45 | 89 | 100 | 100 | 99 | 100 | 100 | 100 | 3 | 91 | 100 | 99 | 100 | 100

Training data for CF and regression-based approaches can then be obtained by
using the performance values as labels. In contrast, for training ranking
approaches, the data is labeled with rankings derived by ordering the
algorithms in a descending order w.r.t. their performance values. Note that
information about the exact performance value itself is lost in ranking
approaches.

We would like to note that the problem underlying this benchmark dataset could
of course be cast as an AC or CASH problem. However, here we make the
assumption that there is no time for costly online evaluations due to the on-
the-fly setting and hence standard AC and CASH methods are not applicable.

Instance Features. For the setting of machine learning, the instances are
classification datasets and associated feature representations are called
meta-features [18]. To derive a feature description of the datasets, we make
use of a specific subclass of meta-features called landmarkers, which are
performance scores of cheap-to-validate algorithms on the respective dataset.
More specifically, we use all 45 landmarkers as provided by OpenML [27], for
which different configurations of the following learning algorithms are
evaluated based on the error rate, area under the (ROC) curve, and Kappa
coefficient: Naive Bayes, One-Nearest Neighbour, Decision Stump, Random Tree,
REPTree and J48. Hence, in total we use 45 features to represent a
classification dataset.

Algorithm Features. The presumably most straight-forward way of representing
an algorithm in terms of a feature vector is to use the values of its
hyperparameters. Thus, we can describe each individual algorithm by a vector
of their hyperparameter-values. Based on this, the general feature description
is obtained by concatenation of the vectors. As already mentioned, the neural
network-based dyad ranking approach implicitly learns a more sophisticated
joint feature map. Due to the way in which we generated the set of candidate
algorithms $\mathcal{A}$, we can compress the vector sharing features for
algorithms of the same type. Additionally, we augment the vector by a single
categorical feature denoting the type of algorithm. Given any candidate
algorithm, its feature representation is obtained by setting the type of
algorithm indicator feature to its type, each element of the vector
corresponding to one of its hyperparameters to the specific value, and other
entries to $0$. Categorical parameteters, i.e. features, are one-hot encoded
yielding a total of 153 features to represent an algorithm.

###  5.2 Baselines

To better relate the performance of the different approaches to each other and
to the problem itself, we employ various baselines. While RandomRank assigns
ranks to algorithms simply at random, AvgPerformance first averages the
observed performance values for each candidate algorithm and predicts the
ranking according to these average performances. $k$-NN LR retrieves the $k$
nearest neighbors from the training data, averages the performances and
predicts the ranking which is induced by the average performances. Since
AvgRank is commonly used as another baseline in the standard AS setting, we
note that we omit this baseline on purpose. This is because meaningful average
ranks of algorithms are difficult to compute in the XAS setting, where the
number of algorithms evaluated, and hence the length of the rankings of
algorithms, vary from dataset to dataset.

###  5.3 Experimental Setup

In the following experiments, we investigate the performance of the different
approaches and baselines in the setting of XAS for the example of the proposed
benchmark dataset as described in Section 5.1.

We conduct a 10-fold cross validation to divide the dataset into 9 folds of
known and 1 fold of unknown instances. From the resulting set of known
performance values, we then draw a sample of 25, 50, or 125 pairs of
algorithms for every instance under the constraint that the performances of
the two algorithms is not identical. Thus, a maximum fill degree of 4%, 8%
respectively 20% of the performance matrix is used for training, as algorithms
may occur more than once in the sampled pairs. The sparse number of training
examples is motivated by the large number of algorithms in the XAS setting.
The assumption that performance values are only available for a small subset
of the algorithms is clearly plausible here. Throughout the experiments, we
ensure that all approaches are provided the same instances for training and
testing, and that the label information is at least based on the same
performance values.

In the experiments, we compare various models with each other. This includes
two versions of Alors, namely Alors (REGR) and Alors (NDCG) optimizing for a
regression respectively ranking loss. Furthermore, we consider a state-of-the-
art regression approach learning a RandomForest regression model per algorithm
(PAReg). Note that for those algorithms with no training data at all, we make
PAReg predict a performance of $0$, as recommending such an algorithm does not
seem reasonable. Lastly, we consider two approaches leveraging a dyadic
feature representation, internally fitting either a RandomForest for
regression (DFReg) or a feed-forward neural network for ranking (DR). For both
dyadic approaches, the simple concatenation of instance and algorithm features
is used as a feature map. In contrast to the other methods, the ranking model
is only provided the information which algorithm of a sampled pair performs
better, as opposed to the exact performance value that is given to other
methods. A summary of the type of features and label information used by the
different approaches/baselines is given on the left side of Table 3.

Table 3: Overview of the data provided to the approaches and their
applicability to the considered scenarios.

|  Approach | $f_{I}$ | $f_{A}$ | Label |  |  Approach | $f_{I}$ | $f_{A}$ | Label
---|---|---|---|---|---|---|---|---|---
approaches |  Alors (REGR) | ✓ | ✗ | $m$ |  baselines |  RandomRank | ✗ | ✗ |
Alors (NDCG) | ✓ | ✗ | $m$ |  AvgPrfm | ✗ | ✗ | $m$
PAReg | ✓ | ✗ | $m$ |  AvgRank | ✗ | ✗ | $\pi$
DFReg | ✓ | ✓ | $m$ |  $k$-NN LR | ✓ | ✗ | $m$
DR | ✓ | ✓ | $\pi$ |  |  |  |  |

The test performance of the approaches is evaluated by sampling $10$
algorithms for every (unknown) instance to test for. The comparison is done
with respect to different metrics detailed further below, and the outlined
sampling evaluation routine is repeated $100$ times.

Statistical significance w.r.t performance differences between the best method
and any other method is determined by a Wilcoxon rank sum test with a
threshold of 0.05 for the p-value. Significant improvements of the best method
over another one is indicated by $\bullet$.

Experiments were run on nodes with two Intel Xeon Gold “Skylake” 6148 with 20
cores each and 192 GB RAM.

###  5.4 Performance Metrics

On the test data, we compute the following performance metrics measuring
desirable properties of XAS approaches.

regret@$k$ is the difference between the performance value of the best
algorithm within the predicted top-$k$ of algorithms and the actual best
algorithm. The domain of regret@$k$ is $[0,1]$, where 0 is the optimum meaning
no regret.

NDCG$@k$ is a position-dependent ranking measure (normalized discounted
cumulative gain) to measure how well the ranking of the top-$k$ algorithms can
be predicted. It is defined as

$\text{NDCG}@k(\pi,\pi^{*})=\frac{\text{DCG}@k(\pi)}{\text{DCG}@k(\pi^{*})}=\dfrac{\left(\sum\limits_{n=1}^{k}\frac{2^{m(i,\pi_{n})-1}}{\log(n+2)}\right)}{\left(\sum\limits_{n=1}^{k}\frac{2^{m(i,\pi^{*}_{n})-1}}{\log(n+2)}\right)}\,,$

where $i$ is a (fixed) instance, $\pi$ is a ranking and $\pi^{*}$ the optimal
ranking, and $\pi_{n}$ gives the algorithm on rank $n$ in ranking $\pi$. The
NDCG emphasizes correctly assigned ranks at higher positions with an
exponentially decaying importance. NDCG ranges in $[0,1]$, where $1$ is the
optimal value.

Kendall’s $\tau$ is a rank correlation measure. Given two rankings (over the
same set of elements) $\pi$ and $\pi^{\prime}$, it is defined as

$\tau(\pi,\pi^{\prime})=\frac{C-D}{\sqrt{(C+D+T_{\pi})\cdot(C+D+T_{\pi^{\prime}})}}$
(11)

where $C$/$D$ is the number of so-called concordant/discordant pairs in the
two rankings, and $T_{\pi}$/$T_{\pi^{\prime}}$ is the number of ties in
$\pi$/$\pi^{\prime}$. Two elements are called a concordant/discordant pair if
their order within the two rankings is identical/different, and tied if they
are on the same rank. Intuitively, this measure determines on how many pairs
the two rankings coincide. It takes values in $[-1,1]$, where $0$ means
uncorrelated, $-1$ inversely, and $1$ perfectly correlated.

###  5.5 Results

The results of the experiments are shown in Table 4. It is clear from the
table that the methods for standard algorithm selection tend to fail
especially in the scenarios with only few algorithm performance values per
instance. This includes the approach of building a distinct regression model
for each algorithm (PAReg) as well as for the collaborative filtering approach
Alors, independently of the loss optimized for, even though the NDCG variant
has a slight edge over the regression one. Moreover, Alors even fails to
improve over simple baselines, such as AvgPerformance and $k$-NN LR. With an
increasing number of training examples, PAReg improves over the baselines and
also performs better than Alors, but never yields the best performance for any
of the considered settings or metrics.

Table 4: Results for the performance metrics Kendall’ tau ($\tau$), NDCG@k
(N@3, N@5), and regret@$k$ (R@1, R@3) for varying number of performance value
pairs used for training. The best performing approach is highlighted in bold,
the second best is underlined, and significant improvements of the best
approach over others is denoted by $\bullet$.

Approach | 4% fill rate / 25 performance value pairs | 8% fill rate / 50 performance value pairs
---|---|---
$\tau$ | N@3 | N@5 | R@1 | R@3 | $\tau$ | N@3 | N@5 | R@1 | R@3
PAReg | 0.1712 $\bullet$ | 0.9352 $\bullet$ | 0.9433 $\bullet$ | 0.0601 $\bullet$ | 0.0185 $\bullet$ | 0.2537 $\bullet$ | 0.9453 | 0.9594 | 0.0493 | 0.0136
Alors (NDCG) | 0.0504 $\bullet$ | 0.9205 $\bullet$ | 0.9223 $\bullet$ | 0.0686 $\bullet$ | 0.0225 | 0.0472 $\bullet$ | 0.9155 $\bullet$ | 0.9164 $\bullet$ | 0.0614 $\bullet$ | 0.0208
Alors (REGR) | 0.0303 $\bullet$ | 0.9117 $\bullet$ | 0.9191 $\bullet$ | 0.0794 $\bullet$ | 0.0190 $\bullet$ | 0.0807 $\bullet$ | 0.9172 $\bullet$ | 0.9304 $\bullet$ | 0.0754 $\bullet$ | 0.0285 $\bullet$
DR |  0.3445 |  0.9523 |  0.9604 | 0.0381 |  0.0089 |  0.3950 |  0.9584 |  0.9685 |  0.0322 |  0.0087
DFReg |  0.3819 |  0.9564 |  0.9652 |  0.0302 |  0.0079 |  0.3692 |  0.9573 |  0.9661 |  0.0300 |  0.0123
RandomRank | -0.0038 $\bullet$ | 0.8933 $\bullet$ | 0.9105 $\bullet$ | 0.0878 $\bullet$ | 0.0272 $\bullet$ | -0.0038 $\bullet$ | 0.8933 $\bullet$ | 0.9105 $\bullet$ | 0.0878 $\bullet$ | 0.0272 $\bullet$
AvgPerformance | 0.1384 $\bullet$ | 0.9388 $\bullet$ | 0.9433 $\bullet$ |  0.0337 | 0.0090 | 0.2083 $\bullet$ | 0.9355 $\bullet$ | 0.9508 $\bullet$ | 0.0493 $\bullet$ | 0.0199 $\bullet$
1-NN LR | 0.1227 $\bullet$ | 0.9290 $\bullet$ | 0.9310 $\bullet$ | 0.0733 $\bullet$ | 0.0230 $\bullet$ | 0.1059 $\bullet$ | 0.9246 $\bullet$ | 0.9296 $\bullet$ | 0.0564 $\bullet$ | 0.0209
2-NN LR | 0.1303 $\bullet$ | 0.9278 $\bullet$ | 0.9310 $\bullet$ | 0.0642 $\bullet$ | 0.0193 $\bullet$ | 0.0874 $\bullet$ | 0.9269 $\bullet$ | 0.9343 $\bullet$ | 0.0541 $\bullet$ | 0.0206
Approach | 20% fill rate / 125 performance value pairs |  |  |  |  |
$\tau$ | N@3 | N@5 | R@1 | R@3 |  |  |  |  |
PAReg | 0.3003 $\bullet$ | 0.9525 | 0.9632 | 0.0395 | 0.0107 |  |  |  |  |
Alors (NDCG) | 0.0540 $\bullet$ | 0.9220 $\bullet$ | 0.9242 $\bullet$ | 0.0542 $\bullet$ | 0.0228 $\bullet$ |  |  |  |  |
Alors (REGR) | 0.1039 $\bullet$ | 0.9160 $\bullet$ | 0.9329 $\bullet$ | 0.0604 $\bullet$ | 0.0222 $\bullet$ |  |  |  |  |
DR |  0.4507 |  0.9696 |  0.9715 |  0.0241 |  0.0055 |  |  |  |  |
DFReg |  0.4264 |  0.9629 |  0.9720 |  0.0292 |  0.0071 |  |  |  |  |
RandomRank | -0.0038 $\bullet$ | 0.8933 $\bullet$ | 0.9105 $\bullet$ | 0.0878 $\bullet$ | 0.0272 $\bullet$ |  |  |  |  |
AvgPerformance | 0.2541 $\bullet$ | 0.9437 $\bullet$ | 0.9536 $\bullet$ | 0.0523 $\bullet$ | 0.0084 |  |  |  |  |
1-NN LR | 0.1152 $\bullet$ | 0.9245 $\bullet$ | 0.9318 $\bullet$ | 0.0594 $\bullet$ | 0.0249 $\bullet$ |  |  |  |  |
2-NN LR | 0.1142 $\bullet$ | 0.9292 $\bullet$ | 0.9350 $\bullet$ | 0.0412 | 0.0176 $\bullet$ |  |  |  |  |

In contrast to this, the proposed dyadic feature approaches clearly improve
over both the methods for the standard AS setting and the considered baselines
for all the metrics. Interestingly, DFReg performs best for the setting with
only 25 performance value pairs, while DR has an edge over DFReg for the other
two settings. Still, the differences between the dyadic feature approaches are
never significant, whereas significant improvements can be achieved in
comparison to the baselines and the other AS approaches.

Moreover, our study demonstrates the heterogeneity of the benchmark dataset.
As described in [22], a relevant measure for heterogeneity is the per-instance
potential for improvement over a solution that is static across instances,
i.e., what is often called the single best algorithm or solver (SBS). In this
case study, the SBS is represented by the AvgPerformance baseline, which is
always worse than the oracle with respect to all measures and in particular
the regret@$k$ measures. Hence, as the superior performances of our approach
compared to the AvgPerformance demonstrate, the benchmark dataset offers a
potential for per-instance algorithm selection.

The results of our study show that models with strong generalization
performance can be obtained despite the small number of training examples.
Moreover, the results suggest that there is a need for the development of
specific methods addressing the characteristics of the XAS setting. This
concerns the large number of different candidate algorithms as well as the
sparsity of the training data.

##  6 Related Work

As most closely related work, we subsequently highlight several AS approaches
to learning latent utility functions. For an up-to-date survey, we refer to
[12].

A prominent example of a method learning a regression-based latent utility
function is [32], which features an empirical hardness model per algorithm for
estimating the runtime of an algorithm, i.e., its performance, for a given
instance based on a ridge regression approach in the setting of SAT solver
selection. Similarly, [13] learn per-algorithm hardness models using
statistical (non-)linear regression models for algorithms solving the winner
determination problem. Depending on whether a given SAT instance is presumably
satisfiable or not, conditional runtime prediction models are learned in [9]
using ridge linear regression.

In [5], a label-ranking-based AS approach for selecting collaborative
filtering algorithms in the context of recommender systems is presented
leveraging nearest neighbor and random forest label rankers.

Similar to our work, [19] leverages algorithm features in the form of a binary
vector indicating which algorithm is considered to learn a probabilistic
ranking model considering up to tens of algorithms. AS was first modeled as a
CF problem in [23], using a probabilistic matrix factorization technique to
select algorithms for the constraint solving problem. Assuming a complete
performance matrix for training, low-rank latent factors are learned in [14]
using singular value decomposition to obtain a selector en par with the
oracle. Lastly, in [26] a decision-theoretic approach is proposed leveraging
survival analysis to explicitly acknowledge timeouts of algorithms in the
learning process.

##  7 Conclusion

In this paper, we introduced the extreme algorithm selection (XAS) setting and
investigated the scalability of various algorithm selection approaches in this
setting. To this end, we defined a benchmark based on the OpenML CC-18
benchmark suite for classification and a set of more than 1,200 candidate
algorithms. Furthermore, we proposed the use of dyadic approaches,
specifically dyad ranking, taking into account feature representations of both
problem instances (datasets) and algorithms, which allows them to work on very
few training data. In an extensive evaluation, we found that the approaches
exploiting dyadic feature representations perform particularly well according
to various metrics on the proposed benchmark and outperform other state-of-
the-art AS approaches developed for the standard AS setting.

The currently employed algorithm features allow for solving the cold start
problem only to a limited extent, i.e., only algorithms featuring known
hyperparameters can be considered as new candidate algorithms. Investigating
features to describe completely new algorithms is a key requirement for the
approaches considered in this paper, and therefore an important direction for
future work.

## Acknowledgements

This work was supported by the German Research Foundation (DFG) within the
Collaborative Research Center “On-The-Fly Computing” (SFB 901/3 project no.
160364472) and by the Paderborn Center for Parallel Computing (PC2) through
computing time.

## References

  * [1] Bengio, S., Dembczynski, K., Joachims, T., Kloft, M., Varma, M.: Extreme classification (dagstuhl seminar 18291). Dagstuhl Reports 8(7), 62–80 (2018)
  * [2] Bischl, B., Kerschke, P., Kotthoff, L., Lindauer, M., Malitsky, Y., Fréchette, A., Hoos, H.H., Hutter, F., Leyton-Brown, K., Tierney, K., et al.: Aslib: A benchmark library for algorithm selection. Artificial Intelligence 237 (2016)
  * [3] Cao, Z., Qin, T., Liu, T., Tsai, M., Li, H.: Learning to rank: from pairwise approach to listwise approach. In: ICML (2007)
  * [4] Cheng, W., Hüllermeier, E., Dembczynski, K.J.: Label ranking methods based on the plackett-luce model. In: ICML (2010)
  * [5] Cunha, T., Soares, C., de Carvalho, A.C.: A label ranking approach for selecting rankings of collaborative filtering algorithms. In: SAC (2018)
  * [6] Feurer, M., Klein, A., Eggensperger, K., Springenberg, J., Blum, M., Hutter, F.: Efficient and robust automated machine learning. In: NIPS (2015)
  * [7] Frank, E., Hall, M., Witten, I.: The WEKA Workbench. Data Mining (2016)
  * [8] Goldberg, D., Nichols, D., Oki, B.M., Terry, D.: Using collaborative filtering to weave an information tapestry. Communications of the ACM 35(12) (1992)
  * [9] Haim, S., Walsh, T.: Restart strategy selection using machine learning techniques. In: SAT (2009)
  * [10] Happe, M., Meyer auf der Heide, F., Kling, P., Platzner, M., Plessl, C.: On-the-fly computing: A novel paradigm for individualized it services. In: ISORC (2013)
  * [11] Hutter, F., Hoos, H.H., Leyton-Brown, K.: Sequential model-based optimization for general algorithm configuration. In: LION (2011)
  * [12] Kerschke, P., Hoos, H.H., Neumann, F., Trautmann, H.: Automated algorithm selection: Survey and perspectives. ECJ 27(1) (2019)
  * [13] Leyton-Brown, K., Nudelman, E., Shoham, Y.: Learning the empirical hardness of optimization problems: The case of combinatorial auctions. In: CP (2002)
  * [14] Malitsky, Y., O’Sullivan, B.: Latent features for algorithm selection. In: SOCS (2014)
  * [15] Mısır, M., Sebag, M.: Alors: An algorithm recommender system. Artificial Intelligence 244 (2017)
  * [16] Mohr, F., Wever, M., Hüllermeier, E.: ML-Plan: Automated machine learning via hierarchical planning. Machine Learning 107(8-10) (2018)
  * [17] Mohr, F., Wever, M., Tornede, A., Hüllermeier, E.: From automated to on-the-fly machine learning. In: INFORMATIK (2019)
  * [18] Nguyen, P., Hilario, M., Kalousis, A.: Using meta-mining to support data mining workflow planning and optimization. JAIR 51 (2014)
  * [19] Oentaryo, R.J., Handoko, S.D., Lau, H.C.: Algorithm selection via ranking. In: AAAI (2015)
  * [20] Rice, J.R.: The algorithm selection problem. In: Advances in computers, vol. 15. Elsevier (1976)
  * [21] Schäfer, D., Hüllermeier, E.: Dyad ranking using plackett-luce models based on joint feature representations. Machine Learning 107(5) (2018)
  * [22] Schneider, M., Hoos, H.H.: Quantifying homogeneity of instance sets for algorithm configuration. In: LION (2012)
  * [23] Stern, D.H., Samulowitz, H., Herbrich, R., Graepel, T., Pulina, L., Tacchella, A.: Collaborative expert portfolio management. In: AAAI (2010)
  * [24] Thornton, C., Hutter, F., Hoos, H.H., Leyton-Brown, K.: Auto-weka: Combined selection and hyperparameter optimization of classification algorithms. In: SIGKDD (2013)
  * [25] Tornede, A., Wever, M., Hüllermeier, E.: Algorithm selection as recommendation: From collaborative filtering to dyad ranking. In: CI Workshop, Dortmund (2019)
  * [26] Tornede, A., Wever, M., Werner, S., Mohr, F., Hüllermeier, E.: Run2Survive: A decision-theoretic approach to algorithm selection based on survival analysis. In: ACML (2020)
  * [27] Vanschoren, J., van Rijn, J.N., Bischl, B., Torgo, L.: Openml: Networked science in machine learning. SIGKDD Explorations 15(2) (2013)
  * [28] Vembu, S., Gärtner, T.: Label ranking algorithms: A survey. In: Preference Learning. Springer (2010)
  * [29] Wang, W., Zheng, V.W., Yu, H., Miao, C.: A survey of zero-shot learning: Settings, methods, and applications. ACM TIST (2019)
  * [30] Wang, Y., Wang, L., Li, Y., He, D., Chen, W., Liu, T.: A theoretical analysis of ndcg ranking measures. In: COLT (2013)
  * [31] Weimer, M., Karatzoglou, A., Le, Q.V., Smola, A.J.: Cofi rank-maximum margin matrix factorization for collaborative ranking. In: NIPS (2008)
  * [32] Xu, L., Hutter, F., Hoos, H.H., Leyton-Brown, K.: Satzilla: portfolio-based algorithm selection for sat. JAIR 32 (2008)


# Supplementary Information
Ferroelectric Semimetals with $\alpha$-Bi/SnSe van der Waals heterostructures
and its Topological Currents

D. J. P. de Sousa  Department of Electrical and Computer Engineering,
University of Minnesota, Minneapolis, Minnesota 55455, USA  Seungjun Lee
Department of Electrical and Computer Engineering, University of Minnesota,
Minneapolis, Minnesota 55455, USA  Qiangsheng Lu  Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee
37831, USA  Rob G. Moore  Materials Science and Technology Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA  Matthew Brahlek
Materials Science and Technology Division, Oak Ridge National Laboratory, Oak
Ridge, Tennessee 37831, USA  J-.P. Wang  Department of Electrical and Computer
Engineering, University of Minnesota, Minneapolis, Minnesota 55455, USA
Department of Physics, University of Minnesota, Minneapolis, Minnesota 55455,
USA  Guang Bian<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Missouri, Columbia, Missouri 65211, USA  Tony Low<EMAIL_ADDRESS>Department of Electrical and Computer Engineering, University of Minnesota,
Minneapolis, Minnesota 55455, USA  Department of Physics, University of
Minnesota, Minneapolis, Minnesota 55455, USA

## Symmetries

Let $D(g)$ be the representation for the symmetry operator $g$. Then, $g$ is a
symmetry of the system if $D(g)H(\textbf{k})D^{-1}(g)=H(g\textbf{k})$. The
symmetry action on the basis is

$\displaystyle g\hat{c}_{\textbf{k}\alpha
s}^{\dagger}g^{-1}=\displaystyle\sum_{s^{\prime}}\hat{c}_{g\textbf{k},\alpha^{\prime}s^{\prime}}^{\dagger}[D(g)]_{\alpha^{\prime}s^{\prime},\alpha
s}.$ (S1)

The symmetries of the low energy $p_{z}$ states of monolayer $\alpha$-Bi are
two-fold rotations in relation to the z and x axis, $C_{2z}$ and $C_{2x}$
respectively, inversion $\mathcal{P}$ and time-reversal $\mathcal{T}$. Their
representations are

$\displaystyle D(C_{2x})=i\sigma_{x}s_{x},$ $\displaystyle
D(C_{2z})=-i\sigma_{x}s_{z},$ $\displaystyle
D(\mathcal{P})=-i\sigma_{x}s_{0},$ $\displaystyle
D(\mathcal{T})=-i\sigma_{0}s_{y}\mathcal{K},$ (S2)

where $\mathcal{K}$ is the complex conjugation operator. The combined
$D(C_{2z}C_{2y})=-i\sigma_{0}s_{y}$ symmetry corresponds to a $y$ mirror plane
$\mathcal{M}_{y}$. Here, the $\sigma$ pauli matrices operator in orbital space
and $s$ in spin space. The basis is $\Psi=[\phi_{1}\ \ \phi_{2}]^{\dagger}$,
where $\phi_{j}=[\phi_{1}^{\uparrow}\ \ \phi_{1}^{\downarrow}]^{T}$, with
$\phi_{1}^{s}=\phi_{p_{z}}^{A,s}+\phi_{p_{z}}^{D,s}$ and
$\phi_{2}^{s}=\phi_{p_{z}}^{B,s}+\phi_{p_{z}}^{C,s}$, where
$\phi_{p_{z}}^{l,s}$ corresponds to the amplitude of finding electron in
sublattice $l=A,B,C,D$ of the monolayer $\alpha$-Bi unit cell, which are
identified in Fig. S1(b).

Figure S1: (a) Monolayer $\alpha$-Bi crystal structure. The sublattices of
$\alpha$-Bi monolayer are highlighted in the right panel. (b) The basis
consists of combinations of pairs of $p_{z}$ state amplitudes in sublattices A
and D, and sublattices B and C. This is due to their symmetric relations such
as in black-phosphorus monolayers . (c) and (d) illustrate the sublattice
symmetry breaking due to proximity with the local polarization of SnSe.

Proximity with the in-plane polarization of SnSe breaks the symmetry between
the sublattices A (D) and B (C). This is illustrated in Figs. (S1) (c) and (d)
for the two ferroelectric configurations. Hence, proximity with SnSe breaks
the two-fold rotation around the x axis, $C_{2x}$, and inversion,
$\mathcal{P}$.

## Insulator-to-Semimetal transition

We describe here the insulator-to-semimetal transition induced by proximity
effects. To simplify our analysis, we assume
$u_{0}=\eta_{x}=\eta_{y}=\lambda_{R}=0$, without loss of generality. Hence,
the continuum Hamiltonian becomes

$\displaystyle
H(\textbf{k})=h_{x}(\textbf{k})\sigma_{x}s_{0}+h_{y}(\textbf{k})\sigma_{y}s_{0}+h_{\textrm{SOC}}(\textbf{k})\sigma_{z}s_{z}+\Delta\sigma_{z}s_{0}$
(S3)

where $h_{x}(\textbf{k})=\delta_{0}+\gamma_{x}k_{x}^{2}+\gamma_{y}k_{y}^{2}$,
$h_{y}(\textbf{k})=\xi k_{x}$ and
$h_{\textrm{SOC}}(\textbf{k})=\lambda_{I}k_{y}$, with $\lambda_{I}$ being the
intrinsic spin-orbit coupling parameter that we assume to be positive,
$\lambda_{I}>0$.It is convenient to work in a basis containing symmetric and
antisymmetric superpositions of all sublattice amplitudes. This is achieved by
employing the unitary transformation $U=(\sigma_{x}+\sigma_{z})s_{0}/2$, which
mixes the sublattice amplitudes as: $\bar{\Psi}=U\Psi=[\bar{\phi}_{1}\ \
\bar{\phi}_{2}]^{\dagger}$, with
$\bar{\phi}_{1}=(\phi_{A}+\phi_{D}+\phi_{B}+\phi_{C})/2$ and
$\bar{\phi}_{2}=(\phi_{A}+\phi_{D}-\phi_{B}-\phi_{C})/2$. The transformed
Hamiltonian, $\bar{H}(\textbf{k})=U^{\dagger}H(\textbf{k})U$, becomes,

$\displaystyle\bar{H}(\textbf{k})=h_{x}(\textbf{k})\sigma_{z}s_{0}-h_{y}(\textbf{k})\sigma_{y}s_{0}+\sigma_{x}(h_{\textrm{SOC}}(\textbf{k})s_{z}+\Delta
s_{0}).$ (S4)

The gap closure takes place along the $k_{x}=0$ region of momentum space.
Hence, we solve for the eigenvalues of

$\displaystyle\bar{H}(k_{y})=\displaystyle\left(\begin{tabular}[]{cc}$(\delta_{0}+\gamma_{y}k_{y}^{2})s_{0}$&$\Delta
s_{0}+\lambda_{I}k_{y}s_{z}$\\\ $\Delta
s_{0}+\lambda_{I}k_{y}s_{z}$&$-(\delta_{0}+\gamma_{y}k_{y}^{2})s_{0}$\end{tabular}\right).$
(S7)

To this end, we explicitly write the Schrödinger equation,
$\bar{H}(k_{y})\bar{\Psi}=\epsilon(k_{y})\bar{\Psi}$, in terms of its
components

$\displaystyle\displaystyle(\delta_{0}+\gamma_{y}k_{y}^{2}-\epsilon)\bar{\phi}_{1}+(\Delta
s_{0}+\lambda_{I}k_{y}s_{z})\bar{\phi}_{2}=$ $\displaystyle 0,$
$\displaystyle\displaystyle(\delta_{0}+\gamma_{y}k_{y}^{2}+\epsilon)\bar{\phi}_{2}-(\Delta
s_{0}+\lambda_{I}k_{y}s_{z})\bar{\phi}_{1}=$ $\displaystyle 0.$ (S8)

The above pair of equations can be easily decoupled into

$\displaystyle\\{[(\delta_{0}+\gamma_{y}k_{y}^{2})^{2}-\epsilon^{2}]s_{0}+(\Delta
s_{0}+\lambda_{I}k_{y}s_{z})^{2}\\}\bar{\phi}_{1}=0,$ (S9)

which can be used to obtain the spin components of the symmetric combination
of the sublattice amplitudes. The energy eigenvalues corresponding to the two
distinct spin states are easily obtained by requiring the existence of non-
trivial solutions for the components of $\bar{\phi}_{1}$. The results is

$\displaystyle\epsilon_{ns}^{\eta}(k_{y})=n\sqrt{(\delta_{0}+\gamma_{y}k_{y}^{2})^{2}+\Delta^{2}+\lambda_{I}^{2}k_{y}^{2}+2s\lambda_{I}\Delta
k_{y}},$ (S10)

where $n=\pm 1$ and $s=\pm 1$ are the band and spin indexes, respectively.
Note that in the absence of spin-orbit coupling and sublattice symmetry
breaking, i.e., $\lambda_{I}=\Delta=0$, a gap closure takes place if
$\operatorname{sign}(\delta_{0}\gamma_{y})=-1$ at
$\Lambda_{y}^{\eta}=\eta\sqrt{-\delta_{0}/\gamma_{y}}$, where $\eta=\pm$ being
the valley index. The band crossings are protected by a mirror
symmetry…distinct symmetry eigenvalues etc.

Figure S2: Insulator-to-semimetal transition; Evolution of the spin-resolved
bands of the $p_{z}$ Hamiltonian of $\alpha$-Bi with the proximity-induced
sublattice symmetry breaking parameter $\Delta$. We considered the parameter
values $u_{0}=0$ eV, $\eta_{x}=0.1$ eV Å 2, $\eta_{y}=0.5$ eV Å 2,
$\delta_{0}=-0.05$ eV, $\gamma_{x}=0.13$ eV Å 2, $\gamma_{y}=1.925$ eV Å 2,
$\lambda_{I}=-0.2$ eV Å , $\lambda_{R}=0$ eV , with (a) $\Delta=0$ eV, (b)
$\Delta=-0.016$ eV, (c) $\Delta=-0.032$ eV $=\Delta_{c}$ and (d)
$\Delta=-0.048$ eV.

In the presence of intrinsic spin-orbit coupling and polarization-induced
sublattice symmetry breaking terms, i.e., $\lambda_{I}\neq 0,\Delta\neq 0$, it
is still possible to attain a gapless phase with band crossings at the valleys
$\Lambda_{y}^{\eta}$. In this scenario, the energies become

$\displaystyle\epsilon_{ns}^{\eta}(k_{y}=\eta\Lambda_{y})=n\sqrt{\Delta^{2}+\lambda_{I}^{2}\Lambda_{y}^{2}+2s\eta\lambda_{I}\Delta\Lambda_{y}},$
(S11)

where $\Lambda_{y}=\sqrt{-\delta_{0}/\gamma_{y}}$, with an energy gap given by

$\displaystyle E_{g}=2|\Delta+s\eta\lambda_{I}\Lambda_{y}|.$ (S12)

Note that a gap closure can be attained for a particular choice of $\Delta$.
We define the quantity $\Delta_{c}=\lambda_{I}\sqrt{-\delta_{0}/\gamma_{y}}$,
which will be very useful in the following. Hence, an insulator-to-semimetal
transition takes place with $\Delta$. This is shown in Fig. (S2), where we
have considered a lattice version of the symmetry derived Hamiltonian
presented in the main text. As $\Delta$ is increased from $0$ eV to a finite
value, the previously degenerated opposite spin bands split in momentum due to
inversion symmetry breaking. The spin degeneracy breaking is accompanied by a
decrease in the energy gap of opposite spin bands around $\pm\Lambda_{y}$,
which eventually vanishes at the critical value $\Delta_{c}$, as shown in Fig.
S2(c). An energy gap reopens for $\Delta>\Delta_{c}$.

Next, we examine two distinct situations: I) $\Delta>0$ and II) $\Delta<0$,
corresponding to the two distinct polarization states.

### (I) The $\Delta>0$ case

The gap closes, i.e., $E_{g}=0$, at the critical value $\Delta=\Delta_{c}$
whenever $\eta s=-1$. The first possibility, with $\eta=+1$ and $s=-1$,
corresponds to a gap closure of the spin down bands at valley located at
$k_{y}=+\Lambda_{y}$. The second possibility, with $\eta=-1$ and $s=+1$,
corresponds to a gap closure of the spin up bands at the valley located at
$k_{y}=-\Lambda_{y}$. Note that the spin down (up) bands remain gapped at
$k_{y}=-\Lambda_{y}$ ($k_{y}=+\Lambda_{y}$). The gaps are
$E_{g}=4\lambda_{I}\sqrt{-\delta_{0}/\gamma_{y}}$. Therefore, the system
becomes semimetallic with two fully spin-polarized valleys separated along the
$k_{y}$ direction. The opposite spin polarization of the two valleys is a
manifestation of time-reversal symmetry.

### (II) The $\Delta<0$ case

A gap closure is attained at the critical value $\Delta=-\Delta_{c}$ whenever
$\eta s=+1$. The first possibility, with $\eta=+1$ and $s=+1$, corresponds to
a gap closure of the spin up bands at valley located at $k_{y}=+\Lambda_{y}$.
The second possibility, with $\eta=-1$ and $s=-1$, corresponds to a gap
closure of the spin down bands at the valley located at $k_{y}=-\Lambda_{y}$.
The spin up (down) bands remain gapped at $k_{y}=-\Lambda_{y}$
($k_{y}=+\Lambda_{y}$), with
$E_{g}=4\lambda_{I}\sqrt{-\delta_{0}/\gamma_{y}}$. Hence, the system is a
semimetal at $\Delta=-\Delta_{c}$ with two fully spin-polarized valleys
separated along the $k_{y}$ direction.

Figure S3: Locking between valley, spin and the sign of $\Delta$. Upon
switching the ferroelectric polarization of SnSe layer, and consequently the
sign of $\Delta$, the spin flavor in a given valley is reversed. We considered
the parameter values $u_{0}=0$ eV, $\eta_{x}=0.1$ eV Å 2, $\eta_{y}=0.5$ eV Å
2, $\delta_{0}=-0.05$ eV, $\gamma_{x}=0.13$ eV Å 2, $\gamma_{y}=1.925$ eV Å 2,
$\lambda_{I}=-0.2$ eV Å , $\lambda_{R}=0$ eV , with (a) $\Delta=-0.032$ eV,
(b) $\Delta=0.032$ eV.

The two situations described above are shown in Fig. S3. We conclude that the
only effect of changing the sign of $\Delta$ when $|\Delta|=\Delta_{c}$ is to
change the spin polarizations of the two valleys. Hence, the valley spin
polarization, and consequently, the spin/Valley Hall effect are coupled to the
polarization states of the SnSe monolayer underneath $\alpha$-Bi.

## Effective pseudo-Weyl Hamiltonian

From the previous analysis, one should expect Weyl-like Hamiltonians
describing the electron states around each valley. Because the crossing bands
are fully spin polarized, we proceed by extracting the Hamiltonians for each
spin sector. To this end, we start by rewriting Hamiltonian (S3) in the “spin-
polarized” basis $\tilde{\Psi}=[\phi_{1}^{\uparrow}\ \ \phi_{2}^{\uparrow}\ \
\phi_{1}^{\downarrow}\ \ \phi_{2}^{\downarrow}]^{T}$. This is achieved through
the unitary transformation

$\displaystyle\tilde{U}=\displaystyle\left(\begin{tabular}[]{cccc}1&0&0&0\\\
0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&1\end{tabular}\right).$ (S17)

The transformed Hamiltonian,
$\tilde{H}(\textbf{k})=\tilde{U}^{\dagger}H(\textbf{k})\tilde{U}$, is block
diagonal in spin-space,
$\tilde{H}(\textbf{k})=\operatorname{diag}[H^{\uparrow}(\textbf{k})\ \
H^{\downarrow}(\textbf{k})]$ , with

$\displaystyle
H^{s}(\textbf{k})=\displaystyle\left(\begin{tabular}[]{cc}$\Delta+sh_{\textrm{SOC}}(\textbf{k})$&$h_{x}(\textbf{k})-ih_{y}(\textbf{k})$\\\
$h_{x}(\textbf{k})+ih_{y}(\textbf{k})$&$-\Delta-
sh_{\textrm{SOC}}(\textbf{k})$\end{tabular}\right),$ (S20)

where $s=\pm 1$ refers to the two spin states. Note that Hamiltonian (S20)
assumes the compact form
$H^{s}(\textbf{k})=\textbf{h}^{s}(\textbf{k})\cdot\bm{\sigma}$, where
$\textbf{h}^{s}(\textbf{k})=(h_{x}^{s}(\textbf{k}),h_{y}^{s}(\textbf{k}),h_{z}^{s}(\textbf{k}))$
with $h_{x}^{s}(\textbf{k})=h_{x}(\textbf{k})$,
$h_{y}^{s}(\textbf{k})=h_{y}(\textbf{k})$ and
$h_{z}^{s}(\textbf{k})=\Delta+sh_{\textrm{SOC}}(\textbf{k})$. Later, we are
going to utilize this fact to investigate the spin Hall response. If
$h_{0}(\textbf{k})\neq 0$, we obtain
$H^{s}(\textbf{k})=h_{0}(\textbf{k})\sigma_{0}+\textbf{h}^{s}(\textbf{k})\cdot\bm{\sigma}$.

Expanding Eq. (S20) around the valley $\eta\Lambda_{y}$, i.e, assuming
$\textbf{k}=\textbf{q}+\eta(0,\Lambda_{y})$ to linear order in q with
$|\textbf{q}|/\Lambda_{y}\ll 1$, we find

$\displaystyle H_{W}^{\eta
s}(\textbf{q})=(\Delta+s\eta\lambda_{I}\Lambda_{y})\sigma_{z}+\hbar
v_{x}q_{x}\sigma_{y}+\eta\hbar v_{y}q_{y}\sigma_{x},$ (S21)

which is an anistropic version of a pseudo-Weyl Hamitonian, with pseudospin
(sublattice) degree of freedom replacing the real spin, with a spin-valley
dependent mass term. The velocities are $v_{x}=\xi/\hbar$ and
$v_{y}=2\sqrt{-\delta_{0}\gamma_{y}}/\hbar$. By fitting the first principles
band structure, we obtain that the velocities are $v_{x}\approx 2.38\times
10^{5}$ m/s and $v_{y}\approx 4.44\times 10^{5}$ m/s, encoding information
about the structural anisotropy.

## Topological Phase Transition Induced by in-Plane Ferroelectricity

We present here the details of the quantum spin Hall state transitions induced
by fluctuations in $\Delta$ away from $\Delta_{c}$. The antisymmetric matrix
$\langle u_{n}(\textbf{k})|\mathcal{T}|u_{m}(\textbf{k})\rangle$ for the
occupied states $|u_{n,m}(\textbf{k})\rangle$ has a single element $\langle
u_{1}(\textbf{k})|\mathcal{T}|u_{2}(\textbf{k})\rangle$ which corresponds to
its Pfaffian in our four-bands model. Transitions between topologically
trivial and non-trivial insulating phases are accompanied by the absence or
presence, respectively, of zeros of the Pfaffian in the half Brillouin zone .
Thus, we proceed by explicitly calculating the Pfaffian $P(\textbf{k})=\langle
u_{1}(\textbf{k})|\mathcal{T}|u_{2}(\textbf{k})\rangle$ over half Brillouin
zone.

The normalized eigenstates of
$H^{s}(\textbf{k})=h_{0}(\textbf{k})\sigma_{0}+\textbf{h}^{s}(\textbf{k})\cdot\bm{\sigma}$
are

$\displaystyle|u^{s}_{n}(\textbf{k})\rangle=\displaystyle\frac{1}{\sqrt{2h_{s}(\textbf{k})}}\left(\begin{tabular}[]{c}$\sqrt{h_{s}(\textbf{k})+nm_{s}(\textbf{k})}e^{-i\phi(\textbf{k})/2}$\\\
$n\sqrt{h_{s}(\textbf{k})-nm_{s}(\textbf{k})}e^{i\phi(\textbf{k})/2}$\end{tabular}\right),$
(S24)

where $n=\pm 1$ is the band index, $s=\uparrow,\downarrow$ is the spin index,
$h_{s}(\textbf{k})=\sqrt{h_{x}(\textbf{k})^{2}+h_{y}(\textbf{k})^{2}+m_{s}(\textbf{k})^{2}}$,
$m_{s}(\textbf{k})=\Delta+sh_{SOC}(\textbf{k})$ and
$\phi(\textbf{k})=\arctan(h_{y}(\textbf{k})/h_{x}(\textbf{k}))$. The
associated energies are
$\epsilon_{n}^{s}(\textbf{k})=h_{0}(\textbf{k})+nh_{s}(\textbf{k})$. The
states with the lowest energies $\epsilon_{1}=\epsilon_{-}^{\uparrow}$ and
$\epsilon_{2}=\epsilon_{-}^{\downarrow}$ are

$\displaystyle|u_{1}(\textbf{k})\rangle=\displaystyle\left(\begin{tabular}[]{c}$\mathcal{A}_{-}^{\uparrow}(\textbf{k})e^{-i\phi(\textbf{k})/2}$\\\
-$\mathcal{A}_{+}^{\uparrow}(\textbf{k})e^{i\phi(\textbf{k})/2}$\\\ 0\\\
0\end{tabular}\right),\ \
|u_{2}(\textbf{k})\rangle=\displaystyle\left(\begin{tabular}[]{c}0\\\ 0\\\
$\mathcal{A}_{-}^{\downarrow}(\textbf{k})e^{-i\phi(\textbf{k})/2}$\\\
-$\mathcal{A}_{+}^{\downarrow}(\textbf{k})e^{i\phi(\textbf{k})/2}$\end{tabular}\right),$
(S33)

respectively. Here,
$\mathcal{A}_{n}^{s}(\textbf{k})=\sqrt{[h_{s}(\textbf{k})+nm_{s}(\textbf{k})]/2h_{s}(\textbf{k})}$.
The time-reversal operator in this basis is

$\displaystyle\mathcal{T}=\displaystyle\left(\begin{tabular}[]{cccc}0&0&-1&0\\\
0&0&0&-1\\\ 1&0&0&0\\\ 0&1&0&0\end{tabular}\right)\mathcal{K}.$ (S38)

A direct computation gives

$\displaystyle
P(\textbf{k})=-\mathcal{A}_{-}^{\uparrow}(\textbf{k})\mathcal{A}_{-}^{\downarrow}(\textbf{k})e^{i\phi(\textbf{k})}-\mathcal{A}_{+}^{\uparrow}(\textbf{k})\mathcal{A}_{+}^{\downarrow}(\textbf{k})e^{-i\phi(\textbf{k})}.$
(S39)

Figure S40 shows the evolution of the Pfaffian with $\Delta$ over the half
Brillouin zone $k_{y}>0$. For $\Delta=0$ eV, the Pfaffian zeros occur along a
ring aroung the $\Gamma$ point. At $\Delta=-16$ meV, a single Pfaffian zero
occur at $k_{x}=0$, indicating that the system is a quantum spin Hall
insulator; Time-reversal symmetry topologically protects such zero from
disappearing. For $\Delta=-48$ meV, there are no Pfaffian zeros and the system
is in a trivial insulating phase. The transition to trivial insulator occurs
at $\Delta=\Delta_{c}$ through the pseudo-Weyl phase, as described in previous
sections.

Figure S4: Evolution of the Pfaffian of $\langle
u_{n}(\textbf{k})|\mathcal{T}|u_{m}(\textbf{k})\rangle$ over half of Brillouin
zone for (a) $\Delta=0$ meV, (b) $\Delta=-16$ meV and (c) $\Delta=-48$ meV.

The presence or absence of a Pfaffian zero is associated with a change of the
sign of the Dirac mass term. To show this is in fact the case, we first assume
$\Delta\neq 0$. Because the Pfaffian zero is confined at $k_{x}=0$ in this
case, $h_{y}(\textbf{k})=0$ and
$h_{x}(\textbf{k})=\delta_{0}+\gamma_{y}k_{y}^{2}$. Further simplification can
be attained by focusing on $k_{y}\approx\Lambda_{y}$, which enables us to
assume $h_{x}(\textbf{k})\approx 0$. Within these assumptions
$\mathcal{A}_{n}^{s}(\textbf{k})\approx\sqrt{[|m_{s}(\textbf{k})|+nm_{s}(\textbf{k})]/2|m_{s}(\textbf{k})|}$.
Note that $|z|\pm z=[1\pm\operatorname{sign}(z)]|z|$ for arbitrary $z$.
Therefore,
$\mathcal{A}_{n}^{s}\approx\sqrt{[1+n\operatorname{sign}(m_{s})]/2}$. If
$m_{\uparrow}<0$, as in our case, $\mathcal{A}_{+}^{\uparrow}=0$ and

$\displaystyle|P(\textbf{k})|=\sqrt{(1-\operatorname{sign}(m_{\downarrow}))/2},$
(S40)

showing that a mass sign change results in transitions between finite and
vanishing $|P(\textbf{k})|$.

We conclude that fluctuations of $\Delta$ around the critical value
$\Delta_{c}$ will transition the system to topologically trivial or non-
trivial insulating states. Because, the polarization magnitude in two-
dimensional SnSe-type ferroelectrics can be tuned by perpendicular electric
fields or strain , we argue that the topological phase transitions can be
fully controlled in $\alpha$-Bi/SnSe van der Waals heterostructures.

## Topological surface states

In this section, we study the $\alpha$-Bi/SnSe bilayer surface states and
their dependence on the ferroelectric polarization of SnSe from the
perspective of the lattice model. Free standing $\alpha$-Bi is a time-reversal
symmetric topological insulator and, therefore, possesses topologically robust
helical surface states whose energies lie within its bulk band gap. From the
previous section, we anticipate an unique interplay between topological
surface states and ferroelectricity across the insulator-to-semimetal
transition, which raises the question about how the helical surface states
evolve.

Figure S5: Evolution of topological surface states along the zigzag edges with
the polarization-induced symmetry breaking parameters $\Delta$ across the
insulator-to-semimetal transition. Local spin density of states intensity map
for the upper edge in panels (a)-(d) and for the lower edge in panels (e)-(h).
The bottom sketches summarize the information about the evolution of helical
surface states with the SnSe polarization.

To investigate the evolution of the topological surface modes, we can
discretize the bulk Hamiltonian along the $x$ or $y$ direction, for armchair
and zigzag edges. This gives

$\displaystyle\mathcal{\hat{H}}=\sum_{j}[H_{S}(\textbf{k}_{||})c_{\textbf{k}_{||},j}^{\dagger}c_{\textbf{k}_{||},j}+\hat{W}_{\rm{c}}(\textbf{k}_{||})c_{\textbf{k}_{||},j+1}^{\dagger}c_{\textbf{k}_{||},j}+\hat{W}_{\rm{c}}^{\dagger}(\textbf{k}_{||})c_{\textbf{k}_{||},j-1}^{\dagger}c_{\textbf{k}_{||},j}],$
(S41)

from where the “onsite” and “hopping” Hamiltonians for the discretized system
can be obtained as $H_{S}$ and $W_{c}$, respectively. In addition,
$\textbf{k}_{||}=k_{x}$ or $\textbf{k}_{||}=k_{y}$ is momentum along the
direction where the system is still translationally periodic, $j$ is the
principal layer (PL) index and $c_{\textbf{k}_{||},j}$
($c_{\textbf{k}_{||},j}^{\dagger}$) annihilates (creates) a particle with
momentum $\textbf{k}_{||}$ at the $j$-th PL. A semi-infinite lead is described
by a tridiagonal block matrix Hamiltonian, given by

$\displaystyle H=\left(\begin{array}[]{cccc}\ddots&\ddots&&\\\
\ddots&H_{\rm{S}}&\hat{W}_{\rm{c}}&\\\
&\hat{W}_{\rm{c}}^{\dagger}&H_{\rm{S}}&\hat{W}_{\rm{c}}\\\
&&\hat{W}_{\rm{c}}^{\dagger}&H_{\rm{S}}\\\ \end{array}\right),$ (S46)

The surface spin density of states $S_{z}(\textbf{k}_{||},\epsilon)$, is

$\displaystyle
S_{\alpha}(\textbf{k}_{||},\epsilon)=-\displaystyle\frac{1}{\pi}\textrm{Im}\textrm{Tr}[g_{\textrm{S}}(\textbf{k}_{||},\epsilon+i\eta)\sigma_{0}s_{\alpha}],$
(S47)

where $\eta=10^{-3}$ eV is the broadening parameter assumed in all
calculations.

Figure S5 shows the evolution of the helical surface states, along both zigzag
edges referred to as upper and lower edges, with the polarization-induced
symmetry breaking parameters $\Delta$ across the insulator-to-semimetal
transition. At $\Delta=0$ eV, upper and lower edge states are related by a
$\mathcal{M}_{x}$ mirror plane, as seen from Figs. S5(a) and (e). The system
is in a quantum spin Hall state with dissipationless couterpropagating spins,
as shown in the corresponding skectch. The $\Delta=0$ eV corresponds to the
free standing $\alpha$-Bi case.

At $\Delta=-0.016$ eV, the $\mathcal{M}_{x}$ symmetry is broken due to the
finite bulk polarization. Hence, upper and lower edge states cannot be related
by such a symmetry and their helical surface state dispersion now differ
qualitatively, as shown in Figs. S5(b) and (f). We interpret such an
asymmetric feature as being a consequence of the development of surface
charges on the zigzag edges due to the finite bulk polarization. Here, the
surface charges can only perturb the surface states, but cannot destroy them
since a topological phase transition has not taken place.

Figure S6: Evolution of topological surface states along the armchair edges
with the polarization-induced symmetry breaking parameters $\Delta$ across the
insulator-to-semimetal transition. Local spin density of states intensity map
for the upper and bottom edges in panels (a)-(d). In this case, the two edges
are related by the mirror symmetry $\mathcal{M}_{y}$. The bottom sketches
summarize the information about the evolution of helical surface states with
the SnSe polarization.

At $\Delta=-0.032$ eV, the critical transition value $\Delta_{c}$, a gap
closure is attained and the system transitions to a semimetallic phase. Note
that the upper zigzag surface still host helical edge states, Figs. S5(c),
while the lower zigzag edge does not, Figs. S5(g). Here, the surface charges
on the lower edge are enough to dissolve the edge states into the bulk, which
is accompanied by the merging of the helical modes to the bulk bands, as shown
Figs. S5(g).

At $\Delta=-0.048$ eV, or higher values, the energy gap reopens again and the
system is now topologically trivial with no helical states. The edge modes of
the upper zigzag edge remain attached to conduction bands while those of the
lower edge have been successfully become bulk valence states.

Figure S6 shows the evolution of the helical surface states, along both
armchair edges, with the polarization-induced symmetry breaking parameters
$\Delta$ across the insulator-to-semimetal transition. Helical surface states
are revealed at $\Delta=0$ eV, consistent with the quantum spin Hall
insulating state of free standing $\alpha$-Bi. At $\Delta\neq 0$ eV, we
observe that the surface states of both armchair edges are identically the
same, irrespective of the magnitude of $\Delta$. This is due to the fact that
the polarization does not break the mirror plane $\mathcal{M}_{y}$ and, hence,
the two armchair edges modes are related to each other. We only plot the local
spin density of states for one of the armchair edges in Figure S6, which is
referred to as the left edge and is illustrated in the corresponding sketch.
As $\Delta$ is increased in magnitude across the transition point, the both
edge modes are dissolved into the bulk.

Note that the situation here is quite distinct from the zigzag edge case,
where helical surface state still survive in one of the edges (upper). The
interplay between polarization and topological surface states enables one to
selectively control in which edge the helical states propagate by means of
electric field switching.

## Berry curvature and spin Hall Response

Here, we examine the spin Hall effect and its dependence on the in-plane
polarization. Because the spin Hamiltonians are decoupled, it is also possible
to write the velocity of a state $n$ at valley $\eta$ for a given spin $s$. To
linear order in electric field, E, the velocity is

$\displaystyle\textbf{v}_{n}^{\eta
s}(\textbf{k})=\displaystyle\frac{1}{\hbar}\frac{\partial\epsilon_{ns}^{\eta}(\textbf{k})}{\partial\textbf{k}}-\frac{e}{\hbar}\textbf{E}\times\bm{\Omega}_{n}^{\eta
s}(\textbf{k}),$ (S48)

where ${\Omega}_{n}^{\eta s}(\textbf{k})$ is the Berry curvature. The
components of the Berry curvature vector appearing in the velocity are related
to the components of the Berry curvature tensor, $\Omega_{n,\mu\nu}^{\eta
s}(\textbf{k})$, through $\Omega_{n,\mu\nu}^{\eta
s}(\textbf{k})=\epsilon_{\mu\nu\rho}\Omega_{n,\rho}^{\eta s}(\textbf{k})$,
where $\epsilon_{\mu\nu\rho}$ are the components of the Levi-Cività
antisymmetric tensor.

Equation (S21) can be rewritten as general two-level system Hamiltonian, i.e.,
$H_{W}^{\eta s}(\textbf{q})=\textbf{h}^{\eta s}(\textbf{q})\cdot\bm{\sigma}$,
for which the Berry curvature tensor components are most generally given by

$\displaystyle\Omega_{n,\mu\nu}^{\eta
s}(\textbf{q})=\displaystyle\frac{n}{2[h^{\eta
s}(\textbf{q}))]^{3}}\sum_{klm}\epsilon_{klm}h^{\eta
s}_{k}(\textbf{q})\frac{\partial h^{\eta s}_{l}(\textbf{q})}{\partial
q_{\mu}}\frac{\partial h^{\eta s}_{m}(\textbf{q})}{\partial q_{\nu}},$ (S49)

where $h^{\eta s}(\textbf{q})=\sqrt{h_{x}^{\eta s}(\textbf{q})^{2}+h_{y}^{\eta
s}(\textbf{q})^{2}+h_{z}^{\eta s}(\textbf{q})^{2}}$. The $xy$ component
derived from Hamiltonian (S21) and Eq. (S49) is

$\displaystyle\Omega_{n,xy}^{\eta
s}(\textbf{q})=\displaystyle-n\frac{v_{x}v_{y}}{2}\frac{\eta\Delta+s\lambda_{I}\Lambda_{y}}{[(v_{x}q_{x})^{2}+(v_{y}q_{y})^{2}+(\eta\Delta+s\lambda_{I}\Lambda_{y})^{2}]^{3/2}},$
(S50)

which is fully time-reversal symmetric

$\displaystyle\Omega_{n,xy}^{\eta
s}(\textbf{q})=-\Omega_{n,xy}^{(-\eta)(-s)}(-\textbf{q}),$ (S51)

as one should expect. The explicit expression of the Berry curvature tensor
for the two polarization states are given below

$\displaystyle\Omega_{n,xy}^{\eta
s}(\textbf{q})=\left\\{\begin{tabular}[]{cc}$\displaystyle-n\frac{v_{x}v_{y}\Delta_{c}}{2}\frac{s+\eta}{[(v_{x}q_{x})^{2}+(v_{y}q_{y})^{2}+(\Delta_{c})^{2}(s+\eta)^{2}]^{3/2}}$,&if
$\Delta=+\Delta_{c}$\\\ $\displaystyle
n\frac{v_{x}v_{y}\Delta_{c}}{2}\frac{s-\eta}{[(v_{x}q_{x})^{2}+(v_{y}q_{y})^{2}+(\Delta_{c})^{2}(s-\eta)^{2}]^{3/2}}$,&if
$\Delta=-\Delta_{c}$\end{tabular}\right.,$ (S54)

from which we conclude that
$\Omega_{n,xy}^{-\uparrow}(\textbf{q}\neq\textbf{0})=\Omega_{n,xy}^{+\downarrow}(\textbf{q}\neq\textbf{0})=0$
for $\Delta=+\Delta_{c}$, whereas
$\Omega_{n,xy}^{+\uparrow}(\textbf{q}\neq\textbf{0})=\Omega_{n,xy}^{-\downarrow}(\textbf{q}\neq\textbf{0})=0$
for $\Delta=-\Delta_{c}$, and singular at $\textbf{q}=\textbf{0}$. Note that
these Berry curvature components correspond to the semimetallic bands at the
two valleys. The singular behavior is analogous to that in monolayer graphene,
but with a two-fold degenerate band crossing in the $\alpha$-Bi/SnSe case.

Figure S7: Fermi level Spin Berry curvature of $\alpha$-Bi/SnSe obtained from
first principles. The two peaks around the two valleys are well described by
Eq. (S50), where the linear pseudo-Weyl approximation is valid.

The non-singular Berry curvature components correspond to the spin-polarized
gapped bands. Their magnitudes are

$\displaystyle\Omega_{n,xy}^{+\uparrow}(\textbf{q})=-\Omega_{n,xy}^{-\downarrow}(\textbf{q})=\displaystyle-n\frac{v_{x}v_{y}\Delta_{c}}{[(v_{x}q_{x})^{2}+(v_{y}q_{y})^{2}+(2\Delta_{c})^{2}]^{3/2}},\
\ \Delta=+\Delta_{c}$ (S55)

and

$\displaystyle\Omega_{n,xy}^{+\downarrow}(\textbf{q})=-\Omega_{n,xy}^{-\uparrow}(\textbf{q})=\displaystyle-n\frac{v_{x}v_{y}\Delta_{c}}{[(v_{x}q_{x})^{2}+(v_{y}q_{y})^{2}+(2\Delta_{c})^{2}]^{3/2}},\
\ \Delta=-\Delta_{c}.$ (S56)

Thus, the total Berry curvature at a given valley is the sum of the Berry
curvature of semimetallic and insulating bands. The contributions from the
semimetallic bands is singular, whereas the contributions from the insulating
bands is finite for a larger range of momenta around the valley.

Our first principles calculations, however, reveals a tiny energy gap,
$\delta_{g}\approx 6.2$ meV, in the lowest energy bands of $\alpha$-Bi/SnSe,
implying that the induced symmetry breaking polarization due to the proximity
with SnSe, $\Delta_{\textrm{SnSe}}$, slightly deviates from the critical
value, estimated to be $\Delta_{c}\approx 33.8$ meV. Hence, we identify
$\delta_{g}=4|\Delta_{\textrm{SnSe}}-\Delta_{c}|$, where
$\delta_{g}/4\Delta_{c}\approx 0.046\ll 1$ is the gap ratio between lower and
higher energy bands in $\alpha$-Bi/SnSe. As a consequences, the Berry
curvature components of the bands with lower energy is non-singular and finite
over a certain range in momentum space. Its magnitude is much larger than that
of the gapped bands. Therefore, when analysing the Berry curvature components
of the lower energy bands, we will assume a small energy gap. The perfect
semimetal limit can also be obtained by leaving $\delta_{g}\rightarrow 0$.

Figure S7 shows the fermi level spin Berry curvature obtained from the first
principles calculation. The two Berry curvature hotspots located at the
valleys can be treated as independent at the vicinity of the fermi level, for
which the linear approaximation is valid. This is due to the fact that the
valley index is a good quantum number in this limit. The Berry curvature in
each valid, then, follows Eq. (S50) with two peaks centered at the valleys, as
explicitly shown in Fig. S7(b).

Utilizing Eq. (S48), the current density becomes $J_{x}^{\eta
s}=\sigma_{xy}^{\eta s}E_{y}$, where the conductivity is

$\displaystyle\sigma_{xy}^{\eta
s}(\epsilon_{F})=-\frac{e^{2}}{\hbar}\int\frac{d^{2}\textbf{q}}{(2\pi)^{2}}\sum_{n}f_{n\textbf{q}}^{\eta
s}(\epsilon_{F})\Omega_{n,xy}^{\eta s}(\textbf{q}).$ (S57)

In the following, we will show that the spin Hall conductivity can be reversed
by switching the in-plane polarization of the SnSe layer. To this end, we
begin by defining the spin Hall conductivity as

$\displaystyle\sigma_{xy}^{z}(\epsilon_{F})=\frac{\hbar}{2e}\sum_{\eta
s}s\sigma_{xy}^{\eta s}(\epsilon_{F}).$ (S58)

We focus on the spin response due to states near the Fermi level
$\epsilon_{F}$, for which we have $s=-\eta\operatorname{sign}(\Delta_{c})$.
The spin Hall conductivity becomes

$\displaystyle\sigma_{xy}^{z}(\epsilon_{F})=\frac{e}{2}\operatorname{sign}(\Delta_{c})\int\frac{d^{2}\textbf{q}}{(2\pi)^{2}}\sum_{n\eta
s}\eta f_{n\textbf{q}}^{\eta s}(\epsilon_{F})\Omega_{n,xy}^{\eta
s}(\textbf{q}).$ (S59)

Because of time-reversal symmetry, $f_{n\textbf{q}}^{\eta
s}(\epsilon_{F})=f_{n\textbf{q}}^{(-\eta)(-s)}(\epsilon_{F})$, where
$f_{n\textbf{q}}^{\eta s}(\epsilon_{F})=[1+\exp{(\epsilon_{n}^{\eta
s}(\textbf{q})-\epsilon_{F})/k_{B}T]}^{-1}$ is the Fermi-Dirac occupation
function. Hence, the summation becomes

$\displaystyle\sum_{n\eta s}\eta f_{n\textbf{q}}^{\eta
s}(\epsilon_{F})\Omega_{n,xy}^{\eta
s}(\textbf{q})=\sum_{n}f_{n\textbf{q}}^{+\downarrow}(\epsilon_{F})\Omega_{n,xy}^{+\downarrow}(\textbf{q})-f_{n\textbf{q}}^{-\uparrow}(\epsilon_{F})\Omega_{n,xy}^{-\uparrow}(\textbf{q})=2\sum_{n}f_{n\textbf{q}}^{+\downarrow}(\epsilon_{F})\Omega_{n,xy}^{+\downarrow}(\textbf{q}),$
(S60)

by taking advantage of Eq. (S51) and $\Omega_{n,xy}^{\eta
s}(\textbf{q})=\Omega_{n,xy}^{\eta s}(|\textbf{q}|)$ at small q, i.e., near
the Fermi level. Therefore,

$\displaystyle\sigma_{xy}^{z}(\epsilon_{F})=\operatorname{sign}(\Delta_{c})e\int\frac{d^{2}\textbf{q}}{(2\pi)^{2}}\sum_{n}f_{n\textbf{q}}^{+\downarrow}(\epsilon_{F})\Omega_{n,xy}^{+\downarrow}(\textbf{q}),$
(S61)

which depends on the polarization states of the SnSe layer through
$\operatorname{sign}(\Delta_{c})$. Thus, the spin Hall currents can be
reversed by switching the ferroelectric order parameter of SnSe. The result
derived here is not confined to toy model and, thus, also applied to a more
rigorous first principles calculations as long as the physics at the vicinity
of the Fermi level is considered.

From the explicit valley/spin-resolved Berry curvature derived for the toy
model, it is possible to obtain analytical expressions for the anomalous
current density for electrons at valley $\eta$ and spin $s$. For simplicity,
we examine the zero temperature conductivities in the following, i.e., the
Fermi-Dirac distribution is taken to be $f_{n\textbf{q}}^{\eta
s}(\epsilon_{F})=\Theta(\epsilon_{F}-\epsilon_{ns}^{\eta}(\textbf{q}))$, where
$\Theta$ is the Heaviside function. The integration can be performed very
easily and renders the important results

$\displaystyle\sigma_{xy}^{+[\operatorname{sign}(\Delta)]}(\epsilon_{F})=-\sigma_{xy}^{-[-\operatorname{sign}(\Delta)]}(\epsilon_{F})=\tilde{\sigma}_{xy}(\epsilon_{F}),$
(S62)

where

$\displaystyle\tilde{\sigma}_{xy}(\epsilon_{F})=\left\\{\begin{tabular}[]{cc}$\displaystyle-\frac{e^{2}}{h}\frac{\Delta_{c}}{|\epsilon_{F}|}$,&if
$|\epsilon_{F}|>E_{g}/2$\\\ $\displaystyle-\frac{e^{2}}{h}\frac{1}{2}$,&if
$|\epsilon_{F}|<E_{g}/2$\end{tabular}\right.,$ (S65)

for the gapped bands.

## Screening length

In addressing the feasibility of switching the in-plane polarization with
external electric fields, it is necessary to account for the screening effects
due to itinerant charges. Because the $\alpha$-Bi/SnSe is a semimetal,
however, we expect large screening lengths, as we address here. The Thomas-
Fermi dielectric function, $\varepsilon$, is

$\displaystyle\varepsilon(|\textbf{q}|)=\displaystyle\kappa+\frac{e^{2}}{2\epsilon_{0}|\textbf{q}|}D(\epsilon_{F}),$
(S66)

for the semimetallic bands obatined from Eq. (S21) with $\Delta=\Delta_{c}$
and $\operatorname{sign}(\Delta_{c})+\eta s=0$. The density of states at the
fermi level $\epsilon_{F}$ can be expressed in terms of the carrier
concentration $n$ as

$\displaystyle D(\epsilon_{F})=\displaystyle\frac{1}{\hbar}\sqrt{\frac{n}{2\pi
v_{x}v_{y}}},$ (S67)

where $\kappa$ is the dieletric constant, from where the Thomas-Fermi
screening length is obtained as

$\displaystyle\lambda_{\textrm{TF}}=\displaystyle\frac{2\epsilon_{0}\kappa\hbar}{e^{2}}\sqrt{\frac{2\pi
v_{x}v_{y}}{n}}.$ (S68)

Assuming a dielectric constant of $\kappa=3$ and typical 2D semimetal carrier
concentrations of $n=10^{11}$ cm-2, we obtain $\lambda_{\textrm{TF}}\approx 6$
nm.

## Charge-Spin Interconversion

The charge-spin interconversion coefficients were computed from linear
response theory. We utilize the Kubo formula fashioned after Smrcka-Streda
kubo1; kubo2; kubo3; kubo4:

$\displaystyle\delta
O_{\alpha\beta}^{\gamma}=\displaystyle\int\frac{d\textbf{k}}{(2\pi)^{2}}[\delta
O_{\alpha\beta}^{\gamma,I}(\textbf{k})+\delta
O_{\alpha\beta}^{\gamma,II}(\textbf{k})],$ (S69)

with integrands

$\delta
O_{\alpha\beta}^{{\gamma},I}(\textrm{k})=\displaystyle-\frac{e\hbar}{\pi}{\Gamma}^{2}\sum_{nm}\frac{\operatorname{Re}[\langle
n\textbf{k}|\hat{O}_{\alpha}^{\gamma}|m\textbf{k}\rangle\langle
m\textbf{k}|\hat{v}_{\beta}|n\textbf{k}\rangle]}{[({\epsilon}_{F}-{\epsilon}_{n{\textbf{k}}})^{2}+{\Gamma}^{2}][({\epsilon}_{F}-{\epsilon}_{m{\textbf{k}}})^{2}+{\Gamma}^{2}]},$
(S70a) $\delta
O_{\alpha\beta}^{\gamma,II}(\textrm{k})=\displaystyle-2e\hbar\sum_{n,m\neq
n}(f_{n\textbf{k}}-f_{m\textbf{k}})\frac{\operatorname{Im}[\langle
n\textbf{k}|\hat{O}_{\alpha}^{\gamma}|m\textbf{k}\rangle\langle
m\textbf{k}|\hat{v}_{\beta}|n\textbf{k}\rangle]}{({\epsilon}_{n\textbf{k}}-{\epsilon}_{m\textbf{k}})^{2}+{\Gamma}^{2}},$
(S70b)

valid in the weak disorder limit described by a constant band broadening
$\Gamma$, where $\hat{v}_{\beta}$ is the $\beta=x,y$ component of the velocity
operator, $\hat{O}_{\alpha}^{\gamma}$ is the perturbed physical observable
with spin index $\gamma=x,y,z$ and $|n\textbf{k}\rangle$ is the eigenstate
associated with the band $\epsilon_{n\textbf{k}}$ of the unperturbed system.
The charge and spin Hall responses are obtained through
$\hat{O}_{\alpha}^{\gamma}\rightarrow-e\hat{v}_{\alpha}$,
$\hat{O}_{\alpha}^{\gamma}\rightarrow(2/\hbar)\hat{Q}_{\alpha}^{\gamma}$,
where the spin current operator is defined as
$\hat{Q}_{\alpha}^{\gamma}=(1/2)\\{\hat{s}^{\gamma},\hat{v}_{\alpha}\\}$.

## First principles calculations

Figure S8: Band structure of (a) monolayer $\alpha$-Bi, (b) $\alpha$-Bi/SnSe
van der Waals heterostructure and (c) Bi (100)/2L-SnSe. Horizontal dashed line
demarcates the fermi level.

The first-principles calculations based on density functional theory (DFT)
Kohn1965 was carried by out Vienna ab initio simulation package (VASP)
Kresse1996. The exchange-correlation (XC) functional was treated within the
generalized gradient approximation of Perdew-Burke-Ernzerhof (PBE) Perdew1996
with noncollinear spin polarization. PhysRevB.93.224425 The electronic
wavefunctions were expanded by planewave basis with kinetic energy cutoff of
500 eV. We employed the projector-augmented wave pseudopotentials Blochl1994;
Kresse1999 to describe the valence electrons, and Grimme-D3 van der Waals
correction grimme2010consistent was chosen to describe interlayer interaction.

The experimental values of lattice constants were used for
$\alpha$-Bi,($a_{\textrm{armchair}}=4.84$ Å, $b_{\textrm{zigzag}}=4.49~{}$Å)
lu2023observation and tensile strain was applied to SnSe substrate to make
commensurate $\alpha$-Bi/SnSe heterostructure. The sufficiently large vacuum
region ($>$20 Å) was included in the periodic cell to mimic two-dimensional
layered structure. The crystal structure of monolayer $\alpha$-Bi was fully
realxed while maintaining $Cmce$ symmetry. In the $\alpha$-Bi/SnSe
heterostructure, only SnSe was fully realxed to avoid buckled phase of
$\alpha$-Bi layer. Lu2015

The spontaneous polarization P of $\alpha$-Bi/SnSe heterostructure along the
armchair direction was calculated in the framework of the modern theory of
polarization through Berry phase. spaldin2012beginner In the modern theory of
polarization, P can be rigorously defined by adiabatic Wannier charge center
change along an arbitrary path. Therefore, the absolute value of P of
$\alpha$-Bi/SnSe heterostructure is only numerically available when we have
its centerosymmetric (or in-plane mirror symmetric) counterpart which
guarantees $\textbf{P}=0$. However, due to the geometry of the
heterostructure, it is unfeasible to find centerosymmetric counterpart through
the experimentally reasonable reaction path from $-\textbf{P}$ to P. To detour
this problem, we evaluate P of $\alpha$-Bi/SnSe heterostructure through
following steps:

  * •

Constructing structure of $-\textbf{P}$ using mirror operation. (Fig. 1(a) in
the main manuscript) Strictly speaking, it corresponds to P to $-\textbf{P}$
transition of SnSe and in-plane sliding of $\alpha$-Bi.

  * •

Initializing a total of 11 structures along the path from $-\textbf{P}$ to P
using linear interpolation.

  * •

Finding minimum energy path using the nudged elastic band (NEB) calculation.
In the NEB calculation, only SnSe was fully realxed due to the same reason
described above.

  * •

Evaluating P of each structure. Note that, P of some structures in the
reaction path were unavailable because they were calculated to be metallic due
to $\alpha$-Bi. In this case, we determined P of each structure through linear
interpolation.

  * •

Then, we set P of 6th structure (in-plane mirror, Fig. 1(a) in the main
manuscript) to be zero, and we obtained P of $\alpha$-Bi/SnSe heterostructure.