# Adversarial Counterfactual Environment Model Learning

Xiong-Hui Chen1,3, , Yang Yu1,3,∗, , Zheng-Mao Zhu1, Zhihua Yu1,2, Zhenjun
Chen1,2
Chenghe Wang1, Yinan Wu2, Hongqiu Wu1, Rong-Jun Qin1,3
Ruijin Ding4, Fangsheng Huang2
1 National Key Laboratory of Novel Software Technology, Nanjing University
2 Meituan, 3 Polixir.ai, 4 Tsinghua University
{chenxh, yuy<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>The authors contribute
equallyYang Yu is the corresponding author.

###### Abstract

A good model for action-effect prediction, named environment model, is
important to achieve sample-efficient decision-making policy learning in many
domains like robot control, recommender systems, and patients’ treatment
selection. We can take unlimited trials with such a model to identify the
appropriate actions so that the costs of queries in the real world can be
saved. It requires the model to correctly handle unseen data, also called
counterfactual data. However, standard data fitting techniques do not
automatically achieve such generalization ability and commonly result in
unreliable models. In this work, we introduce counterfactual-query risk
minimization (CQRM) in model learning for generalizing to a counterfactual
dataset queried by a specific target policy. Since the target policies can be
various and unknown in policy learning, we propose an adversarial CQRM
objective in which the model learns on counterfactual data queried by
adversarial policies, and finally derive a tractable solution GALILEO. We also
discover that adversarial CQRM is closely related to the adversarial model
learning, explaining the effectiveness of the latter. We apply GALILEO in
synthetic tasks and a real-world application. The results show that GALILEO
makes accurate predictions on counterfactual data and thus significantly
improves policies in real-world testing.

##  1 Introduction

A good environment model is important for sample-efficient decision-making
policy learning techniques like reinforcement learning (RL) [1, 2, 3, 4, 5].
The agent can take trials with this model to find better policies, then the
costly real-world trial-and-errors can be saved [1, 6, 2] or completely waived
[3, 4, 5]. In this process, the core of the models is to answer queries on
counterfactual data unbiasedly, that is, given states, correctly answer _what
might happen_ if we were to carry out actions _unseen_ in the training data
[7]. Moreover, environment models are also essential for training transferable
meta-policy [2, 8].

The requirement of counterfactual queries makes the environment model learning
being essentially different from standard supervised learning (SL), and model
learning through SL might results in unreliable models: In real-world
applications, the offline data is often collected with selection bias, that
is, for each state, each action might be chosen unfairly (see Fig. 1(a)). When
a dataset is collected with selection bias, the high correlation between
states and actions makes the causality to the next states agnostic, which
finally leads to the poor generalization ability of model prediction on
counterfactual data. In Fig. 1(c), we give an illustration of catastrophic
failures caused by the issue. Generally speaking, it corresponds to a
challenge of training the model in one dataset but testing in another dataset
(i.e., the dataset generated by counterfactual queries) with a different
distribution. This is beyond the capability of SL since it violates the
independent and identically distributed (_i.i.d._) assumption.

(a) An example of data collection with selection bias

(b) response curves on the training data

(c) response curves on counterfactual test data

Figure 1: An example of selection bias and the effects of predictions under
counterfactual queries. Subfigure (a) shows how the data is collected: a ball
locates in a 2D-plane whose position is $(x_{t},y_{t})$ at time $t$. The ball
will move to $(x_{t+1},y_{t+1})$ according to $x_{t+1}=x_{t}+1$ and
$y_{t+1}\sim\mathcal{N}(y_{t}+a_{t},2)$. Here, $a_{t}$ is chosen by a control
policy $a_{t}\sim\mathcal{N}((62.5-y_{t})/15,0.05)$ which tries to keep the
ball near the line $y=62.5$. Subfigure (b) shows the collected training data
(grey dashed line) and the response curves of the two learned models, which
predict the next position of $y$. Both the models discovered the superficial
relation that the corresponding next $y$ will be smaller with a larger action.
However, the truth is not because the larger $a_{t}$ causes smaller $y_{t+1}$,
but the policy selects a small $a_{t}$ when $y_{t}$ is large. When we conduct
counterfactual queries by fixing $y_{t}$ and reassigning action $a_{t}$ with
counterfactual action $a_{t}+\Delta a$, where $\Delta a\in[-1,1]$ is a
variation of action value, the model of SL will exploit the correlation and
give opposite responses, while GALILEO’s predictions are closer to the ground-
truths. The result is in Subfigure (c), where the darker a region is, the more
samples are fallen in.

In offline model-based RL, the problem is called distribution shift [6, 7, 2]
which has received great attention. However, previous algorithms do not handle
the model learning challenge directly but propose techniques to suppress
policy sampling and learning in risky regions [6, 9]. Although these
algorithms have made great progress in offline policy optimization in many
tasks, so far, how to learn a better environment model in this scenario has
rarely been discussed.

On the other hand, in decision making based on individual treatment effects
(ITEs) estimation [10, 11, 12, 13], learning a model with better
generalization ability in counterfactual data has been studied for the problem
like patients’ treatment selection. In this problem, we are interested in
knowing the effect differences of treatments under a uniform policy for each
individual, called ITEs. With an model to estimate ITEs, a policy is deduced
by selecting the treatments with maximal positive ITEs for each individual.
Since the ITEs are defined as the effect difference of treatments under a
uniform policy, ITEs can be estimated unbiasedly based on the model learned on
a balanced dataset with the independent and identically distribution as the
dataset collected by a uniform policy. The balanced dataset can be derived by
different techniques like propensity-score model [14, 15] or representation
function [16, 17, 18]. However, only having an unbiased estimation of a
uniform policy is not enough for general policy learning process since we have
to query numerous different policies’ feedback besides the uniform policy.
Moreover, these methods only focus on one-step instead of multi-step
environment model learning.

In this work, we propose an objective called Adversarial Counterfactual-Query
Risk Minimization (ACRM). ACRM introduces an adversarial policy to maximize
the prediction error of the learned model on counterfactual data queried by
the policy, then learns a model to be as small prediction risks as possible
under counterfactual data queried by adversarial policies. Second, we give an
approximation of the counterfactual data distribution queried by the optimal
adversarial policy under any given environment model. Then we use the
variational representation to give a tractable solution to model the
approximated data distribution. As a result, we derive a practical approach
named Generative Adversarial offLIne counterfactuaL Environment mOdel learning
(GALILEO) for ACRM. GALILEO is compatible with one-step and multi-step
environment model learning and we give two instantiations for the two
scenarios respectively. We also discover that ACRM is closely related to the
generative-adversarial model learning, which learns a model in an adversarial
counterfactual dataset constructed by adaptively reshaping the distribution of
offline data, explaining the effectiveness of the latter.

Experiments are conducted in two synthetic environments and two realistic
environments. The results in the synthetic environments show that GALILEO can
reconstruct correct responses for counterfactual queries. The evaluation
results in two realistic environments also demonstrate that GALILEO has better
ability in counterfactual query compared with baselines. We finally search for
a policy based on the learned model in a real-world online platform. The
policy significantly improves performance in concerned business indicators.

##  2 Related Work

Generative-adversarial-nets is a useful tool in many domains [19, 20, 21, 22].
We give related adversarial algorithms for model learning in the following and
leave other related work in Appx. E. GANTIE [23] uses a generator to fill
counterfactual outcomes for each data pair and a discriminator to judge the
source (treatment group or control group) of the filled data pair. The
generator is trained to minimize the output of the discriminator. GANITE is
trained until the discriminator cannot determine which of the components is
the factual outcome. [24] propose SCIGAN to extend GANITE to continuous
treatment effect estimation (a.k.a., dosage-response estimation) via a
hierarchical discriminator architecture. In real-world applications,
environment model learning based on Generative Adversarial Imitation Learning
(GAIL) [22] has also been adopted for sequential decision-making problems.
GAIL can be regarded as an extension of GAN in multi-step environments. It is
first proposed for policy imitation [22], which uses the imitated policy to
generate trajectories by interacting with the environment. The policy is
learned with the trajectories through RL which maximizes the cumulative
rewards given by the discriminator. [3, 25, 4] use GAIL for environment model
learning by regarding the environment model as the generator and the behavior
policy as the “environment” in standard GAIL. Although these studies tried to
generate counterfactual data through adversarial generative methods, there was
no proof nor solid evidence that the counterfactual data matching the ground
truth can be so generated. Our study shows that these algorithms are partial
implementations of GALILEO to achieve ACRM.

##  3 Background and Notations

###  3.1 Individualized treatment effects estimation

We first introduce individualized treatment effects (ITEs) estimation in the
two-treatment scenario [26]. ITEs are defined as
$\textit{ITE}(x):=\mathbb{E}[M^{*}(y|x,1)|A=1,X=x]-\mathbb{E}[M^{*}(y|x,0)|A=0,X=x]$,
where $y$ is the feedback of the environment $M^{*}(y|x,a)$, $X$ denotes the
state vector containing pre-treatment covariates (such as age and weight), $A$
denotes the treatment variable which is the action intervening to the state
$X$, and $A$ should be sampled from a uniform policy. In the two-treatment
scenario, $A$ is in $\\{0,1\\}$ where $1$ is the action to intervene and $0$
is the action to do nothing. A correct ITEs estimation should be done in
Randomized Controlled Trials (RCT) in which we have the same probability of
samples of $A=1$ and $A=0$ for each $X$. Here we use lowercase $x,a$ and $y$
to denote samples of random variables $X,A$ and $Y$, and use
$\mathcal{X},\mathcal{A}$ and $\mathcal{Y}$ to denote space of the samples. In
the real world, RCT is expensive or even infeasible [26]. In practice, we
prefer to estimate ITEs under observational studies. In observational studies,
datasets are pre-collected from the real world by a behavior policy such as a
human-expert policy. In this case, a common approach for estimating the ITEs
can be [27]

$\displaystyle\hat{\textit{ITE}}(x_{i})\begin{cases}y^{F}_{i}-M(x_{i},1-a_{i}),a_{i}=1\\\
M(x_{i},1-a_{i})-y^{F}_{i},a_{i}=0\\\ \end{cases}$

in deterministic prediction, where $x_{i}$ and $y^{F}_{i}$ denote the
covariate and factual feedback of the $i$-th sample, and $M\in\mathcal{M}$
denotes an approximated feedback model. $\mathcal{M}$ is space of the model.
In this formulation, the training set is an empirical factual data
distribution $P_{F}=\\{x_{i},a_{i}\\}^{n}_{i}$ and the testing set is an
empirical counterfactual data distribution
$P_{CF}=\\{x_{i},1-a_{i}\\}^{n}_{i}$. If $a$ does not sample from a discrete
uniform policy, i.e., the policy has selection bias, $P_{F}$ and $P_{CF}$ will
be two different distributions, which violate the i.i.d. assumption of
standard supervised learning. In stochastic prediction,
$\hat{\textit{ITE}}(x)=\mathbb{E}[M(y|x,1)]-\mathbb{E}[M(y|x,0)]$ and the
counterfactual distribution for testing is the dataset with action sampling
from a uniform policy.

Individual dose-response curves estimation problem is an extended ITEs
estimation in the continuous-action-space scenario [28, 24]. The goal is to
derive an unbiased model of the potential outcomes under a uniform policy.
Mean Integrated Square Error (MISE) is a commonly used metric to measure the
accuracy of counterfactual queries:
$\mathrm{MISE}=\mathbb{E}\left[\int_{\mathcal{A}}\left(M^{*}(y|x,a)-M(y|x,a)\right)^{2}\mathrm{d}a\right]$.

###  3.2 Counterfactual-Query Risk Minimization

Generally speaking, in ITEs estimation, the extra risks of queries under
counterfactual data are caused by the gap between the policy in training and
testing data distributions. Without further processing, minimizing the
empirical risks cannot guarantee the counterfactual-query risks being
minimized. Assuming that the policy in training data $\mu$ satisfies
$\mu(a|x)>0,\forall a\in\mathcal{A},\forall x\in\mathcal{X}$ (often named
overlap assumption), a classical solution to handle the above Counterfactual-
Query Risk Minimization (CQRM) problem is by re-weighting the empirical risks
through an Inverse Propensity Scoring (IPS) term $\omega$ [29, 30, 31, 32]:

######  Definition 3.1.

The learning objective of counterfactual-query risk minimization through IPS
is formulated as

$\begin{split}\min_{M\in\mathcal{M}}L(M)&=\min_{M\in\mathcal{M}}\mathbb{E}_{x,a,y\sim
p^{\mu}_{M^{*}}}[\omega(x,a)\ell(M(y|x,a),y)],\\\
&s.t.\quad\omega(x,a)=\frac{\beta(a|x)}{\mu(a|x)},\end{split}$ (1)

where $\beta$ and $\mu$ denote the policies in testing and training domains,
and $p^{\mu}_{M^{*}}$ is the joint probability
$p^{\mu}_{M^{*}}(x,a,y):=\rho_{0}(x)\mu(a|x)M^{*}(y|x,a)$ in which
$\rho_{0}(x)$ is the distribution of state. $\mathcal{M}$ is the model space.
$\ell$ is a loss function for model learning.

The $\omega$ is also known as importance sampling (IS) weight, which corrects
the sampling bias. In this objective, $\omega$ is to align the training data
distribution to the testing data. By selecting different $\hat{\omega}$ to
approximate $\omega$ to learn the model $M$, current environment model
learning algorithms for ITEs estimation are fallen into the framework. In
standard supervised learning and some works for ITEs estimation [33, 34],
$\hat{\omega}(x,a)=1$ as the distribution-shift problem is ignored. In [35,
36, 29, 30], $\omega=\frac{1}{\hat{\mu}}$ (i.e., $\beta$ a uniform policy) for
balancing treatment and control group, where $\hat{\mu}$ is an approximation
of behavior policy $\mu$. Note that it is a reasonable weight in ITEs
estimation: ITEs are defined to evaluate the effect of each state between
treatment and control behavior under a uniform policy.

###  3.3 Markov Decision Process

Markov Decision Process (MDP) is a framework that formulates a decision-making
process in a sequential environment [37]. MDP depicts an agent interacting
with a sequential environment through actions. In the first step, states are
sampled from an initial state distribution $x_{0}\sim\rho_{0}(x)$. Then at
each time-step $t\in\\{0,1,2,...\\}$, the agent takes an action
$a_{t}\in\mathcal{A}$ through a policy $\pi(a_{t}|x_{t})\in\Pi$ based on the
state $x_{t}\in\mathcal{X}$, then the agent receives a reward $r_{t}$ from a
reward function $r(x_{t},a_{t})\in\mathbb{R}$ and transits to the next state
$x_{t+1}$, which is given by a transition function
$M^{*}(x_{t+1}|x_{t},a_{t})$ built in the environment. $\Pi$, $\mathcal{X}$,
and $\mathcal{A}$ denote the policy, state, and action spaces, respectively.

##  4 Method

In this section, we first propose a new model-learning objective based on Def.
3.1 for policy optimization tasks in Sec. 4.1; In Sec. 4.2, we derive a
tractable solution to the proposed objective; Finally, we give a practical
implementation in Sec. 4.3.

###  4.1 Problem Formulation

For policy optimization, we require the environment model having
generalization ability in counterfactual queries since we need to query
numerous different policies’ correct feedback from $M$. Refer to the
formulation of counterfactual-query risk minimization through IPS in Def. 3.1,
policy optimization scenario requires $M$ minimizing counterfactual-query
risks under numerous unknown different policies rather than a specific target
policy $\beta$. More specifically, the question is: If $\beta$ is unknown and
can be varied, how should we reduce the risks in counterfactual queries? In
this article, we call the model learning problem in this setting as
counterfactual environment model learning and propose a new objective in the
following to handle the problem.

To be compatible with multi-step environments model learning, we first define
a generalized CQRM through IPS based on Def. 3.1.

######  Definition 4.1.

Given the MDP transition function $M^{*}$, the learning objective of multi-
step generalized counterfactual-query risk minimization through IPS is
formulated as

$\begin{split}\min_{M\in\mathcal{M}}L(M)&=\min_{M\in\mathcal{M}}\mathbb{E}_{x,a,x^{\prime}\sim\rho_{M^{*}}^{\mu}}[\omega(x,a,x^{\prime})\ell_{M}(x,a,x^{\prime})],\\\
&s.t.\quad\omega(x,a,x^{\prime})=\frac{\rho_{M^{*}}^{\beta}(x,a,x^{\prime})}{\rho_{M^{*}}^{\mu}(x,a,x^{\prime})},\end{split}$
(2)

where $\rho_{M^{*}}^{\mu}$ is the training data distribution (collected by
policy $\mu$), $\rho_{M^{*}}^{\beta}$ is the testing data distribution
(collected by policy $\beta$). We define
$\ell_{M}(x,a,x^{\prime}):=\ell(M(x^{\prime}|x,a),x^{\prime})$ for brevity.

In an MDP, given any policy $\pi$,
$\rho_{M^{*}}^{\pi}(x,a,x^{\prime})=\rho_{M^{*}}^{\pi}(x)\pi(a|x)M^{*}(x^{\prime}|x,a)$
where $\rho_{M^{*}}^{\pi}(x)$ denotes the occupancy measure of $x$ for policy
$\pi$, which can be defined as
$\rho_{M^{*}}^{\pi}(x):=(1-\gamma)\mathbb{E}_{x_{0}\sim\rho_{0}}\left[\sum^{\infty}_{t=0}\gamma^{t}{\rm
Pr}(x_{t}=x|x_{0},{M^{*}})\right]$ [37, 22] where ${\rm
Pr}^{\pi}\left[x_{t}=x|x_{0},{M^{*}}\right]$ is the state visitation
probability that the policy $\pi$ visits $x$ at time-step $t$ by executing in
the environment $M^{*}$ and starting at the state $x_{0}$, and
$\gamma\in\left[0,1\right]$ is the discount factor. Here we also define
$\rho^{\pi}_{M^{*}}(x,a):=\rho^{\pi}_{M^{*}}(x)\pi(a|x)$ for simplicity. In
this definition, $\omega$ can be rewritten as:
$\omega(x,a,x^{\prime})=\frac{\rho_{M^{*}}^{\beta}(x)\beta(a|x)M^{*}(x^{\prime}|x,a)}{\rho_{M^{*}}^{\mu}(x)\mu(a|x)M^{*}(x^{\prime}|x,a)}=\frac{\rho_{M^{*}}^{\beta}(x,a)}{\rho_{M^{*}}^{\mu}(x,a)}$.
In one-step environments, for any policy $\pi$,
$\rho_{M^{*}}^{\pi}(x)=\rho_{0}(x)$. Then we have
$\omega(x,a,x^{\prime})=\frac{\rho_{0}(x)\beta(a|x)}{\rho_{0}(x)\mu(a|x)}=\frac{\beta(a|x)}{\mu(a|x)}$,
and the objective is degraded to Eq.(1). Therefore, Def. 3.1 is a special case
of this generalized form.

Since $\beta$ for counterfactual queries is unknown and can be varied in
policy optimization, to reduce the risks in counterfactual queries in this
scenario, we introduce adversarial policies which can induce worst performance
of the model predictions and propose to optimize CQRM under the adversarial
policies. In particular, we propose Adversarial Counterfactual-query Risk
Minimization (ACRM) based on Def. 4.1 to handle this problem.

######  Definition 4.2.

Given the MDP transition function $M^{*}$, the learning objective of
adversarial counterfactual-query risk minimization through IPS is formulated
as

$\begin{split}\hat{M}^{*}&=\min_{M\in\mathcal{M}}\max_{\beta\in\Pi}L(\rho_{M^{*}}^{\beta},M)=\min_{M\in\mathcal{M}}\max_{\beta\in\Pi}\mathbb{E}_{x,a,x^{\prime}\sim\rho_{M^{*}}^{\mu}}[\omega(x,a|\rho_{M^{*}}^{\beta})\ell_{M}(x,a,x^{\prime})],\\\
&s.t.\quad\omega(x,a|\rho_{M^{*}}^{\beta})=\frac{\rho_{M^{*}}^{\beta}(x,a)}{\rho_{M^{*}}^{\mu}(x,a)},\end{split}$
(3)

where the re-weighting term $\omega(x,a|\rho_{M^{*}}^{\beta})$ is conditioned
on the data distribution $\rho_{M^{*}}^{\beta}$ of the adversarial policy
$\beta$. In the following we will ignore $\rho_{M^{*}}^{\beta}$ and use
$\omega(x,a)$ for brevity when the context is clear.

By definition of MDP, $\omega(x,a,x^{\prime})=\omega(x,a)$ since the
transition probability in the ratio can be canceled. Eq.(3) minimizes the
maximum model loss under all counterfactual data distributions
$\rho_{M^{*}}^{\beta}(x,a,x^{\prime}),\beta\in\Pi$ so that to guarantee the
generalization ability for counterfactual data queried by policies in $\Pi$.

$\min_{M\in\mathcal{M}}\max_{\beta\in\Pi}L(\rho_{M^{*}}^{\beta},M)$
$\rho^{\bar{\beta}^{*}}_{M^{*}}$ given $M$ generalized representation
$\bar{\rho}^{\bar{\beta}^{*}}_{M^{*}}$ Thm. 4.5:
$\min_{M\in\mathcal{M}}L(\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}},M)$ tractable
solution Eq.(8) $\kappa$$f$ estimate the distribution of the best-response
policy $\bar{\beta}^{*}$ approximate $\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}$
with variational representation

Figure 2: The overall pipeline to model the tractable solution to ACRM. $f$ is
a generator function defined by $f$-divergence [20]. $\kappa$ is an
intermediary policy introduced in the estimation process.

###  4.2 A Tractable Solution to ACRM

In this section, we propose a tractable solution to optimize Eq.(3). We give
an overall pipeline of the solution modeling in Fig. 2 and the full proof in
Appx. A. In this article, we derive the solution to Eq.(3) based on $\ell_{M}$
being negative log-likelihood loss, then $L(\rho_{M^{*}}^{\beta},M)$ equals
to:

$\displaystyle\mathbb{E}_{x,a\sim\rho_{M^{*}}^{\mu}}\left[\omega(x,a|\rho_{M^{*}}^{\beta})\mathbb{E}_{M^{*}}\left(-\log
M(x^{\prime}|x,a)\right)\right],$

where $\mathbb{E}_{M^{*}}\left[\cdot\right]$ denotes
$\mathbb{E}_{x^{\prime}\sim M^{*}(x^{\prime}|x,a)}\left[\cdot\right]$.

The first problem is how to construct the data distribution
$\rho_{M^{*}}^{\beta^{*}}$ of the best-response policy $\beta^{*}$ in the real
environment $M^{*}$. As it is costly to get extra data from $M^{*}$ in real-
world applications, our solution is to estimate the distribution of
$\rho_{M^{*}}^{\beta^{*}}$ for each $M$ and construct a surrogate objective to
avoid directly querying the real environment $M^{*}$ to optimize $M$. Ideally,
given any $M$, it is obvious that the optimal $\beta$ is the one that makes
$\rho_{M^{*}}^{\beta}(x,a)$ satisfy

$\displaystyle\rho_{M^{*}}^{\beta}(x,a)=\begin{cases}\infty,&(x,a)=\arg\max_{x,a}\mathbb{E}_{M^{*}}\left[-\log
M(x^{\prime}|x,a)\right]\\\ 0,&\textit{otherwise}\end{cases},$

which assigns all densities to the point that has the largest negative log-
likelihood. However, searching for the maximum is impractical, especially in
continuous space. To give a relaxed but tractable solution, we add an $L_{2}$
regularizer to the original objective Eq.(3):

$\begin{split}\min_{M\in\mathcal{M}}\max_{\beta\in\Pi}~{}\bar{L}(\rho_{M^{*}}^{\beta},M)=\min_{M\in\mathcal{M}}\max_{\beta\in\Pi}~{}\mathbb{E}_{x,a\sim\rho_{M^{*}}^{\mu}}\Big{[}\omega(x,a)\mathbb{E}_{M^{*}}\left[-\log
M(x^{\prime}|x,a)\right]\Big{]}-\frac{\alpha}{2}\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2},\end{split}$
(4)

where $\alpha$ denotes the regularization coefficient of
$\rho_{M^{*}}^{\beta}$ and
$\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}=\int_{\mathcal{X},\mathcal{A}}(\rho_{M^{*}}^{\beta}(x,a))^{2}\mathrm{d}a\mathrm{d}x$.
Then we can approximate the optimal distribution
$\rho_{M^{*}}^{\bar{\beta}^{*}}$ via Lemma. 4.3.

######  Lemma 4.3 (Approximated optimal adversarial distribution).

Given any $M$ in $\bar{L}(\rho^{\beta}_{M^{*}},M)$, the distribution of the
best-response policy $\bar{\beta}^{*}$ satisfies:

$\rho_{M^{*}}^{\bar{\beta}^{*}}(x,a)=\frac{1}{\alpha}(D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))+H_{M^{*}}(x,a)),$
(5)

where $D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))$ is the Kullback-Leibler (KL)
divergence between $M^{*}(\cdot|x,a)$ and $M(\cdot|x,a)$, $H_{M^{*}}(x,a)$
denotes the entropy of $M^{*}(\cdot|x,a)$, and $\alpha$ is a constant
normalizer.

Intuitively, $\bar{\beta}^{*}$ make $\rho_{M^{*}}^{\bar{\beta}^{*}}$ have
larger densities on the data where the distribution distance between the
approximation model and the real model (i.e.,
$D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))$) is larger or the stochasticity of the
real model (i.e., $H_{M^{*}}$) is larger. We know that KL divergence is an
instance of a divergence class, called $f$-divergence [20], which measure the
difference between two given probability distributions. We give a generalized
representation of $\rho_{M^{*}}^{\bar{\beta}^{*}}$ by replacing $D_{KL}$ in
Eq.(5) with $f$-divergence and converting the distribution metric from
conditional probabilities to joint probabilities. The latter is for easier
estimation in implementation.

######  Lemma 4.4 (Generalized optimal adversarial distribution of
$\bar{\beta}^{*}$).

Replacing $D_{KL}$ in Lemma 4.3 with $f$-divergence, we have

$\displaystyle\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}(x,a)=$
$\displaystyle\frac{1}{\alpha_{0}(x,a)}\Bigg{(}\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})f(\frac{\rho^{\kappa}_{M}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})})\mathrm{d}x^{\prime}-\rho^{\kappa}_{M^{*}}(x,a)\left(f(\frac{\rho^{\kappa}_{M}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)})-H_{M^{*}}(x,a)\right)\Bigg{)}.$

A new policy $\kappa(a|x)$ is introduced to construct joint probabilities
$\rho^{\kappa}_{M^{*}}(x,a)$ and $\rho^{\kappa}_{M}(x,a)$. $\kappa$ is a
policy satisfying $\forall x\in\mathcal{X},\forall
a\in\mathcal{A},\rho_{M^{*}}^{\kappa}(x,a)>0$ if
$\rho_{M^{*}}^{\bar{\beta}^{*}}(x,a)>0$. Here we call $\kappa$ an intermediary
policy. $f:\mathbb{R}_{+}\rightarrow\mathbb{R}$ is a generator function of
$f$-divergence. $\alpha_{0}(x,a)=\alpha\rho^{\kappa}_{M^{*}}(x,a)$.

$\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}(x,a)$ gives an approximation of the
optimal adversarial data distribution which maximizes the error of the model.
It is constructed on a data distribution $\rho^{\kappa}_{M^{*}}$, which is
collected by an intermediary policy $\kappa$ in the real-world environment
$M^{*}$. The assumption “$\forall x\in\mathcal{X},\forall
a\in\mathcal{A},\rho_{M^{*}}^{\kappa}(x,a)>0$ if
$\rho_{M^{*}}^{\bar{\beta}^{*}}(x,a)>0$” is satisfied if
$\kappa(a|x)>0,\forall x\in\mathcal{X},\forall a\in\mathcal{A}$, while the
assumption of property of $\kappa$ is common in named overlap [26].

######  Theorem 4.5.

Let $\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}$ in Lemma 4.4 as the data
distribution of the best-response policy $\bar{\beta}^{*}$ in Eq.(4) under
model $M_{\theta}$ parameterized by $\theta$, then we can skip the
“$\max_{\beta}[\cdot]$” and directly find the optimal $\theta^{*}$ of
$\min_{\theta}\max_{\beta\in\Pi}\bar{L}(\rho_{M^{*}}^{\beta},M_{\theta})$
(Eq.(4)) through the following objective:

$\displaystyle\theta^{*}=\arg\max_{\theta}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\Bigg{[}\frac{1}{\alpha_{0}(x,a)}\log
M_{\theta}(x^{\prime}|x,a)\underbrace{\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)-f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)+H_{M^{*}}(x,a))\right)}_{W(x,a,x^{\prime})}\Bigg{]},$
(6)

where $\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[\cdot\right]$ denotes
$\mathbb{E}_{x,a,x^{\prime}\sim\rho^{\kappa}_{M^{*}}}\left[\cdot\right]$, $f$
is the generator function in $f$-divergence which satisfies $f^{\prime}(x)\leq
0,\forall x\in\mathcal{X}$. $M_{\bar{\theta}}$ denotes a probability function
with the same parameters as the learned model (i.e., $\bar{\theta}=\theta$)
but the parameter is fixed and only used for sampling.

Thm. 4.5 gives us a surrogate objective of ACRM which avoids learning the
best-response policy $\beta^{*}$ and collecting the dataset running with
$\beta^{*}$ explicitly. ACRM is approximately achieved by using an
intermediary policy $\kappa$ and a pseudo re-weighting module $W$.
Intuitively, $W$ assigns learning propensities for data points with larger
differences between distributions $\rho^{\kappa}_{M_{\bar{\theta}}}$ and
$\rho^{\kappa}_{M^{*}}$. By adjusting the weights, the learning process will
exploit subtle errors in any data point, whatever how many proportions it
contributes, to correct potential generalization errors on counterfactual
data.

In the following, we give a tractable solution to Thm. 4.5.

Note that the part $\rho^{\kappa}_{M^{*}}(x,a)$ in
$\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})$ in Eq.(6) can be cancelled because
$\alpha_{0}(x,a)=\alpha\rho^{\kappa}_{M^{*}}(x,a)$, but we choose to keep it
and ignore $\alpha_{0}(x,a)$. The benefit is that we can estimate
$\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})$ from an empirical data distribution
through data collected by $\kappa$ in $M^{*}$ directly, rather than from a
uniform distribution which is harder to be generated. Although keeping
$\rho^{\kappa}_{M^{*}}(x,a)$ incurs extra bias in theory, the results in our
experiments show that it has not made significant negative effects in
practice. We leave this part of modeling in future work.

The generator function $f$ is still intractable. To solve the problem, first,
we resort to the first-order approximation. Given some
$a\in(1-\xi,1+\xi),\xi>0$, we have

$\displaystyle f(a)\approx f(1)+f^{\prime}(a)(a-1),$ (7)

where $f^{\prime}$ is the first-order derivative of $f$. By Taylor’s formula
and the fact that $f^{\prime}(x)$ of the generator function $f$ is bounded in
$(1-\xi,1+\xi)$, the approximation error is no more than
$\mathcal{O}(\xi^{2})$. Let $a=\frac{p(x)}{q(x)}$ in Eq.(7), the pattern
$f(\frac{p(x)}{q(x)})$ in Eq.(6) can be converted to
$f^{\prime}(\frac{p(x)}{q(x)})$ and estimated through Lemma 4.6:

######  Lemma 4.6 ($f^{\prime}(\frac{p}{q})$ estimation [20]).

By finding the maximized point of $\varphi$, that is,

$\displaystyle\varphi^{*}=\arg\max_{\varphi}\mathbb{E}_{x\sim
p(x)}\left[T_{\varphi}(x)\right]+\mathbb{E}_{x\sim
q(x)}\left[f^{*}(T_{\varphi}(x))\right],$

we have $f^{\prime}(\frac{p(x)}{q(x)})=T_{\varphi^{*}}(x)$, where $p(x)$ and
$q(x)$ are two distributions and $f^{*}$ is Fenchel conjugate of the generator
function $f$ [20, 38].

By computing $f(\cdot)$ with first-order approximation (Eq.(7)) and leveraging
the variational representation $T_{\varphi^{*}}$ to approximate
$f^{\prime}(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}}{\rho^{\kappa}_{M^{*}}})$,
we can optimize the surrogate objective Eq.(6) via

$\begin{split}\hat{\theta}^{*}=\arg&\max_{\theta}~{}\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[A(x,a,x^{\prime})\log{M_{\theta}}(x^{\prime}|x,a)\right]+\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[(H_{M^{*}}(x,a)-A(x,a,x^{\prime}))\log{M_{\theta}}(x^{\prime}|x,a)\right],\\\
s.t.\quad&A(x,a,x^{\prime})=T_{\varphi_{0}^{*}}(x,a,x^{\prime})-T_{\varphi_{1}^{*}}(x,a)\\\
&\varphi_{0}^{*}=\arg\max_{\varphi_{0}}~{}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[T_{\varphi_{0}}(x,a,x^{\prime})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\theta}}}\left[f^{*}(T_{\varphi_{0}}(x,a,x^{\prime}))\right]\\\
&\varphi_{1}^{*}=\arg\max_{\varphi_{1}}~{}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[T_{\varphi_{1}}(x,a)\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\theta}}}\left[f^{*}(T_{\varphi_{1}}(x,a))\right],\end{split}$
(8)

where $\mathbb{E}_{\rho^{\kappa}_{M}}[\cdot]$ is a simplification of
$\mathbb{E}_{x,a,x^{\prime}\sim\rho^{\kappa}_{M}}[\cdot]$ and $\varphi_{0}$
and $\varphi_{1}$ are the parameters of $T_{\varphi_{0}}$ and
$T_{\varphi_{1}}$. Based on Lemma 4.6, we have
$f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\approx
T_{\varphi_{0}^{*}}(x,a,x^{\prime})$ and
$f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)\approx
T_{\varphi_{1}^{*}}(x,a)$. The full proof is in Appx. A. In Eq.(8), the pseudo
re-weighting module $W$ in objective Eq.(6) is split into two terms in the RHS
of the equation. The first terms is a discriminator objective as in the
generative adversarial training by regarding $M_{\theta}$ as the generator,
while second term is a log-likelihood loss re-weighted with
$H_{M^{*}}(x,a)-A(x,a,x^{\prime})$.

###  4.3 Practical Implementation

To give a practical implementation of the solution, three assumptions
introduced in the process of modeling should be handled:
First, $\forall x\in\mathcal{X},f^{\prime}(x)\leq 0$. Based on any specific
$f$-divergence, we can represent $T$ and $f^{*}(T)$ with a discriminator
$D_{\varphi}$. It can be verified that the instance $f(u)=u\log
u-(u+1)\log(u+1)$, $T_{\varphi}(u)=\log D_{\varphi}(u)$, and
$f^{*}(T_{\varphi}(u))=\log(1-D_{\varphi}(u))$ proposed in [20] satisfy this
condition (see Tab. 2 in the Appendix for more details). We select the former
in the implementation.
Second, the approximation of Eq.(7) holds only when $\frac{p(x)}{q(x)}$ is
close to $1$, which might not be satisfied. To handle the problem, we inject a
naive supervised learning loss and replace the second term of the objective
Eq.(8) when the output probability of $D$ is far away from $0.5$
($f^{\prime}(1)=\log 0.5$);
Third, the overlap assumption of $\kappa$. In practice, we need to use the
real-world data to construct the distribution $\rho^{\kappa}_{M^{*}}$ and the
generative data to construct $\rho^{\kappa}_{M_{\bar{\theta}}}$. In the
offline model-learning setting, we only have a real-world dataset
$\mathcal{D}$ collected by the behavior policy $\mu$. We can learn a policy
$\hat{\mu}\approx\mu$ via behavior cloning with the offline data $\mathcal{D}$
[39, 22] and let $\hat{\mu}$ be the intermediary policy $\kappa$. Then we can
regard $\mathcal{D}$ as the empirical data distribution of
$\rho^{\kappa}_{M^{*}}$ and the trajectories collected by $\hat{\mu}$ in the
model $M_{\bar{\theta}}$ as the empirical data distribution of
$\rho^{\kappa}_{M_{\bar{\theta}}}$. But the overlap assumption $\forall
x\in\mathcal{X},\forall a\in\mathcal{A},\mu(a|x)>0$ might not be satisfied. In
behavior cloning, we model $\hat{\mu}$ with a Gaussian distribution and
constrain the lower bound of the variance with a small value
$\epsilon_{\mu}>0$ to keep the assumption holding. Besides, we add small
Gaussian noises $\mathcal{N}(0,\epsilon_{D})$ to the inputs of $D_{\varphi}$
to handle the mismatch between $\rho^{\mu}_{M^{*}}$ and
$\rho^{\hat{\mu}}_{M^{*}}$ due to $\epsilon_{\mu}$.

On the other hand, we notice that the first term in Eq.(8) is similar to the
objective of GAIL [22] by regarding $M_{\theta}$ as the policy to learn and
$\kappa$ as the environment to generate data. For better capability in
sequential environment model learning, here we introduce some practical tricks
inspired by GAIL for model learning [3, 4]: we introduce an MDP for $\kappa$
and $M_{\theta}$, where the reward is defined by the discriminator $D$, i.e.,
$r(x,a,x^{\prime})=\log D(x,a,x^{\prime})$. $M_{\theta}$ is learned to
maximize the cumulative rewards. With advanced policy gradient methods [40,
41], the objective is converted to $\max_{\theta}\left[A(x,a,x^{\prime})\log
M_{\theta}(x,a,x^{\prime})\right]$, where
$A=Q^{\kappa}_{M_{\bar{\theta}}}-V^{\kappa}_{M_{\bar{\theta}}}$,
$Q^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}r(x_{t},a_{t},x_{t+1})\mid(x_{t},a_{t},x_{t+1})=(x,a,x^{\prime}),\kappa,M_{\bar{\theta}}\right]$,
and
$V^{\kappa}_{M_{\bar{\theta}}}(x,a)=\mathbb{E}_{M_{\bar{\theta}}}\left[Q^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})\right]$.
$A$ in Eq.(8) can also be constructed similarly. Although it looks unnecessary
in theory since the one-step optimal model $M_{\theta}$ is the global optimal
model in this setting, the technique is helpful in practice as it makes $A$
more sensitive to the compounding effect of one-step prediction errors: we
would consider the cumulative effects of prediction errors induced by multi-
step of transitions in environments. Therefore, we adopt the trick by defining
$Q^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}\log
D_{\varphi_{0}^{*}}(x_{t},a_{t},x_{t+1})|(x_{t},a_{t},x_{t+1})=(x,a,x^{\prime}),\kappa,M_{\bar{\theta}}\right]$,
$V^{\kappa}_{M_{\bar{\theta}}}(x,a)=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}\log
D_{\varphi_{1}^{*}}(x_{t},a_{t})|(x_{t},a_{t})=(x,a),\kappa,M_{\bar{\theta}}\right]$,
and $A=Q^{\kappa}_{M_{\bar{\theta}}}-V^{\kappa}_{M_{\bar{\theta}}}$.

Based on the above techniques, we propose Generative Adversarial offLIne
counterfactuaL Environment mOdel learning (GALILEO) for sequential environment
model learning. Besides, we give a variant, named GALILEO(0), for one-step
environment model learning by setting $\gamma$ in $Q$ and $V$ to $0$. The
detailed implementation and comparison to previous adversarial methods are in
Appx. D.

##  5 Experiments

In this section, we first conduct experiments in synthetic and benchmark
environments [24] to verify GALILEO on counterfactual queries 111code:
https://github.com/xionghuichen/GALILEO and the compatibility of GALILEO in
sequential and one-step environments. Then we analyze the benefits of the
implementation techniques described in Sec. 4.3 and the problems without them.
Finally, we deploy GALILEO in two complex environments: MuJoCo in Gym [42, 43]
and a real-world food-delivery platform to test the performance of GALILEO in
difficult tasks. The algorithms compared are: (1) Supervised Learning (SL);
(2) Inverse Propensity Weighting (IPW) [44]; (3) eStimating the effects of
Continuous Interventions using GANs (SCIGAN) [24]. We give a detailed
description of the baselines in Appx. F.2.

###  5.1 Test in Synthetic Environment

##### Synthetic sequential environment.

Since there does not exist a task designed for counterfactual queries in
sequential environments, we construct a task, named General Negative Feedback
Control (GNFC), to conduct our experiments, which can represent a classic type
of task with policy having selection bias. Fig. 1(a) is also an example of
GNFC. We give the detailed motivation, the effect of selection bias, and other
details in Appx. F.1.1. We construct tasks on GNFC by adding behavior policies
$\mu$ with different scales of uniform noise $U(-e,e)$ with different
probabilities $p$. In particular, with $e\in\\{1.0,0.2,0.05\\}$ and
$p\in\\{1.0,0.2,0.05\\}$, we construct 9 tasks and name them with the format
of “e*_p*”. For example, e1_p0.2 is the task with behavior policy injecting
with $U(-1,1)$ with 0.2 probability. The results of GNFC tasks are summarized
in Fig. 3(a) and the detailed results can be found in Tab. 8. The results show
that the property of the behavior policy (i.e., $e$ and $p$) dominates the
generalization ability of the baseline algorithms. When $e=0.05$, almost all
of the baselines fail and give a completely opposite dosage-response curve
(see Fig. 4(a) and Appendix. G.2). IPW still perform well when $0.2\leq e\leq
1.0$ but fails when $e=0.05,p<=0.2$. We also found SCIGAN can reach a better
performance than other baselines when $e=0.05,p<=0.2$, but the results in
other tasks are unstable. GALILEO is the only algorithm that is robust to the
selection bias and outputs correct dosage-response curves in all of the tasks.
Based on the experiment, we also indicate that the commonly used overlap
assumption is unreasonable to a certain extent especially in real-world
applications since it is impractical to inject noises into the whole action
space. The problem of overlap assumption being violated should be taken into
consideration otherwise it will be hard to use the algorithm in practice if it
is sensitive to the noise range.

(a) GNFC

(b) TCGA

Figure 3: Illustration of the performance in GNFC and TCGA. The grey bar
denotes the standard error ($\times 0.3$ for brevity) of $3$ random seeds.

##### Synthetic benchmark environment.

(a) e0.05_p0.2 (b) t0_bias20

Figure 4: Illustration of the averaged dosage-response curves.

Previous experiments on counterfactual environment model learning are based on
semi-synthetic data simulation [24, 28]. Since GALILEO is compatible with one-
step environment model learning, we select the same task named TCGA in [24] to
test the variant GALILEO(0). Based on three synthetic dosage-response
functions in TCGA, we construct 9 tasks by choosing different parameters of
selection bias on $\mu$ which is constructed with beta distribution and design
a coefficient $c$ to control the selection bias of the beta distribution. We
name the tasks with the format of “t?_bias?”. For example, t1_bias2 is the
task with the first dosage-response functions and $c=2$. The detail of TCGA is
in Appx. F.1.2. The results of TCGA tasks are summarized in Fig. 3(b) and the
detailed results can be found in Tab. 7 in Appendix. We found the phenomenon
in this experiment is similar to the one in GNFC, which demonstrates the
compatibility of GALILEO to one-step environments. We also found that the
results of IPW are unstable in this experiment. It might be because the
behavior policy is modeled with beta distribution while the propensity score
$\hat{\mu}$ is modeled with Gaussian distribution. Since IPW directly reweight
loss function with $\frac{1}{\hat{\mu}}$, the results are sensitive to the
error on $\hat{\mu}$. GALILEO also models $\hat{\mu}$ with Gaussian
distribution. The results of GALILEO are more stable than IPW since GALILEO
does not re-weight through $\hat{\mu}$ directly.

##### Dosage-response curve visualization.

We plot the averaged dosage-response curves which are constructed by
equidistantly sampling action from the action space and averaging the feedback
of the states in the dataset as the averaged response. We sample $9$ points in
GNFC tasks and $33$ points in TAGC tasks and give parts of the visualization
results in Fig. 4 (all curves can be seen in Appx. G.2). For those tasks where
baselines fails in reconstructing dosage-response curves, GALILEO not only
reaches a better MISE score but reconstructs almost exact responses.

##### Worst-case prediction error.

For the better capacity of environment models for policy optimization, we
propose a new objective ACRM. In theory, GALILEO increases the generalization
ability by focusing on the worst-case samples’ training to achieve ACRM. To
demonstrate the property, we propose a new metric named Mean-Max Square Error
(MMSE):
$\mathbb{E}\left[\max_{a\in\mathcal{A}}\left(M^{*}(x^{\prime}|x,a)-M(x^{\prime}|x,a)\right)^{2}\right]$.
In Appx. G.3, we give the MMSE performance of GALILEO and the baselines in
GNFC.

###  5.2 Ablation Studies

In Sec. 4.3, we introduce several techniques to develop a practical GALILEO
algorithm. Based on task e0.2_p0.05 in GNFC, we give the ablation studies to
investigate the effects of these techniques. We first compare with two
variants that do not handle the assumptions violation problems: (1)
NO_INJECT_NOISE: set $\epsilon_{\mu}$ and $\epsilon_{D}$ to zero, which makes
the overlap assumption is not satisfied; (2) SINGLE_SL: without replacing the
second term in Eq.(8) with standard supervised learning even when the output
probability of $D$ is far away from $0.5$. Besides, we introduced several
tricks inspired by GAIL and give a comparison of these tricks and GAIL: (3)
ONE_STEP: use one-step reward instead of cumulative rewards (i.e., Q and V)
for re-weighting, which is implemented by set $\gamma$ to $0$; (4) SINGE_DIS:
remove $T_{\varphi^{*}_{1}}(x,a)$ and replace it with
$\mathbb{E}_{M_{\theta}}\left[T_{\varphi^{*}_{0}}(x,a,x^{\prime})\right]$,
which is inspired by GAIL that uses value function as a baseline instead of
using another discriminator; (5) PURE_GAIL: remove the second term in Eq.(8).
It can be regarded as a naive adoption of GAIL and a partial implementation of
GALILEO.

Figure 5: Illustration of the ablation studies. The error bars are the
standard error.

We summarize the results in Fig. 5. Based on the results of NO_INJECT_NOISE
and SINGLE_SL, we can see that handling the assumption violation problems is
important and will increase the ability on counterfactual queries. The results
of PURE_GAIL tell us that the partial implementation of GALILEO is not enough
to give stable predictions on counterfactual data; On the other hand, the
result of ONE_STEP also demonstrate that embedding the cumulative error of
one-step prediction is helpful for GALILEO training; Finally, we also found
that SINGLE_DIS nearly has almost no effect on the results. It suggests that,
empirically, we can use
$\mathbb{E}_{M_{\theta}}\left[T_{\varphi^{*}_{0}}(x,a,x^{\prime})\right]$ as a
replacement of $T_{\varphi^{*}_{1}}(x,a)$, which can reduce the computation
costs of the extra discriminator training.

###  5.3 Test in Complex Environments

##### MuJoCo tasks.

MuJoCo is a benchmark task in Gym [42, 43] where we need to control a robot
with specific dynamics to complete some tasks (e.g., standing or running). We
select 3 datasets from D4RL [45] to construct our model learning tasks. We
compare it with a typical transition model learning algorithm used in previous
offline model-based RL algorithm [6, 9], which is a variant of supervised
learning. We name the method OFF-SL. We train models in datasets
Walker2d-medium, Hopper-medium, and HalfCheetah-medium, which are collected by
a behavior policy with 1/3 performance to the expert policy, then we test them
in the corresponding expert dataset.

Table 1: Results on MuJoCo tasks measured with mean square error. We bold the lowest error for each task. task | HalfCheetah | Hopper | Walker2d
---|---|---|---
dataset | training (medium) | test (expert) | training (medium) | test (expert) | training (medium) | test (expert)
GALILEO | 0.12 | 9.26 | 0.0002 | 0.07 | 0.24 | 1.79
OFF-SL | 0.04 | 25.76 | 0.0002 | 0.29 | 0.13 | 2.03

As shown in Tab. 1, although the training errors of both the algorithms have
converged to the same order of magnitude, in all of the tasks, GALILEO reaches
significantly better performances in the expert dataset: the declines of error
are 64%, 76%, and 12%. However, whether GALILEO or OFF-SL, the performance in
expert datasets is as least 5x worse than in the training dataset. The
phenomenon indicates that although GALILEO can make better performance for
counterfactual queries, the risks of using the models are still large and
still challenging to be further solved.

(a) dosage response curve

(b) AUUC

Figure 6: Illustration of the result for BAT.

##### Real-world platform.

We deploy GALILEO in a real-world large-scale food-delivery platform. The goal
of the platform is to balance the demand from take-out food orders and the
supply of delivery men, i.e., helping delivery men finish more orders by
giving reasonable strategies. We focus on a Budget Allocation task to the Time
period (BAT) in the platform (see Appx. F.1.3 for details).

Since there has no oracle environment model for querying, we describe the
results with other metrics. First, we review whether the tendency of the
dosage-response curve is compliant with our commonsense. In this application,
our commonsense is that a larger budget of allowance does not decrease the
supply. As can be seen in Fig. 6(a), the tendency of GALILEO’s response is
compliant but the model of SL gives an opposite direction to the response (All
results are in Fig. 9). If we learn a policy through the model of SL, the
optimal solution is canceling all the allowances, which is obviously incorrect
in practice.

Second, we conduct randomized controlled trials (RCT) in one of the testing
cities. Note that collecting RCT samples is expensive, and should have been
avoided if not for the experiment purpose. Using the RCT samples, we can
evaluate the correctness of the sort order of the model predictions via Area
Under the Uplift Curve (AUUC) [46]. AUUC is computed by sorting the RCT
samples based on the predicted effects in the treatment and control groups.
Then the cumulative treatment effects are computed by scanning the sorted
sample list according to the model predictions. A good model predicts
treatment effects accurately, and sorts samples with larger treatment effects
ahead of the smaller ones, leading to a large AUUC above the random sorting
strategy. Fig. 6(b) shows that GALILEO achieves a larger AUUC above the random
baseline, while the model by supervised learning is even worse than the random
baseline.

We search for the optimal allowance policy via the cross-entropy method
planner [47] or any derivative-free method [48] in the learned model. We
tested the policy online for the supply improvement in $6$ cities for 14 days.
The algorithm compared is a human-expert policy working online, which is also
the behavior policy of the offline datasets. We conduct the online A/B tests
for each of the cities. The results are in Tab. 6 in Appendix. The results
show that the policy learned with GALILEO can make better (the supply
improvements are from 0.14 to 1.63 percentage points) budget allocation than
the behavior policies in all the testing cities. Due to its bad AUUC, the
policy learned in the SL model was not able to be tested online, which could
cause a high real risk.

##  6 Discussion and Future Work

In this work, we discuss the generalization issue of world models. We show
that an extremely simple case where the world model can be totally wrong. We
propose ACRM principle to handle the generalization challenges of the
counterfactual environment model learning. By theoretical modeling, we give a
tractable solution to handle ACRM and propose GALILEO. GALILEO is verified in
synthetic environments, complex robot control tasks, and a real-world
platform, and shows great generalization ability on counterfactual queries.

Giving correct answers to counterfactual queries is important for policy
learning. We hope the work can inspire researchers to develop more powerful
tools for counterfactual environment model learning. The current limitation
lies in that there are several simplifications in the theoretical modeling
process (further discussion is in Appx. B), which can be modeled more
elaborately . Besides, experiments on MuJoCo indicate that these tasks are
still challenging to give correct predictions on counterfactual data. These
should also be further investigated in future work.

## References

  * [1] Stephen James and Edward Johns.  3D simulation for robot arm control with deep Q-learning.  CoRR, abs/1609.03759, 2016.
  * [2] Xiong-Hui Chen, Yang Yu, Qingyang Li, Fan-Ming Luo, Zhiwei Tony Qin, Shang Wenjie, and Jieping Ye.  Offline model-based adaptable policy learning.  In Advances in Neural Information Processing Systems 34 (NeurIPS’21), Virtual Conference, 2021.
  * [3] Jing-Cheng Shi, Yang Yu, Qing Da, Shi-Yong Chen, and Anxiang Zeng.  Virtual-taobao: Virtualizing real-world online retail environment for reinforcement learning.  In The 33th AAAI Conference on Artificial Intelligence, pages 4902–4909, Honolulu, Hawaii, 2019.
  * [4] Wenjie Shang, Yang Yu, Qingyang Li, Zhiwei (Tony) Qin, Yiping Meng, and Jieping Ye.  Environment reconstruction with hidden confounders for reinforcement learning based recommendation.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 566–576, Anchorage, AK, 2019.
  * [5] Wenjie Shang, Qingyang Li, Zhiwei Qin, Yang Yu, Yiping Meng, and Jieping Ye.  Partially observable environment estimation with uplift inference for reinforcement learning based recommendation.  Machine Learning, (9):2603–2640, 2021.
  * [6] Tianhe Yu, Garrett Thomas, Lantao Yu, Stefano Ermon, James Y. Zou, Sergey Levine, Chelsea Finn, and Tengyu Ma.  MOPO: model-based offline policy optimization.  In Advances in Neural Information Processing Systems 33, virtual event, 2020.
  * [7] Sergey Levine, Aviral Kumar, George Tucker, and Justin Fu.  Offline reinforcement learning: Tutorial, review, and perspectives on open problems.  CoRR, abs/2005.01643, 2020.
  * [8] Fan-Ming Luo, Shengyi Jiang, Yang Yu, Zongzhang Zhang, and Yi-Feng Zhang.  Adapting environment sudden changes by learning context sensitive policy.  In Proceedings of the 36th AAAI Conference on Artificial Intelligence (AAAI’22), Virtual Conference, 2022.
  * [9] Rahul Kidambi, Aravind Rajeswaran, Praneeth Netrapalli, and Thorsten Joachims.  MOReL : Model-based offline reinforcement learning.  CoRR, abs/2005.05951, 2020.
  * [10] Guido Imbens.  The role of the propensity score in estimating dose-response functions.  Econometrics eJournal, 1999.
  * [11] Ahmed M. Alaa and Mihaela van der Schaar.  Limits of estimating heterogeneous treatment effects: Guidelines for practical algorithm design.  In Proceedings of the 35th International Conference on Machine Learning, volume 80, pages 129–138, Stockholm, Sweden, 2018.
  * [12] Ahmed M. Alaa and Mihaela van der Schaar.  Bayesian inference of individualized treatment effects using multi-task gaussian processes.  In Advances in Neural Information Processing Systems 30, pages 3424–3432, Long Beach, CA, 2017.
  * [13] Jeffrey A. Smith and Petra E. Todd.  Does matching overcome lalonde’s critique of nonexperimental estimators?  Journal of Econometrics, 125:305–353, 2001.
  * [14] Ahmed M. Alaa, Michael Weisz, and Mihaela van der Schaar.  Deep counterfactual networks with propensity-dropout.  CoRR, abs/1706.05966, 2017.
  * [15] Min Qian and Susan A. Murphy.  Performance guarantees for individualized treatment rules.  Annals of statistics, 39 2:1180–1210, 2011.
  * [16] Uri Shalit, Fredrik D. Johansson, and David A. Sontag.  Estimating individual treatment effect: generalization bounds and algorithms.  In Proceedings of the 34th International Conference on Machine Learning, volume 70, pages 3076–3085, NSW, Australia, 2017.
  * [17] Fredrik D. Johansson, Uri Shalit, and David A. Sontag.  Learning representations for counterfactual inference.  In Proceedings of the 33nd International Conference on Machine Learning, volume 48, pages 3020–3029, New York, NY, 2016.
  * [18] Liuyi Yao, Sheng Li, Yaliang Li, Mengdi Huai, Jing Gao, and Aidong Zhang.  Representation learning for treatment effect estimation from observational data.  In Advances in Neural Information Processing Systems 31, pages 2638–2648, Montréal, Canada, 2018.
  * [19] Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron C. Courville, and Yoshua Bengio.  Generative adversarial nets.  In Advances in Neural Information Processing Systems 27, pages 2672–2680, Montréal, Canada, 2014.
  * [20] Sebastian Nowozin, Botond Cseke, and Ryota Tomioka.  f-GAN: Training generative neural samplers using variational divergence minimization.  In Advances in Neural Information Processing Systems 29, pages 271–279, Barcelona, Spain, 2016.
  * [21] Xiong-Hui Chen, Shengyi Jiang, Feng Xu, Zongzhang Zhang, and Yang Yu.  Cross-modal domain adaptation for cost-efficient visual reinforcement learning.  In Marc’Aurelio Ranzato, Alina Beygelzimer, Yann N. Dauphin, Percy Liang, and Jennifer Wortman Vaughan, editors, Advances in Neural Information Processing Systems 34, pages 12520–12532, virtual event, 2021.
  * [22] Jonathan Ho and Stefano Ermon.  Generative adversarial imitation learning.  In Advances in Neural Information Processing Systems 29, pages 4565–4573, Barcelona, Spain, 2016.
  * [23] Jinsung Yoon, James Jordon, and Mihaela van der Schaar.  GANITE: estimation of individualized treatment effects using generative adversarial nets.  In 6th International Conference on Learning Representations, 2018\.
  * [24] Ioana Bica, James Jordon, and Mihaela van der Schaar.  Estimating the effects of continuous-valued interventions using generative adversarial networks.  In Advances in Neural Information Processing Systems 33, virtual event, 2020.
  * [25] Xinshi Chen, Shuang Li, Hui Li, Shaohua Jiang, Yuan Qi, and Le Song.  Generative adversarial user model for reinforcement learning based recommendation system.  In Proceedings of the 36th International Conference on Machine Learning, pages 1052–1061, Long Beach, CA, 2019.
  * [26] Paul R. Rosenbaum and Donald B. Rubin.  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70:41–55, 1983.
  * [27] Uri Shalit, Fredrik D. Johansson, and David A. Sontag.  Estimating individual treatment effect: generalization bounds and algorithms.  In Proceedings of the 34th International Conference on Machine Learning, volume 70, pages 3076–3085, Sydney, Australia, 2017.
  * [28] Patrick Schwab, Lorenz Linhardt, Stefan Bauer, Joachim M. Buhmann, and Walter Karlen.  Learning counterfactual representations for estimating individual dose-response curves.  In The 34th AAAI Conference on Artificial Intelligence, pages 5612–5619, New York, NY, 2020.
  * [29] Hidetoshi Shimodaira.  Improving predictive inference under covariate shift by weighting the log-likelihood function.  Journal of statistical planning and inference, 90(2):227–244, 2000\.
  * [30] Fredrik D. Johansson, Uri Shalit, Nathan Kallus, and David A. Sontag.  Generalization bounds and representation learning for estimation of potential outcomes and causal effects.  CoRR, abs/2001.07426, 2020.
  * [31] Fredrik D Johansson, Nathan Kallus, Uri Shalit, and David Sontag.  Learning weighted representations for generalization across designs.  arXiv preprint arXiv:1802.08598, 2018.
  * [32] Adith Swaminathan and Thorsten Joachims.  Counterfactual risk minimization: Learning from logged bandit feedback.  In Proceedings of the 32nd International Conference on Machine Learning, volume 37, pages 814–823, Lille, France, 2015.
  * [33] Stefan Wager and Susan Athey.  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association, 113(523):1228–1242, 2018.
  * [34] Jeremy C. Weiss, Finn Kuusisto, Kendrick Boyd, Jie Liu, and David Page.  Machine learning for treatment assignment: Improving individualized risk attribution.  In American Medical Informatics Association Annual Symposium, San Francisco, CA, 2015.
  * [35] Negar Hassanpour and Russell Greiner.  Counterfactual regression with importance sampling weights.  In Proceedings of the 28th International Joint Conference on Artificial Intelligence, pages 5880–5887, Macao, China, 2019.
  * [36] Uri Shalit, Fredrik D. Johansson, and David A. Sontag.  Estimating individual treatment effect: generalization bounds and algorithms.  In Proceedings of the 34th International Conference on Machine Learning, volume 70, pages 3076–3085, Sydney, Australia, 2017.
  * [37] Richard S. Sutton and Andrew G. Barto.  Reinforcement Learning: An Introduction.  2018\.
  * [38] J.B. Hiriart-Urruty and C. Lemaréchal.  Fundamentals of Convex Analysis.  2001\.
  * [39] Dean Pomerleau.  Efficient training of artificial neural networks for autonomous navigation.  Neural Computation, 3(1):88–97, 1991.
  * [40] John Schulman, Sergey Levine, Pieter Abbeel, Michael I. Jordan, and Philipp Moritz.  Trust region policy optimization.  In Proceedings of the 32nd International Conference on Machine Learning, volume 37, pages 1889–1897, Lille, France, 2015.
  * [41] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms.  CoRR, abs/1707.06347, 2017.
  * [42] Emanuel Todorov, Tom Erez, and Yuval Tassa.  MuJoCo: A physics engine for model-based control.  In Proceedings of the 24th IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 5026–5033, Vilamoura, Portugal, 2012.
  * [43] Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba.  Openai gym.  CoRR, abs/1606.01540, 2016.
  * [44] Peter Spirtes.  Introduction to causal inference.  J. Mach. Learn. Res., 11:1643–1662, aug 2010.
  * [45] Justin Fu, Aviral Kumar, Ofir Nachum, George Tucker, and Sergey Levine.  D4RL: datasets for deep data-driven reinforcement learning.  CoRR, abs/2004.07219, 2020.
  * [46] Artem Betlei, Eustache Diemert, and Massih-Reza Amini.  Treatment targeting by AUUC maximization with generalization guarantees.  CoRR, abs/2012.09897, 2020.
  * [47] Danijar Hafner, Timothy P. Lillicrap, Ian Fischer, Ruben Villegas, David Ha, Honglak Lee, and James Davidson.  In Proceedings of the 36th International Conference on Machine Learning, volume 97, pages 2555–2565, Long Beach, CA, 2019.
  * [48] Yu-Ren Liu, Yi-Qi Hu, Hong Qian, Chao Qian, and Yang Yu.  ZOOpt: Toolbox for derivative-free optimization.  SCIENCE CHINA Information Sciences, 2022.
  * [49] Yao Liu, Adith Swaminathan, Alekh Agarwal, and Emma Brunskill.  Provably good batch off-policy reinforcement learning without great exploration.  In Advances in Neural Information Processing Systems 33, virtual event, 2020.
  * [50] Lixin Zou, Long Xia, Pan Du, Zhuo Zhang, Ting Bai, Weidong Liu, Jian-Yun Nie, and Dawei Yin.  Pseudo dyna-q: A reinforcement learning framework for interactive recommendation.  In Proceedings of the 13th ACM International Conference on Web Search and Data Mining, pages 816–824, Houston, TX, 2020.
  * [51] Jin Huang, Harrie Oosterhuis, Maarten de Rijke, and Herke van Hoof.  Keeping dataset biases out of the simulation: A debiased simulator for reinforcement learning based recommender systems.  In Proceedings of the 14th. ACM Conference on Recommender Systems, pages 190–199, virtual event, 2020.
  * [52] Junda Wu, Zhihui Xie, Tong Yu, Qizhi Li, and Shuai Li.  Sim-to-real interactive recommendation via off-dynamics reinforcement learning.  In 2rd Offline Reinforcement Learning Workshop Advances at NeurIPS, 2021.

##  Appendix A Proof of Theoretical Results

In the proof section, we replace the notation of $\mathbb{E}$ with an integral
for brevity. Now we rewrite the original objective
$\bar{L}(\rho^{\beta}_{M^{*}},M)$ as:

$\displaystyle\min_{M\in\mathcal{M}}\max_{\beta\in\Pi}\int_{\mathcal{X},\mathcal{A}}\rho_{M^{*}}^{\mu}(x,a)\omega(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\left(-\log
M(x^{\prime}|x,a)\right)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x-\frac{\alpha}{2}\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2},$
(9)

where $\omega(x,a)=\frac{\rho_{M^{*}}^{\beta}(x,a)}{\rho_{M^{*}}^{\mu}(x,a)}$
and
$\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}=\int_{\mathcal{X},\mathcal{A}}\rho_{M^{*}}^{\beta}(x,a)^{2}\mathrm{d}a\mathrm{d}x$,
which is the squared $l_{2}$-norm. In an MDP, given any policy $\pi$,
$\rho_{M^{*}}^{\pi}(x,a,x^{\prime})=\rho_{M^{*}}^{\pi}(x)\pi(a|x)M^{*}(x^{\prime}|x,a)$
where $\rho_{M^{*}}^{\pi}(x)$ denotes the occupancy measure of $x$ for policy
$\pi$, which can be defined [37, 22] as
$\rho_{M^{*}}^{\pi}(x):=(1-\gamma)\mathbb{E}_{x_{0}\sim\rho_{0}}\left[\sum^{\infty}_{t=0}\gamma^{t}{\rm
Pr}(x_{t}=x|x_{0},{M^{*}})\right]$ where ${\rm
Pr}^{\pi}\left[x_{t}=x|x_{0},{M^{*}}\right]$ is the state visitation
probability that $\pi$ starts at state $x_{0}$ in model $M^{*}$ and receive
$x$ at timestep $t$ and $\gamma\in\left[0,1\right]$ is the discount factor.

We first prove Lemma 4.3 as follows:

##### Proof of Lemma 4.3

For better readability, we first rewrite Lemma 4.3 as follows:

######  Lemma A.1 (Approximated optimal adversarial distribution).

Given any $M$ in $\bar{L}(\rho^{\beta}_{M^{*}},M)$, the distribution of the
best-response policy $\bar{\beta}^{*}$ satisfies:

$\rho_{M^{*}}^{\bar{\beta}^{*}}(x,a)=\frac{1}{\alpha}(D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))+H_{M^{*}}(x,a)),$
(10)

where $D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))$ is the Kullback-Leibler (KL)
divergence between $M^{*}(\cdot|x,a)$ and $M(\cdot|x,a)$, $H_{M^{*}}(x,a)$
denotes the entropy of $M^{*}(\cdot|x,a)$, and $\alpha$ is a constant
normalizer.

###### Proof.

Given a transition function $M$ of an MDP, the distribution of the best-
response policy $\bar{\beta}^{*}$ satisfies:

$\displaystyle\rho_{M^{*}}^{\bar{\beta}^{*}}=$
$\displaystyle\arg\max_{\rho_{M^{*}}^{\beta}}\int_{\mathcal{X},\mathcal{A}}\rho_{M^{*}}^{\mu}(x,a)\omega(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\left(-\log
M(x^{\prime}|x,a)\right)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x-\frac{\alpha}{2}\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}$
$\displaystyle=$
$\displaystyle\arg\max_{\rho_{M^{*}}^{\beta}}\int_{\mathcal{X},\mathcal{A}}\rho_{M^{*}}^{\beta}(x,a)\underbrace{\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\left(-\log
M(x^{\prime}|x,a)\right)\mathrm{d}x^{\prime}}_{g(x,a)}\mathrm{d}a\mathrm{d}x-\frac{\alpha}{2}\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}$
$\displaystyle=$
$\displaystyle\arg\max_{\rho_{M^{*}}^{\beta}}\frac{2}{\alpha}\int_{\mathcal{X},\mathcal{A}}\rho_{M^{*}}^{\beta}(x,a)g(x,a)\mathrm{d}a\mathrm{d}x-\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}$
$\displaystyle=$
$\displaystyle\arg\max_{\rho_{M^{*}}^{\beta}}\frac{2}{\alpha}\int_{\mathcal{X},\mathcal{A}}\rho_{M^{*}}^{\beta}(x,a)g(x,a)\mathrm{d}a\mathrm{d}x-\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}-\frac{\|g(\cdot,\cdot)\|^{2}_{2}}{\alpha^{2}}$
$\displaystyle=$
$\displaystyle\arg\max_{\rho_{M^{*}}^{\beta}}-\left(-2\int_{\mathcal{X},\mathcal{A}}\rho_{M^{*}}^{\beta}(x,a)\frac{g(x,a)}{\alpha}\mathrm{d}a\mathrm{d}x+\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}+\frac{\|g(\cdot,\cdot)\|^{2}_{2}}{\alpha^{2}}\right)$
$\displaystyle=$
$\displaystyle\arg\max_{\rho^{\beta}_{M^{*}}}-\|\rho^{\beta}_{M^{*}}(\cdot,\cdot)-\frac{g(\cdot,\cdot)}{\alpha}\|^{2}_{2}.$

Suppose $\rho^{\beta}_{M^{*}}$ is a map that
$\rho^{\beta}_{M^{*}}:\mathcal{X}\times\mathcal{A}\rightarrow\mathbb{R}$
without constraints, then we have:

$\displaystyle\rho^{\bar{\beta}^{*}}_{M^{*}}(x,a)$
$\displaystyle=\frac{g(x,a)}{\alpha}$
$\displaystyle=\frac{\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log\frac{M^{*}(x^{\prime}|x,a)}{M(x^{\prime}|x,a)}\mathrm{d}x-\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M^{*}(x^{\prime}|x,a)\mathrm{d}x}{\alpha}$
$\displaystyle=\frac{D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))+H_{M^{*}}(x,a)}{\alpha}$
$\displaystyle\propto\Big{(}D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))+H_{M^{*}}(x,a)\Big{)},$
(11)

where $H_{M^{*}}(x,a)=\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)(-\log
M^{*}(x^{\prime}|x,a))\mathrm{d}x^{\prime}$. However, $\rho^{\beta}_{M^{*}}$
is an occupancy measure instead of an arbitrary map function. Assuming the
transition function $M^{*}$ has an upper bound $c$ on the transition
probabilities, that is $0\leq M^{*}(\cdot|x,a)\leq c,\forall
a\in\mathcal{A},\forall x\in\mathcal{X}$, we have
$\rho^{\beta}_{M^{*}}(x,a)\in[0,c],\forall a\in\mathcal{A},\forall
x\in\mathcal{X}$ based on the definition of the occupancy measure. Let
$\alpha$ be a constant value that
$\alpha>\max_{x,a}g(x,a)=\frac{\max_{x,a}(D_{KL}(M^{*},M)+H_{M^{*}})}{c}$,
then we have $\frac{g(x,a)}{\alpha}\in[0,c),\forall a\in\mathcal{A},\forall
x\in\mathcal{X}$ and Eq.(11) holds. ∎

######  Remark A.2.

The optimal solution Eq.(11) relies on
$\rho^{\beta}_{M^{*}}(x,a)\in[0,c],\forall a\in\mathcal{A},\forall
x\in\mathcal{X}$. In some particular $M^{*}$, it is intractable to derive a
$\beta$ that can generate an occupancy specified by $g(x,a)/\alpha$. Consider
the following case: a state $x_{1}$ in $M^{*}$ might be harder to reach than
another state $x_{2}$, e.g., $M^{*}(x_{1}|x,a)<M^{*}(x_{2}|x,a),\forall
x\in\mathcal{X},\forall a\in\mathcal{A}$, then it is impossible to find a
$\beta$ that the occupancy satisfies
$\rho^{\beta}_{M^{*}}(x_{1},a)>\rho^{\beta}_{M^{*}}(x_{2},a)$. In this case,
Eq.(11) can be a sub-optimal solution. Since this work focuses on task-
agnostic solution derivation while the solution to the above problem should
rely on the specific description of $M^{*}$, we leave it as future work.
However, we point out that Eq.(11) is a reasonable re-weighting term even as a
sub-optimum: $\rho_{M^{*}}^{\bar{\beta}^{*}}$ gives larger densities on the
data where the distribution distance between the approximation model and the
real model (i.e., $D_{KL}(M^{*},M)$) is larger or the stochasticity of the
real model (i.e., $H_{M^{*}}$) is larger.

##### Proof of Lemma 4.4

We know that KL divergence is an instance of a divergence class called
$f$-divergence [20], which measures the difference between two probability
distributions. We give a generalized representation of
$\rho_{M^{*}}^{\bar{\beta}^{*}}$ by replacing $D_{KL}$ in Eq.(11) with
$f$-divergence and converting the distribution metric from conditional
probabilities to joint probabilities. The latter is for easier estimating in
practice.

We first introduce the definition of $f$-divergence.

######  Definition A.3 ($f$ -divergence).

Given two distributions $P$ and $Q$, two absolutely continuous density
functions $p$ and $q$ with respect to a base measure $\mathrm{d}x$ defined on
the domain $\mathcal{X}$, we define the $f$-divergence [20],

$\displaystyle
D_{f}(P\|Q)=\int_{\mathcal{X}}q(x)f\left(\frac{p(x)}{q(x)}\right)\mathrm{d}x,$
(12)

where the generator function $f:\mathbb{R}_{+}\rightarrow\mathbb{R}$ is a
convex, lower-semicontinuous function.

Now we prove Lemma A.4 based on Def. A.3. For better readability, we first
rewrite Lemma 4.4 as follows:

######  Lemma A.4 (Generalized optimal adversarial distribution of
$\bar{\beta}^{*}$).

Replacing $D_{KL}$ in Lemma A.1 with $f$-divergence, we have

$\displaystyle\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}=\frac{1}{\alpha_{0}(x,a)}\left(\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})f\left(\frac{\rho^{\kappa}_{M}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\mathrm{d}x^{\prime}-\rho^{\kappa}_{M^{*}}(x,a)\left(f\left(\frac{\rho^{\kappa}_{M}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-H_{M^{*}}(x,a)\right)\right).$

We introduce a new policy $\kappa(a|x)$ to construct joint probabilities
$\rho^{\kappa}_{M^{*}}(x,a)$ and $\rho^{\kappa}_{M}(x,a)$. $\kappa$ is a
policy satisfying $\forall x\in\mathcal{X},\forall
a\in\mathcal{A},\rho_{M^{*}}^{\kappa}(x,a)>0$ if
$\rho_{M^{*}}^{\bar{\beta}^{*}}(x,a)>0$. Here we call $\kappa$ an intermediary
policy. $f:\mathbb{R}_{+}\rightarrow\mathbb{R}$ is a generator function of
$f$-divergence [20]. $\alpha_{0}(x,a)=\alpha\rho^{\kappa}_{M^{*}}(x,a)$.

###### Proof.

$\displaystyle\rho_{M^{*}}^{\bar{\beta}^{*}}(x,a)=$
$\displaystyle\frac{D_{KL}(M^{*}(\cdot|x,a),M(\cdot|x,a))+H_{M^{*}}(x,a)}{\alpha}$
$\displaystyle=$
$\displaystyle\frac{1}{\alpha}\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\left(\log\frac{M^{*}(x^{\prime}|x,a)}{M(x^{\prime}|x,a)}-\log
M^{*}(x^{\prime}|x,a)\right)\mathrm{d}x^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha\rho^{\kappa}_{M^{*}}(x,a)}\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a)M^{*}(x^{\prime}|x,a)\left(\log\frac{M^{*}(x^{\prime}|x,a)}{M(x^{\prime}|x,a)}-\log
M^{*}(x^{\prime}|x,a)\right)\mathrm{d}x^{\prime}$ (13) $\displaystyle=$
$\displaystyle\frac{1}{\alpha\rho^{\kappa}_{M^{*}}(x,a)}\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})\left(\log\frac{\rho_{M^{*}}^{\kappa}(x,a,x^{\prime})}{\rho_{M}^{\kappa}(x,a,x^{\prime})}+\log\frac{\rho^{\kappa}_{M}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}-\log
M^{*}(x^{\prime}|x,a)\right)\mathrm{d}x^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha\rho^{\kappa}_{M^{*}}(x,a)}\Bigg{(}\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})\log\frac{\rho_{M^{*}}^{\kappa}(x,a,x^{\prime})}{\rho_{M}^{\kappa}(x,a,x^{\prime})}\mathrm{d}x^{\prime}-$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)\log\frac{\rho^{\kappa}_{M^{*}}(x,a)}{\rho^{\kappa}_{M}(x,a)}\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\mathrm{d}x^{\prime}-\rho^{\kappa}_{M^{*}}(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M^{*}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\Bigg{)}$ (14)
$\displaystyle\rightarrow$
$\displaystyle\frac{1}{\alpha\rho^{\kappa}_{M^{*}}(x,a)}\Bigg{(}\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})f\left(\frac{\rho^{\kappa}_{M}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\mathrm{d}x^{\prime}-$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)f\left(\frac{\rho^{\kappa}_{M}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)\times
1-\rho^{\kappa}_{M^{*}}(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M^{*}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\Bigg{)}$ (15) $\displaystyle=$
$\displaystyle\frac{1}{\alpha_{0}(x,a)}\left(\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})f\left(\frac{\rho^{\kappa}_{M}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\mathrm{d}x^{\prime}-\rho^{\kappa}_{M^{*}}(x,a)\left(f\left(\frac{\rho^{\kappa}_{M}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-H_{M^{*}}(x,a)\right)\right).$
(16)

Eq.(13) introduces a new policy $\kappa(a|x)$ to construct joint probabilities
$\rho^{\kappa}_{M^{*}}(x,a)$ and $\rho^{\kappa}_{M}(x,a)$. The equality holds
only when the policy satisfies $\forall x\in X,\forall a\in
A,\rho_{M^{*}}^{\kappa}(x,a)>0$ if $\rho_{M^{*}}^{\bar{\beta}^{*}}(x,a)>0$. We
call the policies that satisfy the above condition intermediary policies.
Eq.(14) includes
$\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})\log\frac{\rho_{M^{*}}^{\kappa}(x,a,x^{\prime})}{\rho_{M}^{\kappa}(x,a,x^{\prime})}$
and
$\rho^{\kappa}_{M^{*}}(x,a)\log\frac{\rho_{M^{*}}^{\kappa}(x,a)}{\rho_{M}^{\kappa}(x,a)}$,
which follow the same pattern $q\log\frac{q}{p}$. $q\log\frac{q}{p}$ can be
regarded as the Reverse-KL divergence instantiation of $f$-divergence (See
Reverse-KL divergence of Tab.1 in [20] for more details). We derive Eq.(15) by
replacing the pattern of $q\log\frac{q}{p}$ with $qf(\frac{p}{q})$. Let
$\alpha_{0}(x,a)=\alpha\rho^{\kappa}_{M^{*}}(x,a)$, then we have Eq.(16).
Eq.(16) is a generalized representation of $\rho_{M^{*}}^{\bar{\beta}^{*}}$,
and we use a new notation $\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}$ to denote it:

$\displaystyle\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}=\frac{1}{\alpha_{0}(x,a)}\left(\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})f\left(\frac{\rho^{\kappa}_{M}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\mathrm{d}x^{\prime}-\rho^{\kappa}_{M^{*}}(x,a)\left(f\left(\frac{\rho^{\kappa}_{M}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-H_{M^{*}}(x,a)\right)\right).$

∎

##### Proof of Thm. 4.5 and the tractable solution

We first introduce several useful lemmas for the proof.

######  Lemma A.5.

Rearrangement inequality The rearrangement inequality states that, for two
sequences $a_{1}\geq a_{2}\geq\ldots\geq a_{n}$ and $b_{1}\geq
b_{2}\geq\ldots\geq b_{n}$, the inequalities

$\displaystyle a_{1}b_{1}+a_{2}b_{2}+\cdots+a_{n}b_{n}\geq
a_{1}b_{\pi(1)}+a_{2}b_{\pi(2)}+\cdots+a_{n}b_{\pi(n)}\geq
a_{1}b_{n}+a_{2}b_{n-1}+\cdots+a_{n}b_{1}$

hold, where $\pi(1),\pi(2),\ldots,\pi(n)$ is any permutation of
$1,2,\ldots,n$.

######  Lemma A.6.

For two sequences $a_{1}\geq a_{2}\geq\ldots\geq a_{n}$ and $b_{1}\geq
b_{2}\geq\ldots\geq b_{n}$, the inequalities

$\displaystyle\sum_{i=1}^{n}\frac{1}{n}a_{i}b_{i}\geq\sum_{i=1}^{n}\frac{1}{n}a_{i}\sum\frac{1}{n}b_{i}$

hold.

###### Proof.

By rearrangement inequality, we have

$\displaystyle\sum_{i=1}^{n}a_{i}b_{i}$ $\displaystyle\geq
a_{1}b_{1}+a_{2}b_{2}+\cdots+a_{n}b_{n}$
$\displaystyle\sum_{i=1}^{n}a_{i}b_{i}$ $\displaystyle\geq
a_{1}b_{2}+a_{2}b_{3}+\cdots+a_{n}b_{1}$
$\displaystyle\sum_{i=1}^{n}a_{i}b_{i}$ $\displaystyle\geq
a_{1}b_{3}+a_{2}b_{4}+\cdots+a_{n}b_{2}$ $\displaystyle\vdots$
$\displaystyle\sum_{i=1}^{n}a_{i}b_{i}$ $\displaystyle\geq
a_{1}b_{n}+a_{2}b_{1}+\cdots+a_{n}b_{n-1}$

Then we have

$\displaystyle n\sum_{i=1}^{n}a_{i}b_{i}$
$\displaystyle\geq\sum_{i=1}^{n}a_{i}\sum_{i=1}^{n}b_{i}$
$\displaystyle\sum_{i=1}^{n}\frac{1}{n}a_{i}b_{i}$
$\displaystyle\geq\sum_{i=1}^{n}\frac{1}{n}a_{i}\sum\frac{1}{n}b_{i}$

∎

Now we extend Lemma A.6 into the continuous integral scenario:

######  Lemma A.7.

Given $\mathcal{X}\subset\mathbb{R}$, for two functions
$f:\mathcal{X}\to\mathbb{R}$ and $g:\mathcal{X}\to\mathbb{R}$ that $f(x)\geq
f(y)$ if and only if $g(x)\geq g(y),\ \forall x,y\in\mathcal{X}$, the
inequality

$\displaystyle\int_{\mathcal{X}}p(x)f(x)g(x)\mathrm{d}x\geq\int_{\mathcal{X}}p(x)f(x)\mathrm{d}x\int_{\mathcal{X}}p(x)g(x)\mathrm{d}x$

holds, where $p:\mathcal{X}\to\mathbb{R}$ and $p(x)>0,\forall x\in\mathcal{X}$
and $\int_{\mathcal{X}}p(x)\mathrm{d}x=1$.

###### Proof.

Since $(f(x)-f(y))(g(x)-g(y))\geq 0,\forall x,y\in\mathcal{X}$, we have

$\displaystyle\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)(f(x)-f(y))(g(x)-g(y))\mathrm{d}y\mathrm{d}x\geq
0$
$\displaystyle\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(x)+p(x)p(y)f(y)g(y)-p(x)p(y)f(x)g(y)-p(x)p(y)f(y)g(x)\mathrm{d}y\mathrm{d}x\geq
0$
$\displaystyle\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(x)+p(x)p(y)f(y)g(y)\mathrm{d}y\mathrm{d}x\geq\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(y)+p(x)p(y)f(y)g(x)\mathrm{d}y\mathrm{d}x$
$\displaystyle\int_{x\in\mathcal{X}}\left(\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(x)\mathrm{d}y+\int_{y\in\mathcal{X}}p(x)p(y)f(y)g(y)\mathrm{d}y\right)\mathrm{d}x\geq\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(y)+p(x)p(y)f(y)g(x)\mathrm{d}y\mathrm{d}x$
$\displaystyle\int_{x\in\mathcal{X}}\left(p(x)f(x)g(x)+\int_{y\in\mathcal{X}}p(x)p(y)f(y)g(y)\mathrm{d}y\right)\mathrm{d}x\geq\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(y)+p(x)p(y)f(y)g(x)\mathrm{d}y\mathrm{d}x$
$\displaystyle\int_{x\in\mathcal{X}}p(x)f(x)g(x)\mathrm{d}x+\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(y)g(y)\mathrm{d}y\mathrm{d}x\geq\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(y)+p(x)p(y)f(y)g(x)\mathrm{d}y\mathrm{d}x$
$\displaystyle\int_{x\in\mathcal{X}}p(x)f(x)g(x)\mathrm{d}x+\int_{y\in\mathcal{X}}p(y)f(y)g(y)\mathrm{d}y\geq\int_{x\in\mathcal{X}}\int_{y\in\mathcal{X}}p(x)p(y)f(x)g(y)+p(x)p(y)f(y)g(x)\mathrm{d}y\mathrm{d}x$
$\displaystyle 2$
$\displaystyle\int_{x\in\mathcal{X}}p(x)f(x)g(x)\mathrm{d}x\geq
2\int_{y\in\mathcal{X}}\int_{x\in\mathcal{X}}p(x)p(y)f(x)g(y)\mathrm{d}y\mathrm{d}x$
$\displaystyle 2$
$\displaystyle\int_{x\in\mathcal{X}}p(x)f(x)g(x)\mathrm{d}x\geq
2\int_{x\in\mathcal{X}}p(x)f(x)\mathrm{d}x\int_{x\in\mathcal{X}}p(x)g(x)\mathrm{d}x$
$\displaystyle\int_{x\in\mathcal{X}}p(x)f(x)g(x)\mathrm{d}x\geq\int_{x\in\mathcal{X}}p(x)f(x)\mathrm{d}x\int_{x\in\mathcal{X}}p(x)g(x)\mathrm{d}x$

∎

######  Corollary A.8.

Let $g(\frac{p(x)}{q(x)})=-\log\frac{p(x)}{q(x)}$ where $p(x)>0,\forall
x\in\mathcal{X}$ and $q(x)>0,\forall x\in\mathcal{X}$, for $\upsilon>0$, the
inequality

$\displaystyle\int_{\mathcal{X}}q(x)f(\upsilon\frac{p(x)}{q(x)})g(\frac{p(x)}{q(x)})\mathrm{d}x\geq\int_{\mathcal{X}}q(x)f(\upsilon\frac{p(x)}{q(x)})\mathrm{d}x\int_{\mathcal{X}}q(x)g(\frac{p(x)}{q(x)})\mathrm{d}x,$

holds if $f^{\prime}(x)\leq 0,\forall x\in\mathcal{X}$. It is not always
satisfied for $f$ functions of $f$-divergence. We list a comparison of $f$ on
that condition in Tab. 2.

###### Proof.

$g^{\prime}(x)=-\log x=-\frac{1}{x}<0,\forall x\in\mathcal{X}$. Suppose
$f^{\prime}(x)\leq 0,\forall x\in\mathcal{X}$, we have $f(x)\geq f(y)$ if and
only if $g(x)\geq g(y),\ \forall x,y\in\mathcal{X}$ holds. Thus
$f(\upsilon\frac{p(x)}{q(x)})\geq f(\upsilon\frac{p(y)}{q(y)})$ if and only if
$g(\frac{p(x)}{q(x)})\geq g(\frac{p(y)}{q(y)}),\ \forall x,y\in\mathcal{X}$
holds for all $\upsilon>0$. By defining $F(x)=f(\upsilon\frac{p(x)}{q(x)}))$
and $G(x)=g(\frac{p(x)}{q(x)})$ and using Lemma A.7, we have:

$\displaystyle\int_{\mathcal{X}}q(x)F(x)G(x)\mathrm{d}x\geq\int_{\mathcal{X}}q(x)F(x)\mathrm{d}x\int_{\mathcal{X}}q(x)G(x)\mathrm{d}x.$

Then we know

$\displaystyle\int_{\mathcal{X}}q(x)f(\upsilon\frac{p(x)}{q(x)})g(\frac{p(x)}{q(x)})\mathrm{d}x\geq\int_{\mathcal{X}}q(x)f(\upsilon\frac{p(x)}{q(x)})\mathrm{d}x\int_{\mathcal{X}}q(x)g(\frac{p(x)}{q(x)})\mathrm{d}x$

holds. ∎

######  Lemma A.9 ($f^{\prime}(\frac{p}{q})$ estimation [20]).

By finding the maximized point of $\varphi$, that is,

$\displaystyle\varphi^{*}=\arg\max_{\varphi}\mathbb{E}_{x\sim
p(x)}\left[T_{\varphi}(x)\right]+\mathbb{E}_{x\sim
q(x)}\left[f^{*}(T_{\varphi}(x))\right],$

we have $f^{\prime}(\frac{p(x)}{q(x)})=T_{\varphi^{*}}(x)$. $f^{*}$ is Fenchel
conjugate of $f$ [20, 38].

Table 2: Properties of $f^{\prime}(x)\leq 0,\forall x\in\mathcal{X}$ for $f$-divergences. Name | Generator function $f(x)$ | If $f^{\prime}(x)\leq 0,\forall x\in\mathcal{X}$
---|---|---
Kullback-Leibler | $x\log x$ | False
Reverse KL | $-\log x$ | True
Pearson $\chi^{2}$ | $(x-1)^{2}$ | False
Squared Hellinger | $(\sqrt{x}-1)^{2}$ | False
Jensen-Shannon | $-(x+1)\log\frac{1+x}{2}+x\log x$ | False
GAN | $x\log x-(x+1)\log(x+1)$ | True

Now, we prove Thm. 4.5. For better readability, we first rewrite Thm. 4.5 as
follows:

######  Theorem A.10.

Let $\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}$ be the generalized representation
of $\rho_{M^{*}}^{\bar{\beta}^{*}}$, where $\bar{\beta}^{*}$ is the best-
response policy with respect to $M_{\theta}$ parameterized by $\theta$, then
we can directly find the optimal $\theta^{*}$ of
$\min_{\theta}\max_{\beta\in\Pi}\bar{L}(\rho_{M^{*}}^{\beta},M_{\theta})$
through the following objective:

$\displaystyle\theta^{*}=\arg\max_{\theta}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\Bigg{[}\frac{1}{\alpha_{0}(x,a)}\log
M_{\theta}(x^{\prime}|x,a)\underbrace{\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)-f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)+H_{M^{*}}(x,a)\right)}_{W(x,a,x^{\prime})}\Bigg{]},$

where $\alpha_{0}(x,a)=\alpha\rho^{\kappa}_{M^{*}}(x,a)$,
$\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[\cdot\right]$ denotes
$\mathbb{E}_{x,a,x^{\prime}\sim\rho^{\kappa}_{M^{*}}}\left[\cdot\right]$, $f$
is the generator function in $f$-divergence which satisfies $f^{\prime}(x)\leq
0,\forall x\in\mathcal{X}$, and $\theta$ is the parameters of $M$.
$M_{\bar{\theta}}$ denotes a probability function with the same parameters as
the learned model (i.e., $\bar{\theta}=\theta$) but the parameter is fixed and
only used for sampling.

###### Proof.

Let $\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}$ in Lemma A.4 as the data
distribution of the best-response policy $\bar{\beta}^{*}$ under model
$M_{\theta}$ parameterized by $\theta$, then we can skip the inner
“$\max_{\beta}[\cdot]$” and directly find the optimal $\theta^{*}$ of
$\min_{\theta}\max_{\beta\in\Pi}\bar{L}(\rho_{M^{*}}^{\beta},M_{\theta})$
(Eq.(9)) through
$\min_{\theta}\bar{L}(\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}},M_{\theta})$:

$\displaystyle\min_{\theta}\int_{\mathcal{X},\mathcal{A}}\bar{\rho}_{M^{*}}^{\bar{\beta}^{*}}(x,a)$
$\displaystyle\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\left(-\log
M_{\theta}(x^{\prime}|x,a)\right)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x$
(17)
$\displaystyle=\min_{\theta}\int_{\mathcal{X},\mathcal{A}}\frac{1}{\alpha_{0}(x,a)}\Bigg{(}$
$\displaystyle\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\mathrm{d}x^{\prime}\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)(-\log
M_{\theta}(x^{\prime}|x,a))\mathrm{d}x^{\prime}$
$\displaystyle-\rho^{\kappa}_{M^{*}}(x,a)\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-H_{M^{*}}(x,a)\right)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)(-\log
M_{\theta}(x^{\prime}|x,a))\mathrm{d}x^{\prime}\Bigg{)}\mathrm{d}a\mathrm{d}x$
$\displaystyle=\min_{\theta}\int_{\mathcal{X},\mathcal{A}}\frac{1}{\alpha_{0}(x,a)}\Bigg{(}$
$\displaystyle\int_{\mathcal{X}}\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\mathrm{d}x^{\prime}\left(\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)(-\log\frac{M_{\theta}(x^{\prime}|x,a)}{M^{*}(x^{\prime}|x,a)})\mathrm{d}x^{\prime}+H_{M^{*}}(x,a)\right)$
$\displaystyle-\rho^{\kappa}_{M^{*}}(x,a)\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-H_{M^{*}}(x,a)\right)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)(-\log
M_{\theta}(x^{\prime}|x,a))\mathrm{d}x^{\prime}\Bigg{)}\mathrm{d}a\mathrm{d}x$
$\displaystyle\leq\min_{\theta}\int_{\mathcal{X},\mathcal{A}}\frac{1}{\alpha_{0}(x,a)}\Bigg{(}$
$\displaystyle\underbrace{\rho^{\kappa}_{M^{*}}(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)(-\log\frac{M_{\theta}(x^{\prime}|x,a)}{M^{*}(x^{\prime}|x,a)})\mathrm{d}x^{\prime}}_{\mathrm{based~{}on~{}Corollary}~{}\ref{the:rearrangement-
corollary}}$
$\displaystyle-\rho^{\kappa}_{M^{*}}(x,a)\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-H_{M^{*}}(x,a)\right)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)(-\log
M_{\theta}(x^{\prime}|x,a))\mathrm{d}x^{\prime}\Bigg{)}\mathrm{d}a\mathrm{d}x$
$\displaystyle=\min_{\theta}\int_{\mathcal{X},\mathcal{A}}\frac{1}{\alpha_{0}(x,a)}\Bigg{(}$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)\int_{\mathcal{X}}\left(M^{*}(x^{\prime}|x,a)f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)(-\log
M_{\theta}(x^{\prime}|x,a))\right)\mathrm{d}x^{\prime}$
$\displaystyle-\rho^{\kappa}_{M^{*}}(x,a)\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-H_{M^{*}}(x,a)\right)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)(-\log
M_{\theta}(x^{\prime}|x,a))\mathrm{d}x^{\prime}\Bigg{)}\mathrm{d}a\mathrm{d}x$
$\displaystyle=\max_{\theta}\int_{\mathcal{X},\mathcal{A},\mathcal{X}}\frac{1}{\alpha_{0}(x,a)}$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})\log
M_{\theta}(x^{\prime}|x,a)\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)-f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)+H_{M^{*}}(x,a)\right)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x,$
(18)

where $\bar{\theta}$ is also the parameters of $M$ and $\bar{\theta}=\theta$.
$M_{\bar{\theta}}$ denotes a transition function with the same parameter as
the learned model $M_{\theta}$ but the parameter is fixed instead of back-
propagating. In Eq.(17), $\|\rho_{M^{*}}^{\beta}(\cdot,\cdot)\|^{2}_{2}$ for
Eq.(9) is eliminated as it does not contribute to the gradient of $\theta$.
Assume $f^{\prime}(x)\leq 0,\forall x\in\mathcal{X}$, let
$\upsilon(x,a)=\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M_{\theta^{*}}}(x,a)}>0$,
$p(x^{\prime}|x,a)=M_{\theta}(x^{\prime}|x,a)$ and
$q(x^{\prime}|x,a)=M^{*}(x^{\prime}|x,a)$, the first inequality can be derived
by adopting Corollary A.8 and eliminating the first $H_{M^{*}}$ since it does
not contribute to the gradient of $\theta$.

∎

##### Proof of the tractable solution

Now we are ready to prove the tractable solution:

###### Proof.

Note that the part $\rho^{\kappa}_{M^{*}}(x,a)$ in
$\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})$ in Eq.(18) can be cancelled because
$\alpha_{0}(x,a)=\alpha\rho^{\kappa}_{M^{*}}(x,a)$, but we choose to keep it
and ignore $\alpha_{0}(x,a)$. The benefit is that we can estimate
$\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})$ from an empirical data distribution
through data collected by $\kappa$ in $M^{*}$ directly, rather than from a
uniform distribution which is harder to be generated. Although keeping
$\rho^{\kappa}_{M^{*}}(x,a)$ incurs extra bias in theory, the results in our
experiments show that it has not made significant negative effects in
practice. We leave this part of modeling in future work. In particular, by
ignoring $\alpha_{0}(x,a)$, we have:

$\displaystyle\hat{\theta}^{*}=\arg\max_{\theta}\int_{\mathcal{X},\mathcal{A},\mathcal{X}}$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})\log
M_{\theta}(x^{\prime}|x,a)\left(f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)-f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)+H_{M^{*}}(x,a)\right)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x.$
(19)

The generator function $f$ is still intractable. In the following, we give a
tractable solution to Eq.(19). First, we resort to the first-order
approximation. Given some $a\in(1-\xi,1+\xi),\xi>0$, we have

$\displaystyle f(a)\approx f(1)+f^{\prime}(a)(a-1),$ (20)

where $f^{\prime}$ is the first-order derivative of $f$. By Taylor’s formula
and the fact that $f^{\prime}(x)$ of the generator function $f$ is bounded in
$(1-\xi,1+\xi)$, the approximation error is no more than
$\mathcal{O}(\xi^{2})$. Substituting $a$ with $\frac{p(x)}{q(x)}$ in Eq.(20),
the pattern $f(\frac{p(x)}{q(x)})$ in Eq.(19) can be converted to
$\frac{p(x)}{q(x)}f^{\prime}(\frac{p(x)}{q(x)})-f^{\prime}(\frac{p(x)}{q(x)})+f(1)$,
then we have:

$\displaystyle\hat{\theta}^{*}=\arg\max_{\theta}\int_{\mathcal{X},\mathcal{A}}\Bigg{(}$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}-$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)f\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}+$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)H_{M^{*}}(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\Bigg{)}\mathrm{d}a\mathrm{d}x$
$\displaystyle\approx\arg\max_{\theta}\int_{\mathcal{X},\mathcal{A}}\Bigg{(}$
$\displaystyle\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)\int_{\mathcal{X}}M_{\bar{\theta}}(x^{\prime}|x,a)f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}-$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\left(f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)-f(1)\right)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}-$
$\displaystyle\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}+$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)\left(f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-f(1)\right)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}+$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a)H_{M^{*}}(x,a)\int_{\mathcal{X}}M^{*}(x^{\prime}|x,a)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\Bigg{)}\mathrm{d}a\mathrm{d}x$
$\displaystyle=\arg\max_{\theta}\int_{\mathcal{X},\mathcal{A},\mathcal{X}}$
$\displaystyle\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x)\left(f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)-f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)\right)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x+$
$\displaystyle\int_{\mathcal{X},\mathcal{A},\mathcal{X}}$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})\left(f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)-f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)+H_{M^{*}}(x,a)\right)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x.$

We can estimate
$f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)$
and
$f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)$
through Lemma A.9. In particular,

$\displaystyle\varphi_{0}^{*}=\arg\max_{\varphi_{0}}~{}$
$\displaystyle\mathbb{E}_{x,a,x^{\prime}\sim\rho^{\kappa}_{M^{*}}}\left[T_{\varphi_{0}}(x,a,x^{\prime})\right]+\mathbb{E}_{x,a,x^{\prime}\sim\rho^{\kappa}_{M_{\bar{\theta}}}}\left[f^{*}(T_{\varphi_{0}}(x,a,x^{\prime}))\right]$
$\displaystyle\varphi_{1}^{*}=\arg\max_{\varphi_{1}}~{}$
$\displaystyle\mathbb{E}_{x,a\sim\rho^{\kappa}_{M^{*}}}\left[T_{\varphi_{1}}(x,a)\right]+\mathbb{E}_{x,a\sim\rho^{\kappa}_{M_{\bar{\theta}}}}\left[f^{*}(T_{\varphi_{1}}(x,a))\right],$

then we have
$f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})}{\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})}\right)\approx
T_{\varphi_{0}^{*}}(x,a,x^{\prime})$ and
$f^{\prime}\left(\frac{\rho^{\kappa}_{M_{\bar{\theta}}}(x,a)}{\rho^{\kappa}_{M^{*}}(x,a)}\right)\approx
T_{\varphi_{1}^{*}}(x,a)$. Given $\varphi_{0}^{*}$ and $\varphi_{1}^{*}$, let
$A(x,a,x^{\prime})=T_{\varphi_{0}^{*}}(x,a,x^{\prime})-T_{\varphi_{1}^{*}}(x,a)$,
then we can optimize $\theta$ via:

$\displaystyle\hat{\theta}^{*}=\arg\max_{\theta}\int_{\mathcal{X},\mathcal{A},\mathcal{X}}$
$\displaystyle\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x)\left(T_{\varphi_{0}^{*}}(x,a,x^{\prime})-T_{\varphi_{1}^{*}}(x,a)\right)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x+$
$\displaystyle\int_{\mathcal{X},\mathcal{A},\mathcal{X}}$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})\left(T_{\varphi_{1}^{*}}(x,a)-T_{\varphi_{0}^{*}}(x,a,x^{\prime})+H_{M^{*}}(x,a)\right)\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x$
$\displaystyle=\arg\max_{\theta}\int_{\mathcal{X},\mathcal{A},\mathcal{X}}$
$\displaystyle\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x)A(x,a,x^{\prime})\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x+$
$\displaystyle\int_{\mathcal{X},\mathcal{A},\mathcal{X}}$
$\displaystyle\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})(-A(x,a,x^{\prime})+H_{M^{*}}(x,a))\log
M_{\theta}(x^{\prime}|x,a)\mathrm{d}x^{\prime}\mathrm{d}a\mathrm{d}x.$

∎

##  Appendix B Discussion of the Theoretical Results

We summarize the limitations of current theoretical results and future work as
follows:

  1. 1.

As discussed in Remark A.2, the solution Eq.(11) relies on
$\rho^{\beta}_{M^{*}}(x,a)\in[0,c],\forall a\in\mathcal{A},\forall
x\in\mathcal{X}$. In some particular $M^{*}$, it is intractable to derive a
$\beta$ that can generate an occupancy specified by $g(x,a)/\alpha$. If more
knowledge of $M^{*}$ or $\beta^{*}$ is provided or some mild assumptions can
be made on the properties of $M^{*}$ or $\beta^{*}$, we may model $\rho$ in a
more sophisticated approach to alleviate the above problem.

  2. 2.

In the tractable solution derivation, we ignore the term
$\alpha_{0}(x,a)=\alpha\rho^{\kappa}_{M^{*}}(x,a)$ (See Eq.(19)). The benefit
is that $\rho^{\kappa}_{M^{*}}(x,a,x^{\prime})$ in the tractable solution can
be estimated through offline datasets directly. Although the results in our
experiments show that it does not produce significant negative effects in
these tasks, ignoring $\rho^{\kappa}_{M^{*}}(x,a)$ indeed incurs extra bias in
theory. In future work, techniques for estimating $\rho^{\kappa}_{M^{*}}(x,a)$
[49] can be incorporated to correct the bias.

##  Appendix C Societal Impact

This work studies a method towards counterfactual environment model learning.
Reconstructing an accurate environment of the real world will promote the wide
adoption of decision-making policy optimization methods in real life,
enhancing our daily experience. We are aware that decision-making policy in
some domains like recommendation systems that interact with customers may have
risks of causing price discrimination and misleading customers if
inappropriately used. A promising way to reduce the risk is to introduce
fairness into policy optimization and rules to constrain the actions. We are
involved in and advocating research in such directions. We believe that
business organizations would like to embrace fair systems that can ultimately
bring long-term financial benefits by providing better user experience.

##  Appendix D Details of GALILEO Implementation

###  D.1 Implementation

Based on the specific $f$-divergence, we can represent $T$ and $f^{*}(T)$ with
a discriminator $D_{\varphi}$. It can be verified that $f(u)=u\log
u-(u+1)\log(u+1)$, $T_{\varphi}(u)=\log D_{\varphi}(u)$, and
$f^{*}(T_{\varphi}(u))=\log(1-D_{\varphi}(u))$ proposed in [20] satisfies the
condition $f^{\prime}(x)\leq 0,\forall x\in\mathcal{X}$ (see Tab. 2). We
select the former in the implementation and convert the tractable solution to:

$\displaystyle\hat{\theta}^{*}=\arg\max_{\theta}~{}\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[A(x,a,x^{\prime})\log{M_{\theta}}(x^{\prime}|x,a)\right]+\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[(-A(x,a,x^{\prime})+H_{M^{*}}(x,a))\log{M_{\theta}}(x^{\prime}|x,a)\right]$
$\displaystyle s.t.\quad A(x,a,x^{\prime})=\log
D_{\varphi_{0}^{*}}(x,a,x^{\prime})-\log D_{\varphi_{1}^{*}}(x,a)$
$\displaystyle\quad\quad\quad\varphi_{0}^{*}=\arg\max_{\varphi_{0}}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[\log
D_{\varphi_{0}}(x,a,x^{\prime})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[\log(1-D_{\varphi_{0}}(x,a,x^{\prime}))\right]$
$\displaystyle\quad\quad\quad\varphi_{1}^{*}=\arg\max_{\varphi_{1}}~{}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[\log
D_{\varphi_{1}}(x,a)\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[\log(1-D_{\varphi_{1}}(x,a))\right],$

where $\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}[\cdot]$ is a
simplification of
$\mathbb{E}_{x,a,x^{\prime}\sim\rho^{\kappa}_{M_{\bar{\theta}}}}[\cdot]$.

The approximation of Eq.(20) holds only when $\frac{p(x)}{q(x)}$ is close to
$1$, which might not be satisfied. To handle the problem, we inject a standard
supervised learning loss

$\displaystyle\arg\max_{\theta}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[\log{M_{\theta}}(x^{\prime}|x,a)\right]$
(21)

to replace the second term of the above objective when the output probability
of $D$ is far away from $0.5$ ($f^{\prime}(1)=\log 0.5$).

In the offline model-learning setting, we only have a real-world dataset
$\mathcal{D}$ collected by the behavior policy $\mu$. We learn a policy
$\hat{\mu}\approx\mu$ via behavior cloning with $\mathcal{D}$ [39, 22] and let
$\hat{\mu}$ be the intermediary policy $\kappa$. We regard $\mathcal{D}$ as
the empirical data distribution of $\rho^{\kappa}_{M^{*}}$ and the
trajectories collected by $\hat{\mu}$ in the model $M_{\bar{\theta}}$ as the
empirical data distribution of $\rho^{\kappa}_{M_{\bar{\theta}}}$. But the
assumption $\forall x\in\mathcal{X},\forall a\in\mathcal{A},\mu(a|x)>0$ might
not be satisfied. In behavior cloning, we model $\hat{\mu}$ with a Gaussian
distribution and constrain the lower bound of the variance with a small value
$\epsilon_{\mu}>0$ to keep the assumption holding. Besides, we add small
Gaussian noises $\textbf{u}\sim\mathcal{N}(0,\epsilon_{D})$ to the inputs of
$D_{\varphi}$ to handle the mismatch between $\rho^{\mu}_{M^{*}}$ and
$\rho^{\hat{\mu}}_{M^{*}}$ due to $\epsilon_{\mu}$. In particular, for
$\varphi_{0}$ and $\varphi_{1}$ learning, we have:

$\displaystyle\varphi_{0}^{*}$
$\displaystyle=\arg\max_{\varphi_{0}}~{}\mathbb{E}_{\rho^{\kappa}_{M^{*}},\textbf{u}}\left[\log
D_{\varphi_{0}}(x+u_{x},a+u_{a},x^{\prime}+u_{x^{\prime}})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}},\textbf{u}}\left[\log(1-D_{\varphi_{0}}(x+u_{x},a+u_{a},x^{\prime}+u_{x^{\prime}}))\right]$
$\displaystyle\varphi_{1}^{*}$
$\displaystyle=\arg\max_{\varphi_{1}}~{}\mathbb{E}_{\rho^{\kappa}_{M^{*}},\textbf{u}}\left[\log
D_{\varphi_{1}}(x+u_{x},a+u_{a})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}},\textbf{u}}\left[\log(1-D_{\varphi_{1}}(x+u_{x},a+u_{a}))\right],$

where $\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}},\textbf{u}}[\cdot]$ is a
simplification of
$\mathbb{E}_{x,a,x^{\prime}\sim\rho^{\kappa}_{M_{\bar{\theta}}},\textbf{u}\sim\mathcal{N}(0,\epsilon_{D})}[\cdot]$
and $\textbf{u}=[u_{x},u_{a},u_{x^{\prime}}]$.

As stated in the main body, to consider the cumulative effects of prediction
errors induced by multi-step of transitions in environments, we overwrite the
function $A$ as
$A=Q^{\kappa}_{M_{\bar{\theta}}}-V^{\kappa}_{M_{\bar{\theta}}}$, where
$Q^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}\log
D_{\varphi_{0}^{*}}(x_{t},a_{t},x_{t+1})|(x_{t},a_{t},x_{t+1})=(x,a,x^{\prime}),\kappa,M_{\bar{\theta}}\right]$
and
$V^{\kappa}_{M_{\bar{\theta}}}(x,a)=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}\log
D_{\varphi_{1}^{*}}(x_{t},a_{t})|(x_{t},a_{t})=(x,a),\kappa,M_{\bar{\theta}}\right]$.

By adopting the above implementation techniques, we convert the objective into
the following formulation

$\displaystyle\hat{\theta}^{*}=\arg\max_{\theta}\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[A(x,a,x^{\prime})\log{M_{\theta}}(x^{\prime}|x,a)\right]+\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[(-A(x,a,x^{\prime})+H_{M^{*}}(x,a))\log{M_{\theta}}(x^{\prime}|x,a)\right]$
(22) $\displaystyle s.t.\quad
A(x,a,x^{\prime})=Q^{\kappa}_{M_{\theta}}(x,a,x^{\prime})-V^{\kappa}_{M_{\theta}}(x,a)$
(23) $\displaystyle\quad\quad\quad
Q^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}\log
D_{\varphi_{0}^{*}}(x_{t},a_{t},x_{t+1})|(x_{t},a_{t},x_{t+1})=(x,a,x^{\prime}),\kappa,M_{\bar{\theta}}\right]$
(24) $\displaystyle\quad\quad\quad
V^{\kappa}_{M_{\bar{\theta}}}(x,a)=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}\log
D_{\varphi_{1}^{*}}(x_{t},a_{t})|(x_{t},a_{t})=(x,a),\kappa,M_{\bar{\theta}}\right]$
(25)
$\displaystyle\quad\quad\quad\varphi_{0}^{*}=\arg\max_{\varphi_{0}}\mathbb{E}_{\rho^{\kappa}_{M^{*}},\mathbf{u}}\left[\log
D_{\varphi_{0}}(x+u_{x},a+u_{a},x^{\prime}+u_{x^{\prime}})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}},\mathbf{u}}\left[\log(1-D_{\varphi_{0}}(x+u_{x},a+u_{a},x^{\prime}+u_{x^{\prime}}))\right]$
(26)
$\displaystyle\quad\quad\quad\varphi_{1}^{*}=\arg\max_{\varphi_{1}}\mathbb{E}_{\rho^{\kappa}_{M^{*}},\mathbf{u}}\left[\log
D_{\varphi_{1}}(x+u_{x},a+u_{a})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}},\mathbf{u}}\left[\log(1-D_{\varphi_{1}}(x+u_{x},a+u_{a}))\right].$
(27)

In practice, GALILEO optimizes the first term of Eq.(22) with conservative
policy gradient algorithms (e.g., PPO [41] or TRPO [40]) to avoid unreliable
gradients for model improvements. Eq.(26) and Eq.(27) are optimized with
supervised learning. The second term of Eq.(22) is optimized with supervised
learning with a re-weighting term $-A+H_{M^{*}}$. Since $H_{M^{*}}$ is
unknown, we use $H_{M}$ to estimate it. When the mean output probability of a
batch of data is larger than $0.6$ or small than $0.4$, we replace the second
term of Eq.(22) with a standard supervised learning in Eq.(21). Besides,
unreliable gradients also exist in the process of optimizing the second term
of Eq.(22). In our implementation, we use the scale of policy gradients to
constrain the gradients of the second term of Eq.(22). In particular, we first
compute the $l_{2}$-norm of the gradient of the first term of Eq.(22) via
conservative policy gradient algorithms, named $||g_{\rm pg}||_{2}$. Then we
compute the $l_{2}$-norm of the gradient of the second term of Eq.(22), name
$||g_{\rm sl}||_{2}$. Finally, we rescale the gradients of the second term
$g_{\rm sl}$ by

$\displaystyle g_{\rm sl}\leftarrow g_{\rm sl}\frac{||g_{\rm
pg}||_{2}}{\max\\{||g_{\rm pg}||_{2},||g_{\rm sl}||_{2}\\}}.$

For each iteration, Eq.(22), Eq.(26), and Eq.(27) are trained with certain
steps (See Tab. 5) following the same framework as GAIL.

###  D.2 Connection with Previous Adversarial Algorithms

Standard GAN [19] can be regarded as a partial implementation including the
first term of Eq.(22) and Eq.(26) by degrading them into the one-step
scenario. In the context of GALILEO, the objective of GAN is

$\displaystyle\hat{\theta}^{*}=\arg\max_{\theta}~{}\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[A(x,a,x^{\prime})\log{M_{\theta}}(x^{\prime}|x,a)\right]$
$\displaystyle s.t.\quad A(x,a,x^{\prime})=\log
D_{\varphi^{*}}(x,a,x^{\prime})$
$\displaystyle\quad\quad\quad\varphi*=\arg\max_{\varphi}~{}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[\log
D_{\varphi}(x,a,x^{\prime})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[\log(1-D_{\varphi}(x,a,x^{\prime}))\right]$

In the one-step scenario,
$\rho^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})=\rho_{0}(x)\kappa(a|x)M_{\bar{\theta}}(x^{\prime}|a,x)$.
The term
$\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[A(x,a,x^{\prime})\log{M_{\theta}}(x^{\prime}|x,a)\right]$
can convert to $\mathbb{E}_{\rho^{\kappa}_{M_{\theta}}}\left[\log
D_{\varphi^{*}}(x,a,x^{\prime})\right]$ by replacing the gradient of
$M_{\bar{\theta}}(x^{\prime}|x,a)\nabla_{\theta}\log
M_{\theta}(x^{\prime}|x,a)$ with $\nabla_{\theta}M_{\theta}(x^{\prime}|x,a)$
[37]. Previous algorithms like GANITE [23] and SCIGAN [24] can be regarded as
variants of the above training framework.

The first term of Eq.(22) and Eq.(26) are similar to the objective of GAIL by
regarding $M_{\theta}$ as the “policy” to imitate and $\hat{\mu}$ as the
“environment” to collect data. In the context of GALILEO, the objective of
GAIL is:

$\displaystyle\hat{\theta}^{*}=\arg\max_{\theta}\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[A(x,a,x^{\prime})\log{M_{\theta}}(x^{\prime}|x,a)\right]$
$\displaystyle s.t.\quad
A(x,a,x^{\prime})=Q^{\kappa}_{M_{\theta}}(x,a,x^{\prime})-V^{\kappa}_{M_{\theta}}(x,a)$
$\displaystyle\quad\quad\quad
Q^{\kappa}_{M_{\bar{\theta}}}(x,a,x^{\prime})=\mathbb{E}\left[\sum^{\infty}_{t}\gamma^{t}\log
D_{\varphi^{*}}(x_{t},a_{t},x_{t+1})|(x_{t},a_{t},x_{t+1})=(x,a,x^{\prime}),\kappa,M_{\bar{\theta}}\right]$
$\displaystyle\quad\quad\quad
V^{\kappa}_{M_{\bar{\theta}}}(x,a)=\mathbb{E}_{M_{\bar{\theta}}(x,a)}\left[Q^{\kappa}(x,a,x^{\prime})\right]$
$\displaystyle\quad\quad\quad\varphi^{*}=\arg\max_{\varphi}\mathbb{E}_{\rho^{\kappa}_{M^{*}}}\left[\log
D_{\varphi}(x,a,x^{\prime})\right]+\mathbb{E}_{\rho^{\kappa}_{M_{\bar{\theta}}}}\left[\log(1-D_{\varphi}(x,a,x^{\prime}))\right].$

##  Appendix E Additional Related work

Besides adversarial algorithms for model learning, as we introduce in the main
body, there are several practices to reconstruct an environment in real-world
applications, especially in recommender systems. [50] construct a simulator to
imitate the feedback of customers via a world model. [51] use standard
supervised learning techniques with inverse propensity weighting (IPW) to
handle the selection bias problem of the logged dataset and construct a
debiased simulator; [52] use a real dataset to correct the representation and
reward function of a simulator to improve the fidelity.

##  Appendix F Experiment Details

###  F.1 Settings

####  F.1.1 General Negative Feedback Control (GNFC)

The design of GNFC is inspired by a classic type of scenario that behavior
policies $\mu$ have selection bias and easily lead to counterfactual risks:
For some internet platforms, we would like to allocate budgets to a set of
targets (e.g., customers or cities) to increase the engagement of the targets
in the platforms. Our task is to train a model to predict targets’ feedback on
engagement given targets’ features and allocated budgets.

In these tasks, for better benefits, the online working policy (i.e., the
behavior policy) will tend to cut down the budgets if targets have better
engagement, otherwise, the budgets might be increased. The risk of
counterfactual environment model learning in the task is that: the object with
better historical engagement will be sent to smaller budgets because of the
selection bias of the behavior policies, then the model might exploit this
correlation for learning and get a conclusion that: increasing budgets will
reduce the targets’ engagement, which violates the real causality. We
construct an environment and a behavior policy to mimic the above process. In
particular, the behavior policy $\mu_{GNFC}$ is

$\displaystyle\mu_{GNFC}(x)=\frac{(62.5-{\rm mean}(x))}{15}+\epsilon,$

where $\epsilon$ is a sample noise, which will be discussed later. The
environment includes two parts:

(1) dosage-response function $M_{1}(y|x,a)$:

$\displaystyle M_{1}(y|x,a)=\mathcal{N}({\rm mean}(x)+a,2)$

(2) mapping function $M_{2}(x^{\prime}|x,y)$:

$\displaystyle M_{2}(x^{\prime}|x,a,y)=y-{\rm mean}(x)+x$

The transition function $M^{*}$ is a composite of
$M^{*}(x^{\prime}|x,a)=M_{2}(x^{\prime}|x,a,M_{1}(y|x,a))$. The behavior
policies have selection bias: the actions taken are negatively correlated with
the states, as illustrated in Fig. 7(a) and Fig. 7(b). We control the
difficulty of distinguishing the correct causality of $x$, $a$, and $y$ by
designing different strategies of noise sampling on $\epsilon$. In principle,
with a larger number or more pronounced disturbances, there are more samples
violating the correlation between $x$ and $a$, then more samples can be used
to find the correct causality. Therefore, we can control the difficulty of
counterfactual environment model learning by controlling the strength of
disturbance. In particular, we sample $\epsilon$ from a uniform distribution
$U(-e,e)$ with probability $p$. That is, $\epsilon=0$ with probability $1-p$
and $\epsilon\sim U(-e,e)$ with probability $p$. Then with larger $p$, there
are more samples in the dataset violating the negative correlation (i.e.,
$\mu_{GNFC}$), and with larger $e$, the difference of the feedback will be
more obvious. By selecting different $e$ and $p$, we can construct different
tasks to verify the effectiveness and ability of the counterfactual
environment model learning algorithm.

(a) $y$ in collected trajectories.

(b) $a$ in collected trajectories.

Figure 7: Illustration of information about the collected dataset in GNFC.
Each color of the line denotes one of the collected trajectories. The X-axis
denotes the timestep of a trajectory.

####  F.1.2 The Cancer Genomic Atlas (TCGA)

The Cancer Genomic Atlas (TCGA) is a project that has profiled and analyzed
large numbers of human tumors to discover molecular aberrations at the DNA,
RNA, protein, and epigenetic levels. The resulting rich data provide a
significant opportunity to accelerate our understanding of the molecular basis
of cancer. We obtain features, $\mathbf{x}$, from the TCGA dataset and
consider three continuous treatments as done in SCIGAN [24]. Each treatment,
$a$, is associated with a set of parameters, $\mathbf{v}_{1}$,
$\mathbf{v}_{2}$, $\mathbf{v}_{3}$, that are sampled randomly by sampling a
vector from a standard normal distribution and scaling it with its norm. We
assign interventions by sampling a treatment, $a$, from a beta distribution,
$a\mid\mathbf{x}\sim\operatorname{Beta}\left(\alpha,\beta\right)$. $\alpha\geq
1$ controls the sampling bias and $\beta=\frac{\alpha-1}{a^{*}}+2-\alpha$,
where $a^{*}$ is the optimal treatment. This setting of $\beta$ ensures that
the mode of $\operatorname{Beta}\left(\alpha,\beta\right)$ is $a^{*}$.

The calculation of treatment response and optimal treatment are shown in Table
3.

Table 3: Treatment response used to generate semi-synthetic outcomes for
patient features $\mathbf{x}$. In the experiments, we set $C=10$.

Treatment | Treatment Response | Optimal treatment
---|---|---
1 | $f_{1}(\mathbf{x},a_{1})=C\left(\left(\mathbf{v}_{1}^{1}\right)^{T}\mathbf{x}+12\left(\mathbf{v}_{2}^{1}\right)^{T}\mathbf{x}a_{1}-12\left(\mathbf{v}_{3}^{1}\right)^{T}\mathbf{x}a_{1}^{2}\right)$ | $a_{1}^{*}=\frac{\left(\mathbf{v}_{2}^{1}\right)^{T}\mathbf{x}}{2\left(\mathbf{v}_{3}^{1}\right)^{T}\mathbf{x}}$
2 | $f_{2}(\mathbf{x},a_{2})=C\left(\left(\mathbf{v}_{1}^{2}\right)^{T}\mathbf{x}+\sin\left(\pi\left(\frac{\mathbf{v}_{2}^{2T}\mathbf{x}}{\mathbf{v}_{3}^{2T}\mathbf{x}}\right)a_{2}\right)\right)$ | $a_{2}^{*}=\frac{\left(\mathbf{v}_{3}^{2}\right)^{T}\mathbf{x}}{2\left(\mathbf{v}_{2}^{2}\right)^{T}\mathbf{x}}$
3 |  $f_{3}(\mathbf{x},a_{3})=C\left(\left(\mathbf{v}_{1}^{3}\right)^{T}\mathbf{x}+12a_{3}(a_{3}-b)^{2}\right.$, where $\left.b=0.75\frac{\left(\mathbf{v}_{2}^{3}\right)^{T}\mathbf{x}}{\left(\mathbf{v}_{3}^{3}\right)^{T}\mathbf{x}}\right)$ |  $\frac{3}{b}$ if $b\geq 0.75$, 1 if $b<0.75$

We conduct experiments on three different treatments separately and change the
value of bias $\alpha$ to assess the robustness of different methods to
treatment bias. When the bias of treatment is large, which means $\alpha$ is
large, the training set contains data with a strong bias on treatment so it
would be difficult for models to appropriately predict the treatment responses
out of the distribution of training data.

####  F.1.3 Budget Allocation task to the Time period (BAT)

In the platform, the order in different time periods will be sent to different
allowances to adjust imbalances of supply. For each day, the budget of
allowance is constant. The target of BAT is to find the best budget allocation
policy to improve the overall supply. To solve the problem, we first learn a
model to reconstruct the response of allowance for each delivery man and then
learn a budget allocation policy based on the learned model. In particular,
the input of the model includes the feature of weather, city information,
historical information about the delivery man and so on. The action is the
allowances sent to the delivery man. The model predicts the completed orders
of the time period for each delivery. Then we use a rule-based mapping
function to construct the next time-period states. The behavior policy is
fixed for the long term in this application. So we directly use the data
replay with a small scale of noise (See Tab. 5) to reconstruct the behavior
policy. Although we model the response for each delivery man, for fairness,
the budget allocation policy is just determining the allowance of each time
period and keeps the allowance to each delivery man the same. By adjusting the
budget allocation, we expect to improve the supply of the deliveries. In our
experiment, we evaluate the supply level by computing the averaged five-minute
order-taken rate, that is the percentage of orders picked up within five
minutes.

###  F.2 Baseline Algorithms

The algorithm we compared are: (1) Supervised Learning (SL): training a
environment model to minimize the expectation of prediction error, without
considering the counterfactual risks; (2) inverse propensity weighting (IPW)
[44]: a practical way to balance the selection bias by re-weighting. It can be
regarded as $\omega=\frac{1}{\hat{\mu}}$, where $\hat{\mu}$ is another model
learned to approximate the behavior policy; (3) SCIGAN: a recent proposed
adversarial algorithm for model learning for continuous-valued interventions
[24]. All of the baselines algorithms are implemented with the same capacity
of neural networks (See Tab. 5).

####  F.2.1 Supervised Learning (SL)

As a baseline, we train a multilayer perceptron model to directly predict the
response of different treatments, without considering the counterfactual
risks. We use mean square error to estimate the performance of our model so
that the loss function can be expressed as
$MSE=\frac{1}{n}\sum^{n}_{i=1}\left(y_{i}-\hat{y}_{i}\right)^{2}$, where $n$
is the number of samples, $y$ is the true value of response and $\hat{y}$ is
the predicted response. In practice, we train our SL models using Adam
optimizer and the initial learning rate $3e^{-4}$ on both datasets TCGA and
GNFC. The architecture of the neural networks is listed in Tab. 5.

####  F.2.2 Inverse Propensity Weighting (IPW)

Inverse propensity weighting [44] is an approach where the treatment outcome
model uses sample weights to balance the selection bias by re-weighting. The
weights are defined as the inverse propensity of actually getting the
treatment, which can be expressed as $\frac{1}{\hat{\mu}(a|x)}$, where $x$
stands for the feature vectors in a dataset, $a$ is the corresponding action
and $\hat{\mu}(a|x)$ indicates the action taken probability of $a$ given the
features $x$ within the dataset. $\hat{\mu}$ is learned with standard
supervised learning. The loss function can thus be expressed as
$\frac{1}{n}\sum^{n}_{i=1}\frac{1}{\mu(a_{i}|x_{i})}\left(y_{i}-\hat{y}_{i}\right)^{2}$.
The architecture of the neural networks is listed in Tab. 5.

####  F.2.3 SCIGAN

SCIGAN [24] is a model that uses generative adversarial networks to learn the
data distribution of the counterfactual outcomes and thus generate
individualized dose-response curves. SCIGAN does not place any restrictions on
the form of the treatment-does response functions and is capable of estimating
patient outcomes for multiple treatments, each with an associated dosage
parameter. SCIGAN first trains a generator to generate response curves for
each sample within the training dataset. The learned generator can then be
used to train an inference network using standard supervised methods. For fair
comparison, we increase the number of parameters for the the open-source
version of SCIGAN so that the SCIGAN model can have same order of magnitude of
network parameters as GALILEO. In addition, we also finetune the
hyperparameters (Tab. 4) of the enlarged SCIGAN to realize its full strength.
We set num_dosage_samples 9 and $\lambda=10$.

Table 4: Table of hyper-parameters for SCIGAN.  Parameter | Values
---|---
Number of dosage samples | 3, 5, 7, 9, 11
$\lambda$ | 0.1, 1, 10, 20

###  F.3 Hyper-parameters

We list the hyper-parameter of GALILEO in Tab. 5.

Table 5: Table of hyper-parameters for all of the tasks.

Parameter | GNFC | TAGC | MuJoCo | BAT
---|---|---|---|---
hidden layers of all neural networks | 4 | 4 | 5 | 5
hidden units of all neural networks | 256 | 256 | 512 | 512
collect samples for each time of model update | 5000 | 5000 | 96000 | 96000
batch size of discriminators | 5000 | 5000 | 80000 | 80000
$\epsilon_{D}$ | 0.005 | 0.01 | 0.05 | 0.05
$\epsilon_{\mu}$ | 0.005 | 0.01 | 0.05 | 0.05
times for discriminator update | 2 | 2 | 5 | 5
times for model update | 1 | 1 | 20 | 20
times for supervised learning update | 1 | 1 | 20 | 20
learning rate for supervised learning | 1e-5 | 1e-5 | 1e-5 | 1e-5
$\gamma$ | 0.99 | 0.0 | 0.99 | 0.99
clip_ratio | NAN | NAN | 0.1 | 0.1
max $D_{KL}$ | 0.001 | 0.001 | NAN | NAN
optimization algorithm (the first term of Eq.(22)) | TRPO [40] | TRPO [40] | PPO [41] | PPO [41]

###  F.4 Computation Resources

We use one Tesla V100 PCIe 32GB GPU and a 32-core Intel(R) Xeon(R) Gold 5118
CPU @ 2.30GHz to train all of our model.

##  Appendix G Additional Results

###  G.1 All of the Result Table

We give the result of CNFC in Tab. 8, TCGA in Tab. 7, BAT in Tab. 6 and MuJoCo
in Tab. 1.

###  G.2 Averaged Dosage-Response function

We give the averaged dosage responses for all of the tasks and the algorithms
in Fig. 11 to Fig. 18. We randomly select 20% of the states in the dataset and
equidistantly sample actions from the action space for each sampled state, and
plot the averaged the predicted feedbacks of each action. The real response is
slight different among different figure as the randomly-selected states for
testing is different.

###  G.3 Worst-case Prediction Error

In theory, GALILEO increases the generalization ability by focusing on the
worst-case samples’ training to achieve ACRM. To demonstrate the property, we
propose a new metric named Mean-Max Square Error (MMSE):
$\mathbb{E}\left[\max_{a\in\mathcal{A}}\left(M^{*}(x^{\prime}|x,a)-M(x^{\prime}|x,a)\right)^{2}\right]$
and give the results of MMSE for GNFC in Fig. 9 and for TCGA in Fig. 10.

###  G.4 Detailed results in the MuJoCo Tasks

We select 3 datasets from D4RL [45] to construct our model learning tasks. We
compare it with a typical transition model learning algorithm used in previous
offline model-based RL algorithm [6, 9], which is a variant of supervised
learning. We name the method OFF-SL. We train models in datasets
Walker2d-medium, Hopper-medium, and HalfCheetah-medium, which are collected by
a behavior policy with 1/3 performance to the expert policy, then we test them
in the corresponding expert dataset. We plot the learning curves of GALILEO
and OFF-SL in three MuJoCo tasks in Fig. 8.

(a) halfcheetah-medium (train) (b) hopper-medium (train) (c) walker2d-medium
(train) (d) halfcheetah-expert (test)

(e) hopper-expert (test) (f) walker2d-expert (test)

Figure 8: Illustration of learning curves of the MuJoCo Tasks. The X-axis
record the steps of the environment model update, and the Y-axis is the
corresponding prediction error. The figures with titles ending in “medium”
plot the prediction error in the training data while the titles ending in
“expert” plot the corresponding prediction error in the testing data.

In Fig. 8, we can see that both OFF-SL and GALILEO perform well in the
training datasets. OFF-SL can even reach a bit lower error in halfcheetah and
walker2d. However, when we verify the models through “expert” datasets, which
are collected by other policies, the performance of GALILEO is significantly
more stable and better than OFF-SL. As the training continues, OFF-SL even
gets worse and worse.

###  G.5 Detailed results in the BAT Task

Since there is no oracle environment model for querying, we have to describe
the results with other metrics.

(a) City-A (b) City-B (c) City-C (d) City-D

(e) City-E (f) City-F

Figure 9: Illustration of the dosage response curves in the 6 cities. Although
the ground-truth curves are unknown, it is expected to be monotonically
increasing.

First, we review whether the tendency of the dosage-response curve is
consistent. In this application, with a larger budget of allowance, the supply
will not be decreased. As can be seen in Fig. 9, the tendency of GALILEO’s
response is valid in 6 cities but almost all of the model of SL gives an
opposite direction to the response. If we learn a policy through the model of
SL, the optimal solution is canceling all of the allowances, which is
obviously incorrect in practice.

Finally, we search for the optimal policy via the cross-entropy method planner
[47] based on the learned model. We test the online supply improvement in $6$
cities. The algorithm compared is a human-expert policy, which is also the
behavior policy of the offline datasets. We conduct the online A/B tests for
each of the cities. For each test, we randomly split a city into two
partitions, one is for deploying the optimal policy learned from the GALILEO
model, and the other is as a control group, which keeps the human-expert
policy as before. Before the intervention, we collect $10$ days’ observation
data and compute the averaged five-minute order-taken rates as the baselines
of the treatment and control group, named $b^{t}$ and $b^{c}$ respectively.
Then we start intervention and observe the five-minute order-taken rate in the
following $14$ days for the two groups. The results of the treatment and
control groups are $y^{t}_{i}$ and $y^{c}_{i}$ respectively, where $i$ denotes
the $i$-th day of the deployment. The percentage points of the supply
improvement are computed via difference-in-difference (DID):

$\displaystyle\frac{\sum_{i}^{T}(y^{t}_{i}-b^{t})-(y^{c}_{i}-b^{c})}{T}\times
100,$

where $T$ is the total days of the intervention and $T=14$ in our experiments.

Table 6: Results on BAT. We use City-X to denote the experiments on different cities. “pp” is an abbreviation of percentage points on the supply improvement.  target | City-A | City-B | City-C
---|---|---|---
supply improvement | +1.63pp | +0.79pp | +0.27pp
target | City-D | City-E | City-F
supply improvement | +0.2pp | +0.14pp | +0.41pp

The results are summarized in Tab. 6. The online experiment is conducted in
$14$ days and the results show that the policy learned with GALILEO can make
better (the supply improvements are from 0.14 to 1.63 percentage points)
budget allocation than the behavior policies in all the testing cities. We
give the detailed results which record the supply difference between the
treatment group and the control group in Fig. 10.

(a) City-A (b) City-B (c) City-C (d) City-D

(e) City-E (f) City-F

Figure 10: Illustration of the daily responses in the A/B test in the 6
cities.

Table 7: $\sqrt{MISE}$ results on TCGA. We bold the lowest error for each task. $\pm$ is the standard deviation of three random seeds.  | t0_bias_2.0 | t0_bias_20.0 | t0_bias_50.0
---|---|---|---
GALILEO | 0.34 $\pm$ 0.05 | 0.67 $\pm$ 0.13 | 2.04 $\pm$ 0.12
SL | 0.38 $\pm$ 0.13 | 1.50 $\pm$ 0.31 | 3.06 $\pm$ 0.65
IPW | 6.57 $\pm$ 1.16 | 6.88 $\pm$ 0.30 | 5.84 $\pm$ 0.71
SCIGAN | 0.74 $\pm$ 0.05 | 2.74 $\pm$ 0.35 | 3.19 $\pm$ 0.09
| t1_bias_2.0 | t1_bias_6.0 | t1_bias_8.0
GALILEO | 0.43 $\pm$ 0.05 | 0.25 $\pm$ 0.02 | 0.21 $\pm$ 0.04
SL | 0.47 $\pm$ 0.05 | 1.33 $\pm$ 0.97 | 1.18 $\pm$ 0.73
IPW | 3.67 $\pm$ 2.37 | 0.54 $\pm$ 0.13 | 2.69 $\pm$ 1.17
SCIGAN | 0.45 $\pm$ 0.25 | 1.08 $\pm$ 1.04 | 1.01 $\pm$ 0.77
| t2_bias_2.0 | t2_bias_6.0 | t2_bias_8.0
GALILEO | 1.46 $\pm$ 0.09 | 0.85 $\pm$ 0.04 | 0.46 $\pm$ 0.01
SL | 0.81 $\pm$ 0.14 | 3.74 $\pm$ 2.04 | 3.59 $\pm$ 0.14
IPW | 2.94 $\pm$ 1.59 | 1.24 $\pm$ 0.01 | 0.99 $\pm$ 0.06
SCIGAN | 0.73 $\pm$ 0.15 | 1.20 $\pm$ 0.53 | 2.13 $\pm$ 1.75
Table 8: $\sqrt{MISE}$ results on GNFC. We bold the lowest error for each task. $\pm$ is the standard deviation of three random seeds.  | e1_p1 | e0.2_p1 | e0.05_p1
---|---|---|---
GALILEO | 5.09 $\pm$ 0.06 | 4.73 $\pm$ 0.13 | 5.08 $\pm$ 1.01
SL | 5.15 $\pm$ 0.23 | 4.73 $\pm$ 0.31 | 23.64 $\pm$ 4.86
IPW | 5.22 $\pm$ 0.09 | 5.50 $\pm$ 0.01 | 5.02 $\pm$ 0.07
SCIGAN | 7.05 $\pm$ 0.52 | 6.58 $\pm$ 0.58 | 18.55 $\pm$ 3.50
| e1_p0.2 | e0.2_p0.2 | e0.05_p0.2
GALILEO | 4.88 $\pm$ 0.10 | 6.11 $\pm$ 1.43 | 4.69 $\pm$ 0.02
SL | 5.21 $\pm$ 0.63 | 6.74 $\pm$ 0.15 | 33.52 $\pm$ 1.32
IPW | 5.27 $\pm$ 0.05 | 5.69 $\pm$ 0.00 | 20.23 $\pm$ 0.45
SCIGAN | 16.07 $\pm$ 0.27 | 12.07 $\pm$ 1.93 | 19.27 $\pm$ 10.72
| e1_p0.05 | e0.2_p0.05 | e0.05_p0.05
GALILEO | 4.78 $\pm$ 0.28 | 4.80 $\pm$ 0.08 | 4.88 $\pm$ 0.23
SL | 5.89 $\pm$ 0.88 | 14.25 $\pm$ 3.48 | 37.50 $\pm$ 2.29
IPW | 5.21 $\pm$ 0.01 | 5.52 $\pm$ 0.44 | 31.95 $\pm$ 0.05
SCIGAN | 11.50 $\pm$ 7.76 | 13.05 $\pm$ 4.19 | 25.74 $\pm$ 8.30
Table 9: $\sqrt{MMSE}$ results on GNFC. We bold the lowest error for each task. $\pm$ is the standard deviation of three random seeds.  | e1_p1 | e0.2_p1 | e0.05_p1
---|---|---|---
GALILEO | 3.95 $\pm$ 0.03 | 4.10 $\pm$ 0.01 | 4.03 $\pm$ 0.03
SL | 5.73 $\pm$ 0.33 | 5.80 $\pm$ 0.28 | 18.78 $\pm$ 3.13
IPW | 4.02 $\pm$ 0.05 | 4.15 $\pm$ 0.12 | 22.66 $\pm$ 0.33
SCIGAN | 8.84 $\pm$ 0.54 | 12.62 $\pm$ 2.17 | 24.21 $\pm$ 5.20
| e1_p0.2 | e0.2_p0.2 | e0.05_p0.2
GALILEO | 4.03 $\pm$ 0.10 | 4.31 $\pm$ 0.15 | 4.38 $\pm$ 0.18
SL | 5.87 $\pm$ 0.43 | 7.44 $\pm$ 1.13 | 29.13 $\pm$ 3.44
IPW | 4.12 $\pm$ 0.02 | 6.12 $\pm$ 0.48 | 30.96 $\pm$ 0.17
SCIGAN | 12.87 $\pm$ 3.02 | 14.59 $\pm$ 2.13 | 24.57 $\pm$ 3.00
| e1_p0.05 | e0.2_p0.05 | e0.05_p0.05
GALILEO | 4.79 $\pm$ 0.20 | 4.69 $\pm$ 0.20 | 5.02 $\pm$ 0.07
SL | 6.12 $\pm$ 0.43 | 14.88 $\pm$ 4.41 | 30.81 $\pm$ 1.69
IPW | 13.60 $\pm$ 7.83 | 26.27 $\pm$ 2.67 | 32.55 $\pm$ 0.12
SCIGAN | 9.19 $\pm$ 1.04 | 15.08 $\pm$ 1.26 | 17.52 $\pm$ 0.02
Table 10: $\sqrt{MMSE}$ results on TCGA. We bold the lowest error for each task. $\pm$ is the standard deviation of three random seeds.  | t0_bias_2.0 | t0_bias_20.0 | t0_bias_50.0
---|---|---|---
GALILEO | 1.56 $\pm$ 0.04 | 1.96 $\pm$ 0.53 | 3.16 $\pm$ 0.13
SL | 1.92 $\pm$ 0.67 | 2.31 $\pm$ 0.19 | 5.11 $\pm$ 0.66
IPW | 7.42 $\pm$ 0.46 | 5.36 $\pm$ 0.96 | 5.38 $\pm$ 1.24
SCIGAN | 2.11 $\pm$ 0.47 | 5.23 $\pm$ 0.27 | 5.59 $\pm$ 1.02
| t1_bias_2.0 | t1_bias_6.0 | t1_bias_8.0