# Towards good validation metrics for generative models in offline model-based
optimisation

Christopher Beckham1,2,3, Alexandre Piché1,2,4, David Vazquez1, Christopher
Pal1,2,3,†
1ServiceNow Research, 2Mila, 3Polytechnique Montréal, 4Université de Montréal,
†CIFAR AI Chair

###### Abstract

In this work we propose a principled evaluation framework for model-based
optimisation to measure how well a generative model can extrapolate. We
achieve this by interpreting the training and validation splits as draws from
their respective ‘truncated’ ground truth distributions, where examples in the
validation set contain scores much larger than those in the training set.
Model selection is performed on the validation set for some prescribed
validation metric. A major research question however is in determining what
validation metric correlates best with the expected value of generated
candidates with respect to the ground truth oracle; work towards answering
this question can translate to large economic gains since it is expensive to
evaluate the ground truth oracle in the real world. We compare various
validation metrics for generative adversarial networks using our framework. We
also discuss limitations with our framework with respect to existing datasets
and how progress can be made to mitigate them.111Code:
https://github.com/christopher-beckham/principled-mbo-eval

##  1 Introduction

Figure 1: We want to produce designs $\bm{x}$ that score high according to the
ground truth oracle $y={\color[rgb]{0.67,0.01,0}{f}}(\bm{x})$, but this is
usually prohibitively expensive to compute. We would like to search for
validation metrics that correlate well with the ground truth, yet are cheap to
compute, in order to facilitate faster and more economical generation of novel
designs. In this illustration we imagine a generative model that is able to
sample a complete specification of a chemical compound for drug design.
However, scoring it involves synthesising the compound in a wet lab, which is
expensive and laborious. Being able to reliably sample compounds on large but
unseen $y$’s is the ultimate task in model-based optimisation.

In model-based optimisation (MBO), we wish to learn a model of some unknown
objective function $f:\mathcal{X}\rightarrow\mathcal{Y}$ where $f$ is the
ground truth ‘oracle’, $\bm{x}\in\mathcal{X}$ is some characterisation of an
input and $y\in\mathbb{R}^{+}$ is some score assigned to the input. The larger
the score is, the more desirable $\bm{x}$ is according to some desiderata. In
practice, such a function (a real world process) is often prohibitively
expensive to compute. For instance if $\bm{x}\in\mathcal{X}$ is a
specification of a protein that binds to a specific receptor and we want to
maximise the binding potential, then actually determining this would involve
synthesising and testing it in a wet lab. In other cases, synthesising and
testing a candidate may also be dangerous, for instance components for
vehicles or aircraft. In MBO, we want to learn models that can extrapolate and
generate inputs whose scores are beyond that of what we have seen in our
dataset. Because the most reliable way of validating designs from such a model
also happens to be the most expensive (i.e. evaluating the ground truth
oracle), we are motivated to search for metrics that correlate well with this
function because this can translate to substantial economic savings.

We consider MBO that is _data-driven_ (leverages machine learning) and is
_offline_ , without assuming an active learning loop where the ground truth
can be queried (say, to augment the training set). In offline data-driven MBO
(referred to here as just MBO), a very simple (yet naïve) approach is to
approximate the ground truth oracle $f$ by training a classifier
$f_{\theta}(\bm{x})$ from some dataset $\mathcal{D}$, and exploiting it
through gradient ascent to generate a high-scoring candidate:

$\displaystyle\bm{x}^{*}$
$\displaystyle=\operatorname*{arg\,max}_{\bm{x}}f(\bm{x})\approx\operatorname*{arg\,max}_{\bm{x}}f_{\theta}(\bm{x}),$
(1)
$\displaystyle\implies\bm{x}_{t+1}\leftarrow\bm{x}_{t}+\alpha\nabla_{\bm{x}}f_{\theta}(\bm{x}),\
\text{for }t=\\{1,\dots,T\\},\bm{x}^{*}=\bm{x}_{T}$

The issue here however is that for most problems, this will produce an input
that is either _invalid_ or is poor yet receives a large score from the
approximate oracle (_overestimation_). This is the case when the space of
valid inputs lies on a low-dimensional manifold in a much higher dimension
space (Kumar & Levine, 2020). We assume that invalid inputs will produce a
score $f(\bm{x})=0$, though it is also possible for completely valid inputs to
have a zero score as well. One concern is that the learned oracle $f_{\theta}$
may end up assigning non-zero scores to inputs that are invalid. The extent to
which this is the case would depend on whether the learned oracle $f_{\theta}$
was able to leverage _negative examples_ during training; that is, invalid
inputs whose scores are all equal to zero. While such a procedure could be
used to regularise the learned oracle (Kumar et al., 2021), this issue can be
addressed directly with generative modelling, without the need for negative
examples or mining for adversarial examples (Trabucco et al., 2021).

Assume that we are able to learn a density $p(\bm{x})$. Since it has finite
mass, namely that $\int_{\bm{x}}p(\bm{x})d\bm{x}=1$, by pushing up the
likelihood of only valid inputs, we automatically push down the likelihood of
everything else (invalid inputs). Because we are dealing with additional
variable $y$ which is the score of the input, we can instead consider the
_joint likelihood_ $p(\bm{x},y)$. Through Bayes’ rule, we can expand the joint
likelihood out into $p(\bm{x},y)=p(y|\bm{x})p(\bm{x})$. $p(y|\bm{x})$ can be
thought of as a function of the oracle $f(\bm{x})$; if $f(\bm{x})$ is
stochastic (it usually is), then $p(y|\bm{x})$ will capture the uncertainty
over a range of values of $y$. This however assumes a particular $\bm{x}$ was
observed, regardless of its validity. Conversely, the multiplication by
$p(\bm{x})$ corrects for this based on how likely it is one observes $\bm{x}$,
which is why it is important to be able to measure the joint likelihood of a
sample (Brookes et al., 2019). Therefore, in this work we consider generative
models of the form $p(\bm{x},y)$.

In this work, we propose a principled evaluation framework with a particular
focus on _model selection_ , i.e. how do we select for the model that best
extrapolates? To do this, one must define a _validation metric_ , i.e., some
metric that is a function of the generative model and a held-out validation
set whose examples comprise scores much larger than that of the training set.
We would like to find a validation metric which is _highly correlated with the
ground truth oracle_ , since it is expensive to compute in practice. To
approach this, we present our framework in Section 2, which in turn comprises
the following parts: training (2.1), model selection (2.2), and generation
(2.3). Lastly, in Section 4.1 we discuss some limitations of our framework
with respect to existing datasets, and propose some alternatives.

##  2 Proposed method

Let us assume we have access to $p(\bm{x})$ (or its joint equivalent). Indeed,
$p(\bm{x})$ is measuring a likelihood, and it is likely the case that the
best-scoring input lies in a _low density_ region of the data distribution –
after all, if commonly occuring inputs were high scoring then MBO would not be
needed at all. We do not wish to use likelihood to select for _likely
candidates_ ; instead, we want to use it as a means to measure how well a
generative model can extrapolate. The intuition is that a generative model
that assigns high likelihood to _unseen examples_ whose $y$’s are larger than
what it has seen in training has effectively learned how to _extrapolate_.
This is because the probability of the model sampling such a high-scoring
candidate is directly proportional to the likelihood (density) that it has
assigned to it. Therefore, selecting models based on this notion of
extrapolation means selecting for ones that produce the most _reliable_
samples for $y$’s larger than what it has been trained with. Because we only
have a finite collection of samples $\mathcal{D}$ to work with, we can split
the dataset up in such a way that makes it amenable to measuring
extrapolation. This can be done by constructing training and validation splits
(${\color[rgb]{0,0,0.5}{{\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}}}$
and
${\color[rgb]{0,0.5,0}{{\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}}}$),
where both splits can be seen as samples drawn from their respective
$y$-_truncated_ distributions. In other words, if the ground truth is
${\color[rgb]{0.67,0.01,0}{p}}(\bm{x},y)$ (the distribution whose samples
comprise $\mathcal{D}$), we can define two splits
${\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(\bm{x},y)$ and
${\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)$, where
${\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(\bm{x},y)=\\{(\bm{x},y)\sim{\color[rgb]{0.67,0.01,0}{p}}|y\in[0,\gamma]\\}$
and
${\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)=\\{(\bm{x},y)\sim{\color[rgb]{0.67,0.01,0}{p}}|y\in[\gamma,\infty]\\}$.Therefore,
${\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}$ is a finite collection of
samples from ${\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)$ and
${\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}$ is for
${\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(\bm{x},y)$. This is illustrated in
Figure 2(a).

(a)

(b)

Figure 2: To preface the following, we colour-code most of our equations to
indicate the following: blue to represent the training set
${\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}$ or models trained on the
training set, green to denote models and datasets pertaining to a a held-out
validation set ${\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}$, and
orange to refer to ground truth (usually unknown) or test set
${\color[rgb]{0.67,0.01,0}{\mathcal{D}_{\text{test}}}}$. (2(a)) our evaluation
involves interpreting our dataset of observations $\mathcal{D}$ as coming from
two truncated distributions:
${\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}$ only contains examples
whose $y$’s are in the range $[0,\gamma]$ and conversely
${\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}$ contains those in
$[\gamma,\infty]$. Here, $\gamma$ is a user-specified threshold. (2(b)) Our
framework is split into four parts: training and extrapolation (2.1), model
selection (2.2), and generation (2.3). Shown in light grey is an approximating
oracle ${\color[rgb]{0,0.5,0}{f_{\phi}}}$ which is used for validation, which
is trained on both ${\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}$ and
${\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}$. Since some validation
metrics may not involve the use of the validation oracle, it is shown with
faded/dotted lines.

###  2.1 Training and extrapolation

We assume a conditional generative model
${\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y)$ that has been trained on
${\color[rgb]{0,0,0.5}{{\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}}}$.
If we denote the distribution over $y$’s in the training set as
${\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(y)$ then through Bayes’ rule we can
write the joint likelihood as:
${\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x},y)={\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y){\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(y)$.
This equation essentially says: to generate a sample $(\bm{x},y)$ from this
joint distribution, we first sample
$y\sim{\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(y)$, then we sample
$\bm{x}\sim{\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y)$. Now if we assume
that the generative model has some ability to extrapolate, then we should be
able to replace the prior in this equation with
${\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(y)$, which is the prior distribution
over $y$ for the validation set. We define this as the _extrapolated model_
(Figure 2, _extrapolate_ caption):

$\displaystyle{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x},y)={\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y){\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(y).$
(2)

###  2.2 Model selection

Suppose have a generative model that allows us access to the evaluation to
either a density $p_{\theta}(\bm{x},y)$ or a close approximation of it. This
is possible with all commonly used generative models in deep learning, with
the exception of generative adversarial networks (GANs) (Goodfellow et al.,
2014). To ensure that our model can generate candidates with scores higher
than that observed in the training set (i.e. _extrapolate_), we need to employ
a reliable _model selection_ criteria. Suppose we trained several generative
models $\Theta=\\{\theta_{j}\\}_{j=1}^{M}$ and wanted to select which one was
the best, one simple strategy is to select for the model whose forward KL
divergence with the $[\gamma,\infty]$-_truncated_ ground truth
${\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)$ is the smallest (full proof
is in Section A.2):

$\displaystyle\theta^{*}$
$\displaystyle=\operatorname*{arg\,min}_{\theta\in\Theta}\
\underbrace{\text{KL}\Big{[}{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)\
\|\ {\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x},y)\Big{]}}_{\text{forward
KL}}\approx\operatorname*{arg\,max}_{\theta\in\Theta}\
\frac{1}{N}\sum_{i=1}^{N}\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}({\color[rgb]{0,0.5,0}{\bm{x}_{i}}},{\color[rgb]{0,0.5,0}{y_{i}}}),$
(3)

where samples from the validation set can be seen as our finite number of
samples coming from the truncated ground truth
${\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)$ (Figure 2, model selection
panel). This is the most straightforward way to evaluate extrapolation since
requires only _samples_ from the ground truth, not the ground truth
itself.222Minimising forward KL during _training_ is also known as _maximum
likelihood_ estimation, and this is the case for all commonly-used generative
models, with the exception of GANs.

With the exception of normalising flows, the density evaluated by
${\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}$ in Equation 3 will need to be
approximated. For conditional VAEs (Kingma & Welling, 2013) and diffusion
models (Ho et al., 2020) we would be dealing with the _conditional ELBO_ ,
i.e. $\log
p_{\theta,\gamma}(\bm{x}|y)\geq\mathcal{L}_{\text{C-ELBO}}(\bm{x},y)$. The
_joint ELBO_ can simply be derived by adding the prior $\log p_{\gamma}(y)$:

$\displaystyle\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x},y)=\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x}|y)+\log{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(y)\geq\mathcal{L}_{\text{C-ELBO}}(\bm{x},y)+\log{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(y),$
(4)

The forward KL (or its approximation) is trivial to evaluate for these classes
of models since it is equivalent to the training loss, but instead it is
evaluated with respect to a validation set and not a training set. It is also
theoretically convenient to use the same loss to measure both training and
validation performance. However, it is unclear whether forward KL is the
‘best’ validation metric to use for the aforementioned generative models. We
will define what it means for a metric to be ‘best’ in Section 2.3.

Lastly, we mentioned that GANs do not allow for likelihood evaluation. While
this is unfortunate, evaluation metrics for GANs would allow for one to
compare across different classes of generative model, since there is no
assumption that a likelihood can be computed. Furthermore, they may probe
interesting aspects of generative models that are not succinctly captured with
forward KL (e.g. measuring precision and recall Kynkäänniemi et al. (2019)).
While we defer the reader to Borji (2022) for a comprehensive list, many of
these metrics are based on computing distances between distributions: they
only require a finite number of samples from both the ground truth and
generative model, as well as an external neural network which is either used
to classify inputs (Lopez-Paz & Oquab, 2016) or extract semantically
meaningful features (Heusel et al., 2017). For the remainder of this paper we
will call such metrics _likelihood-free_. A fundamental question we would like
to answer is: what validation metrics correlate the best with the ground
truth?

###  2.3 Generation time

To measure which validation metrics correlate the most with the ground truth,
we first need to devise a strategy for generating the best possible samples
from our extrapolated model. One such strategy is to generate a large set
$\mathcal{S}$ of samples with $p_{\theta,\gamma}(\bm{x}|y=y_{\text{max}})$,
where $y_{\text{max}}$ is the largest value of $y$ observed in the entire
dataset. Another is to simply sample from $p_{0,\gamma}(y)$, which is what we
choose. From $\mathcal{S}$, we can either compute the expected value with
respect to _all_ samples
$\mathbb{E}_{\bm{x}\sim\mathcal{S}}[{\color[rgb]{0.67,0.01,0}{f}}(\bm{x})]$ or
the top $n\%$ scoring ones instead (Trabucco et al., 2022). For each trained
model $\\{\theta_{j}\in\Theta\\}_{j=1}^{M}$ (i.e. each corresponding to a
configuration of hyperparameters), compute the expected value of generated
candidates with respect to the ground truth oracle
$\mathbb{E}_{\bm{x}\sim{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x},y)}[{\color[rgb]{0.67,0.01,0}{f}}(\bm{x})]$.
We can then plot the validation metric against this expected value and compute
some goodness of fit, i.e. $R^{2}$ or Pearson correlation.

For additional discussion on generation strategies, please see Section A.4.

##  3 Related work

##### Design Bench

Design Bench is an evaluation framework by Trabucco et al. (2022) that
facilitates the training and evaluation of MBO algorithms. Design Bench, as of
time of writing, provides _four discrete_ and _four continuous_ datasets.
These datasets can be further categorised based on two attributes: whether a
ground truth oracle exists or not, and whether it is possible for there to
exist out of distribution candidates. One of these datasets for instance is
TFBind8, a dataset consisting of length-8 DNA strings that are transcription
factors for binding. Since the dataset consists of all possible strings, it is
fully realised and the ground truth oracle can simply be a look-up table.
However, since any length-8 DNA string is valid (any combination of
nucleotides is possible), there is no such thing as an invalid candidate.
While there are certainly domains where this is the case, we wish to consider
more difficult problem domains where it is possible for invalid candidates to
exist. This also allows us to justify the use of generative models for MBO,
which naturally learn what is out-of-distribution by learning what is _in-
distribution_. Currently, the only datasets that support both OOD examples and
admit a ground truth oracle are ones based on RL environments, which may
require external dependencies such as Mujoco. One such environment for
instance is Hopper, which involves optimising the weights of a neural network
controller to maximise a reward function.

In terms of evaluation, Design Bench does not prescribe a validation set (only
a training set and test oracle), which we argue is important in order to
address the core question of our work, which is finding validation metrics
that correlate well with the ground truth oracle. While a validation set could
be easily derived from $y$-truncating the training set, this would complicate
evaluation since the training, validation and test splits would all come from
separate truncated distributions. Conversely, our evaluation is set up so that
the validation and test sets are from the _same_ $\gamma$-truncated
distribution.

##### Conditioning by adaptive sampling

CbAS (Brookes et al., 2019), like our proposed method, approaches MBO from a
generative modelling perspective with the aim of avoiding the generation of
invalid inputs. Given some pre-trained ‘prior’ on the input data
$p_{\theta}(\bm{x})$, the authors propose the derivation of the conditional
generative model $p_{\theta}(\bm{x}|y)$ via Bayes’ rule:

$\displaystyle
p_{\theta}(\bm{x}|y)=\frac{p(y|\bm{x})p_{\theta}(\bm{x})}{p_{\theta}(y)}=\frac{p(y|\bm{x})p_{\theta}(\bm{x})}{\int_{\bm{x}}p(y|\bm{x})p_{\theta}(\bm{x})d\bm{x}},$
(5)

where $p(y|\bm{x})$ denotes the oracle in probabilistic form, and is not
required to be differentiable. More generally, the authors use $S$ to denote
some target range of $y$’s, i.e.
$p(S|\bm{x})=\int_{y}p(y|\bm{x})\mathbf{1}_{y\in S}dy$.

$\displaystyle
p_{\theta}(\bm{x}|S)=\frac{p(S|\bm{x})p_{\theta}(\bm{x})}{p_{\theta}(S)}=\frac{p(S|\bm{x})p_{\theta}(\bm{x})}{\int_{\bm{x}}p(S|\bm{x})p_{\theta}(\bm{x})d\bm{x}},$
(6)

Due to the intractability of the denominator term, the authors propose the use
of variational inference to learn a sampling network
${\color[rgb]{0,0,0.5}{q_{\zeta}}}(\bm{x})$ that is as close as possible to
${\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|S)$ as measured by the forward KL
divergence. Here, let us use ${\color[rgb]{0,0,0.5}{p_{\theta}}}(S|\bm{x})$ in
place of $p(S|\bm{x})$:333In their paper the symbol $\phi$ is used, but here
we use $\zeta$ since the former is used to denote the validation oracle.:

$\displaystyle\zeta^{*}$ $\displaystyle=\operatorname*{arg\,min}_{\phi}\
\underbrace{\text{KL}\Big{[}{\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|S)\ \|\
{\color[rgb]{0,0,0.5}{q_{\zeta}}}(\bm{x})\Big{]}}_{\text{forward KL}}$
$\displaystyle=\operatorname*{arg\,min}_{\zeta}\
\text{KL}\Big{[}\frac{{\color[rgb]{0,0,0.5}{p_{\theta}}}(S|\bm{x}){\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x})}{p_{\theta}(S)}\
\|\ {\color[rgb]{0,0,0.5}{q_{\zeta}}}(\bm{x})\Big{]}$
$\displaystyle=\operatorname*{arg\,min}_{\zeta}\
\mathbb{E}_{\bm{x}\sim{\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x})}\Big{[}\underbrace{{\color[rgb]{0,0,0.5}{p_{\theta}}}(S|\bm{x})}_{\text{oracle}}\log{\color[rgb]{0,0,0.5}{q_{\zeta}}}(\bm{x})\Big{]}.$
(7)

In practice, importance sampling must be used since
${\color[rgb]{0,0,0.5}{p_{\theta}}}(S|\bm{x})$ will be vanishingly small. In
terms of evaluation, $\gamma$ is set so that samples from
${\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(\bm{x},y)$ constitute the bottom
20th percentile of values. Here, the oracle
${\color[rgb]{0,0,0.5}{p_{\theta}}}(S|\bm{x})$ is only trained on
${\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}$. Furthermore, during the
variational inference procedure there is no evaluation of Equation 3 on a
held-out validation set, which means one can risk overfitting the search model
${\color[rgb]{0,0,0.5}{q_{\zeta}}}$ to the oracle. Because the expected value
in Equation 3 is over samples from the initial model and not the validation
set, the simplest way to incorporate some notion of validation to use
$p(S|\bm{x})\leftarrow{\color[rgb]{0,0.5,0}{p_{\phi}}}(S|\bm{x})$ for
validation. We can write this out as follows:

$\displaystyle\boxed{\mathcal{M}_{\text{CbAS}}=\mathbb{E}_{\bm{x}\sim{\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x})}\Big{[}{\color[rgb]{0,0.5,0}{p_{\phi}}}(S|\bm{x})\log{\color[rgb]{0,0,0.5}{q_{\zeta}}}(\bm{x})\Big{]}}$
(8)

The main difference between CbAS and our work is that the former requires an
additional fine-tuning step to convert the base model into the extrapolated
model, whereas our extrapolated model is simply induced by switching the
prior. Because the extrapolated model optimisation in CbAS depends on an the
approximate _training oracle_ ${\color[rgb]{0,0,0.5}{p_{\theta}}}(S|\bm{x})$
(Equation 3), some regularisation may be needed to ensure that overfitting
does not occur. Lastly, the extrapolated model does not condition on $y$,
though one could simply define one at inference time via
${\color[rgb]{0,0,0.5}{q_{\zeta}}}(\bm{x}){\color[rgb]{0,0.5,0}{p_{\phi}}}(y|\bm{x})$.

##### Model inversion networks

Generative modelling for MBO was proposed by Kumar & Levine (2020), under the
name _model inversion networks_ (MINs). The name is in reference to the fact
that one can learn the _inverse_ of the oracle
$f_{\theta}^{-1}:\mathcal{Y}\rightarrow\mathcal{X}$, which is a generative
model. It is not clear what model selection criteria was used during training,
though at _generation time_ the authors propose the following equation to
search for high-scoring candidates:444This differs slightly from the
formulation in their paper, since we express the optimisation with soft
constraints, rather than hard constraints.

$\displaystyle y^{*},\bm{z}^{*}=\operatorname*{arg\,max}_{y,\bm{z}}\
y+\epsilon_{1}\underbrace{\log{\color[rgb]{0,0,0.5}{p_{\theta}}}(y|{\color[rgb]{0,0,0.5}{G_{\theta}}}(\bm{z},y))}_{\text{agreement}}+\epsilon_{2}\log
p(\bm{z}),$ (9)

where ${\color[rgb]{0,0,0.5}{G_{\theta}}}(\bm{z},y)$ denotes the GAN but
written explicitly (rather than using sampling notation
$\bm{x}\sim{\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y)$), and $p(\bm{z})$ is
the prior distribution over $\bm{z}$. The agreement term is measuring the log
likelihood of $G_{\theta}(\bm{z},y)$ being classified as $y$ under the
training oracle ${\color[rgb]{0,0,0.5}{f_{\theta}}}$. Note that if we assume
${\color[rgb]{0,0,0.5}{p_{\theta}}}(y|\bm{x})=\mathcal{N}(y;{\color[rgb]{0,0,0.5}{f_{\theta}}}(\bm{x}),\sigma^{2})$,
then maximising the agreement in Equation 9 is equivalent to minimising the
mean squared error (MSE) between the conditioning $y$ and its predicted $y$
with respect to $f_{\theta}$:

$\displaystyle-\log{\color[rgb]{0,0,0.5}{p_{\theta}}}(y|{\color[rgb]{0,0,0.5}{G_{\theta}}}(\bm{z},y))\propto\|{\color[rgb]{0,0,0.5}{f_{\theta}}}({\color[rgb]{0,0,0.5}{G_{\theta}}}(\bm{z},y)-y\|^{2}.$
(10)

Similar to CbAS (Paragraph 3), because there is no mention of a validation
set, we need to write out agreement (Equation 10) in terms of a validation
oracle ${\color[rgb]{0,0.5,0}{f_{\phi}}}$. This indirectly leverages the
validation set because the validation oracle is trained on the validation set
(see Figure 2(a)). Furthermore, the agreement should be defined as an expected
value over the $\gamma$-truncated distribution. Therefore, we formally define
agreement as a validation metric as the following:

$\displaystyle\boxed{\mathcal{M}_{\text{Agr}}=\mathbb{E}_{y\sim{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(y)}\|{\color[rgb]{0,0.5,0}{f_{\phi}}}({\color[rgb]{0,0,0.5}{G_{\theta}}}(\bm{z},y)-y\|^{2}}$
(11)

##### Conservative objective models

Trabucco et al. (2021) propose conservative objective models (COMs), which
regularise the trained oracle by generating adversarial examples via gradient
ascent and penalising their predictions. The authors show theoretically that
their technique is able to lower bound the actual ground truth oracle
prediction on out-of-distribution inputs. Computationally, the method is
appealing since it does not require the training of a generative model. While
this method improves oracle prediction by penalising _overestimation_ , it
does not necessarily address the issue of the oracle scoring an ‘impossible’
candidate. However, Kumar et al. (2021) proposes an additional loss to the
COMs objective, where the scores predicted by the oracle on an ‘invalid set’
(a set of invalid candidates) is forced to be zero. We note that in the case
of density estimation there is no need (in principle) to provide negative
examples because of the constraint that
$\int_{\bm{x},y}p(\bm{x},y)d\bm{x}dy=1$.

###  3.1 Likelihood-free metrics

Likelihood-free metrics are used exclusively in the GAN literature because
there is no straightforward way to compute likelihoods for these class of
models. The search for good metrics is still an active topic of research
(Borji, 2022). Common likelihood-free metrics involve measuring some distance
between distributions in some predefined feature space. For instance, for GANs
trained on natural image datasets the _Fréchet Inception Distance_ (FID)
(Heusel et al., 2017) is used to fit Gaussians to both distributions with
respect to the feature space of an InceptionNet classifier trained on
ImageNet. Since the acronym ‘FID’ specifically refers to a particular
ImageNet-based feature extractor, we will simply say ‘FD’, which generalises
to any feature extractor $f_{h}:\mathcal{X}\rightarrow\mathcal{H}$ (Dowson &
Landau, 1982):

$\displaystyle\mathcal{M}_{\text{FD}}(\bm{X},\tilde{\bm{X}},f_{h})=|\mu(f_{h}(\bm{X}))-\mu(f_{h}(\tilde{\bm{X}}))|\
+$
$\displaystyle\text{Tr}(\Sigma(f_{h}(\bm{X}))+\Sigma(f_{h}(\tilde{\bm{X}}))-2\Sigma(f_{h}(\bm{X}))\Sigma(f_{h}(\tilde{\bm{X}}))^{\frac{1}{2}})$
(12)

where $\mathcal{M}_{\text{FD}}\in\mathbb{R}^{+}$ and lower FD is better. Here,
$\bm{X}\in\mathbb{R}^{N\times p}$ denotes $N$ samples coming from a reference
distribution (i.e. ground truth) and $\tilde{\bm{X}}$ are samples coming from
the generative model. Since we also use neural networks to train the
validation oracle ${\color[rgb]{0,0.5,0}{f_{\phi}}}$, the feature extractor
$f_{h}$ can simply be the hidden representation extracted from a predefined
layer for that network.

We consider $\mathcal{M}_{\text{FD}}$ since it is one of the most commonly
used metrics for generative models, and many implementations are readily
available. However, some caveats are that it requires a very large sample
size, and that it may have a stronger bias towards recall (mode coverage) than
precision (sample quality) (Kynkäänniemi et al., 2019).

##  4 Experiments and Discussion

In this section we work towards answering the following question: which
likelihood-free metrics correlate the best with the ground truth oracle? In
particular, we perform our analysis generative adversarial networks and
compare two validation metrics: Agreement ($\mathcal{M}_{\text{Agr}}$,
Equation 11) and FID ($\mathcal{M}_{\text{FD}}$, Equation 3.1).

##### Dataset

Our codebase is built on top of the Design-Bench (Trabucco et al., 2022)
framework. The dataset we consider is TFBind8 (Barrera et al., 2016), where
candidates are length-8 DNA strings whose score $y\in[0,1]$ measures the
binding activity with a particular human transcription factor. This dataset is
convenient because the ground truth oracle ${\color[rgb]{0.67,0.01,0}{f}}$
exists, and this is because the dataset has full support; in other words, each
possible length-8 DNA string ($4^{8}=65536$ many strings) is contained in the
dataset, so the oracle is simply a lookup table. However, this means there is
no notion of there being an ‘invalid’ candidate, because (in the real world)
any arbitrary permutation of nucleotides is possible. To address this, we
consider a modified version of TFBind8 where, based on some pre-defined
pattern, we artificially select certain DNA strings to be invalid. We consider
a simple strategy where if the DNA string contains more than 3 A’s (adenine),
we flag it as invalid. This means that they are _removed_ from dataset, and
the resulting validation oracle ${\color[rgb]{0,0.5,0}{f}}(\bm{x})$ is not
trained on them. The resulting test oracle is modified so that it will score
${\color[rgb]{0.67,0.01,0}{f}}(\bm{x})=0$ for any of these invalid candidates.

##### Oracle pre-training

The oracle ${\color[rgb]{0,0.5,0}{f_{\phi}}}$ is an MLP comprising of one
hidden layer. This was trained on both the training and validation sets (see
Figure 2), minus a small internal validation set used to tune the number of
hidden layers. No invalid examples (Paragraph 4) were used to train this
oracle.

##### Training

We train a conditional GAN which comprises two networks: the generator
$G_{\theta}(\bm{z},y)$ and discriminator $D(\bm{x},y)$, where $\bm{z}\sim
p(\bm{z})=\mathcal{N}(0,\bm{I})$ is the latent code and $y$ is the score. To
condition on $y$ in the generator we use projection layers that map $y$ to a
high dimensional feature vector that in turn is concatenated with each hidden
layer. Since $\bm{x}\in\\{0,1\\}^{4\times 8}$ is discrete (i.e. represented as
eight one-hot vectors of length 4), the output of the generator must also be
discrete. We address this using the Gumbel-Softmax estimator (Jang et al.,
2016). For the discriminator we use the projection formulation from Miyato &
Koyama (2018). The objective function we use is the non-saturating logistic
loss originally proposed in Goodfellow et al. (2014):

$\displaystyle\min_{D}\mathcal{L}_{D}$
$\displaystyle=-\mathbb{E}_{{\color[rgb]{0,0,0.5}{\bm{x},y}}\sim{\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}}\log\big{[}D({\color[rgb]{0,0,0.5}{\bm{x}}},{\color[rgb]{0,0,0.5}{y}})\big{]}-\mathbb{E}_{\bm{z}\sim
p(\bm{z}),{\color[rgb]{0,0,0.5}{y}}\sim{\color[rgb]{0,0,0.5}{{\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}}}(y)}\log\big{[}1-D(G_{\theta}(\bm{z},{\color[rgb]{0,0,0.5}{y}}),{\color[rgb]{0,0,0.5}{y}})\big{]}$
(13) $\displaystyle\min_{G_{\theta}}\mathcal{L}_{G}$
$\displaystyle=-\mathbb{E}_{\bm{z}\sim
p(\bm{z}),{\color[rgb]{0,0,0.5}{y}}\sim{\color[rgb]{0,0,0.5}{{\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}}}(y)}\log\big{[}D(G_{\theta}(\bm{z},{\color[rgb]{0,0,0.5}{y}}),{\color[rgb]{0,0,0.5}{y}})\big{]},$
(14)

Figure 3: $\mathcal{M}_{\text{FD}}$ vs $\mathcal{M}_{\text{Agr}}$

We consider the following hyperparameters: the number of hidden units per
hidden layer (for a fixed number of two hidden layers) as well as the Gumbel
Softmax temperature $\tau$. For Gumbel-Softmax sampling we use the straight-
through estimator. For each model, over the course of training we train two
versions, each one using either $\mathcal{M}_{\text{FD}}$ or
$\mathcal{M}_{\text{Agr}}$ as the early stopping metric (i.e. retain the
checkpoint corresponding to the smallest value of that metric as seen on the
validation set).

##### Generation

For inference time generation we consider a simple procedure: sample a large
number of candidates with a specific augmented model
${\color[rgb]{0,0,0.5}{p_{\theta,\alpha}}}(\bm{x},y)={\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y)p_{\alpha}(y)$
where $p_{\alpha}=\text{Uniform}(\gamma,1)$, score each candidate with the
validation oracle ${\color[rgb]{0,0.5,0}{f_{\phi}}}$, find the top $M=10$
scoring candidates and compute their average score _with respect to the test
oracle_ ${\color[rgb]{0.67,0.01,0}{f}}$.

(a) $\mathcal{M}_{\text{FD}}$ vs mean score

(b) $\mathcal{M}_{\text{Agr}}$ vs mean score

Figure 4: (4(a)) The smallest $\mathcal{M}_{\text{FD}}$ achieved for each
model plotted against the mean score of the top $M=10$ candidates as per the
inference procedure described in Paragraph 4; (4(b)) same as (4(a)) but for
$\mathcal{M}_{\text{Agr}}$. For both metrics $\mathcal{M}_{\text{Agr}}$ and
$\mathcal{M}_{\text{FD}}$ smaller is better. $\mathcal{M}_{\text{FD}}$ is
negatively and significantly correlated with the mean score of both oracles
($R=-0.57$), whereas this is not the case for $\mathcal{M}_{\text{Agr}}$
($R=0.31$). This indicates that $\mathcal{M}_{\text{FD}}$ is a better
predictor of the oracle score.

##### Results

We find that $\mathcal{M}_{\text{Agr}}$ and $\mathcal{M}_{\text{FD}}$ have a
significant negative correlation (Figure 3) of $R=-0.45$. In Figures 4(a) and
4(b) we plot both likelihood-free metrics against the mean score of their
best-generated candidates (as described in Paragraph 4) and find that
$\mathcal{M}_{\text{FD}}$ has a stronger correlation with producing high-
scoring candidates; since its Pearson correlation is negative ($R=-0.5706$),
as $\mathcal{M}_{\text{FD}}$ gets better (smaller), the mean top-10 score
increases. For $\mathcal{M}_{\text{Agr}}$ this is not the case – it is only
moderately correlated with the mean top-10 score in the reverse direction
($R=0.3061$).

(a) Top three models, ranked by best (smallest) $\mathcal{M}_{\text{FD}}$.

(b) Top three models, ranked by best (smallest) $\mathcal{M}_{\text{Agr}}$.

Figure 5: To underscore some of the issues with the agreement metric, here we
show the conditioning $y$ versus the predictions made by the validation and
test oracle, for
$\\{\bm{x}\sim{\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y)\\}_{j=1}^{K}$,
where $K=32$. The blue region denotes $y$’s that were only seen by the GAN
during training ($y$’s from ${\color[rgb]{0.67,0.01,0}{p_{0,\gamma}}}(y)$),
and orange/green region denotes those from the extrapolated regime
${\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(y)$. Orange denotes _test oracle
predictions_ (${\color[rgb]{0.67,0.01,0}{f}}(\bm{x})$) and green denotes
_validation oracle predictions_ (${\color[rgb]{0,0.5,0}{f_{\phi}}}(\bm{x})$).
We can see that the best models corresponding to $\mathcal{M}_{\text{FD}}$
(5(a)) exhibit monotonically increasing predictions with respect to the
conditioning $y$, and this is not the case for $\mathcal{M}_{\text{Agr}}$
(5(b)), whose predictions are also erratic even in the training region of $y$.
Furthermore, $\mathcal{M}_{\text{FD}}$-based models are able to achieve test
oracle scores very close to the maximum possible attainable value ($y=1.0$).

In Figure 5 we plot the conditioning $y\in\text{linspace}(0,1)$ against the
validation/test oracle predictions for candidates conditionally generated with
that $y$. We call these ‘agreement plots’555Also called ‘paired plots’ in
Brookes et al. (2019).since the sum of squared residuals for each point would
constitute the agreement (with a perfect agreement of zero corresponding to a
diagonal dotted line on each graph). In Figure 5(a) we demonstrate this
amongst the best three models _with respect to_ $\mathcal{M}_{\text{FD}}$, and
in Figure 5(b) we do the same with respect to $\mathcal{M}_{\text{Agr}}$. The
shaded regions denote $\pm 1$ standard deviation from the mean, and the marker
symbols denote the max/worst score for each $y$. We can see that the
predictions shown in Figure 5(a) are more or less monotonically increasing as
a function of $y$, whereas the models with the best agreement exhibit very
erratic or undesirable behaviour (Figure 5(b)).

Even though the mean prediction for any particular $y>\gamma$ lies closer to
the dotted line for the experiments ranked by $\mathcal{M}_{\text{Agr}}$
(Figure 5(b)), the maximum scores level out at a certain score (e.g a little
below 0.8 for the left-most model in Figure 5(b)) whereas for
$\mathcal{M}_{\text{FD}}$ the maximum predictions lie close to 1.0 with
respect to either oracle.

###  4.1 Limitations and Future Work

One dilemma we are faced with in this proposed evaluation setup is the
requirement of the ground truth oracle. As mentioned in the introduction, the
ground truth oracle is expensive to compute and offline MBO assumes we do not
have access to it during training and validation. However, the particular
questions we would like to answer require it so that we can determine with
confidence which validation metrics correlate the best with the ground truth.
Unfortunately, there appear to be very few datasets that contain invalid
candidates _and_ also admit access to the ground truth. This leaves us with
two options: modify ground-truth-containing datasets to support invalid
candidates (like we did with TFBind); design challenging simulation-based
datasets that reflect real world problems; or simply allow for the ground
truth oracle to be replaced with an approximate ‘test’ oracle
${\color[rgb]{0.67,0.01,0}{f_{\psi}}}$. The latter option however is _risky_
because we are now evaluating one approximation
(${\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x}|y)$) with another
(${\color[rgb]{0.67,0.01,0}{f_{\psi}}}(\bm{x})$). One compromise could be to
consider a relaxed scenario where we only assume that ground truth
_verification_ is possible. For instance, consider a simple binary function
${\color[rgb]{0.67,0.01,0}{v}}(\bm{x}):\mathcal{X}\rightarrow\\{0,1\\}$ where
‘1’ denotes ‘valid’ and ‘0‘ denotes ‘invalid’. This verifier could be ‘folded’
into the approximate oracle to produce a new oracle
${\color[rgb]{0.67,0.01,0}{f_{\psi}}}\cdot{\color[rgb]{0.67,0.01,0}{v}}$,
which in turn would assign any invalid candidates a score of zero and
therefore make it more reliable, even if it is still prone to over-scoring bad
candidates that are invalid. One such real-world example is chemistry
retrosynthesis software, which can be used to determine whether a compound is
plausible to synthesise (Liu et al., 2020). This compromise however would
require a slight change to the evaluation pipeline to accommodate a test set
(see Figure 6(b)).

##  5 Conclusion

In this work we propose a principled evaluation framework to measure how well
generative model-based MBO models extrapolate. We would like to leverage this
framework to determine which validation metrics correlate the best with the
ground truth, the latter of which is expensive to compute in practice. The
proposed framework is practical yet simple in the sense that extrapolating a
generative model simply amounts to changing the prior distribution over the
scores. Furthermore, we discussed how existing generative model-based works
can easily be ported to our framework. However, determining what validation
metric correlates most with extrapolative ability requires the ground truth
oracle if we want to be absolutely certain. This is because approximate test
oracles are vulernable to over-scoring bad candidates or even worse, scoring
candidates that are invalid (i.e. have zero density under the ground truth
distribution). As a proof of concept, we compared two validation metrics with
respect to a modified version of TFBind, which contains a ground truth oracle
but is modified to allow for invalid candidates.

Lastly, since finding good validation metrics requires the ground truth oracle
and limits what datasets can be used, we propose a compromise to our
evaluation framework which only assumes that ground truth verification is
possible. This would guard the approximate test oracle against scoring invalid
candidates. This will be considered in a future iteration of this preprint.

## References

  * Asmussen & Glynn (2007) Asmussen, Søren and Glynn, Peter W.  _Stochastic simulation: algorithms and analysis_ , volume 57.  Springer, 2007.
  * Barrera et al. (2016) Barrera, Luis A, Vedenko, Anastasia, Kurland, Jesse V, Rogers, Julia M, Gisselbrecht, Stephen S, Rossin, Elizabeth J, Woodard, Jaie, Mariani, Luca, Kock, Kian Hong, Inukai, Sachi, et al.  Survey of variation in human transcription factors reveals prevalent DNA binding changes.  _Science_ , 351(6280):1450–1454, 2016.
  * Borji (2022) Borji, Ali.  Pros and cons of GAN evaluation measures: New developments.  _Computer Vision and Image Understanding_ , 215:103329, 2022\.
  * Brookes et al. (2019) Brookes, David, Park, Hahnbeom, and Listgarten, Jennifer.  Conditioning by adaptive sampling for robust design.  In _International conference on machine learning_ , pp. 773–782. PMLR, 2019.
  * Dowson & Landau (1982) Dowson, DC and Landau, BV666017.  The Fréchet distance between multivariate normal distributions.  _Journal of multivariate analysis_ , 12(3):450–455, 1982.
  * Goodfellow et al. (2014) Goodfellow, Ian, Pouget-Abadie, Jean, Mirza, Mehdi, Xu, Bing, Warde-Farley, David, Ozair, Sherjil, Courville, Aaron, and Bengio, Yoshua.  Generative adversarial nets.  _Advances in neural information processing systems_ , 27, 2014.
  * Heusel et al. (2017) Heusel, Martin, Ramsauer, Hubert, Unterthiner, Thomas, Nessler, Bernhard, Klambauer, Günter, and Hochreiter, Sepp.  GANs trained by a two time-scale update rule converge to a nash equilibrium.  _CoRR_ , abs/1706.08500, 2017.  URL http://arxiv.org/abs/1706.08500.
  * Ho et al. (2020) Ho, Jonathan, Jain, Ajay, and Abbeel, Pieter.  Denoising diffusion probabilistic models.  _Advances in Neural Information Processing Systems_ , 33:6840–6851, 2020.
  * Jang et al. (2016) Jang, Eric, Gu, Shixiang, and Poole, Ben.  Categorical reparameterization with Gumbel-Softmax.  _arXiv preprint arXiv:1611.01144_ , 2016.
  * Kingma & Welling (2013) Kingma, Diederik P and Welling, Max.  Auto-encoding variational Bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Kumar & Levine (2020) Kumar, Aviral and Levine, Sergey.  Model inversion networks for model-based optimization.  _Advances in Neural Information Processing Systems_ , 33:5126–5137, 2020.
  * Kumar et al. (2021) Kumar, Aviral, Yazdanbakhsh, Amir, Hashemi, Milad, Swersky, Kevin, and Levine, Sergey.  Data-driven offline optimization for architecting hardware accelerators.  _arXiv preprint arXiv:2110.11346_ , 2021.
  * Kynkäänniemi et al. (2019) Kynkäänniemi, Tuomas, Karras, Tero, Laine, Samuli, Lehtinen, Jaakko, and Aila, Timo.  Improved precision and recall metric for assessing generative models.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Liu et al. (2020) Liu, Cheng-Hao, Korablyov, Maksym, Jastrzębski, Stanisław, Włodarczyk-Pruszyński, Paweł, Bengio, Yoshua, and Segler, Marwin HS.  Retrognn: Approximating retrosynthesis by graph neural networks for de novo drug design.  _arXiv preprint arXiv:2011.13042_ , 2020.
  * Lopez-Paz & Oquab (2016) Lopez-Paz, David and Oquab, Maxime.  Revisiting classifier two-sample tests.  _arXiv preprint arXiv:1610.06545_ , 2016.
  * Miyato & Koyama (2018) Miyato, Takeru and Koyama, Masanori.  cGANs with projection discriminator.  _arXiv preprint arXiv:1802.05637_ , 2018.
  * O’Donoghue et al. (2020) O’Donoghue, Brendan, Osband, Ian, and Ionescu, Catalin.  Making sense of reinforcement learning and probabilistic inference.  _arXiv preprint arXiv:2001.00805_ , 2020.
  * Song et al. (2020) Song, Yang, Sohl-Dickstein, Jascha, Kingma, Diederik P, Kumar, Abhishek, Ermon, Stefano, and Poole, Ben.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020.
  * Trabucco et al. (2021) Trabucco, Brandon, Kumar, Aviral, Geng, Xinyang, and Levine, Sergey.  Conservative objective models for effective offline model-based optimization, 2021.
  * Trabucco et al. (2022) Trabucco, Brandon, Geng, Xinyang, Kumar, Aviral, and Levine, Sergey.  Design-Bench: Benchmarks for data-driven offline model-based optimization, 2022.

##  Appendix A Appendix

###  A.1 Acknowledgements

Figure 1 uses clipart from
here666https://lovepik.com/image-401281091/chemical-test-tube.html. The
protein render was lifted from
here777https://www.ebi.ac.uk/chembl/compound_report_card/CHEMBL159/.

###  A.2 Full proof

$\displaystyle\theta^{*}$
$\displaystyle=\operatorname*{arg\,min}_{\theta\in\Theta}\
\underbrace{\text{KL}\Big{[}{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)\
\|\ {\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x},y)\Big{]}}_{\text{forward
KL}}$ $\displaystyle=\operatorname*{arg\,min}_{\theta\in\Theta}\
\mathbb{E}_{{\color[rgb]{0.67,0.01,0}{\bm{x},y}}\sim{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)}\log{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}({\color[rgb]{0.67,0.01,0}{\bm{x}}},{\color[rgb]{0.67,0.01,0}{y}})-\mathbb{E}_{{\color[rgb]{0.67,0.01,0}{\bm{x},y}}\sim{\color[rgb]{0.67,0.01,0}{p}}(\bm{x})}\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}({\color[rgb]{0.67,0.01,0}{\bm{x}}},{\color[rgb]{0.67,0.01,0}{y}})$
$\displaystyle=\operatorname*{arg\,max}_{\theta\in\Theta}\
-\mathbb{E}_{{\color[rgb]{0.67,0.01,0}{\bm{x},y}}\sim{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)}\log{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}({\color[rgb]{0.67,0.01,0}{\bm{x}}},{\color[rgb]{0.67,0.01,0}{y}})+\mathbb{E}_{{\color[rgb]{0.67,0.01,0}{\bm{x},y}}\sim{\color[rgb]{0.67,0.01,0}{p}}(\bm{x})}\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}({\color[rgb]{0.67,0.01,0}{\bm{x}}},{\color[rgb]{0.67,0.01,0}{y}})$
$\displaystyle=\operatorname*{arg\,max}_{\theta\in\Theta}\
\underbrace{\mathbb{H}[{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}]}_{\text{constant}}+\mathbb{E}_{{\color[rgb]{0.67,0.01,0}{\bm{x},y}}\sim{\color[rgb]{0.67,0.01,0}{p_{\gamma}}}(\bm{x},y)}\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}({\color[rgb]{0.67,0.01,0}{\bm{x}}},{\color[rgb]{0.67,0.01,0}{y}})$
$\displaystyle\approx\operatorname*{arg\,max}_{\theta\in\Theta}\
\frac{1}{N}\sum_{i=1}^{N}\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}({\color[rgb]{0,0.5,0}{\bm{x}_{i}}},{\color[rgb]{0,0.5,0}{y_{i}}}),\
\ \ \ \
\\{{\color[rgb]{0,0.5,0}{\bm{x}_{i}}},{\color[rgb]{0,0.5,0}{y_{i}}}\\}_{i=1}^{N}={\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}.$
(15)

###  A.3 Relaxed evaluation procedure

See Figure S6.

(a)

(b)

Figure S6: Under our relaxed evaluation scenario, we no longer assume access
to the ground truth oracle ${\color[rgb]{0.67,0.01,0}{f}}$, but we do assume
access to a ground truth ‘verifier’ ${\color[rgb]{0.67,0.01,0}{v}}(\bm{x})$
which is able to classify whether a candidate is valid or not. This can be
used to induce a more reliable test oracle
${\color[rgb]{0.67,0.01,0}{f}}\approx{\color[rgb]{0.67,0.01,0}{f_{\psi}}}\cdot{\color[rgb]{0.67,0.01,0}{v}}$.
The approximate oracle ${\color[rgb]{0.67,0.01,0}{f_{\psi}}}$ is trained on an
additional split ${\color[rgb]{0.67,0.01,0}{\mathcal{D}_{\text{test}}}}$,
which comes from the same truncated distribution as
${\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}$. The use of a test set
here is necessary since we do not wish to perform model selection on the exact
same dataset that was used to train the approximate test oracle, or else we
may risk reporting a biased estimate of the _true_ oracle performance
${\color[rgb]{0.67,0.01,0}{f}}$. In other words, only
${\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}$ is used for model
selection, while ${\color[rgb]{0.67,0.01,0}{\mathcal{D}_{\text{test}}}}$ is
used to train ${\color[rgb]{0.67,0.01,0}{f_{\psi}}}$. (Note that in practice,
one should use as much data for the test oracle, so we should train it on
${\color[rgb]{0,0,0.5}{\mathcal{D}_{\text{train}}}}$,
${\color[rgb]{0,0.5,0}{\mathcal{D}_{\text{valid}}}}$, and
${\color[rgb]{0.67,0.01,0}{\mathcal{D}_{\text{test}}}}$. This is not shown in
the figure for the sake of readability.)

###  A.4 Model deployment

Once the best model has been found via an appropriate validation metric, one
can train the same type of model on the full dataset $\mathcal{D}$ using the
same hyperparameters as before. Ultimately, we would like to be able to
generate candidates whose $y$’s exceed that of the entire dataset, and at the
same time at plausible according to our generative model. To control how much
we trade-off high likelihood versus high scoring candidates, we can consider
the _exponentially tilted density_ (Asmussen & Glynn, 2007; O’Donoghue et al.,
2020):

$\displaystyle{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x},y)\exp(\eta^{-1}y-\kappa(\eta)),$
(16)

where $\kappa(\eta)$ is a normalisation constant, and smaller $\eta$ puts
larger emphasis on sampling from regions where $y$ is large. Taking the log of
Equation 16, we arrive at:

$\displaystyle\bm{x}^{*},y^{*}=\operatorname*{arg\,max}_{\bm{x},y}\
\log{\color[rgb]{0,0,0.5}{p_{\theta,\gamma}}}(\bm{x},y)+\frac{1}{\eta}y,$ (17)

In practice, it would not be clear what the best $\eta$ should be, but a
reasonable strategy is to consider a range of $\eta$’s, where larger values
encode a higher tolerance for ‘risk’ since these values favour higher scoring
candidates at the cost of likelihood. Note that for VAEs and diffusion models,
$\log p_{\theta}(\bm{x},y)$ will need to be approximated with the ELBO. For
score-based models (Song et al., 2020) only $\nabla_{\bm{x},y}\log
p(\bm{x},y)$ is accessible, but sampling is still straightforward since that
is just the gradient.

###  A.5 Related work

####  A.5.1 CbAS

Figure S7: The training and evaluation of CbAS Brookes et al. (2019) in the
context of our evaluation framework. The extrapolation equation is described
in Equation 3 and involves variational inference to fine-tune
${\color[rgb]{0,0,0.5}{p_{\theta}}}(\bm{x})$ into a search model
${\color[rgb]{0,0,0.5}{q_{\zeta}}}(\bm{x})$. For model selection, the
validation metric is Equation 8, which is just Equation 3 but using the
validation oracle ${\color[rgb]{0,0.5,0}{f_{\phi}}}(\bm{x})$.


# Frustration induced Itinerant Ferromagnetism of Fermions in Optical Lattices

Chengshu Li  Institute for Advanced Study, Tsinghua University, Beijing
100084, China  Ming-Gen He  Hefei National Research Center for Physical
Sciences at the Microscale and School of Physics, University of Science and
Technology of China, Hefei 230026, China  CAS Center for Excellence in Quantum
Information and Quantum Physics, University of Science and Technology of
China, Hefei 230026, China  Chang-Yan Wang  Institute for Advanced Study,
Tsinghua University, Beijing 100084, China  Hui Zhai<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua University, Beijing 100084, China

###### Abstract

When the Fermi Hubbard model was first introduced sixty years ago, one of the
original motivations was to understand correlation effects in itinerant
ferromagnetism. In the past two decades, ultracold Fermi gas in an optical
lattice has been used to study the Fermi Hubbard model. However, the metallic
ferromagnetic correlation was observed only in a recent experiment using
frustrated lattices, and its underlying mechanism is not clear yet. In this
letter, we point out that, under the particle–hole transformation, the single-
particle ground state can exhibit double degeneracy in such a frustrated
lattice. Therefore, the low-energy state exhibits valley degeneracy,
reminiscent of multi-orbit physics in ferromagnetic transition metals. The
local repulsive interaction leads to the valley Hund’s rule, responsible for
the observed ferromagnetism. We generalize this mechanism to distorted
honeycomb lattices and square lattices with flux. This mechanism was first
discussed by Müller-Hartmann in a simpler one-dimension model. However, this
mechanism has not been widely discussed and has not been related to
experimental observations before. Hence, our study not only explains the
experimental findings but also enriches our understanding of itinerant
ferromagnetism.

Itinerant ferromagnetism has been discovered in nature for thousands of years.
However, a complete understanding of its microscopic origin is still
challenging. The Stoner mean-field theory has predicted itinerant
ferromagnetism in metal when the repulsive interaction between fermions
exceeds a critical value Stoner1938 . When such ferromagnetism occurs, the
Pauli exclusion principle increases kinetic energy considerably. Therefore,
the critical interaction strength predicted for Stoner ferromagnetism is
comparable to the Fermi energy. Under such a strong interaction strength, the
correlation effect can no longer be ignored. The competition from other
strongly correlated non-magnetic states usually takes over Stoner
ferromagnetism. Therefore, understanding itinerant ferromagnetism becomes
essentially a strongly correlated problem.

When the Hubbard model was introduced in the middle of the last century
Hubbard1963 ; Gutzwiller1963 ; Kanamori1963 , one major purpose was to
understand the correlation effect in ferromagnetism. Unfortunately, the
consensus on ferromagnetism in the Hubbard model is still limited after many
decades Tasaki1998 ; Tasaki1998b ; Arovas2022 ; Qin2022 . Rigorous results can
only be obtained for exceptional cases such as the Nagaoka ferromagnetism
Thouless1965 ; Nagaoka1966 ; Tasaki1989  and the flat band ferromagnetism
Mielke1991 . The Nagaoka ferromagnetism considers a single hole doping away
from half-filling in the limit of infinite repulsion, and the flat band
ferromagnetism requires fine-tuning of particle hopping to reach a flat band
dispersion. Moreover, the single-band Hubbard model is usually oversimplified
to directly compare with experiments on real materials.

Ultracold atoms in optical lattices directly realize the Hubbard model and
provide a new opportunity to study physics therein Bloch2005 ; Bloch2008 .
However, the temperature of atoms in optical lattices cannot be cooled much
below the kinetic energy, which prevents observing possible low-temperature
orderings, such as fermion pairing, at this stage. In contrast, ferromagnetism
often occurs in a relatively high temperature. It is conceivable that one does
not need to enter an extremely low-temperature regime in order to study
ferromagnetism in optical lattices. Hence, understanding the nature of
ferromagnetism is a potential goal for the current quantum simulations with
optical lattices.

Nevertheless, although the Fermi Hubbard model has been realized with
ultracold atoms in optical lattices for nearly two decades, itinerant
ferromagnetism has not been observed in this system until a very recent
experiment Xu2022 . In this experiment, a novel experimental technology allows
continuously tuning the lattice geometry from a square to a frustrated
triangular lattice. Short-range ferromagnetic correlation has been observed in
the particle doping regime when the lattice geometry is tuned close to
triangular regime. However, a convincing theoretical understanding of the
physical mechanism behind the observed ferromagnetism is still lacking.

Figure 1: (a) The tunable optical lattice realized in a recent experiment. (b)
The ferromagnetic correlation found in the experiment is marked by the shaded
yellow regime in the $n-t^{\prime}$ phase diagram, which is mapped to the
shaded blue area under the particle–hole transformation. (c) Honeycomb lattice
with the nearest neighbor and the horizontal next nearest neighbor hopping.
(d) Square lattice with the nearest neighbor hopping only, but with magnetic
flux in each plaquette. The red lines in (a) and (c) denote the one-
dimensional chain considered by Müller-Hartmann’s paper.

Review of Experimental Setting. This experiment reports an actively phase
stabilized optical lattice that can continuously tune the lattice geometry
from square lattice to triangular lattice. The tight-binding model of this
lattice is shown in Fig. 1(a). This experiment explores spin-$1/2$ fermions in
such lattice, and the model is written as

$\hat{H}=-t\sum\limits_{\langle
ij\rangle\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}-t^{\prime}\sum\limits_{\langle\langle
ij\neswarrow\rangle\rangle}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}+U\sum\limits_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$
(1)

where $\sigma=\uparrow,\downarrow$ denotes two spin components. The hopping
between the nearest neighbor $\langle ij\rangle$ is denoted by $t$, and
between the next nearest neighbor along the dashed line direction
$\langle\langle ij\neswarrow\rangle\rangle$ is denoted by $t^{\prime}$. Both
$t$ and $t^{\prime}$ are positive. The next nearest hopping along another
diagonal direction is negligible. When $t^{\prime}$ increases from zero to
$t^{\prime}=t$, it continuously tunes the lattice geometry from a square
lattice to a triangular lattice. $U$ denotes the interaction strength of on-
site repulsion, and $U/t\approx 9$ in this experiment.

This experiment finds ferromagnetic correlation when $t^{\prime}/t>0.5$ and
when total fermion density $n$ exceeds half-filling $n=1$ and is somewhat
close to $n=1.5$. The regime where ferromagnetism is observed is marked by
shaded yellow area in Fig. 1(b). This experiment has only explored the density
regime $0.5<n<1.5$ and it is not clear whether the ferromagnetic correlation
can also exist when $n$ exceeds $1.5$. The physical origin of this observed
ferromagnetism is also not clear yet. Possible scenarios mentioned include the
existence of van Hove singularity at $n=1.5$ for triangular lattice and the
possible connection to the Nagaoka ferromagnetism Xu2022 ; Merino2006 ;
Demler2022 .

Particle–Hole Transformation. For the benefit of later discussion, we make a
particle–hole transformation
$\hat{c}_{i\sigma}\rightarrow(-1)^{i_{x}+i_{y}}\hat{c}^{\dagger}_{i\sigma}$,
where $i=(i_{x},i_{y})$ is the site label. If $t^{\prime}=0$, this
transformation keeps the form of Hamiltonian Eq. (1) invariant.
$t^{\prime}$-term is what causes frustration in this lattice, and for the same
reason, $t^{\prime}$-term is not invariant under the particle–hole
transformation. Instead, this transformation changes $t^{\prime}$ to
$-t^{\prime}$. This corresponds to inserting a $\pi$-flux in each triangle.
For positive $t^{\prime}$, the system essentially describes particles moving
in a real potential and all hopping terms have negative matrix elements that
satisfies Feynman’s no-node theorem Feynman1998 . Therefore, the single-
particle ground state cannot have degeneracy. However, for negative
$t^{\prime}$, the condition for no-node theorem is no longer satisfied, and
the single-particle ground state can have degeneracy. It turns out that this
degeneracy is crucial for explaining this observed ferromagnetism, as it will
become clear later.

Under this transformation, $n=n_{\uparrow}+n_{\downarrow}$ becomes $2-n$, that
is, the particle doping is mapped to hole doping. Thus, the particle–hole
transformation maps $(t^{\prime},n)$ to $(-t^{\prime},2-n)$ in the phase
diagram shown in Fig. 1(b). Therefore, we can focus solely on the density
regime $0<n\leqslant 1$ but include both positive and negative $t^{\prime}$.
With this mapping, the regime where ferromagnetism is observed is mapped to
the low-density regime with $t^{\prime}<-0.5$, as marked as shaded blue area
in Fig. 1(b).

Numerical Results. We first present our numerical results in Fig. 2 code . All
the calculations below are done for infinite positive $U$. The first
calculation is exact diagonalization for two particles with different system
sizes, reminiscent of different densities. This calculation can cover the low-
density regime up to $n=0.5$. Due to the $SU(2)$ spin rotational symmetry, the
total spin is a good quantum number and all quantum states with the same total
spin are degenerate. We find a level crossing between spin singlet and spin
triplet states, as shown in Fig. 2(a). The spin triplet states have a lower
energy when $t^{\prime}<t^{\prime}_{\text{c}}$, and the value of
$t^{\prime}_{\text{c}}$ depends on density and lies between $-1.0$ and $-0.5$,
as one can see from circles in Fig. 2(c).

Figure 2: (a) Two-body problem in a $4\times 4$ lattice with open boundary
condition. The singlet and triplet state energies as a function of
$t^{\prime}$. (b) $S_{\text{tot}}$ of the ground state calculated by the DMRG
method for $N=12$ particles in a $4\times 8$ strip. (c) The shaded area
schematically shows the possible ferromagnetic regime in the $n-t^{\prime}$
phase diagram. Each point denote $t^{\prime}_{\text{c}}$ between ferromagnetic
and non-magnetic phase at different densities. The numerical calculations
include exact diagonalizations of two-particle problem in different system
sizes and DMRG calculations of finite number of particles on different system
sizes. All DMRG calculations have bond dimension $1400$. In all cases, $t=1$
and $U$ is set to positive infinite.

The second calculation is the density-matrix renormalization group (DMRG)
calculation with finite number of fermions on different strip geometry. We can
calculate $\langle\hat{{\bf S}}_{\text{tot}}^{2}\rangle$ with $\hat{{\bf
S}}_{\text{tot}}=\sum_{i}\hat{{\bf S}}_{i}$ for the ground state, and
$S_{\text{tot}}$ is given by $\langle\hat{{\bf
S}}_{\text{tot}}^{2}\rangle=S_{\text{tot}}(S_{\text{tot}}+1)$. We also find a
transition from $S_{\text{tot}}=N/2$ to $S_{\text{tot}}=0$ around
$t^{\prime}=t^{\prime}_{\text{c}}$, as shown in Fig. 2(b). Fig. 2(c) also
collects $t^{\prime}_{\text{c}}$ obtained by all DMRG calculations with
different fermion number and different system sizes.

A notable feature in Fig. 2(c) is that for all calculations,
$t^{\prime}_{\text{c}}$ approaches $-0.5$ at the low-density limit when
$n\rightarrow 0$. Below we will explain why $t^{\prime}=-0.5$ is a special
point, and the ferromagnetism emerged in this regime will be called the
Müller-Hartmann Mechanism. In this regime, the general trend is that
$t^{\prime}_{\text{c}}$ decreases as $n$ increases.

Figure 3: (a) Single-particle dispersion $\mathcal{E}(k_{x},k_{y})$ of the
Hamiltonian Eq. (1) for three different values of $t^{\prime}$.
$t^{\prime}=-0.2$ (left), $=-0.5$ (middle) and $=-0.8$ (right). (b) The real
and imaginary parts of two valley Wannier wave functions. These two wave
functions are constructed by using the Bloch wave functions around each
minimum of single-particle dispersion, as indicated by the solid circle in the
right dispersion in (a).

Nearby half-filling with $n=1$, it is known that hole doping can result in
Nagaoka ferromagnetism at infinite $U$. Strictly speaking, Nagaoka
ferromagnetism can only be proved for single-hole doping with $t^{\prime}<0$.
However, our numerical results show that when $t^{\prime}=0$, the Nagaoka
ferromagnetism can exist up to $\sim 20\%$ of hole doping, that is, for
$0.8\lesssim n<1$. This is consistent with previous numerical results Liu2012
. Hence, we attribute the ferromagnetism nearby half-filling as Nagaoka
Mechanism. In this regime, $t^{\prime}_{\text{c}}$ also decreases when $n$
decreases from $n=1$, and this density dependence of $t^{\prime}_{\text{c}}$
is opposite to that found in the low-density regime.

Hence, we have identified two different mechanisms of ferromagnetism from our
numerical results, and they emerge from the low-density limit and nearby half-
filling, respectively. They are characterized by different density dependence
of $t^{\prime}_{\text{c}}$. The trend of experimental observed ferromagnetic
boundary, upon the particle–hole transformation, is consistent with the
Müller-Hartmann mechanism.

Figure 4: (a) The ferromagnetic regime indicated by DMRG calculation in the
$n-t^{\prime}$ phase diagram of the distorted honeycomb lattice. The model is
shown in Fig. 1(c), in which $t^{\prime}$ denotes the next nearest hopping
along the horizontal dashed line. The DMRG calculation is performed in a
$3\times 6$ zigzag cylinder, with bond dimension $\chi=1400$. The dashed line
indicates the value of $t^{\prime}$ for the Lifshitz transition. (b) The
single-particle dispersion of the distorted honeycomb lattice, with
$t^{\prime}=-0.1$ (left), $=-0.25$ (middle) and $=-0.4$ (right).

Müller-Hartmann Mechanism. To see why $t^{\prime}=-0.5$ is special, we first
look at the single-particle dispersion which reads

$\displaystyle\mathcal{E}(k_{x},k_{y})$
$\displaystyle=-2t\cos(k_{x})-2t\cos(k_{y})-2t^{\prime}\cos(k_{x}+k_{y})$
$\displaystyle=-4t\cos(k_{+})\cos(k_{-})-2t^{\prime}\cos(2k_{+}),$ (2)

where $k_{\pm}=(k_{x}\pm k_{y})/2$. It is easy to see that the dispersion
minimum occurs at $k_{-}=0$. As shown in Fig. 3(a), when $t^{\prime}/t>-0.5$,
the band dispersion has a unique minimum at $k_{+}=0$. However, when
$t^{\prime}/t<-0.5$, the band dispersion displays two degenerate minima along
$k_{+}$ axes. That is to say, a Lifshitz transition occurs when
$t^{\prime}/t=-0.5$.

The paper by Müller-Hartmann first pointed out that ferromagnetism can occur
in the low-density regime when the band minima display double degeneracy
Muller1995 , followed by few related works in later literatures Pieri1996 .
This paper considered a one-dimensional chain with next nearest hopping, and
the Hamiltonian reads

$\hat{H}=-\sum\limits_{i\sigma}(t\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{i+1,\sigma}+t^{\prime}\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{i+2,\sigma}+\text{h.c.})+U\sum\limits_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}.$
(3)

In this model, the band minima also displays double degeneracy when
$t^{\prime}/t<-1/4$, where Müller-Hartmann argued that ferromagnetism can
occur in the low-density regime. We note that this one-dimensional chain can
be naturally embedded into the two-dimensional lattice we considered, as shown
in Fig. 1(a). In other word, the two-dimensional lattice realized in this cold
atom experiment can be viewed as a two-dimensional generalization of Müller-
Hartmann’s work.

The physical mechanism behind this ferromagnetism can be called the Valley
Hund’s rule valley . Considering the two-particle case, they can be
respectively placed in each minimum of single-particle dispersion. For spin-
triplet states, the spatial wave function is anti-symmetrized, and the
interaction energy automatically vanishes. However, for the spin-singlet
state, one has to apply extra projection to prevent double occupation, which
inevitably increases the kinetic energy.

Using the low-energy Bloch wave functions around each minimum, we can
construct the effective valley Wannier wave functions for each valley, as
shown in Fig. 3(b). These valley Wannier wave functions extend over several
lattice sites. The on-site repulsion between fermions effectively introduces
Hund’s rule between these valley Wannier wave functions. In real materials,
itinerant ferromagnetism usually occurs in transition metals with partially
filled $d$-orbit, and Hund’s rule coupling between these $d$-orbits is crucial
for ferromagnetism. That is to say, the multi-orbital physics plays an
essential role. Müller-Hartmann’s mechanism says that even for a single-band
model, when the single-particle ground state has degeneracy, the valley
degeneracy can lead to an emergent multi-orbital physics, resulting in
itinerant ferromagnetism at the low-density regime. Since Hund’s rule is a
local ferromagnetic coupling, we expect that the short-range ferromagnetic
correlation caused by this mechanism is robust at finite temperatures.
Moreover, since the Hund’s rule coupling is the leading order effect of local
repulsion while the anti-ferromagnetic super-exchange interaction is a second-
order effect of local repulsion, we also expect that this ferromagnetism is
stable with finite interaction strength.

Figure 5: (a) $S_{\text{tot}}$ obtained by DMRG calculation.
$S_{\text{max}}=N/2$. Here we consider square lattice with nearest hopping
only, but with magnetic flux $\phi$. $\phi=0,\pi,2\pi/3,\pi/2$ and $2\pi/5$
from the top to the bottom. The DMRG calculation is performed in a $4\times
10$ strip with bond dimension $\chi=1400$. (b) The single-particle dispersion
around the bottom of the band for different magnetic flux corresponding to
(a).

Generalization to Honeycomb Lattice. To further confirm our understanding of
emergent ferromagnetism in this experiment, we generalize this idea to other
models, where degenerate ground states can also be found. The first
generalization is fermions in a distorted honeycomb lattice. Similar to the
square lattice model, the nearest neighbor hopping is denoted by $t$. We
compress the lattice along the horizontal direction such that we should also
include the next nearest hopping $t^{\prime}$ along the horizontal dash line
direction, as shown Fig. 1(c). This lattice can also be viewed as a set of
$t-t^{\prime}$ chains along $\hat{x}$ direction considered by Müller-
Hartmann’s paper, coupled vertically along $\hat{y}$ direction. We introduce a
similar particle-hole transformation that only changes sign of $t^{\prime}$.
The ground state of this dispersion also shows a Lifshitz transformation from
single to double degeneracy at $t^{\prime}=-0.25$, as shown in Fig. 4(b). The
phase diagram obtained by DMRG calculation is shown in Fig. 4(a). This phase
diagram contains two ferromagnetic regimes, one in the low-density regime
around $t^{\prime}=-0.25$ and the other nearby half-filling for
$t^{\prime}<0$. Clearly, they respectively manifest Müller-Hartmann and
Nagaoka mechanisms Bobrow2018 . Unlike the square lattice case, these two
ferromagnetism mechanisms occur in two disconnected regimes in the phase
diagram, clearly visualizing their difference.

Generalization to Flux Lattice. The second generalization is square lattice
with a magnetic flux in each plaquette. Here we only consider the nearest
hopping in square lattice, and the magnetic flux only introduces a phase in
the hopping and does not cause Zeeman splitting. However, with magnetic flux
$\phi=2\pi/m$ in each plaquette ($m$ is an integer), translation of one
lattice spacing along $\hat{x}$ direction does not commute with one lattice
space translation along $\hat{y}$ direction, and only commutes with
translation of $m$ lattice spacing along $\hat{y}$ direction. This enlarges
the unit cell and leads to $m$-fold degeneracy of dispersion spectrum, as
shown in Fig. 5(b). For a fixed magnetic flux, the only tunable parameter is
density. In Fig. 5(a), we show the total magnetization as a function of
fermion density. It is interesting to note that when $m=1$, the single-
particle dispersion is not be degenerate but the hopping satisfies the
condition for Nagaoka ferromagnetism. Hence, DMRG finds ferromagnetism nearby
half-filling. In contrast, when $m>1$, hopping terms no longer satisfy the
condition for Nagaoka ferromagnetism, but the ground state degeneracy appears.
We find a critical $n_{\text{c}}$ and the system is ferromagnetism when
$n<n_{\text{c}}$. We find that $n_{\text{c}}$ is approximately given by $1/m$,
which corresponds to a half-filled lowest band. Therefore, ferromagnetism
occurs in low-density regime m2 . This clearly shows Nagaoka and Müller-
Hartmann are two different mechanisms.

Conclusion and Discussions. In summary, we analyze a recent experiment that
first discovers itinerant ferromagnetism in ultracold atom realization of the
Fermi Hubbard model. We attribute the mechanism of this ferromagnetism to the
appearance of single-particle ground-state degeneracy. At low density, the
valley degeneracy is reminiscent of multi-orbital physics in transition
metals, and the local repulsive interaction can lead to a valley Hund’s rule,
resulting in ferromagnetic correlations. This mechanism was first discussed by
Müller-Hartmann in a simpler one-dimension model.

We remark that by talking about ferromagnetism, we implicitly assume that the
Hamiltonian of the system should not break time-reversal symmetry. In the
presence of time-reversal symmetry, the single-particle ground state should
not degenerate for normal kinetic energy term due to Feynman’s no-node
theorem. In order for the Müller-Hartmann mechanism to apply, the key is
frustration plus particle-hole transformation. Because of the lattice
frustration, we can effectively inset flux to the hopping terms after a
particle–hole transformation, such that the condition for Feynman’s theorem is
violated. That is why frustration plays an essential role here. Previously,
frustration has been considered an essential ingredient for the emergence of
exotic phases such as spin liquids Balents2010 ; the discussion here shows
that frustration can also be a key ingredient responsible for the emergence of
ferromagnetism.

Acknowledgement. We thank Hong Yao and Yingfei Gu for helpful discussions.
This work is supported by Innovation Program for Quantum Science and
Technology 2021ZD0302005, the Beijing Outstanding Young Scholar Program, the
XPLORER Prize and China Postdoctoral Science Foundation (Grant No.
2022M711868). C.L. is supported by Chinese International Postdoctoral Exchange
Fellowship Program and Shuimu Tsinghua Scholar Program at Tsinghua University.

## References

  * (1) E. C. Stoner, Collective electron ferromagnetism, Proc. Roy. Soc. A. 165, 372 (1938).
  * (2) J. Hubbard, Electron correlations in narrow energy bands, Proc. Roy. Soc. A. 276, 238 (1963).
  * (3) M. C. Gutzwiller, Effect of Correlation on the Ferromagnetism of Transition Metals, Phys. Rev. Lett. 10, 159 (1963).
  * (4) J. Kanamori, Electron Correlation and Ferromagnetism of Transition Metals, Prog. Theor. Phys. 30, 275 (1963).
  * (5) H. Tasaki, The Hubbard model - an introduction and selected rigorous results, J. Phys.: Condens. Matter 10, 4353 (1998).
  * (6) H. Tasaki, From Nagaoka’s Ferromagnetism to Flat-Band Ferromagnetism and Beyond: An Introduction to Ferromagnetism in the Hubbard Model, Prog. Theor. Phys. 99, 489 (1998).
  * (7) D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, The Hubbard model, Annu. Rev. Condens. Matter Phys. 13, 239 (2022).
  * (8) M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard model: A computational perspective, Annu. Rev. Condens. Matter Phys. 13, 275 (2022).
  * (9) D. J. Thouless, Exchange in solid 3He and the Heisenberg Hamiltonian, Proc. Phys. Soc. 86, 893 (1965).
  * (10) Y. Nagaoka, Ferromagnetism in a narrow, almost half-filled $s$ band, Phys. Rev. 147, 392 (1966).
  * (11) H. Tasaki, Extension of Nagaoka’s theorem on the large-U Hubbard model, Phys. Rev. B 40, 9192 (1989).
  * (12) A. Mielke, Ferromagnetism in the Hubbard model on line graphs and further considerations, J. Phys. A: Math. Gen. 24, 3311 (1991).
  * (13) I. Bloch, Ultracold quantum gases in optical lattices, Nat. Phys. 1, 23 (2005).
  * (14) I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (15) M. Xu, L. H. Kendrick, A. Kale, Y. Gang, G. Ji, R. T. Scalettar, M. Lebrat, and M. Greiner, Doping a frustrated Fermi-Hubbard magnet, arXiv:2212.13983.
  * (16) J. Merino, B. J. Powell, and Ross H. McKenzie, Ferromagnetism, paramagnetism, and a Curie-Weiss metal in an electron-doped Hubbard model on a triangular lattice, Phys. Rev. B 73, 235107 (2006).
  * (17) I. Morera, M. Kanász-Nagy, T. Smolenski, L. Ciorciaro, A. Imamoğlu, E. Demler, High-temperature kinetic magnetism in triangular lattices, arXiv:2209.05398.
  * (18) R. P. Feynman, Statistical Mechanics: A Set Of Lectures (CRC Press, 1998) Chapter 11.
  * (19) The code for our numerical calculations are available at https://github.com/chengshul/itinerantFM.
  * (20) L. Liu, H. Yao, E. Berg, S. R. White, and S. A. Kivelson, Phases of the Infinite $U$ Hubbard Model on Square Lattices, Phys. Rev. Lett. 108, 126406 (2012).
  * (21) E. Müller-Hartmann, Ferromagnetism in Hubbard models: Low density route, J. Low Temp. Phys. 99, 349 (1995).
  * (22) P. Pieri, S. Daul, D. Baeriswyl, M. Dzierzawa, and P. Fazekas, Low-density ferromagnetism in the Hubbard model, Phys. Rev. B 54, 9250 (1996).
  * (23) Here valley means two sets of low-energy states nearby degenerate single-particle minima, which is different from valley ferromagnetism discussed in literatures, such as N. J. Zhang, J. Lin, Y. Wang, K. Watanabe, T. Taniguchi, L. Fu, J.I.A. Li, arXiv:2209.12964.
  * (24) E. Bobrow, K. Stubis, and Y. Li, Exact results on itinerant ferromagnetism and the 15-puzzle problem, Phys. Rev. B 98, 180101(R) (2018).
  * (25) Here for $\phi=\pi$ case, the non-magnetic state at $n\sim 0.1-0.3$ is anomalous. We find the results in this density range with $\phi=\pi$ quite sensitive to system geometry, and the finite size effect is quite significant in this density range.
  * (26) L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).