# fAux: Testing Individual Fairness via Gradient Alignment

Giuseppe Castiglione1 Ga Wu1 Christopher Srinivasa1&Simon Prince1
1Borealis AI {giuseppe.castiglione, ga.wu, christopher.srinivasa,
<EMAIL_ADDRESS>
###### Abstract

Machine learning models are vulnerable to biases that result in unfair
treatment of individuals from different populations. Recent work that aims to
test a model’s fairness at the individual level either relies on domain
knowledge to choose metrics, or on input transformations that risk generating
out-of-domain samples. We describe a new approach for testing individual
fairness that does not have either requirement. We propose a novel criterion
for evaluating individual fairness and develop a practical testing method
based on this criterion which we call fAux (pronounced fox). This is based on
comparing the derivatives of the predictions of the model to be tested with
those of an auxiliary model, which predicts the protected variable from the
observed data. We show that the proposed method effectively identifies
discrimination on both synthetic and real-world datasets, and has quantitative
and qualitative advantages over contemporary methods.

##  1 Introduction

Unfair treatment of different populations by machine learning models can
result in undesired social impact Berk et al. (2017); Yucer et al. (2020).
There are three main research challenges associated with this problem. The
first is to identify the source of the bias and understand how this influences
the models Mehrabi et al. (2019); Sun et al. (2020). The second is to modify
the training strategy to prevent unfair predictions Yurochkin et al. (2020);
Ruoss et al. (2020); Yurochkin and Sun (2021). The final challenge, which is
addressed in this paper, is to test the fairness of existing ML models.

To test that a model is fair, we must first agree on what is meant by
‘fairness’. For all definitions, fairness is defined with respect to protected
variables such as race, gender, or age. However, the literature distinguishes
between group fairness (equivalent aggregate treatment of different protected
groups) and individual fairness (equivalent treatment of similar individuals
regardless of their protected group). This paper addresses individual
fairness, but even here there are multiple, potentially conflicting criteria.
For example, a common definition (see Gajane and Pechenizkiy (2017)) is
fairness through unawareness (FTU) in which the model should behave as though
the protected variable is not present. Conversely, Dwork et al. (2012)
proposed fairness through awareness (FTA), which requires that similar
individuals have similar prediction outcomes. Kusner et al. (2017) emphasized
counterfactual fairness (CFF). This takes a data example and synthesizes a
counterfactual example in which the protected variable and its descendants are
changed. It requires that the original predictions and those for the
counterfactual should be similar.

Figure 1: At the heart of our method is a simple idea; if we adjust the model
input so that the predicted protected variable changes, then the model output
should not change. a) Target model predicts $y$ from inputs $\mathbf{x}$. We
wish to test fairness at points 1,2,3,4. b) We construct an auxiliary model
that predicts protected variable $c$ from inputs $\mathbf{x}$. fAux compares
gradients (green arrows) of the two models. Point 2 is unfair because the
target and auxiliary model gradients are large and parallel; the model
prediction changes as the protected variable changes. The other points are
fair since the gradients are orthogonal (point 1) or one or other gradient is
small (points 3,4).

There exist methods to test models with all of these definitions. Agarwal et
al. (2018) and Galhotra et al. (2017) use the FTU definition, John et al.
(2020) and Wachter et al. (2018) concentrate on the FTA definition, and Black
et al. (2019) use the CFF principle. However, each approach has limitations.
Features can act as surrogates of a protected variable, which FTU ignores.
Using FTA needs domain-specific knowledge to define similarity metrics for
inputs and outputs. Using CFF requires building a generative model to produce
counterexamples. Moreover, we show experimentally that methods based on FTA
can exhibit low precision.

In this paper, we introduce fAux (pronounced fox), a new framework for
individual fairness testing which avoids these difficulties. Our contributions
are as follows:

  1. 1.

We propose a simple and general criterion for individual fairness testing: for
a fair model, the derivative of the model prediction with respect to the
protected variable should be small.

  2. 2.

We introduce the concept of an auxiliary model which describes the
relationship between the input features and the protected variable and show
how to use this to help evaluate this criterion. Figure 1 provides the
intuition behind our approach. We specifically focus on unfair treatment that
is created by historical bias in datasets Mehrabi et al. (2019).

  3. 3.

To help evaluate fairness testing, we present a novel synthetic data
generation method that merges multiple real datasets through a probabilistic
graphical model to flexibly simulate realistic data with controllable bias
levels.

##  2 Preliminary Material

In this section, we describe a definition of individual fairness that unites a
number of definitions in the literature. We then analyze the technical
challenges of existing tests, which serve as motivation for our own.
Throughout, we adopt the following notation:

  * •

$X$: Feature (input) variables. When features are observed, we use
$\mathbf{x}$ to represent the feature vector.

  * •

$Y$: Prediction (output) variables. When a label is observed, we use $y$ to
represent the label as a scalar. As we do fairness testing on binary
classification tasks, the prediction $y$ in this paper is a probability.

  * •

$C$: Protected variables (e.g., gender). We use $\mathbf{c}$ to represent the
observed values.

  * •

$\phi$: Distance metric. $\phi_{in}(\cdot,\cdot)$ denotes a metric of input
space, and $\phi_{out}(\cdot,\cdot)$ a metric of output space.

  * •

$f_{tar}$: Target function for fairness testing. This takes features
$\mathbf{x}$ as input and produces a prediction $\hat{y}$.

  * •

$f_{aux}$: Auxiliary model. This takes features $\mathbf{x}$ as input and
produces predictions $\mathbf{c}$ for protected attributes $C$.

###  2.1 Individual Fairness Definitions

Individual fairness describes the tolerable discrepancy of model predictions
at the level of individual data points. There are several characterizations
that operate at this level, and to facilitate their comparison we use the
following definition. Let the observed features $\mathbf{x}$ be generated from
underlying latent variables $\mathbf{z}_{\bot}$ and $\mathbf{z}_{\parallel}$
via a function $f_{g}$:

$\mathbf{x}=f_{g}(\mathbf{z}_{\bot},\mathbf{z}_{\parallel}).$ (1)

Here, $\mathbf{z}_{\bot}$ denotes latent vectors with no correlation with the
protected variables $\mathbf{c}$, and $\mathbf{z}_{\parallel}$ is influenced
by the protected variables $C$ through a function
$\mathbf{z}_{\parallel}=\psi(\mathbf{c})$.

######  Definition 1.

A model $f_{tar}$ is individually fair, if it produces exactly identical
outcomes when given input feature vectors $\mathbf{x}_{i}$ and
$\mathbf{x}_{j}$ which share the same latent vector $\mathbf{z}_{\bot}$:

$f_{tar}(\mathbf{x}_{i})=f_{tar}(\mathbf{x}_{j}),$ (2)

where $\mathbf{x}_{i}=f_{g}(\mathbf{z}_{\bot},\psi(\mathbf{c}))$ and
$\mathbf{x}_{j}=f_{g}(\mathbf{z}_{\bot},\psi(\mathbf{c}^{\prime}))$.

At a high level, the definition states that given a pair of similar data
points $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$, they should receive equal
treatment from a target model. This perspective is shared by the three most
common families of individual fairness tests: Fairness Through Unawareness
(FTU) stipulates that protected variables should not explicitly be used in the
prediction process. “Similar” points are thus points that differ only in the
protected variable, since the latter cannot affect the outputs. Fairness
Through Awareness (FTA) formally define “similar” points using a metric on the
input space Dwork et al. (2012). Finally, Counterfactual Fairness (CFF)
stipulates that the model predictions should not causally depend on protected
variables $C$ Kusner et al. (2017). For a given input $\mathbf{x}_{i}$, the
“similar” individual $\mathbf{x}_{j}$ is the counterfactual constructed by
intervening on $C$.

For these three tests the main technical challenge is the generation of
“similar” pairs, either by searching within a neighbourhood defined by a
metric, or by transformation of the input. We now describe these challenges in
more depth.

####  2.1.1 Technical Challenges in Employing Similarity Metrics

The main challenge to constructing similarity metrics is the task-specific
domain knowledge required Dwork et al. (2012). Some methods simply employ
unweighted $l_{p}$ norms Wachter et al. (2018); John et al. (2020). Others
obtain weights from a linear model trained to predict the protected variable
$c$ from the input $\mathbf{x}$ Ruoss et al. (2020); Yurochkin et al. (2020).
Others still learn metrics from expertly-labelled pairs Mukherjee et al.
(2020).

####  2.1.2 Technical Challenges in Transforming Inputs

With regards to input transformations, the main challenges are stability, and
the risk of generating out-of-distribution samples. Both risks are faced by
approaches that employ adversarial techniques Wachter et al. (2018); Maity et
al. (2021). When the dataset is small, and lower dimensional, Optimal
Transport (OT) can be used to define a mapping between pairs of inputs based
on pairwise-similarity Dwork et al. (2012); Gordaliza et al. (2019). However,
this method scales poorly. To this end, approximate OT methods, based on dual-
formulations Chiappa and Pacchiano (2021) and Generative Adversarial Networks
Black et al. (2019) have also been explored - though again, these present
risks of instability and out-of-distribution samples.

CFF also employs generative models, specifically, causal graphical models.
However, such graphs are rarely available. Moreover, training such models via
unsupervised learning is hard Srivastava et al. (2017), increasing the risk
that the generated inputs are out-of-domain, or have limited coverage.

###  2.2 Mechanisms of Discrimination

We now describe the mechanism of discrimination our test aims to detect. There
are many reasons why models exhibit unfair behaviour, but one of the most
insidious is the mishandling of historical bias Mehrabi et al. (2019). Here,
pre-existing prejudices create spurious correlations between features $X$ and
protected variables $C$111More in-depth discussion and worked examples may be
found in Section Use Cases of the supplementary materials. (Figure 2(a)).
These correlations provide a model with two possible inference paths: one
legitimate, and one that implicitly infers $C$. The latter paths render the
FTU definition unreliable, as models learn to exploit surrogate features even
when $C$ is omitted Barocas et al. (2019). By contrast, depending on the
metric, FTA is sensitive to all variations in the input, and thus may flag
instances where inference was legitimate.

In the next section, we present an approach that can precisely distinguish
between these two paths. It does so scalably and operates within distribution.
Moreover, there is lower overhead, as it requires no domain-expertise, and it
employs only supervised learning techniques.

##  3 Local Fairness Tests with fAux

(a) Data Generation

(b) Fair Model

(c) Unfair Model

Figure 2: Graphical models to show how historical bias causes unfairness. (a)
Generation process of a biased training dataset. Red dashed line denotes that
$Y$ and $C$ may have correlations due to historical bias. (b) Fair models
learn to infer $Y$ while cancelling the impact from protected variables $C$.
Solid arrows show generative dependence. Dashed arrow shows learned inference
mapping. (c) Unfair models infer protected variables to support their
predictions.

We now develop a novel fairness testing method that satisfies Definition 1,
but does not suffer from the limitations described previously. We start by
using the graphical model in Figure 2 to propose a criterion for individual
fairness based on sensitivity analysis. We use this to motivate the Local
Independence Criterion (LIC) which examines whether a model suffers from
historical bias, and show that this satisfies Definition 1. Finally, we
introduce the idea of an Auxiliary Model which is needed to create a practical
test.

###  3.1 Local Independence Criterion

A sufficient condition for a model to violate the individual fairness
criterion of Definition 1 is that its prediction depends on the protected
variable $\mathbf{c}$. We may reveal this dependence locally using the partial
derivative

$\frac{\partial f_{tar}(\mathbf{x})}{\partial\mathbf{c}}\neq 0,$ (3)

which indicates that the prediction is sensitive to a small perturbation of
$\mathbf{c}$. In practice, the machine learning model and data introduce
inevitable noise, thus we relax the above expression with a pre-defined
tolerance $\delta$. We then obtain:

######  Theorem 1.

Let there exist a generative model $f_{g}$ that influences features $X$ with
protected variables $C$, such that
$\mathbf{x}=f_{g}(\mathbf{z}_{\bot},\psi(\mathbf{c}))$. If a machine learning
model $f_{tar}$ violates the Local Independence Criterion (LIC)

$\displaystyle\left|\frac{\partial
f_{tar}(\mathbf{x})}{\partial\mathbf{c}}\right|_{\infty}\leq\delta,$ (4)

by a pre-defined threshold $\delta$, then the model $f_{tar}$ violates the
individual fairness criteria in Definition 1.

A proof is provided in Appendix A.4, but intuitively, the partial derivative
considers the disparate treatment of two infinitesimally-close individuals.

To use the LIC, we need to estimate the derivative $\partial
f_{tar}(\mathbf{x})/\partial\mathbf{c}$, for which the chain rule gives:

$\left|\frac{\partial
f_{tar}(\mathbf{x})}{\partial\mathbf{x}}\frac{\partial\mathbf{x}}{\partial\mathbf{c}}\right|_{\infty}\leq\delta.$
(5)

Unfortunately, the term $\partial\mathbf{x}/\partial\mathbf{c}$ is undefined
without accessing the underlying generative model that maps protected
variables $C$ to the features $X$ and this is rarely available222Note, the
protected variables $C$ are not necessarily continuous as we will model the
mapping through an auxiliary model later. .

###  3.2 Auxiliary Models

In this section we suggest an approximation of
$\partial\mathbf{x}/\partial\mathbf{c}$ that requires neither generative model
nor attempts to model the latent representations $\mathbf{z}\in Z$.

One approach would be to build a model to predict $\mathbf{x}$ from
$\mathbf{c}$ and use the derivative of this model to approximate
$\partial\mathbf{x}/\partial\mathbf{c}$. However, this might be poor as the
number of protected variables is often far smaller than the feature size.
Instead, we consider the inverse problem: we describe the mapping from
features $X$ to the protected variables $C$ using an auxiliary model
$\mathbf{c}=f_{aux}(\mathbf{x})$. We then invert this model in a local
neighbourhood around a given point $\mathbf{x}_{0}$, to approximate the
desired derivative. To this end, given this auxiliary model $f_{aux}$ we apply
the Taylor expansion around $(\mathbf{x}_{0},\mathbf{c}_{0})$:

$\mathbf{c}-f_{aux}(\mathbf{x}_{0})\approx\left(\frac{\partial
f_{aux}(\mathbf{x}_{0})}{\partial\mathbf{x}_{0}}\right)^{\top}(\mathbf{x}-\mathbf{x}_{0})$
(6)

where we replaced $\mathbf{c}_{0}$ with its prediction
$f_{aux}(\mathbf{x}_{0})$. The left-hand side denotes the change in the space
of protected variables. The right hand side is a Jacobian vector product.

We apply the Moore-Penrose pseudo-inverse to find the minimum norm solution
for $\mathbf{x}$:

$f_{aux}^{-1}(\mathbf{c})=\mathbf{x}_{0}+(\mathbf{c}-f_{aux}(\mathbf{x}_{0}))\left(\nabla
f_{aux}^{\top}\nabla f_{aux}\right)^{-1}\nabla f_{aux}^{\top}$ (7)

where we use $\nabla f_{aux}$ to denote $\partial
f_{aux}(\mathbf{x}_{0})/\partial\mathbf{x}_{0}$. This allows us to approximate
$\partial\mathbf{x}/\partial\mathbf{c}$ by

$\frac{\partial\mathbf{x}}{\partial\mathbf{c}}\approx\frac{\partial
f_{aux}^{-1}(\mathbf{c})}{\partial\mathbf{c}}=\left(\nabla
f_{aux}^{\top}\nabla f_{aux}\right)^{-1}\nabla f_{aux}^{\top}.$ (8)

Finally, by combining Equation 8 and the Equation 5, we can approximate the
LIC with:

$\left|\nabla f_{tar}\left(\nabla f_{aux}^{\top}\nabla
f_{aux}\right)^{-1}\nabla f_{aux}^{\top}\right|_{\infty}\leq\delta$ (9)

where we use $\nabla f_{tar}$ to denote $\partial
f_{tar}(\mathbf{x}_{0})/\partial\mathbf{x}_{0}$.

###  3.3 On the Choice of the Auxiliary Model

In this paper, we employ Multilayer Perceptron (MLP) architectures for our
auxiliary models, to minimize the amount of model overhead. Though inverting
such models can potentially yield low-fidelity reconstructions of
$\mathbf{x}$, additional fidelity may require modelling factors of
$\mathbf{x}$ that are independent of $\mathbf{c}$. The LIC avoids the need for
these additional factors since the end goal is estimating the partial
derivative of $\frac{\partial\mathbf{x}}{\partial c}$ only (for more details,
see Appendix A.5). In our experiments, we demonstrate that even such simple
architectures are sufficient to achieve state of the art results.

Beyond their simplicity, another advantage of using MLPs is their flexibility,
as they easily accommodate both real-valued and categorical outputs. This
makes it possible to analyze scenarios in which the protected variable $C$ is
not binary.

###  3.4 Relaxations and Extensions

While the basic fAux described above is sufficient for an individual fairness
test, it depends heavily on the behaviour of $\nabla f_{tar}$, which may be
ill-conditioned. In this section, we introduce several variants of the basic
fAux method.

Normalization of Gradient (fAux+NG): Different features may have very
different valid ranges and so the gradient of either the target or the
auxiliary model could be biased towards a subset of features. To mitigate this
problem, we use an $l_{2}$ normalization of the gradients to give the
criterion:

$\left|\mathrm{norm}(\nabla f_{tar})\;\mathrm{norm}(\nabla
f_{aux}^{\top})\right|_{\infty}\leq\delta.$ (10)

where we removed the inverse term as this normalization is longer needed.

Integrated Gradient (fAux+IG): We substitute the raw gradients $\nabla
f_{tar}$ and $\nabla f_{aux}$, for integrated gradients Sundararajan et al.
(2017) that provide a smoothed gradient signal.

##  4 Experiments

We now evaluate the proposed fAux test to answer the following research
questions (RQs):

  * RQ1: Given the target model $f_{tar}$ trained on synthetic datasets whose ground-truth degrees of bias are known, how well does fAux perform compared to other testing methods?

  * RQ2: Given the target model $f_{tar}$ trained on the real dataset whose ground-truth degree of bias is unknown, can fAux identify discriminatory features?

  * RQ3: How efficient is fAux compared to the existing approaches in terms of inference cost?

  * RQ4: Does fAux have any conditions needed to guarantee reliable test performance? In particular, we want to know how the effect of auxiliary model performance would impact the test performance.

###  4.1 Experimental Setup

####  4.1.1 Candidate Testing Algorithms and Target Model

  * FTA: A local version of FTA, inspired by works such as John et al. (2020). In these approaches, an $\epsilon$-neighbourhood around the input $x_{i}$ is rigorously searched over to bound the output deviation $\phi_{out}(f_{tar}(\textbf{x}_{i}),f_{tar}(\textbf{x}_{j}))$. In the limit that $\epsilon$ goes to zero, this is equivalent to a bound on the $l_{p}$ norm of $\nabla f_{tar}$. We also construct a weighted $l_{p}$-norm using a linear auxiliary model, as is common in the individual fairness literature Ruoss et al. (2020); Yurochkin et al. (2020).

  * Unfair Map: Maity et al. (2021) uses a gradient-flow attack to generate pairs of points that violate FTA. This attack is conducted within a neighbourhood defined by a similarity metric $\phi_{in}$, and to this end we employ the same weighted $l_{p}$ norm used in FTA.

  * FlipTest: Black et al. (2019) approximates CFF by leveraging Wasserstein GANs to generate pairs of inputs.

  * LIC-UB: An upper-bound on the Local Independence Criterion (LIC). For the experiments on synthetic datasets, we can compute the true gradient of the generative model to conduct the LIC check (5). Since error introduced by the approximation is removed, the test performance should achieve its upper-bound.

  * Target Models: Given datasets in the form $D=\\{\cdots(\mathbf{x},y,\mathbf{c})\cdots\\}$, we train multi-layer fully connected networks as the target models $f_{tar}$ with only features $\mathbf{x}$ and label $y$. The target models in our experiments are all classifiers that aim to produce probabilistic predictions $P(y|\mathbf{x})$. However, as previously discussed, an unfair model may infer protected variable $C$, resulting in it implicitly modeling $P(y|\mathbf{x},\mathbf{c})$.

| Synthetic-1 | Synthetic-2 | Synthetic-3 | Synthetic-4 | Synthetic-5 | Synthetic-6 | Synthetic-7 | Synthetic-8
---|---|---|---|---|---|---|---|---
xc/xy |  | magic/
---
backache
| magic/
---
backache
| magic/
---
backache
| magic/
---
backache
| australian/
---
credit
| australian/
---
credit
| australian/
---
credit
| australian/
---
credit
Bias Level | 0.5 | 0.5 | 1 | 1 | 0.5 | 0.5 | 1 | 1
Fusion Approach | outer | concat | outer | concat | outer | concat | outer | concat
FTA |  | 0.307 $\pm$ 0.000
---
| 0.306 $\pm$ 0.000
---
| 0.307 $\pm$ 0.000
---
| 0.307 $\pm$ 0.000
---
| 0.466 $\pm$ 0.000
---
| 0.293 $\pm$ 0.000
---
| 0.355 $\pm$ 0.000
---
| 0.309 $\pm$ 0.000
---
FTA + lin. aux | 0.402 $\pm$ 0.000 | 0.683 $\pm$ 0.003 | 0.451 $\pm$ 0.001 | 0.455 $\pm$ 0.001 | 0.709 $\pm$ 0.000 | 0.658 $\pm$ 0.00 | 0.672 $\pm$ 0.000 | 0.612 $\pm$ 0.000
Unfair Map | 0.407 $\pm$ 0.000 | 0.400 $\pm$ 0.000 | 0.576 $\pm$ 0.000 | 0.820 $\pm$ 0.000 | 0.571 $\pm$ 0.000 | 0.653 $\pm$ 0.000 | 0.687 $\pm$ 0.000 | 0.732 $\pm$ 0.000
FlipTest | 0.598 $\pm$ 0.285 | 0.880 $\pm$ 0.015 | 0.433 $\pm$ 0.259 | 0.633 $\pm$ 0.109 | 0.600 $\pm$ 0.148 | 0.678 $\pm$ 0.056 | 0.701 $\pm$ 0.209 | 0.732 $\pm$ 0.132
fAux |  | 0.332 $\pm$ 0.004
---
| 0.998 $\pm$ 0.001
---
| 0.311 $\pm$ 0.002
---
| 1.000 $\pm$ 0.000
---
| 0.615 $\pm$ 0.032
---
| 0.937 $\pm$ 0.005
---
| 0.564 $\pm$ 0.014
---
| 0.997 $\pm$ 0.001
---
| fAux+NG
---
| 0.876 $\pm$ 0.015
---
| 0.999 $\pm$ 0.001
---
| 0.978 $\pm$ 0.023
---
| 1.000 $\pm$ 0.000
---
| 0.815 $\pm$ 0.015
---
| 0.947 $\pm$ 0.006
---
| 0.910 $\pm$ 0.019
---
| 0.998 $\pm$ 0.001
---
| fAux+IG
---
| 0.717 $\pm$ 0.031
---
| 0.998 $\pm$ 0.000
---
| 0.937 $\pm$ 0.012
---
| 1.000 $\pm$ 0.000
---
| 0.886 $\pm$ 0.015
---
| 0.944 $\pm$ 0.002
---
| 0.979 $\pm$ 0.007
---
| 0.999 $\pm$ 0.001
---
LIC-UB |  | 0.998
---
| 0.999
---
| 1.000
---
| 1.000
---
| 0.965
---
| 0.966
---
| 0.999
---
| 0.999
---

Table 1: Performance Comparison among Individual Fairness Testing Methods on
Synthetic Datasets We report Average-Precision scores with the highest score
in bold font. The confidence interval comes from 10 runs by re-training
auxiliary models. Rows are sorted in order of increasing computational
requirements. The last row corresponds to a theoretical upper bound.

Full implementation details are found in Appendix A.11.

####  4.1.2 Synthetic Datasets with Ground Truth Bias

In producing our synthetic datasets, our goal was to retain the noisy and
nonlinear relationships that are present in real datasets. To this end, we
construct a pipeline which joins real, unbiased datasets together via fusion
operations, based on a intentionally biased data sampling process.
Specifically, given two datasets
$\hat{D}=\\{\cdots(\hat{\mathbf{x}}_{i},\hat{y}_{i})\cdots\\}$ and
$\tilde{D}=\\{\cdots(\tilde{\mathbf{x}}_{j},\tilde{y}_{j})\cdots\\}$, and a
fusion operation $f_{fus}$, we can produce a synthetic dataset $D_{syn}$ such
that
$D_{syn}=\\{\cdots(\mathbf{x},y,c)\cdots\\}=\\{\cdots(f_{fus}(\hat{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{j}),\hat{y}_{i},\tilde{y}_{j})\cdots\\}$,
where $y\stackrel{{\scriptstyle\text{def}}}{{=}}\hat{y}_{i}$,
$c\stackrel{{\scriptstyle\text{def}}}{{=}}\tilde{y}_{j}$, and $f_{fus}$ is a
fusion operation (see below for examples). While this looks simple, the
selection of data indices $i$ and $j$ for fusion is based on the predefined
generative model under the hood. Furthermore, the generative model controls
the degree of bias for the synthetic datasets with hyper-parameters. Thus, the
data generation process reproduces the historical bias described in Section
2.2.

A full description of the synethetic data pipeline (including generative model
specifics) can be found in Appendix A.7. Here, we summarize two key hyper-
parameters of the data generator, which we will use to control ground-truth
data bias and complexity.

  * Bias Level: The bias level controls the level of dependency between $Y$ and $C$ in the range of $[0,1]$. A Higher bias level results in larger correlation between $Y$ and $C$ in the generated dataset.

  * Fusion Function: The fusion function merges feature vectors from the two datasets. We have two variants: Concatenation which stacks features without changing the element values (see Kusner et al. (2017)), and the outer product which blends features perfectly.

####  4.1.3 Evaluating Fairness Tests and Selecting $\delta$

Having access to the ground-truth generative model, we can compute the
individual fairness score (IFS) described in Definition 1 for each generated
synthetic data sample. IFS will serve as the ground-truth label in the
following experiments on synthetic datasets.

In flagging discrimination in practice, it is necessary to set a value for the
threshold $\delta$. This threshold is usually set by regulatory standards that
depend on problem domain or statistics of manual auditing results. Indeed, the
selection of specific thresholds is a subtle question, and in practice,
decision-making may involve multiple thresholds that require further domain-
specific study Corbett-Davies and Goel (2018).

Thus, rather than choose a particular threshold $\delta$, since each candidate
test yields a continuous discrimination score, we may construct a precision-
recall curve. We then compare the tests based on average precision. This
enables us to determine which test is the most reliable across a range of
different thresholds.

####  4.1.4 Real Datasets

We also compare the proposed model with the baseline models on two real
datasets (and one more in the Appendix):

  * Adult Income Dataset Dua and Graff (2017): The protected variable is gender, a binary variable (Female/Male).

  * Bank Marketing Dataset Moro et al. (2014): The protected variable is age (binarized by thresholding at 25).

(a) Adult Dataset

(b) Bank Marketing Dataset Figure 4: Performance Comparison among Individual
Fairness Testing Methods on Real Datasets. For each testing algorithm, we show
the predicted unfairness score for both unfair model (grey box) and fair model
(blue box). Greater difference (between the boxes) shows better performance.
To conserve space, we have only plotted fAux+NG. Plots for the remaining
variants are provided in the Appendix, along with additional plots on a third
dataset.

(a) Bias

(b) Aux. Model Accuracy

(c) Target Model Gradient

Figure 5: Factors impacting fAux performance. We identify the operational
range of fAux by examining PR curves on aggregated runs. (a) Impact from data
bias. (b) Impact from the auxiliary model performance. (c) Impact from the
gradient magnitude of the target model.

###  4.2 Performance on Synthetic Data

Table 1 compares the proposed fAux methods with two state-of-the-art testing
approaches on the synthetic datasets. We make the following observations:

  1. 1.

We note the fAux variants outperform the weighted $l_{p}$ norms across all the
experiments in the Table 1 by a large margin. We attribute this to the more
complex relationships captured by our nonlinear auxiliary model.

  2. 2.

Among all fAux variants, the fAux+NG shows the most promising test
performance. However, fAux+IG becomes the more reliable test as the accuracy
of the auxiliary model decreases. The baseline fAux performs worse than the
other variants. This observation reflects our intuition that controlling the
variance of gradients is critical to maintaining good performance.

  3. 3.

FlipTest demonstrates competitive performance on some datasets but has larger
variability in performance than the proposed fAux models. We attribute this
variability to the diffuculties of training GANs. In particular, we note the
FlipTest performance is much worse on in the outer-product datasets. These are
much higher dimensional, making it more difficult for FlipTest to correctly
model the relations through generative models.

###  4.3 Performance on Real Data

In the absence of ground truth labels for discrimination, we cannot quantify
the precision of a fairness test at the level of individual datapoints.
However, we may still compare their ability to distinguish fair and unfair
models, as is done in Maity et al. (2021); Yurochkin and Sun (2021). We
construct fair models by employing adversarial regularization Adel et al.
(2019); models become spontaneously unfair in the absence of such
regularization. In Figure 4(a), we plot the distribution in the test scores
for fair and unfair models. For both datasets, we note the fAux approach shows
better performance than the existing approaches since it produces observably
higher unfairness scores to the Unfair model than those of the Fair model.
Moreover, with fAux testing, the difference of score distributions between the
Unfair and Fair model is also more distinguishable, which helps the threshold
search using statistic tools.

###  4.4 Exploring Performance

While fAux shows good experimental performance, it is important to explore
potential conditions that may impact its effectiveness. To analyze the
conditions comprehensively, we train target and auxiliary models on the same
synthetic datasets previously considered, but over a wider range of settings
(Appendix A.7.5). We aggregate these results to examine the effectiveness of
the fAux test across different datasets. In particular, we investigated the
following aspects:

  * Sensitivity: The first question is whether fAux can correctly identify models that are only slightly discriminatory. Datasets with a smaller bias parameter (Section 4.1.2) have less correlation between the features and the protected variable, thus we expect that models trained on them will be less discriminatory333We demonstrate this visually in Figure 8 of the Appendix.. Figure 5(a) shows the Precision-Recall (PR) curve of the fAux performance based on the degree of the data bias. We note that when the target model is trained on highly biased data ($bias\geq 0.75$), fAux shows better discrimination detection performance. When the target model is trained on less biased data ($bias\leq 0.5$), the reliability of fAux decreases. In the worst-case ($bias\leq 0.25$), fAux loses its functionality since the PR curve and the diagonal line (the random guess line) have multiple intersections.

  * Auxiliary model accuracy: Figure 5(b) shows the PR curve for different auxiliary model accuracies. Model accuracy varies according to the hyperparameters used, and the difficulty of the classification task. Empirically, we find that auxiliary models with better accuracy lead to better fAux performance. Since the auxiliary models aim to predict the protected variables $C$, their performance impacts how well fAux detects the variance of $C$ given features $X$.

  * Quality of target model: In practice, target model quality is out of scope for tuning fairness testing methods. However, we note that the gradient magnitude $|\nabla f_{tar}|$ impacts the fAux performance as shown in Figure 5(c). When the gradient is extremely large, fAux exhibits worse performance.

##  5 Conclusion

In this paper, we proposed a novel criterion for testing individual fairness
(Theorem 1), which we related to other well-known definitions (summarized in
Definition 1). This new criterion is based on locally testing model
derivatives, and thus avoids the challenge of generating similar pairs of
inputs. In our experiments, we have demonstrated that our new test - named
fAux - outperforms several other state-of-the-art approaches on both synthetic
and real datasets.

## References

  * Adel et al. [2019] Tameem Adel, Isabel Valera, Zoubin Ghahramani, and Adrian Weller.  One-network adversarial fairness.  Proceedings of the AAAI Conference on Artificial Intelligence, 33:2412–2420, 07 2019.
  * Agarwal et al. [2018] Aniya Agarwal, Pranay Lohia, Seema Nagar, Kuntal Dey, and Diptikalyan Saha.  Automated test generation to detect individual discrimination in ai models.  arXiv preprint arXiv:1809.03260, 2018.
  * Barocas et al. [2019] Solon Barocas, Moritz Hardt, and Arvind Narayanan.  Fairness and Machine Learning.  fairmlbook.org, 2019.
  * Bellamy et al. [2018] Rachel K. E. Bellamy, Kuntal Dey, Michael Hind, Samuel C. Hoffman, Stephanie Houde, Kalapriya Kannan, Pranay Lohia, Jacquelyn Martino, Sameep Mehta, Aleksandra Mojsilovic, Seema Nagar, Karthikeyan Natesan Ramamurthy, John Richards, Diptikalyan Saha, Prasanna Sattigeri, Moninder Singh, Kush R. Varshney, and Yunfeng Zhang.  AI Fairness 360: An extensible toolkit for detecting, understanding, and mitigating unwanted algorithmic bias, October 2018.
  * Berk et al. [2017] Richard Berk, Hoda Heidari, Shahin Jabbari, Michael Kearns, and Aaron Roth.  Fairness in criminal justice risk assessments: The state of the art, 2017\.
  * Black et al. [2019] Emily Black, Samuel Yeom, and Matt Fredrikson.  Fliptest: Fairness auditing via optimal transport.  CoRR, 2019.
  * Chiappa and Pacchiano [2021] Silvia Chiappa and Aldo Pacchiano.  Fairness with continuous optimal transport, 2021.
  * Corbett-Davies and Goel [2018] Sam Corbett-Davies and Sharad Goel.  The measure and mismeasure of fairness: A critical review of fair machine learning, 2018.
  * Dua and Graff [2017] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.
  * Dwork et al. [2012] Cynthia Dwork, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Richard Zemel.  Fairness through awareness.  In Proceedings of the 3rd innovations in theoretical computer science conference, pages 214–226, 2012.
  * Gajane and Pechenizkiy [2017] Pratik Gajane and Mykola Pechenizkiy.  On formalizing fairness in prediction with machine learning.  arXiv preprint arXiv:1710.03184, 2017.
  * Galhotra et al. [2017] Sainyam Galhotra, Yuriy Brun, and Alexandra Meliou.  Fairness testing: testing software for discrimination.  In Proceedings of the 2017 11th Joint Meeting on Foundations of Software Engineering, pages 498–510, 2017.
  * Gordaliza et al. [2019] Paula Gordaliza, Eustasio Del Barrio, Gamboa Fabrice, and Jean-Michel Loubes.  Obtaining fairness using optimal transport theory.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 2357–2365. PMLR, 09–15 Jun 2019.
  * Hodges [1958] J. L. Hodges.  The significance probability of the smirnov two-sample test.  Arkiv för Matematik, 3(5):469 – 486, 1958.
  * John et al. [2020] Philips George John, Deepak Vijaykeerthy, and Diptikalyan Saha.  Verifying individual fairness in machine learning models.  In Conference on Uncertainty in Artificial Intelligence, pages 749–758. PMLR, 2020.
  * Kusner et al. [2017] Matt J Kusner, Joshua Loftus, Chris Russell, and Ricardo Silva.  Counterfactual fairness.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 4066–4076. Curran Associates, Inc., 2017.
  * Maity et al. [2021] Subha Maity, Songkai Xue, Mikhail Yurochkin, and Yuekai Sun.  Statistical inference for individual fairness, 2021.
  * Mehrabi et al. [2019] Ninareh Mehrabi, Fred Morstatter, Nripsuta Saxena, Kristina Lerman, and Aram Galstyan.  A survey on bias and fairness in machine learning.  arXiv preprint arXiv:1908.09635, 2019.
  * Moro et al. [2014] Sérgio Moro, P. Cortez, and P. Rita.  A data-driven approach to predict the success of bank telemarketing.  Decis. Support Syst., 62:22–31, 2014.
  * Mukherjee et al. [2020] Debarghya Mukherjee, Mikhail Yurochkin, Moulinath Banerjee, and Yuekai Sun.  Two simple ways to learn individual fairness metrics from data, 2020.
  * Olson et al. [2017] Randal S. Olson, William La Cava, Patryk Orzechowski, Ryan J. Urbanowicz, and Jason H. Moore.  Pmlb: a large benchmark suite for machine learning evaluation and comparison.  BioData Mining, 10(1):36, Dec 2017.
  * Ross [2014] B. C. Ross.  Mutual information between discrete and continuous data sets.  PLoS ONE, 9, 2014.
  * Ruoss et al. [2020] Anian Ruoss, Mislav Balunović, Marc Fischer, and Martin Vechev.  Learning certified individually fair representations, 2020.
  * Srivastava et al. [2017] Akash Srivastava, Lazar Valkov, Chris Russell, Michael U. Gutmann, and Charles Sutton.  VEEGAN: Reducing Mode Collapse in GANs using Implicit Variational Learning.  arXiv e-prints, page arXiv:1705.07761, May 2017.
  * Sun et al. [2020] Wenlong Sun, Olfa Nasraoui, and Patrick Shafto.  Evolution and impact of bias in human and machine learning algorithm interaction.  Plos one, 15(8):e0235502, 2020.
  * Sundararajan et al. [2017] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In International Conference on Machine Learning, pages 3319–3328. PMLR, 2017.
  * Wachter et al. [2018] Sandra Wachter, Brent Mittelstadt, and Chris Russell.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr, 2018.
  * Xu et al. [2019] Lei Xu, Maria Skoularidou, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni.  Modeling tabular data using conditional GAN.  CoRR, abs/1907.00503, 2019.
  * Yucer et al. [2020] Seyma Yucer, Samet Akçay, Noura Al-Moubayed, and Toby P Breckon.  Exploring racial bias within face recognition via per-subject adversarially-enabled data augmentation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pages 18–19, 2020.
  * Yurochkin and Sun [2021] Mikhail Yurochkin and Yuekai Sun.  Sensei: Sensitive set invariance for enforcing individual fairness, 2021\.
  * Yurochkin et al. [2020] Mikhail Yurochkin, Amanda Bower, and Yuekai Sun.  Training individually fair ml models with sensitive subspace robustness, 2020.

##  Appendix A Appendix

###  A.1 Recap the Practical Usage of the Proposed Testing Approach

In order to introduce individual fairness testing into industry as an
important checking point/regulation, an efficient and effective testing
approach that can scale up to applications with millions of customers is
urgently needed. We note that the existing approaches such as counterfactual
fairness, while conceptionally effective, are hard to deploy in practice as
they either consume too much computation resources (that are out of budget) or
require detailed domain knowledge (that are usually too complex to summarize).
Hence, the proposed approach becomes the most practical option. It only needs
to estimate gradient (or feature importance) alignment, which happens to also
be unavoidable when testing for other properties such as adversarial
robustness testing or performing feature importance analysis, etc.

That said, we caution against over reliance on mathematical metrics for
quantifying fairness. In particular, while we show in Table 4 that faux can
flag proxy features, we emphasize this is no substitute for intuition about
the impact of the features in the dataset. We advocate faux as an effective
component of a validation pipeline, but investigators should be mindful of
other sources of bias beyond historical bias e.g. label contamination.

###  A.2 Use Cases

####  A.2.1 Predictive Policing

We consider the use case described in Kusner et al. [2017], where a
hypothetical city government wants to predict crime rates in order to assign
policing resources to different communities. The goal is that, given an
individual, a model may be used to estimate that individual’s predisposition
to violence. A developer trains a model on a dataset obtained by merging
residential information with police records of arrests. However, some
neighborhoods have higher arrest rates due to greater policing there. Because
individuals of different races may congregate in different neighborhoods, this
can lead the model to conclude that members of a particular race are more
likely to break the law.

This bias persists even if the classifier does not take race as an explicit
input, since that information may be inferred from the neighborhood. Thus,
Fairness Through Unawareness (FTU) is not applicable. fAux, however, can
identify such correlations through the gradient of the auxiliary model.
Moreover, it can recognize when these correlations impact model predictions
through gradient alignment.

####  A.2.2 Credit Cards

Consider a hypothetical credit agency that wants to automate the process of
giving credit cards. As inputs, they merge account data, transaction history,
and demographic information to form a high-dimensional, sparse dataset with
millions of entries. They train a model to predict whether the individual will
default on their payments within a certain time window. The model learns that
individuals above a certain income threshold are more likely to make their
payments. However, due to income disparity between men and women, and income
disparity between younger and older individuals, these individuals may be
unfairly treated by the model.

Fairness Through Awareness (FTA) may be deployed to solve this problem, by
going through pairs of individuals, and determining if similar treatment is
offered to similar individuals (as measured by a similarity metric). However,
this has a few shortcomings. Firstly, computing pairwise similarities is
computationally daunting, especially in higher dimensions. Secondly, the
dataset is constructed from a mixture of tabular and time series data, meaning
it incorporates different data domains that have different semantic
interpretation. Finally, the model may be justified in treating two ”similar”
individuals differently if they differ in a business-relevant sensitive
feature. Depending on the similarity metric, FTA may completely miss the
underlying issue, which in this case, is the income disparity, and simply
identify all instances where the model was sensitive.

By contrast, fAux avoids the need for designing similarity metrics, by
inferring which features are relevant for predicting sex/age. By aligning the
feature importances with the the gradients of the target model, fAux can
identify when the model is leveraging this information to make its
predictions, revealing the discrimination.

###  A.3 Discussion: Relation to Adversarial Robustness

Individual fairness testing has strong connections with adversarial robustness
Yurochkin et al. [2020]; John et al. [2020]. The adversarial robustness
criterion near point $\mathbf{x}$

$\forall\mathbf{x}^{\prime}\,s.t\,\left(\phi_{in}(\mathbf{x},\mathbf{x}^{\prime})\leq\epsilon\right)\land\left(\phi_{out}(f(\mathbf{x}),f(\mathbf{x}^{\prime}))<\delta\right)$
(11)

is nearly identical to the criterion of the FTA definition if we reformulate
the FTA criterion as

$\forall\mathbf{x}_{i}\mathbf{x}_{j}\,s.t\,\left(\phi_{in}(\mathbf{x}_{i},\mathbf{x}_{j})\leq\epsilon\right)\land\left(\phi_{out}(f_{tar}(\mathbf{x}_{i}),f_{tar}(\mathbf{x}_{j}))<\delta\right)$
(12)

with additional constraint $\epsilon\geq\delta$. Similarly, the counterfactual
fairness criterion, can be rewritten as

$\displaystyle\forall\mathbf{c}^{\prime}\,s.t\,$ (13)
$\displaystyle\left(\phi_{in}(\mathbf{c},\mathbf{c}^{\prime})\\!\leq\\!\epsilon\right)\\!\land\\!\left(\phi_{out}(P(y|\mathbf{x},\mathbf{c})\\!-\\!P_{C\leftarrow\textit{do}(\mathbf{c}^{\prime})}(y|\mathbf{x},\mathbf{c})))\\!<\\!\delta\right)$

to match the adversarial robustness definition.

The Local Fairness Test proposed in this paper also has strong connections
with adversarial robustness. However, as fAux works on the level of
differentiation, the input distance condition
$\phi_{in}(\mathbf{c},\mathbf{c}+\Delta\mathbf{c})$ becomes trivial since
$\Delta\mathbf{c}$ is negligible. Hence, fAux considers only the output
distance condition:

$\phi_{out}=\left|\nabla f_{tar}\left(\nabla f_{aux}^{\top}\nabla
f_{aux}\right)^{-1}\nabla f_{aux}^{\top}\right|_{\infty},$ (14)

where distance function over the inputs is no longer explicitly provided; the
derivative is the limit of the predictive difference between two near
identical inputs.

However, even though individual fairness testing and adversarial robustness
testing share common properties mentioned above, we note that they are
different tests that aim to reveal different weaknesses of a machine learning
model. Specifically, adversarial robustness testing focuses on detecting
vulnerability of a model with respect to anomalous (or adversarial) inputs,
whereas individual fairness testing pays attention to the sensitivity of
protected attributes on a model.

###  A.4 Proof of Theorem 3.1

The $l_{\infty}$ norm represents a max operation on the most violated
protected attribute c in $\mathbf{c}=\\{c_{1},c_{2}\cdots c_{k}\\}$. To prove
theorem 3.1, we only need to focus on this single protected attribute $c$.
Hence, if a prediction violates the Local Independence Criterion (LIC), we can
express it as:

$\displaystyle\left|\frac{\partial f_{tar}(\mathbf{x})}{\partial
c}\right|>\delta$ (15)

By expanding the left hand side of the inequality with the generative model
$f_{g}$, we note

$\displaystyle\lim_{\Delta{c}\to
0}\frac{f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c)))-f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c+\Delta
c)))}{\Delta{c}}$ (16)

which could be numerically approximated with small $\Delta c$ such that

$\displaystyle\lim_{\Delta{c}\to
0}\frac{f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c)))-f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c+\Delta
c)))}{\Delta{c}}$ (17)
$\displaystyle=\frac{f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c)))-f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c+\Delta
c)))}{\Delta c}+\xi,$

where the $\xi$ is an error term introduced by the approximation.

Therefore, Inequality 15 could be extended as

$\displaystyle\left|f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c)))-f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c+\Delta
c)))\right|$ (18) $\displaystyle>(\delta-\xi)|\Delta c|.$

If we choose $\delta\gg\xi$ and $|\Delta c|$ is positive, we have

$\displaystyle\left|f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c)))-f_{tar}(f_{g}(\mathbf{z}_{\bot},\psi(c+\Delta
c)))\right|>0,$ (19)

which violates the individual fairness definition in Definition 1.

###  A.5 Using Pseudoinverses for Model Inversion

In this section, we provide a more in-depth discussion on the derivation of
Equation 11. The Taylor expansion of Equation 10 is a local linear
approximation of the function $f_{aux}$. The left-hand side of this equation
denotes the change in the score, and is a scalar. The right hand side is a dot
product between two vectors. While we would like to invert this equation to
solve for x, Equation 11 is an indeterminate equation . Nevertheless, using
the Moore-Penrose pseudoinverse, we may still characterize the infinite
solution set for $x$.

To understand this, consider a general linear equation:

$\textbf{y}=A\textbf{x}.$ (20)

Here, $A$ is an $m$ by $n$ matrix, and we want to solve this equation for x.
When $A$ does not have an inverse, solutions may be found to (20), but they
will not be unique. Instead, there will be a family of solutions $x(z)$
parameterized by a vector $z\in\mathbf{R}^{n}$, which will be given by:

$\textbf{x}(\textbf{z})=A^{\dagger}\textbf{y}+(I-A^{\dagger}A)\textbf{z}.$
(21)

Here, $I$ is the $n$ by $n$ identity matrix, and the pseudo-inverse is given
by

$A^{\dagger}=A^{\top}(AA^{\top})^{-1}.$ (22)

Note that equation 21 describes an infinite solution set, insofar as z is
undetermined: the operation of $A$ on the right-hand side of (21) will simply
cause the second term to vanish, so that any vector $z$ may be used.

For the specific case of our local linear approximation of the auxiliary
model, $y=(c-f_{aux}(\textbf{x}_{0})$ is a scalar, $A=\nabla f_{aux}^{\top}$
is a co-vector, and we may solve for $f_{sur}(c)$ as:

$\displaystyle f_{sur}(c)=$ (23)
$\displaystyle\textbf{x}_{0}\\!+\\!(c\\!-\\!f(\textbf{x}_{0}))\frac{\nabla
f_{aux}}{\nabla f_{aux}^{\top}\nabla f_{aux}}\\!+\\!(I\\!-\\!\frac{\nabla
f_{aux}\nabla f_{aux}^{\top}}{\nabla f_{aux}^{\top}\nabla f_{aux}})\textbf{z}$

We illustrate these different terms in Figure 6. Intuitively, the equation
above describes two sources of variation: the first of which is parallel to
$\nabla f_{aux}$, and the second of which is perpendicular to it. The partial
derivative with respect to $C$ is defined by the variance in a function as $C$
is changed, and all other variables are fixed. As a relaxation, we consider
the direction that leads to maximal change in $C$, and minimal change in all
other variables. This corresponds to the minimum norm solution that satisfies
(23).

For a model trained to predict $c$, the gradient should point in the direction
of greatest change in $c$. Therefore, intuitively, the $c$-dependence of the
perpendicular part should be much smaller than the $c$-dependence of the
parallel part. We may thus disregard the perpendicular part, and are left with
the following approximation:

$\frac{\partial\mathbf{x}}{\partial\mathbf{c}}\approx\frac{\nabla
f_{aux}}{\nabla f_{aux}^{\top}\nabla f_{aux}}$ (24)

Figure 6: Inverting a linear function In one dimension, the inverse of a
linear function is uniquely determined by the slope. In higher dimensions, the
slope becomes the magnitude of the gradient. However the inverse is no longer
uniquely determined: the output is constant along directions perpendicular to
the gradient.

###  A.6 Algorithm

Algorithm 1 summarizes the proposed fairness testing approach (fAux) described
in the main paper.

Result: Flag unfair model behaviour

Input: Validation data points
$D=\\{\cdots(\mathbf{x}_{i},\mathbf{c}_{i},y_{i})\cdots\\}$, target model
$f_{tar}$, and threshold $\delta^{\prime}$

Train auxiliary model $\mathbf{c}=f_{aux}(\mathbf{x})$;

for  _each data point $(\mathbf{x},\mathbf{c},y)$ in $D$_ do

Evaluate gradients $\nabla f_{tar}$ and $\nabla f_{aux}$;

Flag unfair behaviour on inputs through $\left|\nabla f_{tar}\left(\nabla
f_{aux}^{\top}\nabla f_{aux}\right)^{-1}\nabla
f_{aux}^{\top}\right|_{\infty}\leq\delta^{\prime}$;

end for

Algorithm 1 Auxiliary Model Test (fAux)

###  A.7 Description of the Synthetic Dataset Pipeline

In Section 2.2 we hypothesized that when datasets contain multiple sources of
variation, models will preferentially exploit the simplest patterns to achieve
their performance objective. In this section, we describe a simple framework
for constructing synthetic datasets that trigger this misbehaviour.
Specifically, we present a scheme for fusing real datasets together, so that
the final dataset contains multiple patterns of variation. This has the
advantage that the final dataset inherits the same types of noisy and
nonlinear relationships that are present in real datasets. We summarize this
procedure in Algorithm 2.

The framework has three main degrees of freedom:

  1. 1.

To simulate a historical bias444In this paper we focus on historical biases,
though the proposed framework is sufficiently general to handle other
mechanisms of bias., we allow the target variable $y$ and the protected
variable $c$ to be correlated while sampling from the dataset. This is
controlled by a joint distribution $P(C,Y)$.

  2. 2.

The fusion function $F$, which combines the datasets together without
distorting the underlying patterns. We consider two extremes: one in which
data is perfectly mixed (outer product), and one in which it is trivially
separable (concatenation).

  3. 3.

The choice of datasets, which controls the difficulty of the learning task. We
measure difficulty according to the best classification score that is obtained
on each individual dataset using an architecture search. We build Conditional
Variational Autoencoders to model these datasets.

In the next few subsections, we provide more details about each of these
components

(a) Unbiased datasets

(b) Fusing synthetic bias

Figure 7: The synthetic bias pipeline. (a) Generative models for unbiased
datasets. In both datasets the observable features $X$ are distributed
according to a target $Y$ and latent variables $Z$. At present there are no
distinct populations in either dataset, and so there is no bias. (b) $X_{1}$
and $X_{2}$ are fused together to form a new dataset $X$, which contains two
distinct patterns of variation. We declare the target of the second dataset
$Y_{2}$ to be a sensitive attribute, labeled with $C$. We introduce a
historical bias using a correlation between $Y$ and $C$. Models trained on the
new dataset may learn to exploit the $C$ dependence of $X$ to predict $Y$,
which corresponds to discriminatory behaviours.

Result: Obtain samples x from a biased dataset

Input : Joint distribution $Pr(c,y)$, latent priors $Pr(\hat{z})$,
$Pr(\tilde{z})$, generative models $\Gamma_{y}$, $\Gamma_{c}$, fusion function
$F$

  1. 1.

Sample $\hat{z}\sim Pr(\hat{z})$ and $\tilde{z}\sim Pr(\tilde{z})$.

  2. 2.

Sample $(y,c)\sim Pr(c,y)$

  3. 3.

Compute $\hat{\mathbf{x}}_{i}=\Gamma_{y}(\hat{z},y)$,
$\tilde{\mathbf{x}}_{j}=\Gamma_{c}(\tilde{z},c)$

  4. 4.

Compute $\textbf{x}=F(\hat{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{j})$

Algorithm 2 Synthetic Data Generation

####  A.7.1 Bias parameter

We can control the amount of bias by specifying a joint distribution $P(C,Y)$.
In our experiments we consider binary $C$ and $Y$, so that $P(C,Y)$ may be
constructed as a 2 by 2 matrix, which involves 4 degrees of freedom. We would
instead like to reduce these to a single degree of freedom, which we term the
bias parameter.

Intuitively, minimal bias occurs when the joint distribution factorizes, so
that:

$P_{min}(C,Y)=Pr(C)Pr(y)$

The amount of bias may then be quantified by considering how much our
specified $P$ deviates from this $P_{min}$, which we do using the entropy:

$H(P)=\sum_{C}\sum_{Y}P(C,Y)\log\left(\frac{P(C,Y)}{P_{min}(C,Y)}\right).$
(25)

We may thus define a $P_{max}$ via:

$P_{max}=\mathrm{argmax}_{P}\;H(P).$ (26)

This optimization is done with the three constraints:

$\displaystyle\sum_{C}\sum_{Y}P_{max}(C,Y)$ $\displaystyle=1$ (27)
$\displaystyle\sum_{Y}P_{max}(C,Y)$ $\displaystyle=Pr(C)$ (28)
$\displaystyle\sum_{C}P_{max}(C,Y)$ $\displaystyle=Pr(Y)$ (29)

Here, the marginals are estimated from the respective datasets by examining
the relative frequencies of $C=1$ and $Y=1$. In addition we have the
positivity constraints:

$P_{max}(C,Y)\geq 0;\forall C,Y.$ (30)

We thus have an optimization problem with 4 degrees of freedom and 3
constraints, and may find $P_{max}$ using a one-dimensional line search. We
then obtain our biased $P$ by linearly interpolating between these two
extremes:

$P_{bias}=(1-bias)P_{min}+(bias)P_{max}$ (31)

We can see the effect of adjusting bias in Figure 8, demonstrating that our
proposed mechanism does indeed provoke unfair treatment in models.

(a) bias = 0

(b) bias = 0.5

(c) bias = 1

Figure 8: Controlling discrimination with tunable bias parameter. Given access
to a generative model, we can measure discrimination by evaluating the
disparity in a model $f$’s outputs over a point x and its counterfactual
$\textbf{x}^{\prime}$. The above histograms demonstrate that datasets with
higher values of the bias parameter provoke higher rates of discrimination in
models trained on them.

####  A.7.2 Fusion Functions

In this section we consider how to fuse datasets together. We consider two
datasets $\hat{D}=\\{\cdots(\hat{\mathbf{x}}_{i},\hat{y}_{i})\cdots\\}$ and
$\tilde{D}=\\{\cdots(\tilde{\mathbf{x}}_{j},\tilde{y}_{j})\cdots\\}$. The
inputs have dimensions $\hat{N}$ and $\tilde{N}$ respectively.

In fusing the datasets, we are slightly constrained: while we aim to mix the
different patterns, they must also be distinct enough that the learning task
can be performed. That is to say, all stochasticity and nonlinearity should
come from the original datasets, whose patterns should be preserved (and thus,
simultaneously learnable) under the fusion function. In our experiments, we
thus consider two extremes of mixing:

  * •

We concatenate the vectors $\hat{\mathbf{x}}_{i}$ and
$\tilde{\mathbf{x}}_{j}$:

$\textbf{x}=\textrm{concat}(\hat{\mathbf{x}}_{i},\tilde{\mathbf{x}}_{j})$ (32)

This corresponds to the setup in Kusner et al. [2017], in which the observable
features $X$ may be partitioned into descendants and non-descendants of the
protected variable $C$555That is to say, given a causal graph, certain
features of $X$ may be children of the protected variable $C$. A simple linear
projection is sufficient to separate these sources of variation:

$\displaystyle\hat{\mathbf{x}}_{i}=\hat{W}\textbf{x}$
$\displaystyle;\;\;\;\hat{W}=\mathrm{concat}(I_{\hat{N}},0)$ (33)
$\displaystyle\tilde{\mathbf{x}}_{j}=\tilde{W}\textbf{x}$
$\displaystyle;\;\;\;\tilde{W}=\mathrm{concat}(0,I_{\tilde{N}})$ (34)

  * •

As another extreme, we consider the case of perfect mixing by taking the outer
product:

$\textbf{x}=\mathrm{vec}(\hat{\mathbf{x}}_{i}\otimes\tilde{\mathbf{x}}_{j}).$
(35)

If $\hat{\mathbf{x}}_{i}$ and $\tilde{\mathbf{x}}_{j}$ each sum to 1, then we
can completely isolate them from x using a different linear projection:

$\displaystyle\hat{\mathbf{x}}_{i}=\hat{W}\textbf{x}$
$\displaystyle;\;\;\;\hat{W}=\frac{1}{\tilde{N}}(I_{\hat{N}}\otimes
1_{\tilde{N}})$ (36) $\displaystyle\tilde{\mathbf{x}}_{j}=\tilde{W}\textbf{x}$
$\displaystyle;\;\;\;\tilde{W}=\frac{1}{\hat{N}}(1_{\hat{N}}\otimes
I_{\tilde{N}})$ (37)

Table 2: Original datasets used in the synthetic pipeline: information on the
datasets from PMLB, along with the hyperparameters of the CVAEs used to model
them. A lower model score corresponds to a more ”difficult” dataset to model.

dataset | # features | # instances |  | # gmm
---
components
hidden size | depth |  | model
---
score
backache | 6 | 180 | 3 | 186 | 6 | 0.545
credit | 7 | 1000 | 7 | 874 | 3 | 0.657
australian | 6 | 690 | 5 | 977 | 5 | 0.721
wdbc | 30 | 569 | 4 | 367 | 1 | 0.815
magic | 10 | 19020 | 2 | 624 | 2 | 0.987

####  A.7.3 Generative Models

For ease of generating counterfactuals, we use Conditional Variational
Autoencoders (CVAEs) to generate samples $X$ from a specified target $Y$. We
discard all categorical features, and retain only those which are real values.
Because the datasets used are tabular, we use a pre-processing strategy
similar to Xu et al. [2019], modelling each feature with a Gaussian Mixture
Model (GMM). This models the univariate feature distribution as a sum of
distinct gaussian modes, with each mode assigning the feature a certain
probability. We may thus transform each feature into a vector of mode
probabilities, and we concatenate these probability vectors to form our inputs
x. Since each probability vector sums to 1, this has the added advantage that
the sum over x will be a fixed integer, namely, the number of features, making
it compatible with the outer product transform described in Section A.7.2.

In addition to the depth and width of the encoder/decoder, we tune the number
of components of the GMM when training our generative models. To evaluate the
final model, we sample a collection of $(x,y)$ pairs, and use them to train a
simple classifer to predict $y$ from $x$. We tune our parameters so as to
maximize the accuracy of this classifier. Intuitively, this results in a
generative model that best preserves the original discriminative relationship
inside a given dataset. Model hyperparameters are shown in Table 2.

####  A.7.4 Datasets Used

In determining our datasets, we aim to find a combination of stochastic and
nonlinear patterns. To quantify the complexity of these patterns, we measure
the maximum accuracy reached by a simple classifier after tuning. Intuitively,
if a model is unable to obtain a high score, that dataset constitutes a more
difficult challenge. We perform this analysis on all of the binary
classification datasets in PMLB Olson et al. [2017], and select those listed
in Table 2, since they represent a spread of complexity. Note that these
datasets are distributed under an MIT License.

Note that, in selecting and training models on these datasets, we retain only
those features which are continuously valued, in order to match our generative
model setup.

####  A.7.5 Specific Settings used in Section 4.4

To generate the plots in Section 4.4 of the main paper, we generate synthetic
datasets according to the following parameter grid:

  * •

We select the bias parameter (Section A.7.1) from the following choices:
$[0,0.25,0.5,0.75,1.0]$.

  * •

We select the fusion function (Section A.7.2 from the set
$[\textrm{concat},\textrm{outer}]$.

  * •

We select all pairs of datasets from Table 2.

We train fully connected neural networks for our auxiliary and target models,
using randomly sampled widths (8-128 units) and depths (1-4 layers).

####  A.7.6 Constructing Ground-Truth Labels and Models

The main utility of using generative models is that we may simultaneously
sample both an instance x and its counterfactual $\textbf{x}^{\prime}$.
Accordingly, by passing x and $\textbf{x}^{\prime}$ through the target model
$f_{tar}$, we may use the $l_{1}$-distance between the model outputs as a
measure of discrimination, which we refer to as the Individual Fairness Score
(IFS):

$\textrm{IFS}=|f_{tar}(\textbf{x}_{i})-f_{tar}(\textbf{x}_{j})|_{1}$ (38)

Figure 8 demonstrates how this CFF score varies as the parameters of our
pipeline are adjusted.

To obtain binary labels, we define a threshold on the CFF score, which we
compute using the variance of the CFF score in the absence of any bias. Above
this threshold, we say that the model is discriminating against the individual
x on the basis of $c$.

Furthermore, knowing the form of the fusion functions as described in Section
A.7.2, we may constrain the input layers of an MLP to separate the different
patterns contained inside the synthetic data. Models which project out the $C$
dependence become perfectly fair target models. Models which project out the
$Y$ dependence become perfectly disentangled auxiliary models.

###  A.8 Sources for the Experiments on Real Datasets

The Adult and Bank datasets were obtained from IBM Fairness 360 Bellamy et al.
[2018], and the SSL dataset was obtained from the Chicago Data Portal, which
may be accessed via this hyperlink.

###  A.9 Full Experiment Results for Real Datasets

Due to the limited space in the main paper, we were unable to show the
complete experiment results. Here, we show the full version of the performance
comparison on real datasets. Figure 10 shows the complete plots on the three
datasets.

(a) Adult Dataset

(b) Bank Dataset

(c) Chicago SSL Dataset

Figure 10: Performance Comparison among Individual Fairness Testing Methods on
Real Datasets. For each testing algorithm, we show the predicted unfairness
score for both unfair model (grey box) and fair model (blue box). Greater
difference (between the boxes) shows better performance.  Table 3: Reliability
of Transparency Reports Here, we compare the candidate testing methods’
ability for flagging surrogate features of protected attributes through
ranking metric NDCG.

| Adult | Chicago SSL | Bank
---|---|---|---
Protected attribute | Sex (Binary) | Race (Binary) | Age (Binarized)
FlipTest | 0.810 $\pm$ 0.150 | 0.840 $\pm$ 0.136 | 0.752 $\pm$ 0.108
fAux | 0.802 $\pm$ 0.058 | 0.796 $\pm$ 0.098 | 0.798 $\pm$ 0.047
fAux+NG | 0.807 $\pm$ 0.105 | 0.884 $\pm$ 0.119 | 0.789 $\pm$ 0.058
fAux+IG | 0.896 $\pm$ 0.037 | 0.958 $\pm$ 0.006 | 0.814 $\pm$ 0.058

Table 4: Resource Consumption Comparison between FlipTest and fAux. The fAux
variants share the same basic architecture and so we use fAux to represent all
variants.

| Resource | FlipTest | fAux |
---|---|---|---|---
Adult | # Parameters | 220 $\pm$ 187 K | 13 $\pm$ 8 K |
Training time | 1012 $\pm$ 532 s | 12 $\pm$ 5 s |
Chicago SSL | # Parameters | 154 $\pm$ 150 K | 7 $\pm$ 6 K |
Training time | 1086 $\pm$ 580 s | 90 $\pm$ 23 s |
Bank | # Parameters | 206 $\pm$ 159 K | 13 $\pm$ 10 K |
Training time | 444 $\pm$ 204 s | 15 $\pm$ 4 s |

###  A.10 Explaining Discrimination With Transparency Reports

As the real datasets do not have corresponding generative models that produce
the data, we cannot compute IFS as the ground-truth label and so direct
performance comparison on individual fairness testing is infeasible. In the
main body of the paper, we opted to compare the scores output by the fairness
test on fair and unfair models. However, another way to evaluate a fairness
test is to examine the explanations it provides for flagging discrimination.
In this section, we examine how reliably fAux can generate such explanations.

In Black et al. [2019] the authors construct transparency reports that rank
features based on their contribution to an unfair decision. To generate scores
for the different tests, we collect all $N_{0}$ instances in the dataset for
which $c=0$, to define a set $S_{0}$. For Fliptest, we use the generator $G$
to compute the following vector:

$\frac{1}{N_{0}}\sum_{\textbf{x}\in S_{0}}(\textbf{x}-G(\textbf{x}))$ (39)

The Fliptest feature scores are then given by the absolute value of this
vector. While fAux does not use an explicit Generator, it does define a
transformation of features through the auxiliary model $f_{aux}$ (23). We can
thus compute transparency reports from the following vector:

$\frac{1}{N_{0}}\sum_{\textbf{x}\in S_{0}}\nabla f_{aux}(\textbf{x})$ (40)

Once again, the feature score is given by the absolute values of this vector.

To quantitatively compare these scores, we use a ranking metric, the
Normalized Discounted Cumulative Gain (NDCG). Our ground truth ranking is
based on an estimate of the mutual information between each feature and the
protected variable $C$. We use a nonparametric entropy-estimator based on Ross
[2014], which is available through scikit-learn. Table 3 then compares the
NDCG scores for the different test rankings. Note that, for categorical
features that are one-hot encoded, we compute an aggregate score for the
feature by taking the mean of the scores of the one-hot components.

We note fAux+IG shows remarkably better performance than FlipTest in terms of
explaining the reason for unfair treatment. In addition, we may compare the
computational efficiency of both approaches. Table 4 shows the computation
resources used for the previous experiment on the real datasets. We note that
the proposed fAux framework is remarkably efficient. It uses $\sim$5% of the
parameters and 10% of the training time of FlipTest on both datasets. In
combination with the results from Table 1 and 3, this highlights the
significant advantages of fAux in terms of both effectiveness and efficiency.

(a) Average KS statistic

(b) Average MSE

Figure 11: Correlation between GAN evaluation metrics and fairness test
performance. Given the high sensitivity of GANs to their hyperparameters, it
is important to find evaluation criteria to tune them. The authors of Fliptest
suggest two, average KS and average MSE, and we compare how these metrics are
related to the average precision of the fliptest. The plots show the effect of
using different values of the transport weight $\lambda$. We find that, while
average KS does correlate well with average precision, average MSE does not.

###  A.11 Hyperparameters for the Fairness Tests

In this section, we review the hyperparameters used for the different fairness
tests for both the synthetic and real datasets. All experiments were
distributed on our internal cluster which contain NVIDIA DGX-1 with 8X NVIDIA
Tesla® V100 32 GB/GPU. The cluster runs Red Hat Enterprise Linux Server
release 7.9, and our experiments employ pytorch version 1.9.0. Unless
otherwise stated, all stochastic algorithms (ex minibatch gradient descent,
parameter initialization, etc) use a random seed of 0.

####  A.11.1 Hyperparameters for Constructing Fair Models

For the real datasets, without knowledge of the fusion-function, we employ
Adel et al. [2019] to train fair models using adversarial regularization.
There were two parameters to tune:

  * •

The regularization strength $\alpha$. This is the initial relative weight
given to the regularizer over the loss function, and this weight exponentially
decays between epochs.

  * •

Adversarial training can sometimes be unstable, and so for every iteration
that we trained the target model, we train the adversarial model for $n_{adv}$
epochs.

Adel et al. [2019] was quick to use on our datasets, and so we tuned these
values manually. In particular, we tuned the regularization strength on our
synthetic datasets since they offered a more accurate assesment of the
tradeoff between accuracy and fairness. We found that an $\alpha$ of 100
yielded good results. We tuned $n_{adv}$ in a similar fashion, and obtained a
value of 3. We used these values when training models on all of the real
datasets.

####  A.11.2 Hyperparameters for fAux

We use MLPs for all auxiliary models, varying the depth and width across
different experimental runs. We use a fixed batch size of 64 across all
experiments, using the ADAM optimizer and a learning rate of 0.001 with early
stopping.

For the real datasets, we additionally employ Adel et al. [2019], so that the
resulting auxiliary model attends to those features that are most strongly
correlated with the protected attribute. We employ the same settings that we
use when constructing the fair models.

####  A.11.3 Hyperparameters for Unfair Map

The Unfair Map Maity et al. [2021] uses a gradient-flow attack restricted to a
sensitive subspace in order to test individual fairness. Tensorflow code for
this paper was publically released through OpenReview666Link here., and our
own (pytorch) implementation is based off this. There are three
hyperparameters to consider:

  * •

The regularization strength, which is responsible for restricting the attacks
to the sensitive subspace.

  * •

The learning rate, which is used to make the gradient updates in the attack.

  * •

The number of steps used in the attack.

In our experiments, we use the same values for these parameters that are
provided in the reference implementation. In addition to these parameters,
this test requires a fair metric, which is constructed from a logistic
regression model. To determine the sensitivity of this approach to the fair
metric used, we train the logistic regression model using scikit-learn using
different random seeds. For the synthetic dataset experiments, we use ten
random seeds, taken from [10, 20, … 100].

####  A.11.4 Hyperparameters for GANs

We use MLPs of varying depth and width for our Generators and Discriminators.
Due to the sensitive nature of GAN training, we tune the batch size, depth,
and number of hidden dimensions, in addition to the following hyperparameters:

  1. 1.

$n_{critic}$, which controls the relative number of training steps between the
Dicriminator and the Generator

  2. 2.

$\lambda$, the weight of the transport cost in the Generator loss function
specific to Black et al. [2019].

Collectively, these parameters account for the architecture, loss function,
and optimization of the GANs, to see which components of the training pipeline
are most influential. As per the original Fliptest paper, we tune the GAN
models using the following metrics:

  1. 1.

The Kolmogorov–Smirnov (KS) two-sample test Hodges [1958] on the marginal
distributions for each feature between the real data x and the generated data
$G(\textbf{x})$. Better GAN models will have a smaller KS-statistic (averaged
across features), as the real and generated distributions will be similar.

  2. 2.

The Mean-Squared Error (MSE) of a linear regression model trained to predict
each observable feature from the remaining features. Better GAN models will
have lower MSE values (averaged across features), as they will have captured
correlation between the features well.

In Figure 11, we examine the correlation of these GAN metrics on the
Synthetic-5 experiment, shown in Table 1. We see that, while a smaller KS
statistic does lead to a higher average precision, the MSE metric is at odds
precision score. In this experiment, the most relevant hyperparameter was the
weight of the transport cost $\lambda$: generally, we find the best
performance is given by very small values (less than $10^{-4}$). We use this
value when testing fair/unfair models on the real datasets.

###  A.12 Scope of Work

In this section, we collect the assumptions and limitations of our approach:

  * •

To assess the impact of modifying the protected variable $C$, we must train an
auxiliary model to predict $C$. This requires the protected variable to be an
explicit variable in the training data. However, because of regulations, it is
not always possible to collect this. In this case, FTA may be a better choice.

  * •

We base our test on definitions of independence. For certain applications,
however, this may be the wrong criteria to use. For instance, in the presence
of label bias, it might make sense to use the protected variable as an input
for biased mitigation. In this case, our technique is not applicable.

  * •

Our test leverages correlations between the input data $X$ and the protected
variable $C$ in order to identify discrimination. Path-dependent
discrimination on the other hand, is out of the scope of this paper.

  * •

The features flagged by fAux as surrogates may have a legitimate business
reason for being used. For example, with the adult dataset, on average women
work fewer hours per week, but it is fair to use this feature to predict
income. fAux will flag features for being correlated with a protected
variable, so that a human validator is needed to interpret the final feature
rankings.

  * •

Our analysis employs a single first order Taylor expansion around an input x.
We thus only consider discrimination that may be revealed by perturbations in
a small neighbourhood around the input. If the target model does not change
appreciably over this neighbourhood (say, for points very far from the
decision boundary), it may be necessary to consider successive first order
expansions. This will be a topic for future investigation.


# Octupole deformation instability in atomic nuclei

Bui Minh Loc<EMAIL_ADDRESS>Center for Exotic Nuclear Studies,
Institute for Basic Science (IBS), Daejeon 34126, Korea.  Nguyen Le Anh
<EMAIL_ADDRESS>Department of Physics, Ho Chi Minh City University of
Education, 280 An Duong Vuong, District 5, Ho Chi Minh City, Vietnam.
Department of Theoretical Physics, Faculty of Physics and Engineering Physics,
University of Science, Ho Chi Minh City, Vietnam  Vietnam National University,
Ho Chi Minh City, Vietnam.  Panagiota Papakonstantinou<EMAIL_ADDRESS>Rare Isotope Science Project, Institute for Basic Science, Daejeon 34047,
Korea  Naftali Auerbach<EMAIL_ADDRESS>School of Physics and
Astronomy, Tel Aviv University, Tel Aviv 69978, Israel

###### Abstract

Recent high-energy heavy-ion collision experiments have revealed that some
atomic nuclei exhibit unusual softness and significant shape fluctuations. In
this Letter, we use the fully self-consistent mean-field theory to identify
all even-even nuclei that are unstable or soft against octupole deformation.
All exceptional cases of enhanced octupole transition strengths in stable
even-even nuclei throughout the nuclide chart are resolved. These results
represent a significant advance in our understanding of the underlying
mechanisms of nuclear octupole deformation and have implications for further
experimental and theoretical studies.

Octupole excitation, Random Phase Approximation, Skyrme force, pairing,
single-particle level

Introduction. The impact of nuclear deformation on the elliptic flow measured
in relativistic heavy-ion collisions was emphasized and formalized in Ref.
[1]. Subsequently, analyses Ref. [2, 3] of the published data from the STAR
Collaboration [4] provide direct evidence of strong octupole correlation in
the ground-state 96Zr in heavy-ion collisions at very high energy. The
question is whether the deformation values in low-energy literature are
consistent with these new data at high energy.

The existence of low-lying octupole ($3^{-}$) vibrations in atomic nuclei, at
excitation energies of only a few MeV, was recognized early on [5]. They are
now among the best-established collective states in nuclear physics and have
been observed in almost all even-even stable nuclei throughout the nuclide
chart. A review of octupole collectivity can be found in Ref. [6], while Ref.
[7] offers a review of nuclear reflection asymmetry in general. Experimental
information on the first $3^{-}$ states of even-even nuclides is compiled in
Refs. [8, 9]. Maximal strength values for the transition from the ground state
to the first $3^{-}$ state are observed at neutron numbers $N=34,56,88$, and
$134$, and atomic (proton) numbers $Z=30,40,62$, and $88$ [10]. The term
“octupole-magic” number has been used to refer to them [7].

As already indicated by the theoretical interpretation in Ref. [5], low-lying
$3^{-}$ collectivity in atomic nuclei is driven by the presence of pairs of
single-particle states with opposite parity, whose angular momenta differ by
3, in the vicinity and on both sides of the Fermi energy. The presence of such
pairs with energy differences much lower than the characteristic shell energy
of 1$\hbar\omega$ becomes possible thanks to spin-orbit splitting. Enhanced
collectivity at lower energies can be expected when the energy difference
between such states becomes especially small as a result of the interplay of
the spin-orbit coupling strength and the other nuclear interaction terms. Such
is the case, for example, of the $2d_{5/2}-1h_{11/2}$ neutron particle-hole
pair in 96Zr. Octupole-magic numbers can be explained in such a way.

Global theoretical analyses of these excitations in even-even nuclei were
given in Refs. [11, 12]. It has been reported that to have a reliable way of
understanding the character and making predictions about octupole states, one
has to go beyond the mean-field approach [13, 12, 14, 15, 16, 17]. However,
several questions remain open even then. Over time, the exceptional
characteristic of some octupole transitions has been revealed experimentally
without a satisfactory explanation. The nucleus 96Zr with $N=56$ is quite
irregular and the various theoretical attempts to study its octupole
excitation are inconsistent with the experiment and with one another, as was
recently discussed in Ref. [18]. Very recently, Ref. [19] reported the
enhancement in the light atomic nucleus 72Se with $Z=34$.

The irregular behavior of the energy of the low-lying $3^{-}$ in 96Zr was
first found and discussed in Ref. [20]. Examining the single-particle
spectrum, one observed that the $2d_{5/2}$-neutron state is fully occupied and
there is close opposite parity $1h_{11/2}$ state unoccupied. The excitation of
each neutron from $2d_{5/2}$ to $1h_{11/2}$, therefore, may form the low-lying
$3^{-}$ state. In contrast, 96Ru has only $2$ neutrons at the $2d_{5/2}$ and
therefore fewer particle-hole states to form the low-lying $3^{-}$ state,
resulting in less collectivity compared to 96Zr. The presence of this
particle-hole pair enhances collectivity and pushes the $3^{-}$ state to lower
energy.

As a collective vibration, the octupole vibration can be studied within the
self-consistent mean-field approach, specifically, the self-consistent random-
phase approximation (RPA), which is derived as the linearized limit of the
time-dependent Hartree-Fock method [21]. This was the approach followed early
on in Ref. [22], where the importance of self-consistency was emphasized.
Self-consistent RPA is a unique tool not only for accurately describing
collective vibrations and connecting their properties to an underlying energy
density functional or effective interaction, but also for diagnosing
instabilities and broken symmetries: For example, if the self-consistent RPA
equations are solved for nuclei assuming a spherical ground state, the
presence of imaginary solutions in, e.g., the quadrupole or octupole channel
will indicate that the ground state is, in fact, quadrupole- or octupole-
deformed, respectively [23, 24, 25]. The spherical ground state is said to
“collapse.” Similarly, softness towards shape fluctuations can be expected
when there are collective vibrations at very low energy.

In this Letter, the instability or softness of certain even-even nuclei
against the variations in the octupole collective variable is indicated
throughout the nuclide chart using the self-consistent mean-field framework.
The enhanced collectivity of octupole excitation around octupole-magic nuclei
is explained and reproduced.

Method. The self-consistent mean-field theory is a powerful tool for studying
the properties of atomic nuclei, such as their shapes, energies, and
excitation spectra [26, 27, 21]. The theoretical foundation of the mean field
is provided by the Hartree-Fock (HF) theory using the Skyrme interaction,
which is one of the most commonly used types of effective interaction. RPA is
used within the framework of the self-consistent mean-field theory to describe
collective motion in atomic nuclei and especially harmonic vibrations.
Therefore, this framework is highly suitable for the study of the first
$3^{-}$ octupole state which is low-lying and strongly collective. Self-
consistency is ensured when the RPA particle-hole interaction is derived from
the same effective interaction used for obtaining the HF ground state.

Computational codes exist for solving the HF and RPA equations and they may
include a pairing interaction for describing open-shell nuclei. The publically
available computational code described in Ref. [28], allows both theorists and
experimentalists to perform computations on a wide range of nuclear
excitations using fully self-consistent Skyrme HF-RPA, assuming spherical
symmetry. All relevant terms of the residual particle-hole interaction are
incorporated into the calculation, including the Coulomb and spin-orbit terms,
the latter being especially relevant for the octupole vibration, as already
discussed. With pairing correlation included, HF is extended to HF-BCS and RPA
to Quasiparticle RPA (QRPA) accordingly [29] with the self-consistency
remaining. We use the above approaches in the present study. The pairing
correlation is included in our calculation. The selected Skyrme forces in the
present work are SkM* [30], SLy4, and SLy5 [31] and, for comparison, the SIII
[32], which was one of the first parameterizations.

The reduced transition strength from the correlated RPA ground state
$|\tilde{0}\rangle$ to the first octupole state $3^{-}_{1}$ is

$B(E3)=|\langle 3^{-}_{1}||\hat{F}_{3}||\tilde{0}\rangle|^{2},$ (1)

where the isoscalar octupole operator is

$\hat{F}_{3,0}(\bm{r})=\sum_{i}^{A}r_{i}^{3}Y_{3,0}(\hat{r}_{i}).$ (2)

All $B(E3)$ are presented in Weisskopf units (W.u.) with
$1~{}\text{W.u.}=0.0594~{}A^{2}e^{2}\text{fm}^{6}$.

The degree of stability or softness of the atomic nucleus against variations
is expressed via the polarizability $\alpha$, which is obtained from the
inverse energy-weighted sum rule of the response function [27, 20]. In the
RPA, it is calculated from the $m_{-1}$ moment

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle 2m_{-1}/A,$ (3)
$\displaystyle m_{-1}$ $\displaystyle=$ $\displaystyle\sum_{n}|\langle
3^{-}_{n}||\hat{F}_{3}||\tilde{0}\rangle|^{2}E_{n}^{-1}.$ (4)

The larger $\alpha$ (W.u./MeV) is, the less stable the system is against
deformations. For example, it is around $1$ W.u./MeV for 208Pb, which is not
considered especially soft.

We have performed calculations for all stable even-even nuclei and we will
discuss a selection of them below. We will compare our results with
experimental values from the NuDat database (https://www.nndc.bnl.gov/nudat/)
for the energies and from Table VII of Ref. [9] for $B(E3)$.

Results and Discussions.

Table 1: The HFBCS-QRPA results for 96Zr are irregular in contrast to the results of 96Ru. The values of $E(3_{1}^{-})$, $B(E3)$, and $\alpha$ are in MeV, W.u., and W.u./MeV, respectively. Force |  96Zr |  96Ru
---|---|---
$E$($3^{-}_{1}$) | $B(E3)$ | $\alpha$ |  $E$($3^{-}_{1}$) | $B(E3)$ | $\alpha$
SIII | collapse | $2.437$ | $22.4$ | $1.2$
SkM* | collapse | $3.065$ | $19.7$ | $0.9$
SLy4 | $0.961$ | $96.4$ | $12.6$ | $3.875$ | $23.1$ | $0.8$
SLy5 | $1.698$ | $51.9$ | $4.0$ | $3.737$ | $22.3$ | $0.8$
Exp. | $1.897$ | $53$ |  | $2.650$ | - |

First, we discuss the difference between 96Zr and 96Ru with respect to the
octupole excitation. The recent STAR measurement showed a significant
difference between 96Zr$+$96Zr and 96Ru$+$96Ru collisions that were explained
using a transport simulation as an octupole deformation of 96Zr in its
intrinsic reference frame and large quadrupole deformation of 96Ru [2]. Our
results for 96Zr and 96Ru with different Skyrme forces are presented in Table
1. The results of 96Ru are not dramatically different from the experiment and
from each other and can be considered well-understood. The 96Ru nucleus has a
small octupole polarizability, $\alpha\approx 1$ as we find, making it hard to
be octupole deformed in, for example, isobaric heavy-ion collisions. In
contrast, the results for 96Zr are irregular. In Ref. [20], the energy of the
$3^{-}_{1}$ state of 96Zr was found to depend very strongly on the single-
particle spectrum that is obtained self-consistently in the framework. When
the gap between occupied and unoccupied single-particle energies of opposite
parity is abnormally small, the excitation energy of the $3^{-}_{1}$ state
gets dramatically low. For example, the gap between $2d_{5/2}$ and $1h_{11/2}$
given by the calculation with SIII force is so small that the $3^{-}_{1}$
state “collapses.”

When the energy of the first $3^{-}$ is imaginary in PRA, we indicate in Table
1 that there is “collapse.” As we assumed a spherically symmetric ground
state, the interpretation is that the respective effective interactions
predict that the ground state is octupole-deformed. In reality, 96Zr has a
$3^{-}_{1}$ state at $1.9$ MeV but our results with the different Skyrme
interactions suggest that this nucleus is soft against octupole deformation
and a small change in the calculation input can even lead to octupole
instability.

The calculation with SLy5 force reproduces both the excitation energy and
transition strength correctly. The octupole collectivity of 96Zr is $52$ W.u.,
and the total octupole polarizability $\alpha$ is $4$ W.u./MeV. The large
enhancement of octupole collectivity of 96Zr observed in reactions is
understood and reproduced by our calculation. The condition for softness is
not unique to 96Zr as there are other nuclei with similar characteristics. Our
HFBCS-QRPA calculation identifies such nuclei in the nuclide chart and well
reproduces the experimental data.

In the case of 96Zr, the pair with strong octupole coupling is
$2d_{5/2}-1h_{11/2}$. From the single-particle spectrum, the other pairs are
$2p_{3/2}-1g_{9/2}$, $2f_{7/2}-1i_{13/2}$, and $2g_{9/2}-1j_{15/2}$. They are
corresponding to the octupole-magic nuclei which are around $34,56,88$, and
$134$ (for the neutron). In addition, the pair $2s_{1/2}-1f_{7/2}$ is also
valid. Therefore, we suggest that the lightest octupole-magic number is
actually $16$. Table 2 shows the results of selected nuclei of interest.

Table 2: The results for selected nuclei with the number of neutrons or protons are around $16,34,56,88$, and $134$ (neutron only). Enhanced octupole transitions are found from light to heavy octupole-magic nuclei, while “collapse” may be obtained for heavy nuclei. Nuclei | Force |  $E$($3^{-}_{1}$) | $B(E3)$ | $\alpha$ |  $E$($3^{-}_{1}$) | $B(E3)$ | $\alpha$
---|---|---|---|---|---|---|---
|  | RPA | QRPA
${}^{32}_{16}$S16 | SkM* | $5.692$ | $15.6$ | $1.0$ | $5.686$ | $15.7$ | $1.0$
| SLy4 | $6.255$ | $18.3$ | $1.0$ | $6.248$ | $18.4$ | $1.0$
| SLy5 | $6.382$ | $18.8$ | $1.0$ | $6.147$ | $20.4$ | $1.0$
| Exp. |  |  |  | $5.006$ | $30$ |
${}^{64}_{30}$Zn34 | SkM* | $1.959$ | $5.8$ | $1.4$ | $3.315$ | $18.5$ | $1.4$
| SLy4 | $3.381$ | $13.3$ | $1.2$ | $4.243$ | $24.7$ | $1.2$
| SLy5 | $3.431$ | $13.6$ | $1.2$ | $4.265$ | $24.8$ | $1.2$
| Exp. |  |  |  | $2.999$ | $20$ |
${}^{72}_{34}$Se38 | SkM* | $1.068$ | $54.8$ | $6.5$ | $1.135$ | $75.3$ | $8.4$
| SLy4 | $2.069$ | $44.8$ | $2.9$ | $2.406$ | $49.3$ | $2.8$
| SLy5 | $1.931$ | $44.2$ | $3.1$ | $2.412$ | $47.1$ | $2.7$
| Exp. |  |  |  | $2.406$ | $32$ |
${}^{90}_{34}$Se56 | SkM* | $1.542$ | $38.6$ | $4.7$ | $1.882$ | $33.6$ | $3.5$
| SLy4 | $1.947$ | $37.1$ | $3.6$ | $2.749$ | $26.8$ | $2.1$
| SLy5 | $1.879$ | $37.7$ | $3.8$ | $2.574$ | $29.6$ | $2.4$
| Exp. |  |  |  | $-$ | $-$ |
${}^{98}_{40}$Zr58 | SkM* | collapse | collapse
| SLy4 | collapse | $0.435$ | $215.6$ | $64.4$
| SLy5 | collapse | $1.332$ | $65.8$ | $6.6$
| Exp. |  |  |  | $1.806$ | $-$ |
${}^{146}_{\phantom{1}56}$Ba90 | SkM* | collapse | collapse
| SLy4 | collapse | $1.695$ | $39.7$ | $2.9$
| SLy5 | collapse | $1.541$ | $45.0$ | $3.4$
| Exp. |  |  |  | $0.821$ | $-$ |
${}^{226}_{\phantom{1}88}$Ra138 | SkM* | collapse | collapse
| SLy4 | collapse | $1.158$ | $59.5$ | $10.3$
| SLy5 | collapse | $1.325$ | $51.4$ | $17.7$
| Exp. |  |  |  | $0.322$ | $54$ |

The nucleus 32S is a double-octupole magic nucleus ($N=Z=16$). The fully
occupied state $2s_{1/2}$ is strongly coupled with the unoccupied state
$1f_{7/2}$. The value of $E(3^{-}_{1})$ is not negative for such light nuclei,
but there is the enhancement of octupole collectivity. Indeed, the
experimental value of $B(E3)$ of 32S is $30$ W.u. making it the strongest
known $B(E3)/A$ value. Reminding that it is $34$ W.u. for 208Pb. Our result is
$20$ W.u with SLy5 force.

Reference [19] reported a recent experiment that showed a notably higher
octupole strength of approximately $32$ W.u. for 72Se, in contrast to 74,76Kr
which had only exhibited the strength of $15$ W.u.. The origin of enhanced
octupole collectivity in 72Se is well-explained in our discussion. The
calculation for 72Se with the SLy5 force gave $E_{x}(3^{-}_{1})=2.41$ MeV
close to the experimental value of $2.40$ MeV, and $B(E3)=47$ W.u.. The
neutron-rich nucleus 90Se that can be produced in a rare-isotope facility
nowadays is also presented in Table 2. In this region, the problem describing
the $B(E3)$ in 64Zn mentioned in Ref. [11] is also resolved as $N=34$ which is
an octupole magic number. It is $24$ W.u. in our result which is not far from
the experimental value of $20$ W.u. (see Table 2).

Of course, the nucleus 98Zr is similar to 96Zr which exhibits a distinct
behavior compared to 90Zr or 96Ru mentioned earlier. Ba and Ra isotopes have
been the subject of numerous studies in the literature (see [6] for a review).
The results for 146Ba and 226Ra are presented for examples. The details of
research on lanthanides and actinides within the QRPA framework will be
addressed in future work.

In Table 2, we also show RPA results without pairing. Pairing plays a role as
it keeps the nucleus less deformed and is essential to reproduce experimental
data, but the presence of pairing does not alter the basic physics we discuss.
It should be noted that previous research on the $3^{-}_{1}$ state using the
Green’s function RPA framework was conducted in Refs.[33, 20], but without
full self-consistency.

The results for octupole-magic nuclei are extremely sensitive to the choice of
Skyrme force. In Tables 1 and 2, the calculations with SLy4 force for 96,98Zr
overestimate the transition strengths $B(E3)$ while the results with SLy5
force are found reasonable, in both energy and strength. The difference
between SLy5 and SLy4 is from the terms which depend on the spin-orbit
densities [31]. It is well-known that the spin-orbit interaction plays a key
point in the single-particle spectrum which, as we saw, largely determines the
octupole magic numbers. A small change in the spin-orbit component makes a
significant change in the result.

Finally, we remark that in our calculations with the SLy5 force, the
“collapse” occurs in 152,154Sm, 154-160Gd, 156-170Dy, 162-170Er. These atomic
nuclei are also octupole-magic nuclei with $Z$ and $N$ around $56$ and $88$,
respectively. The research to find reliable forces or create a new one for the
study of the octupole state and instability of octupole deformation across the
nuclide chart will be discussed in the future.

Conclusions and outlooks. The fully-self consistent mean-field framework
reveals every soft-octupole deformed atomic nucleus. The octupole instability
mode may lead to an octupole shape in the intrinsic frame. Such nuclei can
enhance some unusual phenomena, such as violation of fundamental symmetries
(invariance with respect to coordinate inversion and time reversal) [34]. An
interesting question is whether the evolution of shell structure towards the
drip line leads to new octupole magic numbers and reflection asymmetry in
exotic nuclei. We plan to investigate the question of how well the QRPA
framework would perform in a global context and the interaction of the
octupole with the quadrupole degree of freedom in the future.

Acknowledgment. B. M. L. was supported by the Institute for Basic Science
(IBS-R031-D1). P.P. was supported by the Rare Isotope Science Project of the
Institute for Basic Science funded by the Ministry of Science, ICT and Future
Planning and the National Research Foundation (NRF) of Korea
(2013M7A1A1075764).

## References

  * Giacalone  _et al._ [2021] G. Giacalone, J. Jia, and C. Zhang, Phys. Rev. Lett. 127, 242301 (2021).
  * Zhang and Jia [2022] C. Zhang and J. Jia, Phys. Rev. Lett. 128, 022301 (2022).
  * Nijs and van der Schee [2023] G. Nijs and W. van der Schee, (2023), arXiv:2112.13771 [nucl-th] .
  * STAR Collaboration [2021] STAR Collaboration, (2021), arXiv:2109.00131 [nucl-ex] .
  * Lane and Pendlebury [1960] A. M. Lane and E. D. Pendlebury, Nucl. Phys. 15, 39 (1960).
  * Butler [2016] P. A. Butler, J. Phys. G: Nucl. Part. Phys. 43, 073002 (2016).
  * Butler and Nazarewicz [1996] P. A. Butler and W. Nazarewicz, Rev. Mod. Phys. 68, 349 (1996).
  * Spear [1989] R. H. Spear, At. Data Nucl. Data Tables 42, 55 (1989).
  * Kibédi and Spear [2002] T. Kibédi and R. H. Spear, At. Data Nucl. Data Tables 80, 35 (2002).
  * Spear and Catford [1990] R. H. Spear and W. N. Catford, Phys. Rev. C 41, R1351 (1990).
  * Robledo and Bertsch [2011] L. M. Robledo and G. F. Bertsch, Phys. Rev. C 84, 054302 (2011).
  * Agbemava  _et al._ [2016] S. E. Agbemava, A. V. Afanasjev, and P. Ring, Phys. Rev. C 93, 044304 (2016).
  * Yao  _et al._ [2015] J. M. Yao, E. F. Zhou, and Z. P. Li, Phys. Rev. C 92, 041304 (2015).
  * Bernard  _et al._ [2016] R. N. Bernard, L. M. Robledo, and T. R. Rodríguez, Phys. Rev. C 93, 061302 (2016).
  * Li  _et al._ [2016] Z. P. Li, T. Nikšić, and D. Vretenar, J. Phys. G: Nucl. Part. Phys. 43, 024005 (2016).
  * Nomura  _et al._ [2021a] K. Nomura, L. Lotina, T. Nikšić, and D. Vretenar, Phys. Rev. C 103, 054301 (2021a).
  * Nomura  _et al._ [2021b] K. Nomura, R. Rodríguez-Guzmán, L. M. Robledo, and J. E. García-Ramos, Phys. Rev. C 103, 044311 (2021b).
  * Iskra  _et al._ [2019] Ł. W. Iskra, R. Broda, R. V. F. Janssens, M. Carpenter, B. Fornal, T. Lauritsen, T. Otsuka, T. Togashi, Y. Tsunoda, W. B. Walters, and S. Zhu, Phys. Lett. B 788, 396 (2019).
  * Spieker  _et al._ [2022] M. Spieker, L. A. Riley, P. D. Cottle, K. W. Kemper, D. Bazin, S. Biswas, P. J. Farris, A. Gade, T. Ginter, S. Giraud, J. Li, S. Noji, J. Pereira, M. Smith, D. Weisshaar, and R. G. T. Zegers, Phys. Rev. C 106, 054305 (2022).
  * Abbas  _et al._ [1981] A. Abbas, N. Auerbach, N. Van Giai, and L. Zamick, Nucl. Phys. A 367, 189 (1981).
  * Co’ [2023] G. Co’, Universe 9, 141 (2023).
  * Blaizot [1976] J. P. Blaizot, Phys. Lett. B 60, 435 (1976).
  * Thouless [1960] D. J. Thouless, Nucl. Phys. 21, 225 (1960).
  * Thouless [1961] D. J. Thouless, Nucl. Phys. 22, 78 (1961).
  * Rowe [1968] D. J. Rowe, Rev. Mod. Phys. 40, 153 (1968).
  * Bender  _et al._ [2003] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * Ring and Schuck [1980] P. Ring and P. Schuck,  _The Nuclear Many-Body Problem_ (Springer-Verlag New York, 1980).
  * Colò  _et al._ [2013] G. Colò, L. Cao, N. Van Giai, and L. Capelli, Comput. Phys. Commun. 184, 142 (2013).
  * Colò and Roca-Maza [2021] G. Colò and X. Roca-Maza, arXiv:2102.06562 (2021).
  * Bartel  _et al._ [1982] J. Bartel, P. Quentin, M. Brack, C. Guet, and H.-B. Håkansson, Nucl. Phys. A 386, 79 (1982).
  * Chabanat  _et al._ [1998] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
  * Beiner  _et al._ [1975] M. Beiner, H. Flocard, N. Van Giai, and P. Quentin, Nucl. Phys. A 238, 29 (1975).
  * Bertsch and Tsai [1975] G. F. Bertsch and S. F. Tsai, Phys. Rep. 18, 125 (1975).
  * Auerbach  _et al._ [1996] N. Auerbach, V. V. Flambaum, and V. Spevak, Phys. Rev. Lett. 76, 4316 (1996).