# Grouped approximate control variate estimators

Alex A. Gorodetsky1∗ John D. Jakeman2 Michael S. Eldred2

( 1Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI.
48109
2Optimization and Uncertainty Quantification, Sandia National Laboratories,
Albuquerque, NM, 87123
∗Corresponding author<EMAIL_ADDRESS>
February 2024 )

###### Abstract

This paper analyzes the approximate control variate (ACV) approach to
multifidelity uncertainty quantification in the case where weighted estimators
are combined to form the components of the ACV. The weighted estimators enable
one to precisely group models that share input samples to achieve improved
variance reduction. We demonstrate that this viewpoint yields a generalized
linear estimator that can assign any weight to any sample. This generalization
shows that other linear estimators in the literature, particularly the
multilevel best linear unbiased estimator (ML-BLUE) of Schaden and Ullman in
2020, becomes a specific version of the ACV estimator of Gorodetsky, Geraci,
Jakeman, and Eldred, 2020. Moreover, this connection enables numerous
extensions and insights. For example, we empirically show that having non-
independent groups can yield better variance reduction compared to the
independent groups used by ML-BLUE. Furthermore, we show that such grouped
estimators can use arbitrary weighted estimators, not just the simple Monte
Carlo estimators used in ML-BLUE. Furthermore, the analysis enables the
derivation of ML-BLUE directly from a variance reduction perspective, rather
than a regression perspective.

##  1 Introduction

Multifidelity sampling approaches in uncertainty quantification (UQ) seek to
reduce the mean squared error in the estimation of the statistics of a high-
fidelity quantity of interest. These methods are essential for tractable
model-based analyses that require repeated evaluation of computationally
expensive simulations. Consequently, the last decade has seen tremendous
growth in multifidelity (MF) approaches for reducing error in UQ statistics.
For example, multilevel Monte Carlo (MLMC) [6], multi-index Monte Carlo [8]
(MIMC), multifidelity Monte Carlo (MFMC) [12], and multilevel-multifidelity
Monte Carlo [5] have all emerged as approaches to combine evaluations from
ensembles of model instances, where these instances are typically defined by
varying model forms and/or discretization levels in order to trade accuracy
versus cost.

Many of these existing MF sampling methods were unified by the approximate
control variate (ACV) framework [7]. This framework considers how linear
combinations of statistical estimators can be used to achieve improved
performance. These underlying estimators can be of varying types including
simple Monte Carlo, importance sampling [13], or surrogate based [18]. They
are also typically instantiated with specific sample allocation strategies
across models to ensure that certain estimators are grouped together with
respect to the input samples. Overall, this approach outlined a general
framework for optimal sample allocation with fused estimators, with
flexibility that can be exploited through additional configuration
optimizations such as directed acyclic graph and ensemble selection schemes
[1].

Recently, the Multilevel Best Linear Unbiased Estimator (ML-BLUE) approach
[15, 16] has emerged as a group-based sampling alternative to ACV. ML-BLUE
combines Monte Carlo (MC) estimators from different groupings of models based
on linear regression theory, focusing on the specific sampling context of a
set of independent MC estimators of groups of models. In this context, a
subset of ACV approaches (e.g., ACV-IS with independent samples) is an ML
BLUE. In a broader context (e.g., overlapping non-independent sample sets),
another subset of ACV (e.g. ACV-MF) are outside the class of BLUE. In this
paper, we explore the broader context of the ACV framework — especially the
consideration of sample-grouped estimators — and show that any ML-BLUE
estimator is in fact an approximate control variate, whereas not every ACV is
a ML-BLUE. While ML-BLUE is best for a specific type of model groupings, it
may not be best in the sense of providing the greatest variance reduction
amongst all possible groupings. In this paper, we show how ML-BLUEs can be
reformulated as ACVs and how this insight enables us to show that all
estimators that linearly weight samples from model ensembles are representable
as a linear ACV.

To be more specific, we consider the case of unbiased estimators so that all
of these approaches consider the following problem:

######  Problem 1 (Multifidelity Variance Reduction).

Let $Q$ denote a “high-fidelity” random variable, and let $\hat{Q}(c)$ denote
an unbiased estimate of some statistic of $Q$ achieved with cost $c$. Let $L$
“low-fidelity” random variables be denoted as $Q_{1},\ldots,Q_{L}$. Estimates
of statistics $\hat{Q}_{i}(c_{i})$ of each low fidelity model can also be
obtained at a cost $c_{i}.$ Multifidelity variance reduction strategies seek
to find a new estimator $\hat{Q}^{*}$ for the high-fidelity statistics by
fusing estimators of all models subject to a total cost.111One can also seek
the smallest variance for a fixed cost.

The solution to this problem typically has two components: first one must
determine an estimator ansatz, and, second, one must determine how the costs
are allocated (typically through numbers of samples). These tasks are
interrelated as the best sampling strategy may depend on the estimator ansatz,
and vice versa. Moreover, an “optimal” solution remains open. This paper makes
a significant stride in this direction by developing an estimator that can
represent any linear weighted combination of samples. The optimal choice of
groupings and sample allocations remains open. In particular, the novel
contribution of this paper is the development and analysis of the grouped ACV
(GACV) estimator (Definition 4), which provides a more general way to look at
the ACV where each of the low-fidelity model estimators estimators can be
formed through linear combinations of weighted estimators. Moreover, our
development of GACV estimators enabled the following additional contributions,
in which we:

  1. 1.

Demonstrate that any ML-BLUE defined in [15, 16, 3] and the Ensemble ACV [13]
can be interpreted as specific cases of the GACV (Theorem 3 and (25),
respectively);

  2. 2.

Generalize the idea of independent model groupings to non-independent
groupings and to general statistical estimators (beyond Monte Carlo estimators
of mean and variance) in Theorem 6.

  3. 3.

Show that the ML-BLUE approach can be derived through variance minimization
arguments (Corollary 8).

  4. 4.

Provide empirical validation that non-independent groupings can yield improved
variance reduction (Section 5.2).

Note that while analysis of ML-BLUEs as the solution of a variance reduction
problem has previously been performed in the case of mean and variance
estimation [3], our results are focused on the general case of arbitrary
estimators.

The remainder of this paper is structured as follows. In Section 2, we
describe the connection between ACVs and ML-BLUEs that serves to inspire and
motivate the general approach described in Section 3. The generalized
estimator is then analyzed and its optimal selection of weights is defined in
Section 4. Finally, Section 5 provides a setting for comparisons between the
GACV and the ML-BLUE, and it is shown that the additional generality of the
GACV can produce estimators with smaller variance than ML-BLUE estimators in
certain cases.

###  1.1 Notation

Sets are denoted by caligraphic upper case, such as $\mathcal{W}.$ Exception
is made for ordered sets of real numbers corresponding to vectors and
matrices. Vectors are denoted by bold lowercase, e.g., $\bm{w}$, and matrices
are denoted by bold uppercase, e.g., $\mathbf{W}$. We make specific exception
for zero vectors and matrices where $\mathbf{0}_{L}$ and $\mathbf{0}_{LM}$
denote a vector of size $L$ and a matrix of size $L\times M$, respectively. A
finite set of integers is denoted by $[L]=\\{0,1,\ldots,L\\}.$ The power set
$2^{\mathcal{S}}$ is the set of all possible subsets of $\mathcal{S}$.

The elements of a countable set can be mapped to the integers via an index to
provide an ordering. For example,

$\mathcal{W}=\\{\mathbf{w}^{1},\ldots,\mathbf{w}^{k}\\},\qquad
k\in\mathbb{Z}^{+}$ (1)

denotes a set of $k$ vectors. . While superscripts are used to define an
ordering in a general set, subscripts are used to denote elements of a vector
or matrix — as in $\mathbf{w}_{i}$ or $\mathbf{W}_{ij}$. In contrast to
elements of a vector, the elements of a set may be different objects. For
example, the set in (1) differs from a matrix because each of the vectors may
be of different sizes.

We will make special use of multi-indices to denote subsets of $[L].$ Consider
a subset of $[L]$ with $n$ elements denoted by $\mathcal{S}\in
2^{[L]}\,\backslash\,\emptyset,$ with $|\mathcal{S}|=n.$ A lexigraphical
ordering of this set can be associated with a multi-index vector
$\bm{\lambda}\in[L]^{n}$ with elements $\bm{\lambda}_{i}\in\mathcal{S}$ for
$i=1,\ldots,n,$ and $\bm{\lambda}_{i}<\bm{\lambda}_{j}$ if $i\leq j$. Given
this ordering, an inverse mapping can be established
$\bm{\lambda}^{-1}:\mathcal{S}\to[n]$ from an element to its location in the
multi-index. Finally, associated with this subset is a restriction matrix
$\mathbf{R}$ that determines $\bm{\lambda}$ from the vector $[L]$. That is
$\mathbf{R}\in\mathbb{R}^{n\times(L+1)}$ with elements

$\mathbf{R}_{ij}=\left\\{\begin{array}[]{cc}1&\text{ if }\bm{\lambda}_{i}=j\\\
0&\text{ otherwise}\end{array}\right.,$ (2)

so that $\bm{\lambda}=\mathbf{R}[L].$ Lastly, $\mathbf{1}(\bm{w}\geq 0)$
denotes a masking function which outputs a vector with entries that assume the
value of one when $\bm{w}_{i}\geq 0$ and zero otherwise.

###  1.2 Linear Approximate Control Variates

In this section we describe the linear approximate control variate. We begin
with a description of the classical control variate from which it is derived.
Control variates combine an ensemble of estimators to obtain a new estimator
with reduced variance compared to some baseline. Let the estimators in this
ensemble be $\hat{\mathcal{Q}}=\\{\hat{Q}$,
$\hat{Q}_{1},\ldots,\hat{Q}_{L}\\},$ and let $\hat{Q}$ denote the baseline
estimator whose variance we seek to reduce.

The linear control variate creates a new estimator through the linear
combination

$\hat{Q}^{CV}(\alpha,\alpha_{1},\alpha_{2},\ldots,\alpha_{L},\beta_{1},\ldots,\beta_{L})=\alpha\hat{Q}+\sum_{\ell=1}^{L}\alpha_{\ell}\hat{Q}_{\ell}+\sum_{\ell=1}^{L}\beta_{\ell}\mu_{\ell}$
(3)

where $\mu_{\ell}$, for $\ell=1,\ldots,L,$ are presumed to be the known
expectations of $Q_{1},\ldots,Q_{L}$; and
$\alpha,\alpha_{1},\ldots,\alpha_{L},\beta_{1},\ldots,\beta_{L}$ are real-
valued tunable coefficients.

The control variate estimator $\hat{Q}^{CV}$ can be made unbiased (i.e.
$\mathbb{E}\left[\hat{Q}^{CV}\right]=\mathbb{E}\left[\hat{Q}\right]$) through
two constraints on the coefficients. The first is that the baseline estimator
$\hat{Q}$ is unbiased so that we can set $\alpha=1$. The second is that the
$L$ low-fidelity estimators $\hat{Q}_{i}$ are unbiased so that setting
$\beta_{\ell}=-\alpha_{\ell}$ cancels out the sum in expectation. Under these
constraints the standard CV representation becomes

$\hat{Q}^{CV}(\alpha_{1},\alpha_{2},\ldots,\alpha_{L})=\hat{Q}+\sum_{\ell=1}^{L}\alpha_{\ell}\left(\hat{Q}_{\ell}-\mu_{\ell}\right)$
(4)

The remaining free parameters $\alpha_{1},\ldots,\alpha_{L}$ can be selected
according to the algorithmic strategy, and typically they are chosen to
minimize the overall estimator variance. We refer the reader to [11, 10, 14]
for additional details.

The linear control variate requires knowledge of the means
$\mu_{1},\ldots,\mu_{L}$. The linear approximate control variate was
introduced to circumvent this issue [7]. The approach assumed a simple
heuristic: replace $\mu_{\ell}$, termed the expected control variate mean
(ECVM), by another estimator $\hat{Q}_{\ell}^{\mu}$. To distinguish it from
the first term, we refer to the first term as the correlated mean estimator
(CME) and denote it by $\hat{Q}_{\ell}^{e}$ so that the approximate control
variate is written as [7, Eq 7]

$\displaystyle\hat{Q}^{ACV}(\alpha_{1},\alpha_{2},\ldots,\alpha_{L})$
$\displaystyle=\hat{Q}+\sum_{\ell=1}^{L}\alpha_{\ell}\left(\hat{Q}^{e}_{\ell}-\hat{Q}^{\mu}_{\ell}\right)=\hat{Q}+\sum_{\ell=1}^{L}\alpha_{\ell}\Delta_{i},$
(5)

where $\Delta_{\ell}=\hat{Q}^{e}_{\ell}-\hat{Q}^{\mu}_{\ell}.$ For
$\hat{Q}^{ACV}$ to remain unbiased, it is sufficient to require that either
both $\hat{Q}_{\ell}^{e}$ and $\hat{Q}_{\ell}^{\mu}$ are unbiased or have the
same bias.

Finally, it will often be convenient to treat the baseline estimator in the
same framework as the $L$ other alternatives by rewriting (5) as

$\displaystyle\hat{Q}^{ACV}(\alpha_{1},\alpha_{2},\ldots,\alpha_{L})$
$\displaystyle=\sum_{\ell=0}^{L}\alpha_{\ell}\left(\hat{Q}^{e}_{\ell}-\hat{Q}^{\mu}_{\ell}\right)=\sum_{\ell=0}^{L}\alpha_{\ell}\Delta_{\ell},$
(6)

where $\alpha_{0}=1$, $Q_{0}=Q$, $\hat{Q}^{e}_{0}=\hat{Q}=\Delta_{0}$. Note
that in this case we choose $\hat{Q}^{\mu}_{0}=0$ for notational convenience,
since there is only a single estimator for the high-fidelity model. We also
note that we have not specified any particular form on
$\hat{Q}_{\ell}^{e,\mu}$ – much of the general theory is agnostic to specific
choices, as long as they are unbiased. Again the optimal weights can be chosen
to minimize the variance of this estimator, and the reader is referred to [7].
Furthermore, we note that the underlying estimators can also be chosen by the
user, and need not be restricted to Monte Carlo sampling.

###  1.3 Multilevel BLUE estimators

The ML-BLUE estimator is formed by a linear combination of correlated Monte
Carlo estimators. It is defined via $K$ model groups, where each group
consists of a set of Monte Carlo estimators for $L+1$ models and different
groups consist of estimators formed by independent samples. Specifically the
$k$-th group is a couple ($\mathcal{S}^{k},\mathcal{Z}^{k})$, where the models
in the group are a non-empty element of the powerset $\mathcal{S}^{k}\in
2^{[L]}\,\backslash\,\emptyset$ of size $|\mathcal{S}^{k}|=n_{k}.$ Associated
with this group is the multi-index vector $\bm{\lambda}^{k}\in[L]$ and
restriction matrix $\mathbf{R}^{k}$ as discussed in Section 1.1. The samples
used for the Monte Carlo estimators in this group are
$\mathcal{Z}^{k}=\\{z_{1}^{k},\ldots,z_{m_{k}}^{k}\\}.$ Within each group, the
ML-BLUE estimator constructs the Monte Carlo estimators
$\hat{Q}_{\ell}^{k}=\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Q_{\ell}(z_{i}^{k})$ that
are correlated due to shared samples within the group. Note that
$\mathcal{Z}^{k}$ and $\mathcal{Z}^{j}$ are independent sets so that different
groups are not correlated. Finally, the ML-BLUE estimator weights and combines
these estimators. The weights for each group are given by
$\mathcal{B}=\\{\bm{\beta}^{1},\ldots,\bm{\beta}^{K}\\}$ with
$\bm{\beta}^{k}\in\mathbb{R}^{n_{k}}$, and then the ML-BLUE estimator is given
as [16, Eq 2.6]

$\hat{Q}^{ML-
BLUE}(\mathcal{B})~{}=~{}\sum_{k=1}^{K}\sum_{\ell\in\mathcal{S}^{k}}\bm{\beta}_{\bm{\lambda}^{k^{-1}}(\ell)}^{k}\hat{Q}^{k}_{\ell}~{}=~{}\sum_{k=1}^{K}\sum_{\ell\in\mathcal{S}^{k}}\bm{\beta}_{\bm{\lambda}^{k^{-1}}(\ell)}^{k}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Q_{\ell}(z^{k}_{i}).$
(7)

Note that this estimator is only “best” under an optimal choice of weights,
and the form above with optimal weights corresponds to the solution of a
corresponding least squares regression problem [15, Eq. 2.8]. However, the
theory that follows is valid any choice of weights, including the optimal
weights chosen by the least squares problem cited above. Thus, we are slightly
generalizing the definition of an ML-BLUE to any weighted set of independent
estimators — it should be understood that many of the desirable properties of
the ML-BLUE are in fact due to choosing optimal weights.

It will later be useful to extend each of the coefficient vectors
$\bm{\beta}^{k}$ into a vector of size $L+1$ with zeros filling the elements
that correspond to models not included in group $S^{k}$. These vectors,
$\tilde{\bm{\beta}}^{k}\in\mathbb{R}^{L+1}$ are obtained through
multiplication by the transpose of the corresponding restriction matrices

$\tilde{\bm{\beta}}^{k}=\mathbf{R}^{k^{T}}\bm{\beta}^{k}.$ (8)

Using these zero-filled vectors, the ML-BLUE is equivalently

$\hat{Q}^{ML-
BLUE}(\mathcal{B})~{}=~{}\sum_{k=1}^{K}\sum_{\ell\in\mathcal{S}^{k}}^{L}\tilde{\bm{\beta}}_{\ell}^{k}\hat{Q}^{k}_{\ell}~{}=~{}\sum_{k=1}^{K}\sum_{\ell=0}^{L}\tilde{\bm{\beta}}_{\ell}^{k}\hat{Q}^{k}_{\ell},$
(9)

with the convention that $\hat{Q}^{k}_{\ell}=0$ if $\ell\notin\mathcal{S}^{k}$
in the second equality.

An appropriate choice of coefficients $\bm{\beta}_{\ell}^{k}$ leads to an
estimator that both is unbiased and has reduced variance. We will return to
the optimal set of coefficients in Section 4. For now, we will demonstrate
that, for any choice of coefficients that yields an unbiased estimator, the
MLBLUE is also an ACV with a particular choice of ensemble estimators.

We begin with a requirement needed for unbiasedness.

######  Proposition 1.

Let $\beta_{\ell}=\sum_{k=1}^{K}\bm{\tilde{\beta}}_{\ell}^{k}$ for
$\ell=0,\ldots,L$ be the sum of the weights of each model $\ell$ across the
groups. An unbiased estimator of the form (7) requires

$\left(1-\beta_{0}\right)\mathbb{E}\left[Q_{0}\right]=\sum_{\ell=1}^{L}\mathbb{E}\left[Q_{\ell}\right]\beta_{\ell}$
(10)

###### Proof.

To be unbiased, we require $\mathbb{E}\left[\hat{Q}^{ML-
BLUE}(\beta)\right]=\mathbb{E}[Q_{0}]=\mathbb{E}[Q].$ From definition (9), we
have

$\mathbb{E}\left[\hat{Q}^{ML-
BLUE}(\bm{\beta})\right]=\sum_{k=1}^{K}\sum_{\ell=0}^{L}\tilde{\bm{\beta}}_{\ell}^{k}\mathbb{E}\left[Q_{\ell}\right].$
(11)

We now separate the highest fidelity model from the lower fidelity models and
re-arrange the sum to obtain

$\displaystyle\mathbb{E}\left[\hat{Q}^{MBLUE}(\beta)\right]$
$\displaystyle=\mathbb{E}\left[Q_{0}\right]\sum_{k=1}^{K}\bm{\beta}_{0}^{k}+\sum_{\ell=1}^{L}\mathbb{E}\left[Q_{\ell}\right]\sum_{k=1}^{K}\bm{\beta}_{\ell}^{k}=\mathbb{E}\left[Q_{0}\right]\beta_{0}+\sum_{\ell=1}^{L}\mathbb{E}\left[Q_{\ell}\right]\beta_{\ell}.$
(12)

Thus to achieve an unbiased estimate we require

$\mathbb{E}\left[Q_{0}\right]=\mathbb{E}\left[Q_{0}\right]\beta_{0}+\sum_{\ell=1}^{L}\mathbb{E}\left[Q_{\ell}\right]\beta_{\ell}.$
(13)

Equality necessitates that the coefficients satisfy (10), completing the
proof. ∎

These conditions on the coefficients needed to satisfy the unbiasedness
requirement are clearly dependent on the actual expected values of the
quantities of interest. This makes the model-specific solution uncomputable in
practice because these quantities are not available. Thus we must ask whether
there exists a set of coefficients for which this condition is satisfied
regardless of the choice of models. The linear relationship between
coefficients and the expectations makes the following condition obvious.

######  Proposition 2 (ML-BLUE Unbiasedness).

The ML-BLUE is unbiased if $\beta_{0}=1$ and $\beta_{\ell}=0$ for
$\ell=1,\ldots,L.$

Moreover, all reported optimal coefficients satisfy these conditions, and are
guaranteed to satisfy it due to due to Gauss-Markov-Aitken theorem arguments
[15, 16].

##  2 ML-BLUE estimators are approximate control variates

When originally proposed, ML-BLUE estimators were not recognized as ACV
estimators. However, any ML-BLUE estimator is in fact an approximate control
variate. To prove this, we transform the expression of an ML-BLUE estimator
into the corresponding estimators $\hat{Q}_{\ell}^{e}$, $\hat{Q}_{\ell}^{\mu}$
and weights $\alpha_{\ell}$ that comprise an ACV estimator taking the form of
(5).

######  Theorem 3 (ML-BLUE is an ACV).

For any set of coefficients $\mathcal{B}$, groupings
$\\{\mathcal{S}^{1},\ldots,\mathcal{S}^{k}\\}$ and samples
$\\{\mathcal{Z}^{1},\ldots,\mathcal{Z}^{K}\\}$ for which the ML-BLUE estimator
(Eqs. (7) or (9)) is unbiased, there corresponds an equivalent ACV estimator
(5) defined by weights

$\alpha_{\ell}=\sum_{k=1}^{K}\mathbf{1}(\tilde{\bm{\beta}}_{\ell}^{k}\geq
0)\tilde{\bm{\beta}}_{\ell}^{k}=-\sum_{k=1}^{K}\mathbf{1}(\tilde{\bm{\beta}}_{\ell}^{k}<0)\tilde{\bm{\beta}}_{\ell}^{k},$
(14)

for $\ell=1,\ldots,L,$ and a set of weighted estimators

$\displaystyle\hat{Q}=\sum_{k=1}^{K}\tilde{\bm{\beta}}_{0}^{k}\hat{Q}_{0}^{k}\qquad\hat{Q}_{\ell}^{e}=\sum_{k=1}^{K}\omega_{\ell}^{k,e}\hat{Q}_{\ell}^{k}\qquad\hat{Q}_{\ell}^{\mu}=\sum_{k=1}^{K}\omega_{\ell}^{k,\mu}\hat{Q}_{\ell}^{k},$
(15)

where $\hat{Q}_{\ell}^{k}$ is a Monte Carlo estimator using the samples
$\mathcal{Z}^{k}$. The weights are given by
$\omega_{\ell}^{k,e}=\mathbf{1}(\tilde{\bm{\beta}}_{\ell}^{k}\geq
0)\tilde{\bm{\beta}}_{\ell}^{k}/\alpha_{\ell}$ and
$\omega_{\ell}^{k,\mu}=-\mathbf{1}(\tilde{\bm{\beta}}_{\ell}^{k}<0)\tilde{\bm{\beta}}_{\ell}^{k}/\alpha_{\ell},$
for $\ell=1,\ldots,L$.

###### Proof.

The proof follows from a straightforward grouping of estimators into those
with positive and negative weights. We first rewrite Equation (9) by reversing
the order of the sums, pulling out the high-fidelity model, and separating the
positive and negative coefficients

$\displaystyle\hat{Q}^{ML-BLUE}(\mathcal{B})$
$\displaystyle=\sum_{k=1}^{K}\tilde{\bm{\beta}}_{0}^{k}\hat{Q}_{0}+\sum_{\ell=1}^{L}\left(\sum_{k=1}^{K}\tilde{\bm{\beta}}_{\ell}^{k}\hat{Q}^{k}_{\ell}\right)$
(16)
$\displaystyle=\sum_{k=1}^{K}\tilde{\bm{\beta}}_{0}^{k}\hat{Q}_{0}+\sum_{\ell=1}^{L}\left(\sum_{k=1}^{K}\mathbf{1}(\tilde{\bm{\beta}}_{\ell}^{k}\geq
0)\tilde{\bm{\beta}}_{\ell}^{k}\hat{Q}^{k}_{\ell}+\mathbf{1}(\tilde{\bm{\beta}}_{\ell}^{k}<0)\tilde{\bm{\beta}}_{\ell}^{k}\hat{Q}^{k}_{\ell}\right)$
(17)
$\displaystyle=\sum_{k=1}^{K}\tilde{\bm{\beta}}_{0}^{k}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Q_{0}(z^{k}_{i})+\sum_{\ell=1}^{L}\alpha_{\ell}\left(\sum_{k=1}^{K}\omega_{\ell}^{k,e}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Q_{\ell}(z^{k}_{i})-\sum_{k=1}^{K}\omega_{\ell}^{k,\mu}\frac{1}{m_{k}}\sum_{i=1}^{m_{k}}Q_{\ell}(z^{k}_{i})\right)$
(18)

where the final equality used (14) and the weight definitions. Substituting
the definitions (15) we obtain

$\displaystyle\hat{Q}^{ML-
BLUE}(\mathcal{B})=\hat{Q}+\sum_{\ell=1}^{L}\alpha_{\ell}\left(\hat{Q}_{\ell}^{e}-\hat{Q}_{\ell}^{\mu}\right).$
(19)

The right hand side is equivalent to an ACV estimator, if each of the
components can be shown to be unbiased. Note that the weights clearly satisfy
the following properties $\sum_{k=1}^{K}\omega_{\ell}^{k,1}=1$,
$\omega_{\ell}^{k,0}>0$, $\sum_{k=1}^{K}\omega_{\ell}^{k,2}=1$ and
$\omega_{\ell}^{k,2}>0$. Consequently,

$\mathbb{E}\left[\hat{Q}\right]=\sum_{k=1}^{K}\tilde{\bm{\beta}}_{0}^{k}\mathbb{E}\left[Q_{0}\right]=\mathbb{E}\left[Q_{0}\right]\qquad\mathbb{E}\left[\hat{Q}_{\ell}^{e}\right]=\sum_{k=1}^{K}\omega_{\ell}^{k,e}\mathbb{E}\left[Q_{\ell}\right]=\mathbb{E}\left[Q_{\ell}\right]\qquad\mathbb{E}\left[\hat{Q}_{\ell}^{\mu}\right]=\sum_{k=1}^{K}\omega_{\ell}^{k,\mu}\mathbb{E}\left[Q_{\ell}\right]=\mathbb{E}\left[Q_{\ell}\right].$

The last equality of the first expression for $\mathbb{E}[\hat{Q}]$ arises due
to Proposition 2. The other expressions follow from the summation properties
of the weights. ∎

##  3 Generalized estimators

The connection between ML-BLUEs and ACVs suggests two generalizations: (1) we
can use general estimators within a grouping strategy, and (2) group
estimations may benefit from relaxing the independence requirement and reusing
samples across correlated input sets. Note that this last benefit has already
been observed within the class of ACV estimators, e.g., the ACV-MF estimator
[7, Def. 3] can outperform the ACV-IS estimator [15, Sec 5.4] when exchanging
sample independence for sample reuse222Note that the ACV-IS estimator
described in [15, Sec. 5.4] is slightly different from ACV-IS in [7, Def. 2]
because the ACV paper reuses a small portion of samples for the ECVM. The
authors believe that the version in [15] is more in the spirit of the
“independent” aspect of the ACV-IS, and we follow their work here. However,
empirical results in both papers suggest similar behavior in both cases..

Using the idea of groupings and the connection we established between ML-BLUE
and ACV, we now provide several novel estimators based on what we call the
grouped ACV ansatz

$\displaystyle\hat{Q}^{GACV}$
$\displaystyle=\sum_{k=1}^{K}\omega_{0}^{k,e}\hat{Q}_{0}^{k}+\sum_{\ell=1}^{L}\alpha_{\ell}\left(\sum_{k=1}^{K}\omega_{\ell}^{k,e}\hat{Q}_{\ell}^{k,e}-\sum_{k=1}^{K}\omega_{\ell}^{k,\mu}\hat{Q}_{\ell}^{k,\mu}\right).$
(20)

Like ML-BLUE, there are $K$ groups, but now there is no requirement for each
set of group samples to be independent. Nor is there any requirement for the
individual estimators to arise from Monte Carlo. Moreover, both the weights on
the $e$ and $\mu$ terms may be nonzero for the same $\ell$ and $k$ if
$\hat{Q}_{\ell}^{k,e}$ and $\hat{Q}_{\ell}^{k,\mu}$ do not use identical
samples.

To make (20) easier to work with, first note that (20) is an unbiased
estimator if $\hat{Q}_{0}^{k}$ and $\hat{Q}_{\ell}^{k,(e,\mu)}$ are unbiased
(or the latter have the same bias) for $\ell=1,\ldots,L$ while
$\sum_{k=1}^{K}\omega_{\ell}^{k,e}=\sum_{k=1}^{K}\omega_{\ell}^{k,\mu}$ and
$\sum_{k=1}^{K}\omega_{0}^{k,\mu}=1.$ Second, note that the control variate
weight $\alpha_{\ell}$ can be embedded into the coefficients $\omega$, with no
change in constraints, so that we can write

$\displaystyle\hat{Q}^{GACV}$
$\displaystyle=\sum_{k=1}^{K}\omega_{0}^{k,e}\hat{Q}_{0}^{k}+\sum_{\ell=1}^{L}\left(\sum_{k=1}^{K}\omega_{\ell}^{k,e}\hat{Q}_{\ell}^{k,e}-\sum_{k=1}^{K}\omega_{\ell}^{k,\mu}\hat{Q}_{\ell}^{k,\mu}\right),$
(21)

This expression can be generalized further by noting that we implicitly
assumed $\hat{Q}_{\ell}^{k,e}\neq\hat{Q}_{\ell}^{k,\mu}$, since if they were
equal, these two terms would be combined. Consequently, the second term can be
treated as another set of $K$ groups resulting in the following expression
that contains a total of $2K$ groups

$\displaystyle\hat{Q}^{GACV}~{}$
$\displaystyle=~{}\sum_{k=1}^{K}\omega_{0}^{k,e}\hat{Q}_{0}^{k}+\sum_{k=1}^{K}\sum_{\ell=1}^{L}\omega_{\ell}^{k,e}\hat{Q}_{\ell}^{k,e}+\sum_{k=K+1}^{2K}\sum_{\ell=1}^{L}\omega_{\ell}^{k,\mu}\hat{Q}_{\ell}^{k,\mu}~{}=~{}\sum_{k=1}^{K}\omega_{0}^{k,e}\hat{Q}_{0}^{k}+\sum_{k=1}^{2K}\sum_{\ell=1}^{L}\tilde{\bm{\beta}}_{\ell}^{k}\hat{Q}_{\ell}^{k},$
(22)

where the last equality combines the two double summations of the previous
equality by setting the weights
$\tilde{\bm{\beta}}^{k}_{\ell}=\omega_{\ell}^{k,e}$ if $k\leq K$ and
$\omega_{\ell}^{k,\mu}$ otherwise, and by defining the estimator
$\hat{Q}_{\ell}^{k}$ as $\hat{Q}_{\ell}^{k,e}$ if $k\leq K$ and
$\hat{Q}_{\ell}^{k,\mu}$ otherwise. The following definition extends this idea
developed for two separate groupings to a more general strategy that can be
used with $K$ groupings.

######  Definition 4 (Generalized Linear Grouped ACV estimator).

The generalized grouped ACV estimator is defined by $K$ model groupings
consisting of subsets of estimators of $L+1$ models. The $k$-th group is
defined by a linear combination of estimators for models defined in
$\mathcal{S}^{k}\in 2^{[L]}\,\backslash\,\emptyset$ of size
$|\mathcal{S}^{k}|=n_{k}$ according to

$\displaystyle\hat{Q}^{GACV}(\mathcal{B})~{}=~{}\sum_{k=1}^{K}\sum_{\ell\in\mathcal{S}^{k}}\bm{\beta}_{\bm{\lambda}^{k^{-1}}(\ell)}^{k}\hat{Q}_{\ell}^{k}~{}=~{}\sum_{k=1}^{K}\sum_{\ell=0}^{L}\tilde{\bm{\beta}}_{\ell}^{k}\hat{Q}_{\ell}^{k},$
(23)

where $\mathcal{B}=\\{\bm{\beta}^{1},\ldots,\bm{\beta}^{K}\\}$ with
$\bm{\beta}^{k}\in\mathbb{R}^{n_{k}}$. Furthermore, $\bm{\lambda}^{k}$ and
$\mathbf{R}^{k}$ are the multi-index and restriction operators associated with
each $\mathcal{S}^{k}$, as discussed in Section 1.1. Finally,
$\tilde{\bm{\beta}}^{k}$ is defined as in Equation (8).

Note this estimator appears identical to (9), with this exception that
$\hat{Q}_{\ell}^{k}$ need not be a simple Monte Carlo estimator. For example,
it could be an importance sampling estimator or a more complex estimator
derived as some function obtained by first constructing a surrogate model. A
second, and more critical, difference is that the estimators across groups
need not be independent. Indeed each group can even use an identical set of
models, but with different estimators. Nevertheless, the linear structure
yields the same requirements on the weights as the ML-BLUE estimator to ensure
unbiasedness. Additionally, the estimator can be written as an ACV with
essentially the same argument as in Section 2. Note that each estimator
$\hat{Q}_{\ell}^{k}$ can itself be potentially optimized; for example, one
could use an ensemble estimator that serves to connect ACV and ML-BLUE
estimators where the weights of this ensemble can be optimized. Moreover, the
sample allocation can be optimized for each estimator and the ensemble
allocation can be optimized as well. Prior to analyzing this estimator, we
provide some connections to existing approaches.

Note that the GACV is in fact the most general linear ansatz for any control-
variate type variance reduction scheme. If we simply set the estimator of each
group to be a particular evaluation of a model
$\hat{Q}_{\ell}^{k}=Q_{\ell}^{k}(z_{1}^{k})$, then we allow any linear
combinations of this model. Thus, the GACV can be viewed as simply combining
individual evaluations of each of the $L+1$ models – regardless of whether
they are evaluated at the same inputs or not, and regardless of the
combinations of models evaluated at a particular input.

As just one other example available in the literature, consider the Ensemble
ACV estimator [13]. The ensemble ACV is a variant of the proposed GACV
estimator with a specific choice of weights. This estimator is obtained by
first rewriting (23) as

$\displaystyle\hat{Q}^{GACV}$
$\displaystyle=\sum_{k=1}^{K}\left(\omega_{0}^{k,e}\hat{Q}_{0}^{k}+\sum_{\ell=1}^{L}\left(\omega_{\ell}^{k,e}\hat{Q}_{\ell}^{k,e}-\omega_{\ell}^{k,\mu}\hat{Q}_{\ell}^{k,\mu}\right)\right),$
(24)

and then making the following assumptions: (1) each of the $K$ groups involves
the same models, but independent estimators; (2) the high-fidelity estimators
are equally weighted across groups, $\omega_{0}^{k,e}=\frac{1}{K}$; and (3)
the weights of a given model across all groups is identical,
$\omega_{\ell}^{k,e}=\omega_{\ell}^{k,\mu}=\frac{\bm{\alpha}_{\ell}}{K}.$
Under these assumptions we have

$\displaystyle\hat{Q}^{\text{ensemble}-ACV}(\bm{\alpha})$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\left(\hat{Q}_{0}^{k}+\sum_{\ell=1}^{L}\bm{\alpha}_{\ell}\left(\hat{Q}_{\ell}^{k,e}-\hat{Q}_{\ell}^{k,\mu}\right)\right),$
(25)

precisely the ensemble ACV estimator given by [13, Eq. 34]. We refer to that
paper for analysis of the specific estimator, and how the weights
$\bm{\alpha}$ can be chosen.

##  4 Analysis

In this section, we analyze the generalized estimator (23) by computing its
variance, and then defining an optimal selection of weights. First we note
that, for the estimator to be unbiased, we again require

$\sum_{k=1}^{K}\tilde{\bm{\beta}}_{0}^{k}=1\quad\textrm{ and
}\quad\sum_{k=1}^{K}\tilde{\bm{\beta}}_{\ell}^{k}=0,\quad\text{ for
}\ell=1,\ldots,L$ (26)

under the assumption that the estimator for each model $\ell$,
$\hat{Q}_{\ell}^{k}$, has the same bias for each $k$.

Next let us define some covariances between and within groups. First the
covariance matrix of all estimators in group $k$ is denoted by

$\mathbf{C}^{k}=\mathbb{V}\left[\bm{\hat{Q}}^{k}\right]=\begin{bmatrix}\mathbb{V}\left[\hat{Q}_{\bm{\lambda}^{k}_{1}}\right]&\mathbb{C}\text{ov}\left[\hat{Q}_{\bm{\lambda}^{k}_{1}},\hat{Q}_{\bm{\lambda}^{k}_{2}}\right]&\cdots&\mathbb{C}\text{ov}\left[\hat{Q}_{\bm{\lambda}^{k}_{1}},\hat{Q}_{\bm{\lambda}^{k}_{n_{k}}}\right]\\\
&\ddots&&\vdots\\\ &&\ddots&\\\
{\text{Sym}}&&&\mathbb{V}\left[\hat{Q}_{\bm{\lambda}^{k}_{n_{k}}}\right]\end{bmatrix}\in\mathbb{R}^{n_{k}\times
n_{k}}.$ (27)

Next, the covariance between estimators in groups $k$ and $k^{\prime}$ is
given by

$\mathbf{C}^{kk^{\prime}}=\begin{bmatrix}\mathbb{C}\text{ov}\left[\hat{Q}_{\lambda_{1}^{k}}^{k},\hat{Q}_{\lambda_{1}^{k^{\prime}}}^{k^{\prime}}\right]&\cdots&\mathbb{C}\text{ov}\left[\hat{Q}_{\lambda_{1}^{k}}^{k},\hat{Q}_{\lambda_{n_{k^{\prime}}}^{k^{\prime}}}^{k^{\prime}}\right]\\\
\vdots&\ddots&\vdots\\\
\mathbb{C}\text{ov}\left[\hat{Q}_{\lambda_{n_{k}}^{k}}^{k},\hat{Q}_{\lambda_{1}^{k^{\prime}}}^{k^{\prime}}\right]&\cdots&\mathbb{C}\text{ov}\left[\hat{Q}_{\lambda_{n_{k}}^{k}}^{k},\hat{Q}_{\lambda_{n_{k^{\prime}}}^{k^{\prime}}}^{k^{\prime}}\right]\end{bmatrix}\in\mathbb{R}^{n_{k}\times
n_{k^{\prime}}}.$ (28)

Together, these covariance matrices can be written as a single matrix

$\mathbf{C}=\begin{bmatrix}\mathbf{C}^{1}&\mathbf{C}^{12}&\cdots&&\mathbf{C}^{1K}\\\
&\mathbf{C}^{2}&\mathbf{C}^{23}&\cdots&\vdots\\\ &&\ddots&&\\\
&&&\mathbf{C}^{K-1}&\mathbf{C}^{(K-1)K}\\\
\text{Sym}&&&&\mathbf{C}^{K}\end{bmatrix}$ (29)

Note that, to obtain an optimal set of weights, we will assume that the
individual model estimators of each group are structured so that the
covariance among them is invertible.

The variance of the estimator can be written as a function of these matrices.
This is summarized by the following result.

######  Proposition 5 (Variance of the GACV).

The variance of the GACV estimator in Definition 4 is given by

$\mathbb{V}\left[\hat{Q}^{GACV}(\mathcal{B})\right]=\bm{\beta}^{T}\mathbf{C}\bm{\beta},\qquad\textrm{
where }\qquad\bm{\beta}=\begin{bmatrix}\bm{\beta}^{1}\\\ \bm{\beta}^{2}\\\
\vdots\\\ \bm{\beta}^{K}\end{bmatrix},$ (30)

and $\mathbf{C}$ is given by (29).

###### Proof.

The proof simply applies the property of covariance of a sum of random
variables.

$\displaystyle\mathbb{V}\left[\hat{Q}^{GACV}(\mathcal{B})\right]$
$\displaystyle=\sum_{k=1}^{K}\sum_{\ell\in\mathcal{S}^{k}}\sum_{k^{\prime}=1}^{K}\sum_{\ell^{\prime}\in\mathcal{S}^{k^{\prime}}}\tilde{\bm{\beta}}_{\ell}^{k}\tilde{\bm{\beta}}_{\ell^{\prime}}^{k^{\prime}}\mathbb{C}\text{ov}\left[\hat{Q}_{\ell}^{k},\hat{Q}_{\ell^{\prime}}^{k^{\prime}}\right].$
(31)

Vectorizing this expression yields the stated result. ∎

We now seek to find the optimal weights $\mathcal{B}^{*}$ that are defined to
have minimum variance and zero bias

$\mathcal{B}^{*}=\operatorname*{arg\,min}_{\mathcal{B}}\mathbb{V}\left[\hat{Q}^{GACV}(\mathcal{B})\right]\quad\text{subject
to}\quad\sum_{k=1}^{K}\tilde{\bm{\beta}}_{0}^{k}=1\quad\textrm{ and
}\quad\sum_{k=1}^{K}\tilde{\bm{\beta}}_{\ell}^{k}=0,\quad\text{ for
}\ell=1,\ldots,L.$ (32)

The solution to this problem is given by the following theorem.

######  Theorem 6 (Optimal weights of the GACV).

The solution to (32) is given by the stacked weight vector

$\bm{\beta}^{*}=\mathbf{C}^{-1}\mathbf{R}^{T}\left(\mathbf{R}\mathbf{C}^{-1}\mathbf{R}^{T}\right)^{-1}\bm{e}^{0},$
(33)

where $\bm{e}^{0}\in\mathbb{R}^{L+1}$ is the standard basis having the first
element equal to one and the rest equal to zero, $\mathbf{C}$ is given in
(29), and

$\mathbf{R}=\begin{bmatrix}\mathbf{R}^{1^{T}}&\cdots&\mathbf{R}^{k^{T}}\end{bmatrix}\in\mathbb{R}^{L+1\times\prod_{k=1}^{K}n_{k}}$
(34)

is a horizontal concatenation of all the restriction matrices. Furthermore,
the corresponding estimator variance is

$\mathbb{V}\left[\hat{Q}^{GACV}(\mathcal{B}^{*})\right]=(\bm{e}^{0})^{T}\left(\mathbf{R}\mathbf{C}^{-1}\mathbf{R}^{T}\right)^{-1}\bm{e}^{0}$
(35)

###### Proof.

Using Proposition 5, the Lagrangian form of (32) becomes

$L(\mathcal{B},\bm{\gamma})=\bm{\beta}^{T}\mathbf{C}\bm{\beta}+2\bm{\gamma}^{T}\left(\mathbf{R}\bm{\beta}-\bm{e}^{0}\right).$
(36)

With setting the Lagrangian derivative to zero,

$\nabla_{\mathcal{B}}L(\mathcal{B},\bm{\gamma})=2\mathbf{C}\bm{\beta}+2\mathbf{R}^{T}\bm{\gamma}=0,$
(37)

we obtain $\bm{\beta}=-\mathbf{C}^{-1}\mathbf{R}^{T}\bm{\gamma}.$ Plugging
this expression into the Lagrangian, we can then seek to maximize the dual

$g(\bm{\gamma})=\bm{\gamma}\mathbf{R}\mathbf{C}^{-1}\mathbf{R}^{T}\bm{\gamma}-2\bm{\gamma}^{T}\mathbf{R}\mathbf{C}^{-1}\mathbf{R}\bm{\gamma}-2\bm{\gamma}^{T}\bm{e}^{0}=-\bm{\gamma}\mathbf{R}\mathbf{C}^{-1}\mathbf{R}^{T}\bm{\gamma}-2\bm{\gamma}^{T}\bm{e}^{0}$
(38)

The maximum occurs at
$\bm{\gamma}=-\mathbf{R}\mathbf{C}^{-1}\mathbf{R}^{T}\bm{e}^{0}.$ Using this
equation in the expression for the weights, we obtain the stated result. ∎

We can specialize this result to the case of independent groups by realizing
that this change causes $\mathbf{C}$ to be block diagonal. Plugging this
structure into Theorem 6 yields the following corollary.

######  Corollary 7 (Optimal weights under independent groups).

The optimal weights for a GACV estimator that uses $K$ independent groups of
estimators is given by

$\bm{\beta}^{k}=\mathbf{C}^{k^{-1}}\mathbf{R}^{k}\left(\sum_{k=1}^{K}\mathbf{R}^{k^{T}}\mathbf{C}^{k^{-1}}\mathbf{R}^{k}\right)^{-1}\bm{e}^{0},$
(39)

and the resulting variance is

$\mathbb{V}\left[\hat{Q}^{GACV}(\mathcal{B}^{*})\right]=\bm{e}^{0^{T}}\left(\sum_{k=1}^{K}\mathbf{R}^{k^{T}}\mathbf{C}^{k^{-1}}\mathbf{R}^{k}\right)^{-1}\bm{e}^{0}.$
(40)

Note that while these results look similar to the ML-BLUE optimal weights and
variance, they are actually a generalized case because they handle arbitrary
estimators. The precise ML-BLUE estimator is recovered by assuming Monte Carlo
estimators of the expectation. In this case, the covariance of each group of
estimators $\mathbf{C}^{k}$ is directly related to the covariance of the
underlying random variables $\mathbf{\hat{C}}^{k}$ in the group through the
simple relation $\mathbf{C}^{k}=\frac{1}{m^{k}}\mathbf{\hat{C}}^{k}$. Plugging
this relation into (39) and (40) recovers the corresponding solution to the
least-squares problem described in [15, Eq. 2.7 and 2.8].

######  Corollary 8 (Optimal ML-BLUE estimator).

Let each estimator within a group be a Monte Carlo estimator for the
expectation using shared samples across the group
$\hat{Q}_{\ell}^{k}=\frac{1}{m_{k}}\sum Q_{\ell}(z_{i}^{k}).$ Furthermore, let
the covariance of all models within a group be $\mathbf{\hat{C}}^{k}$ and let
each of the $K$ groups use independent samples. Then the optimal GACV
estimator is given by

$\bm{\beta}^{k}=m^{k}\mathbf{\hat{C}}^{k^{-1}}\mathbf{R}^{k}\left(\sum_{k=1}^{K}m^{k}\mathbf{R}^{k^{T}}\mathbf{\hat{C}}^{k^{-1}}\mathbf{R}^{k}\right)^{-1}\bm{e}^{0},$
(41)

and the resulting variance is

$\mathbb{V}\left[\hat{Q}^{GACV}(\mathcal{B}^{*})\right]=\bm{e}^{0^{T}}\left(\sum_{k=1}^{K}m^{k}\mathbf{R}^{k^{T}}\mathbf{\hat{C}}^{k^{-1}}\mathbf{R}^{k}\right)^{-1}\bm{e}^{0}.$
(42)

##  5 Numerical experiments

In this section, we turn to an empirical investigation of the benefits of non-
independent group ACV estimators, as compared to the special case of an
independent group ML-BLUE estimator. Our initial aim is to simply demonstrate
existence of cases where non-independent groupings are beneficial. The
ramifications of this investigation suggest that (1) significant open
questions remain regarding how best to distribute resources in a linear
multifidelity variance reduction scheme; (2) the concept of a BLUE must be
understood to refer to the best of a certain sub-class of linear models, not a
general statement pertaining to all possible estimators formed by linear
combinations of model evaluations.

We begin this investigation in Section 5.1 by comparing two classical ACV
estimators, ACV-IS and ACV-MF, the former of which is an ML-BLUE.
Specifically, we show that there exists combinations of correlations and model
costs for which removing the assumption of grouping independence required by
ML BLUE results in smaller estimator variance. Next, in Section 5.2, we
demonstrate how model groupings of an ML-BLUE can be converted into a non-
independent grouping estimator, and that again there exists settings where the
non-independent GACV has better variance reduction.

###  5.1 ACV-IS vs. ACV-MF

As an initial example, we demonstrate the benefit of moving beyond the
independent groupings used by ML-BLUE and employing groups with shared (non-
independent) samples by comparing the ACV-IS and ACV-MF estimators. As shown
in [15], for cases with more than two models, the ACV-IS estimator is a BLUE
while the ACV-MF estimator is not. Both of these estimators are special cases
of the GACV.

The ACV-IS discussed in [15] consists of the groups
$\mathcal{S}=\\{[L]\\}\cup_{i=1}^{L}\\{\\{i\\}\\},$ and takes the form

$\hat{Q}^{ACV-
IS}=\frac{1}{m^{1}}\sum_{i=1}^{m^{1}}Q(z^{1}_{i})+\sum_{\ell=1}^{L}\alpha_{\ell}\left(\frac{1}{m^{1}}\sum_{i=1}^{m^{1}}Q_{\ell}(z^{1}_{i})-\frac{1}{m^{\ell}}\sum_{i=1}^{m^{\ell}}Q_{\ell}(z^{\ell}_{i})\right),\quad\text{
where }\mathcal{Z}^{i}\bigcap\mathcal{Z}^{j}=\emptyset\text{ when }i\neq j.$
(43)

ACV-MF [7] consists of the same groups, but the second term of the control
variate re-uses the input-samples that were used by the higher fidelity
models:

$\hat{Q}^{ACV-
MF}=\frac{1}{m^{1}}\sum_{i=1}^{N}Q(z^{1}_{i})+\sum_{\ell=1}^{L}\alpha_{\ell}\left(\frac{1}{m^{1}}\sum_{i=1}^{N}Q_{\ell}(z^{1}_{i})-\frac{1}{{m}^{\ell}}\sum_{i=1}^{{m}^{\ell}}Q_{\ell}(z^{\ell}_{i})\right),\quad\text{
where }\mathcal{Z}^{i}\subset\mathcal{Z}^{j}\text{ when }i<j.$ (44)

Because the samples are reused, the estimators in each group are not
independent and therefore ACV-MF is not a BLUE based on the definition in
[15]. Note also that the number of samples in the groups for ACV-MF is higher
than in ACV-IS, even though the overall cost is the same.

We perform several experiments on the three model setting where both the ACV-
IS and ACV-MF sets are

$\mathcal{S}=\\{\\{{0,1,2\\},\\{1\\},\\{2}\\}\\}.$

In this setting, we consider an equal cost allocation across the two
estimators. Thus, we let $n$ denote the total number of high-fidelity
evaluations, $m_{1}$ denote the total number of evaluations of the first low-
fidelity model $Q_{1}$, and $m_{2}$ denote the total number of evaluations of
the second low-fidelity model $Q_{2}$. This setup implies the following number
of samples per group: for ACV-IS, group 1 has $m^{1}=n$, group two has a
unique set $m^{2}=m_{1}-n$, and group three has a unique set $m^{3}=m_{2}-n$;
and for ACV-MF, group 1 has $m^{1}=n$, group two has all available mid-
fidelity samples $m^{2}=m_{1}$, and group three has all available low-fidelity
samples $m^{3}=m_{2}$.

We compare the variance of these estimators in Figure 1 for $n=5$, a range of
$m_{1}$ and $m_{2}$, and two different correlation settings.333The variance of
each model does not affect the ratio of estimator variances. The left panel
corresponds to a correlation structure where the correlation between $Q$ and
$Q_{2}$ is lower than the correlation between $Q_{1}$ and $Q_{3}$, whereas the
right panel has the opposite structure. The qualitative behavior between the
two cases is different. In the left panel, the benefit of ACV-MF over ACV-IS
is maximized for smaller numbers of $Q_{1}$ samples and larger numbers of
$Q_{2}$ samples. The red lines in this plot indicate unity contours which
correspond to transitions between estimator preference. Thus there exist a
large regime where ACV-MF has lower variance than ACV-IS for this case. The
right panel indicates that this behavior is highly dependent on the
correlation structure. Here, ACV-IS has lower variance over all sample
allocations, but its benefit decreases as the number of samples of $Q_{1}$
increases.

Overall, this example demonstrates that just because an estimator is a BLUE
does not mean that it makes the best use of all available resources. Clearly,
there exist cases where reusing samples across groups yields better variance
reduction for a given budget. In the next section, we show more generally how
a GACV estimator that reuses samples can be constructed automatically from
certain common group and sample allocation schemes used in ML-BLUE, and that
this construction can yield variance reduction.

(a) $\rho_{01}=0.95,\rho_{02}=0.8,\rho_{12}=0.9.$

(b) $\rho_{01}=0.95,\rho_{02}=0.93,\rho_{12}=0.9.$

Figure 1: Ratio of estimator variance $\mathbb{V}\left[\hat{Q}^{ACV-
IS}\right]/\mathbb{V}\left[\hat{Q}^{ACV-MF}\right]$ for fixed cost using a
three model setting for two choices of correlation structures as a function of
number of evaluations of each model. The $y$ axis is the number of additional
samples used for $Q_{2}$ over $Q_{1}$. Red lines indicate a variance ratio of
one, values higher than one indicate the ACV-MF estimator obtains greater
variance reduction. These results indicate that ACV-MF estimator can out-
perform the ACV-IS estimator (which is BLUE), depending on correlation and
computational cost of models.

###  5.2 Converting an arbitrary ML-BLUE estimator into a nested GACV
estimator

In this section we first provide an algorithm that converts certain ML-BLUE
sample allocations into GACV sample allocations. We then show numerical
experiments that indicate that this GACV estimator can outperform an optimal
ML-BLUE estimator.

####  5.2.1 Nested sample grouped estimator

In situations with well-correlated approximate models, the optimal groups for
the ML-BLUE estimator often follow a recursive structure. For example, in the
case where all combination of models in a group are allowed, the optimal
structure often becomes fully nested so that $\mathcal{S}^{1}=[L]$, and
$\mathcal{S}^{k}=\\{k-1,\ldots,L\\}$ for $k=2,\ldots,L+1$ [15, 16]. Because it
may be computationally expensive to search over the space of full models, the
SAOB-M estimator was introduced in the same work to limit each group to at
most $1\leq M\leq L+1$ models. In this case, there still exists a specific
type of nested structure $\mathcal{S}^{1}=[M-1]$, and
$\mathcal{S}^{k}=\\{k-1,\ldots,\min(k+M-2,L)\\}$ for $k=2,\ldots,L+1$ — see
e.g., [16, Fig. 6]. Note that SAOB-L would correspond to the fully nested
estimator. A second viewpoint of the same structure is that the $\ell$-th
model belongs to groups
$\mathcal{S}^{\max(\ell-M,1)},\ldots,\mathcal{S}^{\ell+1}.$ As an example, for
a SAOB-3 estimator for $L=4$, the groups are

$\mathcal{S}^{1}=\\{0,1,2\\},\quad\mathcal{S}^{2}=\\{1,2,3\\},\quad\mathcal{S}^{3}=\\{2,3,4\\},\quad\mathcal{S}^{4}=\\{3,4\\},\quad\mathcal{S}^{5}=\\{4\\}.$
(45)

In this section, we provide a non-independent GACV estimator that has these
identical groupings, but re-distributes the sample allocation amongst non-
independent groups to achieve an estimator of equivalent cost. Consider a
(potentially optimal) sample allocation scheme for an ML-BLUE estimator that
suggests $m^{k}$ samples for group $\mathcal{S}^{k}.$ Our nested estimator
will determine new sample allocations $\hat{m}^{k}$ that respects a certain
nested sample structure. Specifically, if $\hat{Z}^{k}$ are the samples used
for group $k$, then we construct GACV estimators that use
$\mathcal{Z}^{k+1}\subset\mathcal{Z}^{k}$. Moreover, each of these groupings
will now have $\hat{m}^{k}=\left|\mathcal{Z}^{k}\right|$ samples.

To enable a valid comparison, we must ensure that the total number of
evaluations of each model is consistent between the ML-BLUE estimator and the
GACV estimator with equivalent model groups but a nested sample structure that
induces dependencies. In other words, we construct a mapping
$\hat{m}^{k}=f_{k}(m^{1},\ldots,m^{k})$ such that the total number of
evaluations of each model are equal. With this goal, let $n^{\ell}$ define the
total number of samples of model $\ell$ such that for a ML-BLUE allocation

$n^{\ell}=\sum_{k=1,\
\ell\in\mathcal{S}^{k}}^{L+1}m^{k}=\sum_{k=\max(\ell-M+2,1)}^{\ell+1}m^{k}$
(46)

and for the associated GACV allocation

$n^{\ell}=\max\\{\hat{m}^{k}:\ \text{s.t.}\
\ell\in\mathcal{S}^{k}\\}=\hat{m}^{\ell+1},$ (47)

for $\ell=0,\ldots,L.$ The last equality for the GACV allocation exploits the
nested property of the input sets for that estimator. We begin with
$\hat{m}^{1}=m^{1}.$ Then Equating (46) with (47) we obtain

$\hat{m}^{\ell+1}=\sum_{k=\max(\ell-M+2,1)}^{\ell+1}m^{k},\quad\ell=1,\ldots,L.$
(48)

Table 1 shows a concrete example for a prototypical example of a 5 model case
with a maximum grouping of three models $(M=3)$.

Table 1: Example SAOB-3 allocation for a 5 model case ($L=4$). The ML-BLUE allocation is [5,5,5,7,18] for groups $\mathcal{S}^{1}$ through $\mathcal{S}^{5}$. The number of evaluations of each model in each group is shown for both the ML-blue allocation and the equivalent GACV allocation using (48). Each element shows the ML-BLUE (MLB) allocation and the GACV allocation (GACV) for each model in each group. The bottom row shows that the sum of the number evaluations is equal using (46) and (47). Group Model | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---
| MLB | GACV | MLB | GACV | MLB | GACV | MLB | GACV | MLB | GACV
$\mathcal{S}^{1}$ | 5 | 5 | 5 | 5 | 5 | 5 | 0 | 0 | 0 | 0
$\mathcal{S}^{2}$ | 0 | 0 | 5 | 10 | 5 | 10 | 5 | 10 | 0 | 0
$\mathcal{S}^{3}$ | 0 | 0 | 0 | 0 | 5 | 15 | 5 | 15 | 5 | 15
$\mathcal{S}^{4}$ | 0 | 0 | 0 | 0 | 0 | 0 | 7 | 17 | 7 | 17
$\mathcal{S}^{5}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 18 | 30
Number of evals $n^{\ell}$ | 5 | 5 | 10 | 10 | 15 | 15 | 17 | 17 | 30 | 30

To summarize, an allocation for a SAOB-M grouping can be converted into a
concrete GACV nested sample estimator according to Algorithm 1. Thus, this
algorithm can serve as a drop-in replacement for ML-BLUE estimators with this
structure.

Algorithm 1 Nested-sample GACV Monte Carlo estimator conversion from SAOB-M
estimator

1:Input

2: $M$ Maximum size of model group in SAOB-M estimator

3: $m^{k}$ for $k=1,\ldots,L+1$ Sample allocations for each SAOB-Mgroup

4:Output

5: $\hat{Q}_{\ell}^{k}$ for $\ell=0,\ldots,L$ and $k=1,\ldots,L+1$ Nested-
sample estimators within the GACV (23)

6: $\hat{m}^{k}$ for $k=1,\ldots,L+1$ Number of samples in equivalent cost
nested GACV estimator

7:Compute the number of evaluations of each group

8:$\hat{m}^{1}=m^{1}$

9:for $\ell=1$ to $L$ do

10: $\hat{m}^{\ell+1}\leftarrow$ Equation (48)

11:end for

12:$\bm{z}=\\{z_{1},\ldots,z_{\hat{m}^{L+1}}\\}$ $\triangleright$ Generate all
inputs to be used

13:Compute the estimators in each group

14:for $\ell=0$ to $L$ do

15: for $k=1$ to $L+1$ do

16: if $\ell\in\mathcal{S}^{k}$ then

17:
$\hat{Q}_{\ell}^{k}=\frac{1}{\hat{m}^{k}}\sum_{a=1}^{\hat{m}^{k}}Q_{\ell}(z_{a})$

18: end if

19: end for

20:end for

The variance of this estimator is given by the following corollary.

######  Corollary 9 (Nested-sample GACV estimator from independent-sample ML
BLUE).

Let $\mathbf{\hat{C}}^{k}$, $k=1,\ldots,K$ denote the covariance between the
models within group $k$ and $\mathbf{\hat{C}}^{kk^{\prime}}$ denote the
covariance of models between models of group $k$ and $k^{\prime}$. Then the
optimal GACV estimator computed using Algorithm 1 is given by optimal weights
(33) with corresponding variance (35), where the covariance matrix is

$\mathbf{C}=\begin{bmatrix}\frac{1}{\hat{m}^{1}}\mathbf{\hat{C}}^{1}&\frac{1}{\max(\hat{m}^{1},\hat{m}^{2})}\mathbf{\hat{C}}^{12}&\cdots&&\frac{1}{\max(\hat{m}^{1},\hat{m}^{K})}\mathbf{\hat{C}}^{1K}\\\
&\frac{1}{\hat{m}^{2}}\mathbf{\hat{C}}^{2}&\frac{1}{\max(\hat{m}^{2},\hat{m}^{3})}\mathbf{\hat{C}}^{23}&\cdots&\vdots\\\
&&\ddots&&\\\
&&&\frac{1}{\hat{m}^{K-1}}\max{\hat{C}}^{K-1}&\frac{1}{\max(\hat{m}^{K-1},\hat{m}^{K})}\mathbf{\hat{C}}^{(K-1)K}\\\
\text{Sym}&&&&\frac{1}{\hat{m}^{K}}\mathbf{\hat{C}}^{K}\end{bmatrix}$ (49)

###### Proof.

The Monte Carlo sampling strategy implies that the covariance of the
estimators are directly related to the covariance of the underlying models.
Specifically, The covariances in $\mathbf{C}^{kk^{\prime}}$ become

$\mathbb{C}\text{ov}\left[\hat{Q}_{\lambda_{i}^{k}}^{k},\hat{Q}_{\lambda_{j}^{k^{\prime}}}^{k^{\prime}}\right]=\mathbb{C}\text{ov}\left[\frac{1}{\hat{m}^{k}}\sum_{a=1}^{\hat{m}^{k}}Q_{\lambda_{i}^{k}}(z_{a}),\frac{1}{\hat{m}^{k^{\prime}}}\sum_{a=1}^{\hat{m}^{k^{\prime}}}Q_{\lambda_{j}^{k^{\prime}}}(z_{a})\right]$
(50)

Assume without loss of generality that $\hat{m}^{k}\leq\hat{m}^{k^{\prime}}$.
Then we have

$\displaystyle\mathbb{C}\text{ov}\left[\hat{Q}_{\lambda_{i}^{k}}^{k},\hat{Q}_{\lambda_{j}^{k^{\prime}}}^{k^{\prime}}\right]$
$\displaystyle=\mathbb{C}\text{ov}\left[\frac{1}{\hat{m}^{k}}\sum_{a=1}^{\hat{m}^{k}}Q_{\lambda_{i}^{k}}(z_{a}),\frac{1}{\hat{m}^{k^{\prime}}}\sum_{a=1}^{\hat{m}^{k}}Q_{\lambda_{j}^{k^{\prime}}}(z_{a})+\frac{1}{\hat{m}^{k^{\prime}}}\sum_{a=\hat{m}^{k}+1}^{\hat{m}^{k^{\prime}}}Q_{\lambda_{j}^{k^{\prime}}}(z_{a})\right]$
(51)
$\displaystyle=\frac{1}{\hat{m}^{k}\hat{m}^{k^{\prime}}}\mathbb{C}\text{ov}\left[\sum_{a=1}^{\hat{m}^{k}}Q_{\lambda_{i}^{k}}(z_{a}),\sum_{a=1}^{\hat{m}^{k}}Q_{\lambda_{j}^{k^{\prime}}}\right]=\frac{1}{\max(\hat{m}^{k},\hat{m}^{k^{\prime}})}\mathbb{C}\text{ov}\left[Q_{\lambda_{i}^{k}},Q_{\lambda_{j}^{k}}\right],$
(52)

while the covariances in $\mathbf{C}^{k}$ are the same as those in Corollary
8. In this case we would have
$\mathbf{C}^{kk^{\prime}}=\frac{1}{\max(\hat{m}^{k},\hat{m}^{k^{\prime}})}\mathbf{\hat{C}}^{kk^{\prime}}.$
Combining these results yields the stated covariance. ∎

####  5.2.2 Numerical demonstration

We now provide an empirical comparison between a SAOB-M estimator and the
derived nested-sample GACV estimator. Because the relative performance of the
ML BLUE and nested group ACV estimators depends on the covariance between
models and the relative computational cost of each model, we compared the
variance of these two estimators for 1000 randomly generated covariance
matrices $\mathbf{\hat{C}}$ and cost vectors $\bm{w}$. The covariance matrices
were generating using the scipy.stats [17] function random_correlations with
arguments that generated correlations roughly between 0.8 and 1. Moreover, we
fixed the cost of the high-fidelity model to 1, and randomly drew the
computational costs of each low fidelity model I.I.D. from the log uniform
distribution over the interval $[0.01,1]$. The covariance matrices were
ordered such that the covariance between the high-fidelity model the low-
fidelity models decreases with the column index of the matrix. This ordering
was enforced to reflect the hierarchical nature of the model groups we
considered. The randomly generated costs were also sorted from highest to
lowest for the same reason.

For each model setting (covariance and cost) sampled above, we first perform
an optimal sample allocation for ML-BLUE, ensure that it assigns samples to
each of the model groups, and then compute a derived nested-sample GACV
estimator. After this procedure, we compare the theoretical estimator variance
for these allocations.

To be specific, the optimal ML-BLUE sample allocation is obtained by solving a
semi-definite-program (SDP) [2]

$\displaystyle t^{*},\bm{m}^{*}=\arg\min_{t,\bm{m}\geq
0}\quad\text{s.t}\quad\begin{bmatrix}\sum_{k=1}^{K}m^{k}\mathbf{R}^{k^{T}}\mathbf{\hat{C}}^{k^{-1}}\mathbf{R}^{k}&\bm{e}^{0}\\\
\bm{e}^{0^{T}}&t\end{bmatrix}\succeq 0,\qquad\bm{m}^{T}\bm{w}\leq
W,\qquad\bm{m}^{T}\bm{h}\geq 1,$ (53)

where $\bm{m}=[m^{1},\ldots,m^{K}]$ is the sample allocation, $W$ is the total
allowable computational cost,
$\bm{w}=[w_{0},w_{1},\ldots,w_{L}]\in\mathbb{R}^{L+1}$ is the cost of each
model; and $\bm{h}\in\mathbb{R}^{K}$ is a vector with $v_{k}=1$ if the high-
fidelity model is in $\mathcal{S}^{k}$, and $v_{k}=0$ otherwise.

The following numerical results were generated using the PyApprox package [9].
Figure 2 displays a histogram of the variance ratio of the ML-BLUE SAOB
estimator compared to the derived nested GACV estimator for various $(L,M)$
combinations. Note that $(L,L+1)$ would refer to a fully recursive grouping.
The main insight from this graph is that there exists a set of problem
conditions under which non-independent groupings with sample reuse perform
better than independent groupings without sample reuse. Under the experimental
conditions explored here, this occurs a majority of the time, with seemingly
greater benefit of GACV as $M$ is lowered compared to $L$. Not visible on this
plot, as it is hidden behind the line at unity, is the fully recursive case
with $(L,M)=(4,5)$ which yields identical performance between the ML-BLUE and
nested GACV estimators. Further theoretical investigation is needed to
identify what properties of the covariance matrix in Corollary 9 lead to this
behavior.

Finally, we see that there are also cases where the ratio of performance
between GACV and ML-BLUE reaches a factor of two or more, suggesting that it
is worthwhile to check the variance of both estimators (based on a pilot
sample or other covariance estimate) to ensure that the best estimator is used
in practical applications.

Figure 2: Histogram of the ratio of the ML-BLUE SAOB-M estimator variance to
the GACV estimator variance for several $(L,M)$ pairs. The black line denotes
equal performance between estimator where experiments to the right of the line
indicate that the variance reduction of GACV is better. The GACV often
exhibits better variance reduction, and this improvement is typically larger
when $M$ is lower compared to $L$. The fully recursive case $(4,5)$ results in
a variance ratio of 1.

##  6 Conclusion

In this paper, we have extended the approximate control variate estimator to
allow for model groupings. We have shown that this estimator encompasses all
known ML-BLUE estimators, and is able to represent any linear combination of
model samples. Using this relationship, we have derived a nested-sample GACV
estimator from any SAOB-M ML-BLUE estimator and empirically demonstrated that
having non-independent sample groupings can yield improvements in variance
reduction.

The discovery of the proposed estimator has opened a number of questions. A
number of areas for future work are clearly evident to fully realize the
potential of this class of estimators. Future work should address the issue of
optimal sample allocation and estimator design to avoid having to perform
model selection for choosing a hierarchical estimator based on a pilot set of
covariances. Moreover, new optimization approaches should be developed to
enable less restricted sampling schemes than for models with independent
groupings. New schemes should be devised to simultaneously allocation samples
and discover optimal groupings. Finally, extensions to multiple outputs and
additional statistics [4] should be made for the GACV estimator.

##  7 Acknowledgements

Alex A. Gorodetsky was supported in part by Sandia National Laboratories and,
in part by an NSF CAREER Award CMMI-2238913. John D. Jakeman was supported by
the US Department of Energy’s Office of Advanced Scientific Computing Research
program. Michael S. Eldred was supported by the National Nuclear Security
Administration’s Accelerated Strategic Computing program.

Sandia National Laboratories is a multi-mission laboratory managed and
operated by National Technology & Engineering Solutions of Sandia, LLC
(NTESS), a wholly owned subsidiary of Honeywell International Inc., for the
U.S. Department of Energy’s National Nuclear Security Administration
(DOE/NNSA) under contract DE-NA0003525. This written work is authored by an
employee of NTESS. The employee, not NTESS, owns the right, title and interest
in and to the written work and is responsible for its contents. Any subjective
views or opinions that might be expressed in the written work do not
necessarily represent the views of the U.S. Government. The publisher
acknowledges that the U.S. Government retains a non-exclusive, paid-up,
irrevocable, world-wide license to publish or reproduce the published form of
this written work or allow others to do so, for U.S. Government purposes. The
DOE will provide public access to results of federally sponsored research in
accordance with the DOE Public Access Plan.

## References

  * [1] Geoffrey F Bomarito, Patrick E Leser, James E Warner, and William P Leser.  On the optimization of approximate control variates with parametrically defined estimators.  Journal of Computational Physics, 451:110882, 2022.
  * [2] M. Croci, K.E. Willcox, and S.J. Wright.  Multi-output multilevel best linear unbiased estimators via semidefinite programming.  Computer Methods in Applied Mechanics and Engineering, 413:116130, August 2023.
  * [3] Mayeul Destouches, Paul Mycek, and Selime Gürol.  Multivariate extensions of the multilevel best linear unbiased estimator for ensemble-variational data assimilation.  arXiv preprint arXiv:2306.07017, 2023.
  * [4] Thomas O Dixon, James E Warner, Geoffrey F Bomarito, and Alex A Gorodetsky.  Covariance expressions for multi-fidelity sampling with multi-output, multi-statistic estimators: Application to approximate control variates.  arXiv preprint arXiv:2310.00125, 2023.
  * [5] Gianluca Geraci, Michael S Eldred, and Gianluca Iaccarino.  A multifidelity multilevel monte carlo method for uncertainty propagation in aerospace applications.  In 19th AIAA non-deterministic approaches conference, page 1951, 2017.
  * [6] Michael B Giles.  Multilevel monte carlo path simulation.  Operations research, 56(3):607–617, 2008.
  * [7] Alex A Gorodetsky, Gianluca Geraci, Michael S Eldred, and John D Jakeman.  A generalized approximate control variate framework for multifidelity uncertainty quantification.  Journal of Computational Physics, 408:109257, 2020.
  * [8] Abdul-Lateef Haji-Ali, Fabio Nobile, and Raúl Tempone.  Multi-index monte carlo: when sparsity meets sampling.  Numerische Mathematik, 132:767–806, 2016.
  * [9] J.D. Jakeman.  Pyapprox: A software package for sensitivity analysis, bayesian inference, optimal experimental design, and multi-fidelity uncertainty quantification and surrogate modeling.  Environmental Modelling & Software, 170:105825, 2023.
  * [10] Stephen S Lavenberg, Thomas L Moeller, and Peter D Welch.  Statistical results on control variables with application to queueing network simulation.  Operations Research, 30(1):182–202, 1982.
  * [11] Stephen S Lavenberg and Peter D Welch.  A perspective on the use of control variables to increase the efficiency of monte carlo simulations.  Management Science, 27(3):322–335, 1981.
  * [12] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger.  Optimal model management for multifidelity monte carlo estimation.  SIAM Journal on Scientific Computing, 38(5):A3163–A3194, 2016.
  * [13] Trung Pham and Alex A Gorodetsky.  Ensemble approximate control variate estimators: Applications to multifidelity importance sampling.  SIAM/ASA Journal on Uncertainty Quantification, 10(3):1250–1292, 2022.
  * [14] Reuven Y Rubinstein and Ruth Marcus.  Efficiency of multivariate control variates in monte carlo simulation.  Operations Research, 33(3):661–677, 1985.
  * [15] Daniel Schaden and Elisabeth Ullmann.  On multilevel best linear unbiased estimators.  SIAM/ASA Journal on Uncertainty Quantification, 8(2):601–635, 2020\.
  * [16] Daniel Schaden and Elisabeth Ullmann.  Asymptotic analysis of multilevel best linear unbiased estimators.  SIAM/ASA Journal on Uncertainty Quantification, 9(3):953–978, 2021\.
  * [17] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [18] Hang Yang, Yuji Fujii, Kon-Well Wang, and Alex A Gorodetsky.  Control variate polynomial chaos: Optimal fusion of sampling and surrogates for multifidelity uncertainty quantification.  International Journal for Uncertainty Quantification, 13(3), 2023\.


# Does BERT Learn as Humans Perceive?
Understanding Linguistic Styles through Lexica

Shirley Anugrah Hayati  Dongyeop Kang Lyle Ungar
University of Pennsylvania  Georgia Institute of Technology University of
Minnesota
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>research conducted at
the University of Pennsylvania

###### Abstract

People convey their intention and attitude through linguistic styles of the
text that they write. In this study, we investigate lexicon usages across
styles throughout two lenses: human perception and machine word importance,
since words differ in the strength of the stylistic cues that they provide. To
collect labels of human perception, we curate a new dataset, Hummingbird, on
top of benchmarking style datasets. We have crowd workers highlight the
representative words in the text that makes them think the text has the
following styles: politeness, sentiment, offensiveness, and five emotion
types. We then compare these human word labels with word importance derived
from a popular fine-tuned style classifier like BERT. Our results show that
the BERT often finds content words not relevant to the target style as
important words used in style prediction, but humans do not perceive the same
way even though for some styles (e.g., positive sentiment and joy) human- and
machine-identified words share significant overlap for some styles.111Our
dataset and code are available at https://github.com/sweetpeach/hummingbird

##  1 Introduction

To express their interpersonal goal and attitude, people often use different
styles in their communication. The style of a text can be as important as its
literal meaning for effective communication Hovy (1987). NLP researchers have
built many models to identify different styles in text, including politeness
Danescu-Niculescu-Mizil et al. (2013), emotion Alm et al. (2005); Mohammad et
al. (2018), and sentiment Socher et al. (2013). Recently, transformer-based
Vaswani et al. (2017) pretrained language models, such as BERT Devlin et al.
(2019), have achieved impressive performance on many NLP tasks, including
stylistic studies. However, explaining what these deep learning models learn
remains a challenge. Thus, there is a growing effort to understand how these
models behave Rogers et al. (2021); Rajagopal et al. (2021).

In this work, we attempt to understand style variation through the contrasting
words identified by humans and BERT as determining a style. Given the
subjective nature of styles, we are interested in capturing human’s inherent
perception of stylistic cues in the text and compare this with the BERT’s
“perception”. Specifically, we investigate the extent to which BERT’s word
importance, as estimated using Shapley value-based attribution scores
Mudrakarta et al. (2018), aligns with human perception in stylistic text
classification.

Figure 1:  Both humans and BERT models label the sentence (a) as “polite”,
whereas in sentence (b), the humans label it as “anger” but BERT does not.
Pink highlight: high human perception score. Blue: BERT’s important words.
Purple: the word is seen as a strong cue by both human and BERT. The darker
color means that the score for human perception or machine word importance is
higher. Best seen in color.

When humans identify styles in a text, specific words play an important role
in recognizing the style, such as hedges for identifying politeness Danescu-
Niculescu-Mizil et al. (2013). We call such words stylistic cues. For example,
in Figure 1(a), humans perceive the words “understand,” “like,” and “accept”
as strong stylistic cues for politeness. But does the BERT model learn the
same words as indicative? It turns out that although the model learns that the
word “accept” is an important feature for classifying the text as polite, it
disagrees with humans for “understand” and “like” by identifying these words
as signals for impoliteness. This leads to a concern that lexical explanation
from BERT could be unreliable and motivates us to investigate more deeply into
the lexical cues used by humans and BERT. Since styles overlap significantly
Kang and Hovy (2021), we cover multiple styles: politeness, sentiment,
offensiveness, anger, disgust, fear, joy, and sadness.

Our contributions are as follows:

  * •

This is the first comparative study to examine stylistic lexical cues from
human perception and BERT. To characterize their discrepancy, we developed a
dataset, called Hummingbird, where crowd-workers relabeled benchmarking
datasets for style classification tasks.

  * •

We found that human and BERT cues are quite different; BERT pays more
attention to content words, and word-level human labels provide more accurate
multi-style correlations than sentence-level machine predictions.

  * •

Our work differs from previous works which have generated stylistic lexica
from manually-curated seed words or thesauri Davidson et al. (2017); Mohammad
and Turney (2010); Instead, in our work, the full text is given to annotators,
providing more context to the selection of the cue words.

##  2 Collection of Human and BERT’s Importance Scores on Stylistic Words

While there are many datasets with stylistic labels, to the best of our
knowledge, there is no available dataset of stylistic texts with human labels
on the individual words that drive the human perception. Therefore, on top of
existing benchmark style datasets, we we develop Hummingbird, a new dataset
with human-identified stylistic words in those stylistic sentences .

### Dataset

We use the following datasets for our style classification tasks and as a
starting point for collecting human perception scores on lexical level:
StanfordPoliteness Danescu-Niculescu-Mizil et al. (2013), Sentiment Treebank
Socher et al. (2013), Davidson et al. (2017)’s tweet dataset for
offensiveness, and SemEval 2018’s dataset for emotion classification Mohammad
et al. (2018).

Style | Label | Interannotator | F1
---|---|---|---
Distribution | Agreement
Politeness | 22.8%(+) / 41.2%(-) | 62.8 | 69.4
Sentiment | 24.6%(+) / 53.6%(-) | 71.1 | 96.5
Offensiveness | 33.6% | 75.7 | 98.0
Anger | 35.0% | 73.5 | 82.0
Disgust | 41.6% | 71.2 | 80.7
Fear | 16.4 % | 76.1 | 84.6
Joy | 22.6% | 82.7 | 86.5
Sadness | 26.4% | 72.4 | 78.2
Table 1: Dataset statistics: % of stylistic texts labeled by annotators in 500
texts. (+) refers to polite or positive, (-) refers to impolite or negative.
Inter-annotator agreement: percent agreement scores for two or more people. F1
score is the performance of BERT models on the existing test sets.

Politeness | Positive Sentiment | Joy
---|---|---
$\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  $\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  $\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$
lovely | hilarious | disappointed | delightful | deep | shocking | excited | moved | movies∗
delightful | thank | scenes∗ | lovely | thanks | scare | love | share | managing
loving | moved | suffers | smart | fun | move | entertaining | performances | referring
smart | good | hi# | solid | deftly | absolutely | great | congrats | documentary
trouble | clear | optimism | excited | best | wow# | perfect | smile | baseball∗

Table 2: Top 5 words where humans and BERT agree or disagree.
$\uparrow$$\uparrow$: both human and BERT agree. $\uparrow$: high human
perception score but low BERT’s importance score. $\uparrow$: high BERT’s
importance score but low human perception score. BERT-only agreement includes
more content words (∗) or interjections (#) than human-only agreement.

### Human Perception Scores

To collect human perception scores, we first pick 500 stylistically-diverse
texts from the four style datasets by the following method. First, we fine-
tune BERT on the training sets of the exiting datasets using the original
train/dev/test splits. The models’ performance is shown in Table 1. We then
run each model on every development set. For example, we run a sentiment
classifier on our emotion dataset. From this, we obtain the probability score
from the model for predicting each style.

To encourage that the chosen texts exhibit diverse styles, we sort them based
on their probability scores and compute the standard deviation of these scores
across the eight styles, following Kang and Hovy (2021). We then select the 50
most polite texts, 50 most impolite texts, 50 positive texts, 50 negative
texts, 100 offensive texts, and 200 emotional texts (40 from each emotion
style), resulting in total 500 texts from the four different style datasets.

We hired 622 workers to annotate them with human perception on
Prolific222https://www.prolific.co/ from November to December 2020. We
required the workers to be in the United States and payed them an average of
$9.6/hour. Each worker was asked what styles they perceive each of the texts
to exhibit. If they think the text has certain styles, workers then highlight
the words in the text which they believe make them think the text had those
styles (Pink highlights in Figure 1).

Three workers label the same pair of sentence and style, and we take majority
voting for the style labels.333See Appendix for original dataset details.
Crowd-workers obtained an average percentage agreement of 73.2% on majority
labeling, which is a substantial agreement, for text-level as shown in Table 1
and an average percentage agreement of 27.7% for word-level agreement.

Then, for a word $w_{i}$ in a text $t=w_{1}..w_{N}$, the human perception
score is defined as:

$H(w_{i})=\frac{\sum_{j=1}^{\\#\textrm{annotators}}h_{j}(w_{i})}{\\#\textrm{annotators}}$
(1)

where $h_{j}\in{-1,0,1}$ is the score given by the $j^{\textrm{th}}$
annotator. Each annotator’s label will contribute a score of either 1 for a
word that is perceived as a positive cue, -1 for a negative cue, and otherwise
0 (neutral or no emotion).

### BERT’s Importance Scores

To obtain the word importance (attribution) scores from BERT, we first trained
BERT-based models, yielding with F1 scores in Table 1. We then use the popular
technique of layered integrated gradients Mudrakarta et al. (2018) provided by
Captum Kokhlikyan et al. (2020). This technique is a variant of integrated
gradients, an interpretability algorithm that attributes an importance score
to each input feature by approximating the integral of the gradients of the
model’s output with respect to the inputs along a straight line from given
baselines to the inputs Sundararajan et al. (2017).

Since BERT could tokenize a word $w$ into several word pieces, the importance
of a word is an average of the scores of the word pieces $x$ that make up the
word. For an input of word pieces $x$, if we have a function
$F\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n}\rightarrow[0,1]$ as a
neural network, and an input $x=(x_{1},...,x_{n})\in\mathbb{R}^{n}$, an
attribution of the prediction at input $x$ relative to a baseline input
$x^{\prime}$ is a vector
$A_{F}=(x,x^{\prime})=(a_{1},...,a_{n})\in\mathbb{R}^{n}$ where $a_{i}$ is the
attribution of $x_{i}$ to the prediction $F(x)$. We use the default setting of
Captum for the baseline input $x^{\prime}$ which is zero scalar. Finally, we
obtain [-1,1] attribution score for each token like the blue highlights in
Figure 1.

Figure 2:  Pearson’s $r$ between human and BERT for the eight styles
($p<0.001$).

##  3 Human-BERT Agreement through Lexical Analysis

We study how similar human perception and BERT’s word importance are, within
each style (intra-style) and across styles (multi-styles).

### Intra-stylistic Analyses

We measure the correlation between human perception of stylistic words and
BERT’s word importance, by computing the Pearson’s $r$ for them across all
words in the vocabulary, as shown in Figure 2. Naïve refers to our baseline
which is that we simply count word frequencies in the stylistic text. For
example, if the style is positive sentiment, for a word $w$, we computed how
many times $w$ appears for sentences labeled as “positive”. We calculated the
Pearson’s r between this word count and the sentences’ styles across all
sentences. This Pearson’s r score is the baseline score of the word importance
for word $w$.

We find that BERT’s word importances correlate more highly with human
judgements than this baseline; neither BERT nor humans rely purely on co-
occurrence frequencies. Some styles are easier to identify by both human and
BERT, such as joy and sentiment with Pearson $r$=0.288 and 0.273. The yellow
bar suggests that human-BERT agreement is higher when the word appears more
often, especially for offensiveness (0.088 vs. 0.224).

Figure 3:  Simple classification using top-N human and BERT features for all
the eight styles. Best seen in color.

We now look into which words BERT and humans agree and disagree on. Table 2
shows such words selected based on the difference of the word ranks of the
human perception score and those from BERT’s word importance. To include only
highly stylistic words, words are selected only if their scores are greater
than a threshold of 0.3. When humans and BERT agree ($\uparrow$$\uparrow$),
they attend to words that are clearly associated with the styles (e.g joy,
positive) and are general (“lovely”, “delightful”, “excited”).

In contrast, BERT often finds words that suggest contexts in which the
sentiment is likely to occur. For example, top-5 words from BERT-only
agreement ($\uparrow$) contain more content words such as “scenes” for
politeness and “movies” and “basebell” for joy than those from human-only
agreement ($\uparrow$). In particular, we see that for politeness and positive
sentiment, BERT pays more attention to interjections (e.g., “hi”, “wow”) than
humans. For offensiveness and fear in Table 4 in the Appendix, humans perceive
hashtags as important cues but BERT does not. Interestingly, humans perceive a
seemingly positive word, “charming,” as offensive while BERT does not, perhaps
missing sarcasm. These content words or words irrelevant to the target style
are mostly learned due to the biased training dataset, leading to inaccurate
prediction by the machine.

Then, we evaluate the impact of important words perceived by human and BERT in
the existing test set using a simple occurrence-based classification method.
From the ranked word list by their human perception score and BERT’s word
importance scores, we label a text as having the target style, if at least one
word in the test sentence exists in the top-N word list. For this study, we
only select words that appear three times or more in the dataset.

In Figure 3, human’s word list outperforms BERT’s for most styles, even with
this small size of annotations compared to the large size of original datasets
used for training the BERT model. Interestingly, for some negative styles
(e.g., impoliteness, negative sentiment, fear), BERT’s word list performs
better. We observe that words from offensive dataset (mostly swear words) are
more consistently labeled as impolite and negative by human annotators.
However, these words are not often seen in the original politeness and
sentiment datasets. It explains why features from BERT models which are
trained on the original, large datasets get higher F1 score. As for fear, we
found that content words, such as "facebook" and "theatre", appear in the test
data. Here we see that BERT relies on content words (topic-related words) to
help predict the style, which is fragile to out-of-domain samples.

### Multi-stylistic Analyses

As we extend our analyses on multi-style correlation from a lexical viewpoint,
we found that humans and machines give similar correlations among the styles.
For instance, joy, positive sentiment, and politeness are all positively
correlated, as are anger, disgust, and offensiveness (Figure 4). However, the
multi-style correlation strength is greater for human perceptions than for
machine importance.

Figure 4:  Pearson’s $r$ word correlation matrix across styles. The upper
triangle (blue and red) represents human perception scores, while the lower
triangle (green and brown) represents machine word importances.

The weaker correlation across styles for machines is confirmed in Figure 5,
which presents a lower-dimensional visualization for the stylistic
representation of each word. Stylistic words are more clustered in human
perception, while for BERT, the separation between highly stylistic words and
non-stylistic words is less clear. Figure 5 also shows the geometric closeness
across the style clusters, giving extra information beyond the pairwise
correlations in Figure 4. In human scores, styles cluster into two extremes:
politeness, positive sentiment, and joy to the left, and anger, negative
sentiment, offensiveness, and impoliteness to the right, with disgust, fear,
and sadness, between them. This leads to more accurate style correlation
analysis than machine-based analysis on the text level Kang and Hovy (2021).

Figure 5:  T-SNE Van der Maaten and Hinton (2008) for human (top) and machine
(bottom). Each word is represented as a vector of its perception score for the
styles in this order: politeness, sentiment, offensiveness, anger, disgust,
fear, joy, and sadness.

##  4 Conclusion

We showed that BERT’s word importances for style prediction, as calculated
using integrated gradients, correspond only very loosely with the word
importances given by human annotators. These differences likely result from
several factors: 1) Word-importances computed for words which appear rarely in
the text tend to be noisy. 2) BERT, as a contextual pretrained model, take
more context into account for deciding the style of the text while human
intuitively choose the most obvious “stylistic” words to judge the style of
the text. 3) Styles are subjective matter, so human annotators may have
different perception toward the style of a sentence.

Future Directions This work also provides a public dataset as a first step for
researchers to further investigate these issues. We plan to scale up our data
collection in their size and style types including other higher-level of
styles such as sarcasm and humor. We also explore a possibility of informing
BERT to pay more attention on human-annotated lexica.

Limitations We acknowledge that while the inter-annotator agreement for the
sentence-level style is quite high, there is a huge variation for the word-
level agreement. As a caveat, the annotators could be unreliable. We do find
that annotators label different words as being important than those that drive
BERT predictions. Note that we do not claim that BERT is “wrong” and humans
are “always reliable”; only that they are different. BERT’s important words
can help the model predict correctly, but they are perceived as stylistic
features as humans do. Studying this difference is our major goal of this
paper. We believe that if a word is perceived as “stylistic” by the majority
of people, this word can be regarded as an important cue for the model.
Learning this variability of human perception on styles could be an
interesting future work using Hummingbird.

##  5 Ethical Considerations

A full analysis of style, such as politeness or expression of anger, depends
upon the context of the utterance: who is saying it to whom in what situation.
Such analysis is beyond the scope of this work, which looks only at how the
style of the utterance is perceived without context by a small number of crowd
workers. Methods such as we have used here should be extended to look at the
more subtle contextual interpretations of style and, eventually, at the ways
in which perceived styles may differ from intended styles.

Many people have (correctly) drawn attention to the role that (mis)perceptions
of style can foster gender or racial discrimination Kang and Hovy (2021).
Closer attention to the words which drive style perception is an important
first step towards addressing such problems.

Commercial platforms such as Crystal, Grammarly, and Textio offer "style
checkers". Such software would benefit from analyses that extend the work
presented here, in that they could compare the words that human editors
suggest indicate a given style to the words that NLP methods select as most
important for recognizing different styles. Such comparisons, particularly
when contextualized, should allow construction of better software to help
writers control the effect their writing has on the people reading it.

## Acknowledgments

We would like to thank Garrick Sherman for helping with the server setup
during data collection and the anonymous reviewers for their thoughtful
comments.

## References

  * Alm et al. (2005) Cecilia Ovesdotter Alm, Dan Roth, and Richard Sproat. 2005.  Emotions from text: machine learning for text-based emotion prediction.  In _Proceedings of human language technology conference and conference on empirical methods in natural language processing_ , pages 579–586.
  * Danescu-Niculescu-Mizil et al. (2013) Cristian Danescu-Niculescu-Mizil, Moritz Sudhof, Dan Jurafsky, Jure Leskovec, and Christopher Potts. 2013.  A computational approach to politeness with application to social factors.  In _Proceedings of the 51st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 250–259, Sofia, Bulgaria. Association for Computational Linguistics.
  * Davidson et al. (2017) Thomas Davidson, Dana Warmsley, Michael Macy, and Ingmar Weber. 2017.  Automated hate speech detection and the problem of offensive language.  In _Proceedings of the International AAAI Conference on Web and Social Media_ , volume 11.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Hovy (1987) Eduard Hovy. 1987.  Generating natural language under pragmatic constraints.  _Journal of Pragmatics_ , 11(6):689–719.
  * Kang and Hovy (2021) Dongyeop Kang and Eduard Hovy. 2021.  Style is not a single variable: Case studies for cross-stylistic language understanding.  Association for Computational Linguistics.
  * Kokhlikyan et al. (2020) Narine Kokhlikyan, Vivek Miglani, Miguel Martin, Edward Wang, Bilal Alsallakh, Jonathan Reynolds, Alexander Melnikov, Natalia Kliushkina, Carlos Araya, Siqi Yan, et al. 2020.  Captum: A unified and generic model interpretability library for pytorch.  _arXiv preprint arXiv:2009.07896_.
  * Loshchilov and Hutter (2017) Ilya Loshchilov and Frank Hutter. 2017.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_.
  * Mohammad et al. (2018) Saif Mohammad, Felipe Bravo-Marquez, Mohammad Salameh, and Svetlana Kiritchenko. 2018.  SemEval-2018 task 1: Affect in tweets.  In _Proceedings of The 12th International Workshop on Semantic Evaluation_ , pages 1–17, New Orleans, Louisiana. Association for Computational Linguistics.
  * Mohammad and Turney (2010) Saif Mohammad and Peter Turney. 2010.  Emotions evoked by common words and phrases: Using Mechanical Turk to create an emotion lexicon.  In _Proceedings of the NAACL HLT 2010 Workshop on Computational Approaches to Analysis and Generation of Emotion in Text_ , pages 26–34, Los Angeles, CA. Association for Computational Linguistics.
  * Mudrakarta et al. (2018) Pramod Kaushik Mudrakarta, Ankur Taly, Mukund Sundararajan, and Kedar Dhamdhere. 2018.  Did the model understand the question? In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1896–1906, Melbourne, Australia. Association for Computational Linguistics.
  * Rajagopal et al. (2021) Dheeraj Rajagopal, Vidhisha Balachandran, Eduard Hovy, and Yulia Tsvetkov. 2021\.  Selfexplain: A self-explaining architecture for neural text classifiers.  _arXiv preprint arXiv:2103.12279_.
  * Rogers et al. (2021) Anna Rogers, Olga Kovaleva, and Anna Rumshisky. 2021.  A primer in bertology: What we know about how bert works.  _Transactions of the Association for Computational Linguistics_ , 8:842–866.
  * Socher et al. (2013) Richard Socher, Alex Perelygin, Jean Wu, Jason Chuang, Christopher D. Manning, Andrew Ng, and Christopher Potts. 2013.  Recursive deep models for semantic compositionality over a sentiment treebank.  In _Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing_ , pages 1631–1642, Seattle, Washington, USA. Association for Computational Linguistics.
  * Sundararajan et al. (2017) Mukund Sundararajan, Ankur Taly, and Qiqi Yan. 2017.  Axiomatic attribution for deep networks.  In _International Conference on Machine Learning_ , pages 3319–3328. PMLR.
  * Van der Maaten and Hinton (2008) Laurens Van der Maaten and Geoffrey Hinton. 2008.  Visualizing data using t-sne.  _Journal of machine learning research_ , 9(11).
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in Neural Information Processing Systems_ , volume 30, pages 5998–6008. NIPS.

##  Appendix A Existing Datasets for Style Classification

We use existing style datasets from StanfordPoliteness Danescu-Niculescu-Mizil
et al. (2013) for politeness, Sentiment TreeBank Socher et al. (2013) for
sentiment, Davidson et al. (2017)’s dataset for offensiveness, and SemEval
2018 Task 1: Affect in Tweets for emotion classification Mohammad et al.
(2018). We convert non-binary labels or scores to binary labels to standardize
the multi-style analysis, resulting in eight styles. Table 3 shows their
dataset sizes and train/dev/test splits.

Styles | Train | Dev | Test
---|---|---|---
Politeness | 9,859 | 530 | 567
Sentiment | 236,077 | 1,045 | 2,126
Offensiveness | 22,277 | 1,251 | 1,255
Anger | 6,839 | 887 | 3,260
Disgust | 6,839 | 887 | 3,260
Fear | 6,839 | 887 | 3,260
Joy | 6,839 | 887 | 3,260
Sadness | 6,839 | 887 | 3,260
Table 3: Dataset Statistics

StanfordPoliteness is collected from StackExchange and Wikipedia requests. The
labels are continous values of [-2, 2] so we convert it to binary labels of
“polite” and “impolite” by converting all values greater than 0 as polite and
the rest are impolite. Sentiment TreeBank dataset consists of movie review
texts, and we only use the coarse label of “positive” and “negative” labels
for training. Davidson et al. (2017) collected their data from Twitter, and we
only consider “offensive” and “none” labels. SemEval 2018 dataset is collected
from tweets and it has total 11 emotions for the same 1̃0.9k instances: anger,
anticipation, disgust, fear, joy, love, optimism, pessimism, sadness,
surprise, and trust. We select anger, disgust, fear, joy, and sadness since
these emotions have the highest F1-score compared to the rest. Each emotion
has two labels: “anger” or “not anger”, “disgust” or “not disgust”, and so on.

##  Appendix B Training Configuration

We use the lower-cased BERT-base model with 12 hidden layers, 12 attention
heads, hidden size 768, for training our style classifiers on GeForce GTX
TITAN X GPU. Drop-out rate is 0.1, learning rate is $2\times 10^{-5}$, and the
optimizer is AdamW Loshchilov and Hutter (2017). Vocabulary size is 30,522 and
max position embeddings is 512. Training ran for 3 epochs, and each epoch took
around 4 minutes.

Figure 6:  Instruction page for crowd workers.  Figure 7:  Annotation page for
crowd workers.  Figure 8:  Demographic survey for crowd workers.

Politeness | Positive Sentiment | Joy
---|---|---
$\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  $\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  $\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$
lovely | hilarious | disappointed | delightful | deep | shocking | excited | moved | movies
delightful | thank | scenes | lovely | thanks | scare | love | share | managing
loving | moved | suffers | smart | fun | move | entertaining | performances | referring
smart | good | hi | solid | deftly | absolutely | great | congrats | documentary
trouble | clear | optimism | excited | best | wow | perfect | smile | baseball
happy | share | weather | hilarious | high | optimism | loving | example | audience
charming | friend | sounds | great | pretty | news | amazing | morning | scenes
compellling | rest | genre | how’s | friday | journey | happy | how’s | fan
serious | chance | grief | #comfort | smile | nice | delightful | among | dream
Offensiveness | Anger | Fear
$\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  $\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  $\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$
negro | charming | suffers | awful | ignore | negro | despair | rage | shocking
pussy | haunting | used | insult | works | goat | loss | tired | haunting
filth | murderous | ma | bitches | throw | exam | horrific | #heartbreaking | childish
bitch | wee | those | nasty | dumb | closed | creepy | sucks | insult
shit | #porn#android… | got | #horrific | mental | baby | smashed | mental | journey
smashed | holy | made | bitch | bust | trump | #horrific | spectacle | hate
fat | worst | bad | fuck | broke | sex | feeling | #sad | doors
nasty | #dreadful | men | smashed | rage | gotta | nervous | imma | delayed
ass | cheating | ho | fucking | scare |  | dreadful | war | movies
dreadful | tf | sex | discussing | murderous |  | murderous | depression | nightmare
Disgust | Sadness |  |  |
$\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  $\uparrow$$\uparrow$ |  $\uparrow$ |  $\uparrow$ |  |  |
insult | shocking | trouble | crying | sucks | smashed |  |  |
fuck | #terrified | business | #depressing | #horrific | closed |  |  |
terrible | failing | referring | miss | absolutely | often |  |  |
sucks | #saddened | trump | #disappointment | jail | ways |  |  |
damn | #nervous | correct | sad | flight | decline |  |  |
utterly | horror | negro | shocking | #hatred | creepy |  |  |
#horrific | mad | ruins | despair | nasty | insult |  |  |
hate | cold | church | crash | ignore | clumsy |  |  |
pussy | broke | #terrifying | #sad | negro | white |  |  |
bitches | bitter | exam | ruins | gives | emotional |  |  |

Table 4: Top 10 words where humans and BERT agree and disagree for all the
eight styles. We only select words that appear $>=$ 2\. $\uparrow$$\uparrow$:
both human and BERT agree. $\uparrow$: high human perception score but low
word importance score. $\uparrow$: high word importance score but low human
perception score.

##  Appendix C Annotation Interface

For each text-style pair (total: 500 texts × 8 styles = 4,000 pairs), we ask
three different annotators to select the style label for text and highlight
463 the words which make them think the text has that style with instructions
shown in Figure 6. To guarantee that the workers are serious with this task,
we provide a screening practice session which resembles the exact task but
with a text that is very obvious to be annotated as in Figure 7. The real task
interface is also the same as Figure 7. Figure 8 displays an interface where
we also ask the worker’s demographic profile.

##  Appendix D Important Words Perceived by Humans and the Machine

Table 4 shows top twenty words where humans and BERT agree and disagree for
all styles.