Tackling Algorithmic Bias in Neural-Network Classifiers using Wasserstein-2 Regularization

TL;DR

This paper introduces a Wasserstein-2 regularization method to reduce algorithmic bias in neural network classifiers, applicable to large datasets and compatible with standard training procedures.

Contribution

A novel Wasserstein-2 based regularization technique for neural networks that mitigates bias without altering architecture and scales efficiently to large datasets.

Findings

Effective bias reduction demonstrated on Adult census, MNIST, CelebA datasets.

Method is architecture-agnostic and compatible with standard SGD training.

Regularization improves fairness metrics without sacrificing accuracy.

Abstract

The increasingly common use of neural network classifiers in industrial and social applications of image analysis has allowed impressive progress these last years. Such methods are however sensitive to algorithmic bias, i.e. to an under- or an over-representation of positive predictions or to higher prediction errors in specific subgroups of images. We then introduce in this paper a new method to temper the algorithmic bias in Neural-Network based classifiers. Our method is Neural-Network architecture agnostic and scales well to massive training sets of images. It indeed only overloads the loss function with a Wasserstein-2 based regularization term for which we back-propagate the impact of specific output predictions using a new model, based on the Gateaux derivatives of the predictions distribution. This model is algorithmically reasonable and makes it possible to use our regularized loss with standard stochastic gradient-descent strategies. Its good behavior is assessed on the reference Adult census, MNIST, CelebA datasets.

Tables

Table 1. Table 1: Adult dataset with sensitive variable Gender

	Acc	DI	$O_{1,0}$	$O_{0,0}$	$O_{1,1}$	$O_{0,1}$
LR	0.84	0.47	0.18	0.05	0.71	0.80
NN	0.82	0.37	0.20	0.05	0.73	0.64
ZFN	0.79	0.64	0.10	0.07	0.42	0.43
ZFA	0.66	0.95	0.39	0.42	0.81	0.83
LRrW	0.65	0.94	0.40	0.44	0.85	0.86
NNrW	0.78	0.68	0.25	0.20	0.85	0.83

Table 2. Table 2: Impact of the regularization on the generalization properties of two trained classifiers using the experimental protocol (SR) of Section  4.2 .

	Train			Test			Test
	DI	GP0	GP1	DI	GP0	GP1	TP0	TN0	TP1	TN1
$\lambda =0.$	0.87	1.00	1.00	0.86	0.92	0.96	0.83	0.97	0.96	0.97
$\lambda =2.0e-5$	0.98	0.94	0.96	0.96	0.93	0.95	0.91	0.95	0.93	0.97

Equations

$$DI(g_{\theta },X,S)=\frac{\min {\left[\mathbb{P}(g_{\theta }(X)=1|S=s)\right]}_{s\in \{0,1\}}}{\max {\left[\mathbb{P}(g_{\theta }(X)=1|S=s)\right]}_{s\in \{0,1\}}}.$$

$$O_{s,y}=\mathbb{P}(\hat{Y}=1|S=s,Y=y).$$

$$\mathcal{R}(\theta ):=\mathbb{E}[\ell (f_{\theta }(X),Y)],$$

$${\mathcal{W}}_2^2({\mu }_{\theta ,0},{\mu }_{\theta ,1})={\int }_0^1{\left({{\mathcal{H}}_0}^{-1}(\tau )-{{\mathcal{H}}_1}^{-1}(\tau )\right)}^{2}d\tau .$$

$$\hat{\theta }=\underset{\theta \in \Theta }{\mathrm{argmin}}\left\{\mathcal{R}(\theta )+\lambda {\mathcal{W}}_2^2({\mu }_{\theta ,0},{\mu }_{\theta ,1})\right\},$$

$$D{\mathcal{W}}_2^2(\mu ,\nu )(\alpha ,\beta ):=\underset{t\searrow 0}{\lim }\frac{{\mathcal{W}}_2^2(\mu +t\alpha ,\nu +t\beta )-{\mathcal{W}}_2^2(\mu ,\nu )}{t}$$

$$D{\mathcal{W}}_2^2(\mu ,\nu )(\alpha ,\beta )=2\int {\int }_0^x(t-F_{\nu }^{-1}(F_{\mu }(t)))dtd\alpha (x)+2\int {\int }_0^x(t-F_{\mu }^{-1}(F_{\nu }(t)))dtd\beta (x),$$

$$\int {\int }_0^x(t-F_{\nu }^{-1}(F_{\mu }(t)))dtd\alpha (x)=\int {\int }_0^x\frac{t-F_{\nu }^{-1}(F_{\mu }(t))}{u(t)}d\mu (t)d\alpha (x),$$

$$\int {\int }_0^x\frac{t-F_{\nu }^{-1}(F_{\mu }(t))}{u(t)}d\mu (t)d{\alpha }_{ϵ}^i(x)={\int }_{f_{\theta }(x_i)-ϵ}^{f_{\theta }(x_i)+ϵ}{\int }_0^x\frac{t-F_{\nu }^{-1}(F_{\mu }(t))}{u(t)}d\mu (t)d{\alpha }_{ϵ}^i(x)$$

$$\frac{1}{2ϵ}\left[{\int }_{f_{\theta }(x_i)-ϵ}^{f_{\theta }(x_i)+ϵ}\underset{\Phi (x)}{\underbrace{\left({\int }_0^x\frac{t-F_{\nu }^{-1}(F_{\mu }(t))}{u(t)}d\mu (t)\right)}}2ϵd{\alpha }_{ϵ}^i(x)\right],$$

$$\frac{1}{2ϵ}\left[{\int }_0^{f_{\theta }(x_i)+ϵ}\Phi (x)dx-{\int }_{f_{\theta }(x_i)-ϵ}^0\Phi (x)dx\right],$$

$$\frac{f_{\theta }(x_i)-F_{\nu }^{-1}(F_{\mu }(f_{\theta }(x_i)))}{u(f_{\theta }(x_i))}.$$

$$\underset{ϵ\searrow 0}{\lim }D{\mathcal{W}}_2^2(\mu ,\nu )({\alpha }_{ϵ}^i,0)=\frac{f_{\theta }(x_i)-{\mathcal{H}}_1^{-1}({\mathcal{H}}_0(f_{\theta }(x_i)))}{{\mathcal{H}}_0^{\prime }(f_{\theta }(x_i))}$$

$$L:\mathcal{P}(\mathbb{R}):$$

$$\mu$$

$$DL(\mu ;\beta )=\int ld\beta .$$

$$L(\mu )=L({\mu }_{\theta ,0})P(S=0)+L({\mu }_{\theta ,1})P(S=1),$$

$$DL({\mu }_{\theta };{\alpha }_{ϵ}^i)$$

$$\frac{1}{\#B}\sum _{i\in B}(\nabla l)(f(x_i)),$$

$${\mu }_{\theta }⟼\lambda {\mathcal{W}}_2^2({\mu }_{\theta ,0},{\mu }_{\theta ,1}).$$

$${\eta }^j=\underset{i}{\min }(f_{\theta }(x_i))+j{\Delta }_{\eta },j=1,\dots ,J_{\eta },$$

$${\mathcal{W}}_2^2({\mu }_{\theta ,0}^n,{\mu }_{\theta ,1}^n)={\int }_0^1{\left(H_0^{-1}(\tau )-H_1^{-1}(\tau )\right)}^{2}d\tau ,$$

$$\frac{{\Delta }_{\tau }(f_{\theta }(x_i)-co{r}_{|1-s_i|}(f_{\theta }(x_i)))}{(H_{s_i}^{j_i+1}-H_{s_i}^{j_i})},$$

$${\Delta }_{\tau }\left[\text{mathds}{1}_{s_i=0}\frac{f_{\theta }(x_i)-co{r}_1(f_{\theta }(x_i))}{n_0\left(H_0^{j_i+1}-H_0^{j_i}\right)}-\text{mathds}{1}_{s_i=1}\frac{co{r}_0(f_{\theta }(x_i))-f_{\theta }(x_i)}{n_1\left(H_1^{j_i+1}-H_1^{j_i}\right)}\right],$$

$$DMSE(g_{\theta },X,S)=\frac{\min {\left[MSE(g_{\theta },X,Y|S=s)\right]}_{s\in \{0,1\}}}{\max {\left[MSE(g_{\theta },X,Y|S=s)\right]}_{s\in \{0,1\}}},$$

$$MSE(g_{\theta },X,Y|S=s)$$

$$cMSE(f_{\theta },X,Y|S=s)=\frac{1}{n_s}\sum _{i=1,s_i=s}^n(f_{\theta }(x_i)-y_i{)}^2,$$

$${\tilde{H}}_s^j=n_s^{-1}\sum _{i=1,s_i=s}^n\text{mathbbm}{1}_{(f_{\theta }(x_i)-y_i{)}^2<{\eta }^j},$$

$$\hat{\theta }=\underset{\theta \in \Theta }{\mathrm{argmin}}\left\{R(\theta )+\lambda W_2^2({\tilde{\mu }}_{\theta ,0}^n,{\tilde{\mu }}_{\theta ,1}^n)\right\}.$$

$$2{\Delta }_{\tau }[\text{mathds}{1}_{s_i=0}\frac{(f_{\theta }(x_i)-y_i{)}^2-co{r}_1\left((f_{\theta }(x_i)-y_i{)}^2\right)}{n_0\left({\tilde{H}}_0^{j_i+1}-{\tilde{H}}_0^{j_i}\right){\left(f_{\theta }(x_i)-y_i\right)}^{-1}}$$

Full text

K****eywords Algorithmic bias  $\cdot$ Image Classification  $\cdot$ Neural-Networks  $\cdot$ Regularization

1 Introduction

1.1 Algorithmic bias in Machine Learning

Recent neural-network models have become extremely popular for a large variety of applications in image analysis. They have indeed made it possible to strongly improve the accuracy of classic image analysis strategies for the detection of specific shapes, as originally shown in [1, 2], or other applications of image analysis. A huge effort has then been made to design neural network architectures with desirable properties [3, 4, 5] and to efficiently optimize their parameters [6, 7]. How powerful they might be, these models may however suffer from algorithmic bias. This bias impacts the predictions, generally in the detriment of a given subset of observations. For instance, the prediction errors may be higher in a distinct images subset compared than in the other images. The portion of predictions of a given label may also be particularly low in a given subgroup of images with no valid reason. The effect of this bias in the society has been studied recently in various papers [8, 9, 10]. In this paper, we focus on a new strategy to tackle it when detecting specific image features using neural networks.

A popular illustration of the effects of algorithmic bias on images was shown in [9], where the authors measured that different publicly available commercial face recognition online services achieved their lowest accuracy on females. For instance, Face++ had a 99.3% accuracy to recognize male politicians and only a 78.7% for female politicians, so about 95.9% of the errors were obtained on females. In the general case, the literature dealing with algorithmic bias in Machine Learning considers that the bias comes from the impact of a non-informative variable denoted sensitive (or protective) variable. This variable is the gender in the former example. It splits the observations into two subsets for which the Machine Learning predictions have different behaviors, implying that they are biased with respect to this variable. These biases are widely discussed in the Fair learning community, which studies their effects in the society. Because of the increasingly massive use of machine learning-based applications in the everyday life, justice, or resource allocations, this community has become extremely dynamic, as shown by the strong emergence of a conference like ACM FAT111https://www.fatml.org/*. The Computer Vision community has also drawn a significant attention to such questions recently with for instance the contribution of [11], the organization of the CVPR 2018 workshop Computer vision meets fairness, and the CVPR 2019 workshop Fairness Accountability Transparency and Ethics in Computer Vision.

To make obvious that these biases are common in Machine-Learning based applications of image analysis, we detail hereafter three major causes of algorithmic bias in Machine-Learning: (Cause 1) Machine learning algorithms are first meant to automatically take accurate and efficient decisions that mimic human expertise, based on reference datasets that are potentially biased with respect to one or several protective variables. An algorithm is then likely to reproduce the errors contained in a biased training database, although there is generally no intention of doing so. A biased training set, may be due to specific persons, who label the training data with little experience. It may also be due to a reference algorithm which labels automatically the data with a bias on a specific type of data. In the fair-learning community, it is often considered as being due to societal biases in the sampled populations or to systematic biases made by the persons who took the reference decisions. Note that these errors are considered as systematic in the sensitive sub-group, at least in law. (Cause 2) The algorithmic bias may also be due to the fact that a classifier overfits or underfits a specific sub-group of training observations and adequately fits the other observations. This leads to different generalization error properties in the sensitive subgroup compared with the other observations. (Cause 3) A more subtle issue, which is at the heart of [10] is finally related to the influence of confounding variables: the prediction accuracy can be penalized in a subset of observations because they share common properties with another subset of observations having a different output. These common properties can be related to an explicit variable in the data or to unobserved latent variables. In neural-network based applications with images as inputs, the confounding variables will generally be a subset of pixels with unforeseen specific properties. This was made popular by the Husky or Wolf example [12] where huskies are classified as wolves when they are represented in images with snow because the wolves of the training set are generally represented as surrounded by snow.

Tackling the algorithmic bias in machine learning is then an important and ambiguous task. This is particularly true when using deep neural-network models on images for which the decision rules are humanly impossible to interpret in the general case. In this paper, we then propose a new solution to tackle the algorithm bias issue in neural-network based classification of image. The key advantages of our strategy are that it is Neural-Network model agnostic and that it scales particularly well to massive training sets of images. Two variants of our method are presented, each of them being pertinent in specific contexts. The first variant helps the classifiers to predict similar frequencies of positive outputs in two groups of data (e.g. in pictures representing males or females). It is then related to the notion of statistical parity. The second variant favor similar error rates in the subgroups and is therefore related to the complementary notion of Equalized odds [13]. This second variant is of particular interest in industrial frameworks, where a similar accuracy has to be guaranteed in different contexts. Hence removing the bias is an important task when dealing with the generalization of predictions with different conditions modeled by the variable $S$. It appears to be an important tool for domain adaptation.

2 Bibliography

2.1 Standard measures of the algorithmic bias

Different indices are commonly used to measure the algorithmic bias. We refer to [13, 14, 15] for recent reviews of these measures. In order to introduce them properly, we first denote $(x_{i},s_{i},y_{i})_{i=1,\dots,n}$ the training observations, where $x_{i}\in\mathbb{R}^{p}$ is an input image with $p$ pixels or voxels. Although we will treat 2D RGB images in our tests the $x_{i}$ may be in any dimension (2D, 3D, 2D+t, $\ldots$) and contain various amount of channels. In our paper, the output prediction $x_{i}\in\{0,1\}$ related to $x_{i}$ is supposed to be binary for the sake of simplicity. We refer to [16, 17] for discussions about the multiple protective attributes and regression cases. The protective variable $S\in\{0,1\}$ indicates whether observation $i$ is in a group which may be subject to algorithmic bias or not. It is common to use $s_{i}=0$ for the group which may be penalized and $s_{i}=1$ for the other data. The parameters $\theta$ of a binary classifier $g_{\theta}$ are trained using the observations $(x_{i},s_{i},y_{i})_{i=1,\dots,n}$. The trained classifier is then used to predict the outputs $y_{i}$ of new input observations $x_{i}$, with $i>n$, so a prediction $\hat{y}_{i}=g_{\theta}(x_{i})\in\{0,1\}$. Note that for each observation $i$, the neural network classifier returns in practice a score $f_{\theta}(x_{i})\in[0,1]$ which measures how $g_{\theta}(x_{i})$ likely is to be equal to $1$. In general, $g_{\theta}(x_{i})=\mathds{1}_{f_{\theta}(x_{i})>\frac{1}{2}}$.

The most standard measure of algorithmic bias is the so-called Statistical Parity, which is often quantified in the fair learning literature using the Disparate Impact (DI). This was been introduced in the US legislation in 1971222https://www.govinfo.gov/content/pkg/CFR-2017-title29-vol4/xml/CFR-2017-title29-vol4-part1607.xml and measures how the variable $Y$ depends on $S$.

More specifically, $DI(Y,S)$ is defined as the ratio between $\mathbb{P}(Y=1|S=0)$ and $\mathbb{P}(Y=1|S=1)$, where we suppose that $S=0$ is the group which may be discriminated with respect to $Y$. The smaller this index, the stronger the discrimination over the group $S=0$. A threshold $\tau_{0}=0.8$ is commonly used to judge whether the discrimination level of an algorithm is acceptable or not [18, 19, 20, 21]. This fairness criterion can be straightforwardly extended to the outcome of an algorithm by using $\hat{Y}=g_{\theta}(X)$ instead of the true variable $Y$, and by making no hypothesis on the potentially discriminated group:

[TABLE]

A concern with the notion of statistical parity is that it implies that the decisions should be statistically the same for all subgroups of $S$, while in many situations, the end-user is more interested to detect whether the prediction algorithm has the same behavior for the observations out of each subgroup. In particular, the statistical parity does not take into account the false positive and false negative predictions, and more generally the prediction errors. As discussed in [13], the notions of equality of odds and opportunity additionally use the true predictions $Y$ and may then be more suitable than the statistical parity when a similar prediction accuracy (i.e. $Acc=\mathbb{P}(\hat{Y}=Y)$) is desired in the subgroups $S=0$ and $S=1$. We use the following notation to discuss these notions

[TABLE]

In the binary case, the classifier $g_{\theta}$ gives an equal opportunity with respect to the protective attribute $S$ when $O_{0,1}$ equals $O_{1,1}$, which means that the true positive prediction rate is the same in $S=0$ and $S=1$. Equalized odds with respect to $S$ are also satisfied if $O_{0,0}$ is additionally equal to $O_{1,0}$, which means that the false positive prediction rate is also the same $S=0$ and $S=1$. As the exact equality between two empirical probabilities does not necessarily make sense in practice, the predictions of a binary classifier will be considered as fair if $O_{0,y}$ reasonably close to $O_{1,y}$, for the true outputs $y=0$ and $y=1$.

We finally want to emphasize that a purely random binary classifier has a $DI$ close to $1$, and values of $O_{0,0}$ and $O_{0,1}$ similar to $O_{1,0}$ and $O_{1,1}$, respectively. Its predictions are indeed independent of $S$. Despite of these nice properties, its accuracy is obviously almost equal to $0.5$, so such a classifier is useless in practice. We then consider in this paper that a fair binary classifier has only an interest if its accuracy is reasonably close to $1$.

2.2 Tackling the algorithmic bias

In all generality, the notion of fairness in Machine Learning in [8, 22, 23, 24], consists in modeling the algorithmic bias as the independence (or the conditional independence) between the output $Y$ of an algorithm and a variable $S$ (or a group of variables). Exhaustive bibliographies dealing with how to tackle this bias can be found in [24, 25, 21, 19, 26, 27, 13]. Achieving this independence is indeed a difficult task which is obtained in the literature by favoring various fairness criteria that convey a part of the independence between $Y$ and $S$. In a sense, the algorithmic bias criteria used in the literature enable to weaken the notion of independence to achieve a good trade-off between fairness and prediction accuracy. This idea was recently formalized as a price for fairness in [28]. Among the most popular fairness criteria, the Disparate Impact, the Equality of opportunity and the Equality of Odds [13] presented in Section 2.1 are indirectly based on the  probability distribution of the outputs $f_{\theta}(X)$ with respect to a protected variable $S$. Their interpretation is clear and they are, in the authors opinion, powerful and straightforwardly interpretable tools to detect unfair decision rules.

An important concern with these measures is however that they are not smooth, which we believe is necessary to blend mathematical theory with algorithmic practice, in particular when using gold-standard gradient descent based optimization methods: For instance, if two observations $i$ and $j$ have similar non-binarized outputs $|f_{\theta}(x_{i})-f_{\theta}(x_{j})|<\epsilon$ (with $\epsilon$ small) but different binarized outputs $\hat{y_{i}}=g_{\theta}(x_{i})=\mathds{1}_{f_{\theta}(x_{i})>0.5}=0$ and $\hat{y_{j}}=g_{\theta}(x_{j})=\mathds{1}_{f_{\theta}(x_{j})>0.5}=1$ these fairness criteria do not take into account the fact that the outputs $f_{\theta}$ are nearly the same. They won’t therefore directly favor observation $i$ instead of another observation to make the predictions more fair. As a result, this can make it hard to numerically find a good balance between fair and accurate predictions.

The covariance of the outputs $Y$ with respect to $S$ has also been investigated for instance in [19] to directly measure linear independence. The variability of the loss function with respect to the protected variable has also been considered for similar purposes in [23]. These two strategies have again these concerns dealing with the fact that they only measure the fairness based on $Y$. This clearly justifies for us the use of fairness measures taking into account the whole distribution, as the Wasserstein Distance, in particular when learning fair decision rules and not only when detecting unfair decision rules. Another point view is given by considering the distance between conditional distributions. General divergences, Kullback-Leibler divergence or total variation have been considered for Hilbert-Schmidt dependency in [29, 30, 31]. Contrary to a transport based strategy like the Wasserstein distance, these measure are however not symmetric and cannot therefore be considered as metrics. Moreover the Wasserstein metric will make it possible to compare different measures, even if the intersection of their supports is null, which can be the case in real-life applications.

2.3 Overview of the contribution

We therefore propose to use the Wasserstein metric when training Neural Networks, which is smoothly defined and does not degenerate on empirical distributions given by finite samples. This distance has been already used in previous works dealing with fairness either to repair data or to build fair algorithms [24, 18].  Wasserstein distance appears in this framework as a smooth criterion to assess the sensitivity w.r.t to the protected variable. As a matter of fact, it can be considered as a distance between the quantiles of score function of the predictor for the two groups (S=0 and S=1). The distance between the quantiles of these two distributions therefore acts as a level of fairness measuring whether the spread of the scores is homogeneously spread whatever the values of the protected attribute, hence acting as a  sensitivity index of the predicted values $f_{\theta}$ with respect to $S$. Remark that, in [24], which is the closest work to our contribution, the authors specifically used Wasserstein-1 to post-process the output of a model to favor fair predictions. In our work, we instead use a more feasible criterion based on Wasserstein-2 regularization, as motivated in [21] to impose fairness constraints. In some technical aspects, this also distinguishes our method from the optimization models used in Wasserstein GANS [32, 33], which are also based on Wasserstein-1 regularization. Computing for a neural network the gradient of Wasserstein-2 distances with respect to machine learning model parameters is however an important technical lock to address, which is at the heart of Section 3. We finally remark that relationships between Wasserstein distances and usual disparate impact are further discussed in [34].

We additionally focus on how to train neural network classifiers with fairness constraints for images since neural networks are particularly flexible models and can treat huge volumes of data. Very little work has however been done so far to ensure fair decisions with neural networks [35, 36, 37]. The existing literature also does not explicitly explain how to compute the gradients of the loss terms ensuring fair decisions. It additionally does not use Wasserstein-based regularization, which was recently shown in [24] as an interesting solution for this purpose. In addition to the use Wasserstein-2 regularization instead of Wasserstein-1 regularization, our methodology differs from [24] in different aspects. We first make our predictions using neural-networks and not logistic-regression in order to train more flexible decision rules. By using the standard back-propagation formalism, we therefore do not directly compute the Wasserstein distance gradients with respect to the model parameters but with respect to the model outputs. As the Wasserstein distances are based on probability distribution, we then formalize and rigorously model a new interpretation of the impact of specific neural-network outputs on the continuous distribution of its outputs. We then show how to discretize this model with a numerical strategy that scales well to large datasets. We finally explicitly develop how to handle the estimation of these gradients when using mini-batch training, which is extremely common for neural networks.

3 Methodology

3.1 Main notations

We recall that $(x_{i},s_{i},y_{i})_{i=1,\dots,n}$ are the training observations, where $x_{i}\in\mathbb{R}^{p}$ and $y_{i}\in\{0,1\}$ are the input and output observations, respectively. Remark here that we denote by $x_{i}$, resp. $y_{i}$, i.i.d observations of the random variable $X$, resp. $Y$. The protective variable $S\in\{0,1\}$ indicates whether observation $i$ is in a sensitive group ($s_{i}=0$) or not ($s_{i}=1$). A classifier $\hat{g}_{\theta}$ with parameters $\theta$ is trained so that its predictions $\hat{y_{i}}=g_{\theta}(x_{i})\in\{0,1\}$ are, as often as possible, equal to the output observations $y_{i}$. Importantly, the binary predictions $g_{\theta}(x_{i})$ are defined as equal to $\mathds{1}_{f_{\theta}(x_{i})>0.5}$, where the continuous score $f_{\theta}\in[0,1]$ is the actual neural network output. The problem can then be solved by minimizing a risk

[TABLE]

or empirically $R(\theta)=\frac{1}{n}\sum_{i=1}^{n}\ell(f_{\theta}(x_{i}),y_{i})$, where the loss function $\ell$ represents the price paid for inaccuracy of predictions [38].

3.2 Binary classification using Neural Networks

The optimization of the parameters $\theta$ over a compact set $\Theta$ is typically achieved using standard stochastic gradient descent [39, 40] or its variants [41, 6] which require less computations than standard gradient descent and enable to explore more efficiently the parameters space. At each iteration of the stochastic gradient descent, the average gradient of $R(\theta)$, and then $\mathcal{R}(\theta)$, is approximated on a subset $B$ of several observations. This subset is denoted a mini-batch and contains an amount of $\#B$ observations.

Training neural-networks is specific in the sense that the parameters $\theta$ are indirectly optimized based on the derivatives of $R(\theta)$ with respect to network outputs $f_{\theta}(x_{i})$, $i\in B$. In our tests, we defined $\ell(f_{\theta}(x_{i}),y_{i})$ is equal to $\left(f_{\theta}(x_{i})-y_{i}\right)^{2}$. The empirical derivatives are then computed using $2\left(f_{\theta}(x_{i})-y_{i}\right)$. These derivatives are then back-propagated in the neural network and the parameters $\theta$ are updated using the stochastic gradient descent approach [42, 38]. Many modern tools such as TensorFlow, Keras or PyTorch, make it simple to implement such training strategies based on automatic differentiation. Note that although we used the mean squared error loss in our paper, the regularization method proposed in this paper can be directly used with other losses (e.g. Binary Cross Entropy) as long as their derivative with respect to $f_{\theta}(x_{i})$ can be computed.

3.3 Wasserstein-2 based regularization to favor low discriminate impacts

We denote by $\mu_{\theta,0}$ and $\mu_{\theta,1}$ the output distributions of $f_{\theta}(X)$ for observations in the groups $S=0$ and $S=1$, respectively, and denote by $h_{0}$ and $h_{1}$ their densities. Our regularization strategy consists in ensuring that the Wasserstein-2 distance (or Kantorovich-Rubinstein metric) between the distributions of $\mu_{\theta,0}$ and $\mu_{\theta,1}$ remains small compared with $\mathcal{R}(\theta)$. Their corresponding cumulative distribution functions are $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$. The Wasserstein-2 distance between the two conditional distributions is defined as

[TABLE]

where ${\mathcal{H}_{s}}^{-1}$ is the inverse of the cumulative distribution function $\mathcal{H}_{s}$. Using the distance Eq. (4) as a regularizer ensures that the distributions of $f_{\theta}(X)$ for the observations in the groups $S=0$ and $S=1$ remain reasonably close to each-other. The training problem in a continuous setting is then:

[TABLE]

where $\lambda$ is a weight giving more or less influence to the regularization term compared with the prediction accuracy.

3.4 Fast estimation of Wasserstein-2 pseudo-gradients in a batch

We have seen in Section 3.2 that neural-network parameters are typically trained using the average gradient of the empirical risk $\mathcal{R}(\theta)$ with respect to all parameters on mini-batch observations. Our goal here is to minimize Eq. (5) when we observe empirical versions of the distributions $\mu_{\theta,0}$ and $\mu_{\theta,1}$ in the regularization term. Standard automatic differentiation, which is ubiquitously used to train deep neural networks, cannot indeed be used to minimize Eq. (5) since the Wasserstein constraint Eq. (4) is applied on the distributions $\mu_{\theta,0}$ and $\mu_{\theta,1}$ and not on a subset of observations. In this section, we then explain how we compute, and which sense we give, to the gradient of Eq. (4) with respect to a finite subset of neural-network predictions $f_{\theta}(x_{i})$, $i\in B$. To achieve this, we model the regularization term of Eq. (5) as a function of the data distribution, denoted as $\mu_{\theta}$ and introduce a notion of derivative of $\mu_{\theta}\mapsto\Phi(\mu_{\theta})$ with respect to the distributional changes. Discretization of our methodology is then developed in Section 3.4.4.

Note that similar arguments have for instance been used in [43] in the discrete case or [44] in the semi-discrete case. A similar notion of derivative in a geometric sense was also developed in the Proposition 7.17. of [45]. None of these contributions however specifically modeled how to backpropagate the impact of probability density variations.

3.4.1 Modeling the impact of probability density variations on the transport cost

In all generality, let $\mu,\nu\in\mathcal{P}(\mathbb{R})$ be two probability measures over $\mathbb{R}$. Consider as well $(\alpha,\beta)\in\mathcal{M}(\mathbb{R})\times\mathcal{M}(\mathbb{R})$, where $\mathcal{M}(\mathbb{R})$ denotes the space of regular measures over $\mathbb{R}$. Inspired by [46, 47], we can remark that the transport cost is Gâteaux differentiable in the direction $(\alpha,\beta)$, with derivative $D\mathcal{W}^{2}_{2}(\mu,\nu)(\alpha,\beta)$ if the limit

[TABLE]

exits and is finite. Note that the value $\mathcal{W}^{2}_{2}(\mu+t\alpha,\nu+t\beta)$ only makes sense when $\alpha(\mathbb{R})=\beta(\mathbb{R})$, where $\alpha(\mathbb{R})$ and $\beta(\mathbb{R})$ denote the total masses of the measures and are not necessarily positive (see e.g. [48], chapter 6). Otherwise, the transport would not be defined any longer. Remark that $\mu+t\alpha$ and $\nu+t\beta$ should not take negate values too, as they represent probability measures.

Proposition 1

Let $\mu,\nu\in\mathcal{P}(\mathbb{R})$ with density w.r.t Lebesgue measure and supported on respective connected sets and $(\alpha,\beta)\in\mathcal{M}(\mathbb{R})\times\mathcal{M}(\mathbb{R})$ be such that

$\alpha(\mathbb{R})=\beta(\mathbb{R})$, 2. 2.

$\operatorname{Supp}(\alpha)$* and $\operatorname{Supp}(\beta)$ are both compact sets,* 3. 3.

$\operatorname{Supp}(\alpha)\subset\operatorname{Supp}(\mu)$* and $\operatorname{Supp}(\beta)\subset\operatorname{Supp}(\nu)$.*

Then the Transport cost is Gâteaux differentiable in the direction $(\alpha,\beta)$ and it satisfies

[TABLE]

where $F_{\mu}$ and $F_{\nu}$ are the cumulative distribution functions of $\mu$ and $\nu$.

A proof of Proposition 1 is given in Appendix A.

In our application, our goal is to quantify the impact of variations of $f_{\theta}(x_{i})$ on either $\mu_{\theta,0}$ or $\mu_{\theta,1}$, depending on the label of $s_{i}$. Suppose now that we are in the case $s_{i}=0$. In order to later model a perturbation of only one of the two distributions, we consider now the direction $(\alpha,0)$, where $\alpha(\mathbb{R})=0$, so that $\mu+t\alpha$ remains a probability on $\mathbb{R}$. Let $u$ be the density of $\mu$. As $d\mu(t)=u(t)d(t)$, we can write that

[TABLE]

in the case where the right side is finite. This will be the case in our specific application, as the predictions $f_{\theta}(x_{i})$ are supposed to be finite.

3.4.2 Influence of the perturbation of an output $f_{\theta}(x_{i})$ on the Wasserstein-2 distance

We now specify Eq. (7) to quantify the influence of a perturbation of the output prediction $f_{\theta}(x_{i})$, which is the very first step of a back-propagation algorithm when training a neural network. To model model a local perturbation around $f_{\theta}(x_{i})$ using a measure, we first consider $\epsilon>0$ and set $\alpha_{\epsilon}^{i}=\frac{1}{2\epsilon}\mathbb{I}_{[0,f(x_{i})+\epsilon]}-\frac{1}{2\epsilon}\mathbb{I}_{[f(x_{i})-\epsilon,0[}$. Hence $\alpha_{\epsilon}^{i}$ is supported in $[f_{\theta}(x_{i})-\epsilon,f_{\theta}(x_{i})+\epsilon]$, such that $\alpha_{\epsilon}^{i}(\mathbb{R})=0$ and is anti-symmetric with respect to $f_{\theta}(x_{i})$. We then have

[TABLE]

This equation can be re-written as:

[TABLE]

which is equal to

[TABLE]

and tends towards $\Phi^{\prime}(f_{\theta}(x_{i}))$ when $\epsilon\rightarrow 0$. In this case, Eq. (8) is then approximated by:

[TABLE]

We then have the following proposition:

Proposition 2

Let $\mu=\mu_{\theta,0}$ and $\nu=\mu_{\theta,1}$, then we have that

[TABLE]

This represents the influence of a perturbation made by the output prediction $f_{\theta}(x_{i})$ if $s_{i}=0$. Note that the cumulative distribution functions $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ will straightforwardly be swapped in this equation if $s_{i}=1$.

3.4.3 Impact of several observations in a mini-batch $B$

The following lemma allows to extend the previous result to the mini-batch case.

Lemma 1

let $\mu\in\mathcal{P}(\mathbb{R})$ be a probability measure and $l$ be a real measurable function, then the function

[TABLE]

is Gâteaux differentiable for every direction $\beta\in\mathcal{M}(\mathbb{R})$ such that $\int ld\beta<\infty$ with derivative

[TABLE]

A proof of Lemma 1 is given in Appendix A. If $\mu_{\theta}$ depends on the protected attribute $S$, it is therefore easy to see that

[TABLE]

then we can use the previous result with the direction $\alpha^{i}_{\epsilon}$ defined before and deduce that

[TABLE]

As a consequence, we can define the mean influence of the batch $B$ as the mean of the previous limit, which yields:

[TABLE]

Together with Proposition 1, this proves the Gâteaux differentiability of

[TABLE]

when using batch training.

3.4.4 Influence of an output $f_{\theta}(x_{i})$ in the discrete case

We now extend the methodology of Section 3.4.2 to the discrete case, where the number of observations $f_{\theta}(x_{i})$ is finite and the cumulative distributions are only known on discrete grids of values. We first denote $\mu_{\theta,s}^{n}$ the approximation of $\mu_{\theta,s}$ made using $n$ observations. For the group $S=s$, we denote $H_{s}(\eta)$ the discrete approximation of $\mathcal{H}_{s}(\eta)$. It represents the empirical portion of predictions $f_{\theta}(X)$ in group $s$ having a value lower or equal to $\eta\in[0,1]$, i.e. $H_{s}(\eta)=n_{s}^{-1}\sum_{i=1}^{n}\mathbbm{1}_{f_{\theta}(x_{i})<\eta\textrm{ and }s_{i}=s}$, where $n_{s}$ is the number of observations in group $s$. As illustrated in Fig. 1, the $H_{s}$ are pre-computed on discrete grids of values:

[TABLE]

where $\Delta_{\eta}=J_{\eta}^{-1}(\max_{i}(f_{\theta}(x_{i}))-\min_{i}(f_{\theta}(x_{i})))$ and $J_{\eta}$ is the number of discretization steps. We denote $H_{s}^{j}=H_{s}(\eta^{j})$. The regularization function Eq. (4) is then estimated by using the empirical distributions

[TABLE]

where $H_{s}^{-1}$ is the inverse of the empirical function $H_{s}$, i.e. $H_{s}^{-1}(\tau)$ is the $\tau$’th quantile of the observed values $f_{\theta}(x_{i})$ for $s_{i}=s$. As the functions $H_{s}(\eta)$ are discontinuous and are only known for the values $\eta^{j}$ of Eq. (10), they are inverted using linear interpolations of the $H_{s}^{j}$. The integration of Eq. (11) is finally performed on a discrete grid of $J_{\tau}$ values $\tau$ with a discretization step-size $\Delta_{\tau}=1/J_{\tau}$. We denote $W_{2}^{2}(\mu^{n}_{\theta,0},\mu^{n}_{\theta,1})$ this approximation.

In this discrete context, the derivative of $H_{s}$ will be approximated by $\left(H_{s}^{j_{i}+1}-H_{s}^{j_{i}}\right)/\Delta_{\tau}$ using finite differences, where $s_{i}=s$ and $j_{i}$ is defined such that $\eta^{j_{i}}\leq f_{\theta}(x_{i})<\eta^{j_{i}+1}$. Interestingly, this discrete derivative will never be equal to zero for a value of $f_{\theta}(x_{i})$ as long as this prediction was used to compute the $H_{s}^{j}$. As illustrated Fig. 2, we also denote $cor_{s}(f_{\theta}(x_{i}))$ the value corresponding to $f_{\theta}(x_{i})$ in the group $S=|1-s|$ when integrating Eq. (11), i.e. $cor_{s}(f_{\theta}(x_{i}))=H_{s}^{-1}(H_{|1-s|}(f_{\theta}(x_{i})))$.

Then the empirical counterpart of Eq. (9) is

[TABLE]

for $s_{i}\in\{0,1\}$. As a consequence, the variation of $\mu_{\theta}\mapsto W_{2}^{2}(\mu_{\theta,0}^{n},\mu_{\theta,1}^{n})$ due to the prediction $f_{\theta}(x_{i})$, with $i\in B$, can be back-propagated using the expression

[TABLE]

where we added the terms $n_{0}$ and $n_{1}$ to Eq. (12) with $s_{i}=0$ and $s_{i}=1$ in order to have a balanced influence of the classes $S_{0}$ and $S_{1}$ when computing the average gradients in a mini-batch.

3.4.5 Algorithmic pertinence of Eq. (13) on large training sets

For a group $s$, we denoted the discrete cumulative distribution $H_{s}=\{H_{s}^{j}\}_{j=1,\ldots,J}$ such as $H_{s}^{j}$ is equal to $n_{s}^{-1}\sum_{i=1}^{n}\mathbbm{1}_{f_{\theta}(x_{i})<\eta^{j}\textrm{ and }s_{i}=s}$, where $n_{s}$ is the number of observations in group $s$. These quantiles can be computed once for all before estimating all derivatives in a batch $B$. Computing $H_{0}$ and $H_{1}$ using all observations has an algorithmic cost $o(n)$. They can also be computed on reasonably large random subsamples of $n_{sub}$ observations if $n$ is extremely large, so the algorithmic cost can be reduced to $o(n_{sub})$. This makes this pre-computation tractable in all cases. For a given observation $i\in B$, we also denoted $j_{i}$ the index such that $\eta^{j_{i}}\leq f_{\theta}(x_{i})<\eta^{j_{i}+1}$. This index can be found for a very reasonable algorithmic cost of $o(\log_{2}(J))$ using a divide and conquer approach. Finally that the $cor_{s}$ can be computed using linear interpolation, which is again algorithmically reasonable.

As a consequence, Eq. (13) can be efficiently computed in a batch as the $\eta^{j}$ are known, the $H_{{s_{i}}}^{j}$ are pre-computed once for all in each batch, the search for the $j_{i}$ is algorithmically reasonable, and the linear interpolations to compute the $cor_{s}(f_{\theta}(x_{i}))$ is algorithmically cheap.

Remark that in order to make fine and stable the empirical estimation of $\mu_{\theta,0}$ and $\mu_{\theta,1}$, more observations than only those of $B$, but still far less than all $n$ observations, can be used. For instance, in the results obtained Section 4.3, the batches contained 2000 observations and the empirical output distributions were computed using 8000 additional observations. Importantly, it was not required to store any information for the automatic differentiation when computing these additional observations (specifically, the torch.no_grad mode was used in PyTorch). This strongly limited the memory footprint of these predictions.

3.5 Favoring similar error rates

3.5.1 Introduction

An important variant of Sections 3.3 and 3.4 consists in favoring decision rules leading to similar true positive and/or true negative rates for $S=1$ and $S=0$, and not only a similar portion of true decisions. This is related to the notion of equality of opportunities [13], which was discussed in Section 2.1.

In order to favor similar true positive predictions only, the method of Section 3.4 can be simply extended by only computing the cumulative distributions on the observations for which the true prediction $y_{i}$ is $1$. The same idea holds for the true negative predictions by using the observations for which $y_{i}=0$. Favoring a similar error rate in the groups $S=0$ and $S=1$ is slightly more complex as described hereafter.

We first use a non-binary definition of equalized odds conditions (presented in Section 2.1), as in Eq. (1), and merge these two conditions on the observations $Y=1$ and $Y=0$. To do so, we will use the following index denoted by the Disparate Mean Squared Error (DMSE):

[TABLE]

where the mean squared error of the classifier is empirically defined as:

[TABLE]

where $n_{s}=\sum_{i=1}^{n}\mathds{1}_{s_{i}=s}$ and we recall that $g_{\theta}(X_{i})=\mathds{1}_{f_{\theta}(X_{i})>0.5}$. It is straightforward to show that this index is equal to $1$ if the equalized odds conditions are respected and that the more disproportionate the error rate between $S=0$ and $S=1$ the closer this index to [math].

Motivated by the same optimization concerns as in Section 3.3, we first propose to use a slightly modified definition of $MSE$ which is continuous with respect to the parameters $\theta$

[TABLE]

and to use this definition in the DMSE index, Eq. (14). Our main contribution is then to propose to minimize this index by using the Wasserstein-2 distance between the densities of squared errors obtained in the groups $S=0$ and $S=1$.

3.5.2 Back-propagating the impact of an observation

We denote by $\tilde{\mu}_{\theta,0}$ and $\tilde{\mu}_{\theta,0}$ the densities of the squared error rates $(f_{\theta}(X)-Y)^{2}$ in the groups $S=0$ and $S=1$, respectively, and $\tilde{\mu}_{\theta,0}^{n}$ and $\tilde{\mu}_{\theta,0}^{n}$ their discrete counterparts. For a group $s$, we then denote the discrete cumulative distribution $\tilde{H}_{s}=\{\tilde{H}_{s}^{j}\}_{j=1,\ldots,J}$ such as

[TABLE]

where $n_{s}$ is the number of observations in group $s$, and the $\{\eta^{j}\}_{j=1,\ldots,J}$ are regularly sampled thresholds between [math] and $1$ (as in Section 3.4). Instead of solving Eq. (5), we solve here

[TABLE]

We develop in Appendix B how to back-propagate efficiently the impact of an observation $f_{\theta}(x_{i})$ on the regularization term $W_{2}^{2}(\tilde{\mu}^{n}_{\theta,0},\tilde{\mu}^{n}_{\theta,1})$. This leads to the fact that a variation of $\tilde{\mu}_{\theta}\mapsto W_{2}^{2}(\tilde{\mu}_{\theta,0},\tilde{\mu}_{\theta,1})$ modeling this impact can be back-propagated using the expression

[TABLE]

where the $cor_{s}$ have the same meaning as in Eq. (13) but are computed on the $\tilde{H}_{s}$ instead of the $H_{s}$. Computing this equation has a very similar algorithmic cost as the estimation step of Eq. (13). The Wasserstein-2 regularization is therefore computationally reasonable on the squared error. Remark that other extensions of our method are possible. In particular, we extend our method in Section C to the Wasserstein-1 penalty measure and the Logistic regression case.

3.6 Training procedure

Our batch training procedure is summarized Alg. 1 for the disparate impact case, i.e. to keep similar output distributions $Y$ in the groups $S=0$ and $S=1$. The procedure is the same for the squared error case, except that the $\{H_{s},\mu_{\theta,s}\}$ are replaced by the $\{\tilde{H}_{s},\tilde{\mu}_{\theta,s}\}$ (see Eq. (17)) and Eq. (13) is replaced by Eq. (3.5.2). It is important to remark that we implemented this training procedure in PyTorch333https://pytorch.org/ by only writing a specific autograd.Function for our regularized loss term.

4 Results

We assess in this section different aspects of the proposed method. In Section 4.1, we compare it to two other methods on the Adult Census dataset. Although this is a tabular dataset, it has become the gold standard dataset to assess the level of fairness of new classification strategies. This justifies its use in the beginning of the results section. We then evaluate in Section 4.2 the influence of the weight $\lambda$ on the other standard MNIST image dataset. In order to discuss in depth the results of this section, we simulated two kinds of bias in the training dataset with a controlled level of bias. We finally compare in Section 4.3 the two proposed regularization alternatives on the large CelebA image dataset. In this section, we denote Reg. Prediction the strategy of Section 3.3 which favor similar predictions in $S=0$ and $S=1$, and we denote Reg. Error the strategy of Section 3.5 which favor similar predictions errors in the two groups.

4.1 Adult census dataset

4.1.1 Dataset

In order compare the proposed method with different alternatives, we used the Adult Census dataset444https://archive.ics.uci.edu/ml/datasets/adult. It contains $n=45222$ subjects and $p=14$ input variables. The binary output variable $Y$ indicates whether a subject’s incomes is above (positive class, so $Y=1$) or below (negative class, so $Y=0$) 50K USD. We also consider variable Gender as sensitive. Note that the authors have extensively studied this dataset in [10]. After discussing that the training set clearly contains more males than females with $Y=1$, they have made clear that naive bias correction techniques are inefficient to train decision rules leading to similar rates of predictions with $\hat{Y}=1$. This is because of the strong influence of confounding variables, which corresponds to the Cause 3 of algorithmic bias in Section 1.1. Note that as discussed in [10], the training dataset also contains about two times more males than females and that the males have more frequently $Y=1$ than the females, leading potentially to the Cause 2 of bias.

4.1.2 Experimental protocol

The goal in these tests is to reach as much as possible a good balance between the prediction accuracy and the statistical parity of positive outputs w.r.t the gender. We then denote NNrW and NN neural-network based classification strategies with and without the regularization method of Section 3.3. We first compared these neural-network classifiers to a logistic regression classifier with Wasserstein-1 regularization (LRrW), as in [24]. Specifically, we mimicked [24] by using the technique of Appendix C.2 with the regularization of Appendix C.1. We also tested different variants of the reference method of [19], with constraints that explicitly favor high discriminate impacts (ZFA) or low rates of false negative predictions (ZFN). We finally used the non-regularized Logistic-Regression of Scikit-Learn555https://scikit-learn.org/stable/modules/linear_model.html (LR), which is the baseline method of [19].

Default parameters were used for the reference methods of LR, ZFA and ZFN and 300 iterations were used for the gradient descents LRrW. The strategies NNrW, and NN use the same elementary network architecture. It is made of three fully-connected hidden layers with Rectified Linear Units (ReLU) activation functions. A sigmoid activation function is then used in the output layer to allow binary classification. Optimization was made using Adam method [6] with the default parameters on PyTorch. An amount of 100 epochs and a batch size of 50 observations were used. The method of Appendix C.3 was also used to tune $\lambda$ in NNrW. Finally, similar pre-treatment as those of [10] were made on the data to transform the categorical variables into quantitative variables.

All scores are given on test data after randomly splitting all available data into training and test datasets, with 75% and 25% of the observations each, respectively. To measure the classification quality and fairness we considered the accuracy of the predictions (formally defined Section 2.1), the disparate impact Eq. (1) and the empirical probabilities Eq. (2), which were derived from the equalized odds principle of [13] as $O_{s,y}=\mathbb{P}(\hat{Y}=1|S=s,Y=y)$.

4.1.3 Results

Results are given in Table 1. The results of LRrW are first similar to those of ZFA and can be considered as the most fair of all, both from the statistical parity ($DI\approx 1$) and the equalized odds ($O_{1,0}\approx O_{0,0}$ and $O_{1,1}\approx O_{0,1}$) perspectives. Their accuracy is however clearly lower than $0.8$, so they cannot be considered as reasonably accurate. The baseline methods LR and NN, which have no fairness constraints, are now the most accurate ones with $Acc>0.82$. Their DI is however particularly low as they favor much more false positive predictions in $S=1$ than $S=0$ i.e. $O_{1,0}\gg O_{0,0}$. The methods leading to a good balance between accuracy and fairness are finally ZFN and NNrW. Although ZFA predicts well the negative outputs ($O_{1,0}$ and $O_{0,0}$ lower than $0.11$) it does not capture properly how to predict true positive outputs ($O_{1,1}$ and $O_{0,1}$ lower than $0.44$, meaning more than 56% error). This can be explained by the fact that only 24% of the training observations have positive outputs. For a similar accuracy and DI, the strategy NNrW appears as more balanced when comparing to the true and false positive rates. The true positive predictions ($O_{1,1}$ and $O_{0,1}$) have less than 17% of errors and the false positive predictions ($1-O_{1,0}$ and $1-O_{0,0}$) have less than 26% of errors. The NNrW therefore appears as the one leading to the best balance between fairness and accuracy.

4.2 MNIST dataset

4.2.1 Unbiased data

We then tested our method on a semi-synthetic image dataset, where we explicitly controlled the bias in the training data. To do so, we used 20000 training images and 8000 test images out of the classic MNIST dataset666http://yann.lecun.com/exdb/mnist/. Each image $x_{i}$ has $28\times 28$ pixels and represents a handwritten digit. We considered in our experimental protocol that the images $x_{i}$ representing digits strictly lower than 5 (i.e. 0, 1, 2, 3 or 4) have an output $y_{i}=0$ and the other images (i.e. representing 5, 6, 7, 8 or 9) have an output $y_{i}=1$. Note that the frequency of each digit was relatively stable in these datasets. Each digit was indeed represented between 1889 and 2232 times in the training set and between 757 and 905 times in the test set. We then randomly draw in the training set and the test set a label $s_{i}\in\{0,1\}$ for each observation. A Bernoulli law with $p=0.5$ was used to draw $s_{i}$, so there is roughly the same amount of observations with $s_{i}=0$ and $s_{i}=1$. All observations with $s_{i}=0$ were then rotated by 180 degrees.

In order to automatically predict the outputs $y_{i}$ using the images $x_{i}$, we used a basic Convolutional Neural Network (CNN) architecture with three stacks of Convolution/MaxPooling layers. We then trained this classifier using 50 epochs and a batch size of 200 observations. The results were standard for this dataset with $99.8\%$ and $97.4\%$ good predictions on the training and test sets, respectively. On the test set, the true positive rate (TP) was $0.963$ and the false positive rate (FP) was $0.0144$.

4.2.2 Adding a bias in the training data

We now impaired the training data in the group $S=0$ only in order to control how the trained classifiers will be biased. The test set remained unchanged. Two strategies were used to do so: Semi-Random (SR) We randomly set $y_{i}=0$ to $65\%$ of the observations representing the digit 7 in the group $S=0$, Separate-Treatment (ST) We first coarsely trained our classifier on the training set of Section 4.2.1 with 2 epochs and we then set $y_{i}=0$ to the $30\%$ of observations with $s_{i}=0$ and originally $y_{i}=1$ that had the lowest predictions. Remark that about $60\%$ of the transformed observations were true positives, so among the digits higher or equal to 5 in the group $S=0$ of the training set, the $18\%$ being considered as the most similar to 0, 1, 2, 3 or 4 were mislabelled.

In can be remarked that the two kind of biases simulated here corresponds to the Cause 1 of algorithmic bias in Section 1.1. As our goal is to correct a disparate impact in two groups of images which are supposed to have the same distribution of output predictions, we use the regularization energy of Section 3.3. The regularization therefore favors similar predictions.

4.2.3 Results

The main results obtained using the experimental protocols SR and ST are presented in Fig. 3. In both cases, it is clear that increasing the influence of the regularization term (i.e. $\lambda$) leads to gradually increased DIs, until this influence is too high ($\lambda>2.0e-5$ in our tests). In this case, the trained decision rules start to significantly loose their predictive power, so they are useless. We can also remark that the portion of good predictions in the groups $S=0$ and $S=1$ are particularly similar when the DI is close to $1$. This is particularly true when using the protocol protocol SR.

Our nicest results here is that the gain of fairness we obtained by using regularized predictions came with an improved predictive power for SR and no loss of predictive power for ST. The regularization then improved the generalization properties of the trained neural networks in these experiments.

We further explain this phenomenon by discussing the detailed results of Table 2. It represents the DI, and the portion of good predictions in the groups $S=0$ and $S=1$ (GP0 and GP1) obtained on the training and the test sets using two training strategies: experimental protocol SR with $\lambda=0$ (no regularization) and $\lambda=2.0e-5$ (good level of regularization according to Fig. 3). True positive (TP) and true negative (TN) rates are also given on the test set in the groups $S=0$ and $S=1$.

As expected, the trained neural network with no regularization is very accurate on the training set and its DI reflects the simulated bias in the training set. This DI is stable when generalizing on the test set but the predictions are less accurate. When using our regularization strategy, the DI becomes very close to 1 on the training set and the predictions remain reasonably accurate. Contrary to the non-regularized cases, both the prediction accuracy and the DI are however stable when generalizing to the test set. This leads to a similar accuracy as when using no regularization but a clearly improved DI. Note that the gain of fairness is mainly due to an improved true positive rate in the group $S=0$.

4.3 CelebA dataset

4.3.1 Experimental protocol

We now present more advanced results obtained on the Large-scale CelebFaces Attributes (CelebA) Dataset777http://mmlab.ie.cuhk.edu.hk/projects/CelebA.html [49]. It contains more than 200K celebrity images, each with 40 binary attribute annotations. Contrary to what we studied in Section 4.2, the images of the CelebA dataset cover large pose variations and background clutter. This makes this dataset far more complex to study than the MNIST dataset. The binary attributes are for instance Eyeglasses, PaleSkin, Smiling, Young, Male and Attractive.

Interestingly, it is relatively simple to train a classifier detecting some objective attributes, such as Eyeglasses, using a modern Neural Network (NN) architecture. For instance, we reached in the test set more than $98\%$ accuracy when detecting whether the represented celebrities weared eyeglasses. Other attributes, such as Attractive, are more complex to handle because they are subjective. As shown later in our results, it is however relatively simple to have more than $84\%$ of good predictions when predicting whether a person in the test set is considered as attractive. In practice, this suggests that the persons who labeled the data were relatively coherent when choosing who was attractive or not. The main issue with the fact that Attractive is subjective however comes to the fact that it may be influenced by other undesirable attributes. In this section, we then use our regularization strategy to limit the impact of the attribute Young when predicting whether someone is considered as attractive or not, while preserving as much as possible a good prediction accuracy.

Remark that the algorithmic bias we will tackle in this section may be due to different causes identified in Section 1.1. It is obviously related to the Cause 1 because of the subjectivity of the labeling process. It is also likely to be due to the confounding variables of Cause 3, at least in hidden latent spaces. For instance, wearing eyeglasses is less frequent for young persons than older ones, making potentially the young persons with eyeglasses assimilated to older persons. In a similar vein, there are about two $3.5$ times more persons being considered as young than older ones, leading potentially to the Cause 2 of bias. By using this dataset, which is representative of what can be used in industrial applications, we have then high chances to learn biased decision rules.

We used the ResNet-18 convolutional Neural Network (NN) architecture [4] to predict the attribute Attractive based on the CelebA images. The sensitive variable $S$ used for the regularization was also the attribute Young. Note that the original NN architecture was unchanged when implementing our strategy, as our regularization method is fully encoded in the loss function. We used to two proposed regularization strategies to favor similar predictions (see Section 3.3) and similar prediction errors (see Section 3.5) in the groups $S=1$ (young) and $S=0$ (not young). In both cases, different classifiers were trained by using $\lambda=0$ (no regularization) and gradually increasing $\lambda$s with values ranging between [math] and $1e-2$.

Additional results giving the convergence of the training algorithm for different values of $\lambda$ are also presented in Appendix D.

4.3.2 Results

The main results obtained on the CelebA dataset are presented in Fig. 4. As in Fig. 3, it represents the disparate impacts (DI) the mean squared errors (MSE) and the portion of good predictions for the images in the groups $S=0$ and $S=1$ (GP0 and GP1) obtained on the test set. To efficiently compare the gain of fairness with the loss of prediction accuracy, we also present a point cloud representing the MSE and the minimal value between $GP0/GP1$ and $GP1/GP0$ of each trained classifier.

The global behavior of our strategy is similar to what we observed in Section 4.2 on the MNIST dataset. When increasing $\lambda$ until about $5.5e-3$, the results are more and more fair, i.e. the DI is higher and higher, and GP1 is closer and closer to GPO. It appears that a DI of 0.6 can be reached with little loss of predictive power in this test, and that $GP0$ can be very close to $GP1$ under this constraint. It can also be remarked that for $\lambda<5.5e-3$, favoring similar predictions (see Section 3.3) or similar prediction errors (see Section 3.5) in the groups $S=0$ and $S=1$ had a relatively similar effect. When using the method of Section 3.3 with $\lambda>5e-3$, the predictions are however unstable: They are sometimes similar to those obtained using $\lambda\approx 4.5e-3$, sometimes very fair ($DI>0.8$) but inaccurate ($MSE>0.2$), and sometimes clearly unfair ($DI\approx 0.3$) and inaccurate ($MSE>0.2$). When using the method of Section 3.5 with $\lambda>6e-3$, the predictions first appear as similar to those obtained without regularization. Then, they become more unstable for $\lambda>8.5e-3$, as for the results represented with $\lambda>8.5e-3$. Another point to emphasize is that although a DI of about 0.6 can be reached with little loss of predictive power, DIs above 0.8 can only be reached with strong loss of predictive power. This is interestingly not the case for the difference between the portion of good predictions GP0 and GP1 (which takes into account the true prediction $\hat{Y}$). The ratio $GP0/GP1$ is indeed very close to 1 when favoring similar prediction errors with $\lambda>5e-3$. Importantly, the regularization finally allowed to get clearly more fair results with little loss of predictive power for weights $\lambda\in[3.5e-3,5.5e-3]$.

It is now interesting to discuss what was the practical impact of the regularization algorithm in these tests. We recall that $Y$ represents the subjective attribute Attractive, which represents who is attractive for the persons who labelled the data and that $S$ is the attribute Young. We then measured the scores differences $f_{\bar{\theta}}(x_{i})-f_{\tilde{\theta}}(x_{i})$, where the parameters $\bar{\theta}$ were obtained by favoring similar error rates with $\lambda=4.e-3$ and the parameters $\tilde{\theta}$ were obtained without regularization. The three highest positive and negative differences between $f_{\bar{\theta}}(x_{i})$ and $f_{\tilde{\theta}}(x_{i})$ are represented in Fig. 5 This figure makes clear that the regularization with $S$ encoding the attribute Young made more attractive the persons with eyeglasses and relatively dark tans, and made less attractive the women with very pale skins and no eyeglasses. In addition to the fact that the young images are positively correlated with attractive and female in the training set, this can be explained by different reasons: The frequency of persons with eyeglasses is 4.71 times more frequent in old persons than young ones; Males are 2.01 time more frequent in old persons than young ones; Pale skins are 1.8 times more frequent in young persons than old ones. In all cases, this suggests that the regularization has lead to coherent results with regard to the distribution of the attributes in the training set.

5 Discussion

In this paper, we have proposed a new method to temper the algorithmic bias in Neural-Network based classifiers. The first key advantage of this strategy is that it can be integrated to any kind of Neural-Network architecture, as it only overloads the loss term when training optimal decision rules. This additionally makes it straightforward to integrate it to existing deep learning solution. As demonstrated on the CelebA dataset, its second key advantage is that it scales particularly well to large image datasets, which are increasingly ubiquitous in industrial applications of artificial intelligence. In terms of methodology, the central idea of this work was to define two alternatives of a fairness penalty term which can be naturally used in stochastic gradient descent based optimization algorithms. The first alternative favors similar prediction outputs in two groups of data and the second one favors similar prediction errors in these groups. The main technical lock we had to address was to approximate the gradient of this penalty term in the real-life context where the size of mini-batches can be relatively small. Our results have shown the good properties our regularization strategy.

Future work will first consist in extending our strategy to the non-binary classification case, which should be straightforward by using one-hot-encoding representations of the outputs. A more ambitious extension, with potentially numerous applications, would be to address the regression case too. We may also explore the use of this strategy on different kinds of data, such as those treated in Natural Language Processing, as it can be used on any kind of data. In order to make the optimization process more stable when the algorithmic bias is particularly complex, a promising strategy would be finally to work on the latent spaces of the neural-network hidden layers and not the outputs directly.

Appendix A Proof of Gâteaux Differentiability

A.1 Proof of Proposition 1

Proof 1

Recall from the duality of the transport problem that it can be rewritten as the optimization problem

[TABLE]

where the optimization is over the set of convex functions. Following [50], we then denote here $f$ a convex function and $f^{*}$ its conjugate. Let $f_{0}$ and $f_{t}$ be the solutions for $\mathcal{W}^{2}_{2}(\mu,\nu)$ and $\mathcal{W}^{2}_{2}(\mu+t\alpha,\nu+t\beta)$. We then have

[TABLE]

and

[TABLE]

Note that, following [51, 52], the convergence $\mu+t\alpha\xrightarrow{w}\mu$ implies that $f_{t}\rightarrow f_{0}$ uniformly on the compact sets of the support of $\mu$ if the support is connected. As a consequence, if $\operatorname{Supp}(\alpha)\subset\operatorname{Supp}(\mu)$ and $\operatorname{Supp}(\beta)\subset\operatorname{Supp}(\nu)$, then under the assumption of compact supports we have

[TABLE]

Recall from [53] (Section 2.2.) that in the real line case the potential of the transport, solution of Eq. (20), is a primitive of the transport map, which is the map $F^{-1}_{\mu}(F_{\nu})$. Then the conclusion becomes straightforward.

A.2 Proof of Lemma 1

Proof 2

Set $t>0$ and $\beta$, then compute

[TABLE]

In consequence we have that

[TABLE]

Appendix B Extension to error rates Eq. (3.5.2)

By using the same strategy as in Proposition 1 with $(f_{\theta}(x_{i})-y_{i})^{2}$ and the $\tilde{H}_{g}$ instead of $f_{\theta}(x_{i})$ and the $H_{g}$, we can compute:

[TABLE]

.

This equation specifically expresses how to back-propagate the impact of a squared error $(f_{\theta}(x_{i})-y_{i})^{2}$, and not an output observation $f_{\theta}(x_{i})$, on $W_{2}^{2}(\tilde{\mu}^{n}_{\theta,0},\tilde{\mu}^{n}_{\theta,1})$. From Lemma 1, we know that the chain rule applies on the Gâteaux differentiability which is used to back-propagate the impact of $f_{\theta}(x_{i})$ on the distributions. As

[TABLE]

we can then simply deduce that the impact of an observation $f_{\theta}(x_{i})$ on $W_{2}^{2}(\tilde{\mu}^{n}_{\theta,0},\tilde{\mu}^{n}_{\theta,1})$ can be back-propagated using Eq. (3.5.2):

[TABLE]

Appendix C Extensions of the method of Section 3.4

C.1 Wasserstein-1 distances

Our approach can be straightforwardly extended to approximate Wasserstein-1 distances. By using the same reasoning as in Section 3.4, it can be shown that the impact of an observation $f_{\theta}(x_{i})$ on $W_{1}(\tilde{\mu}^{n}_{\theta,0},\tilde{\mu}^{n}_{\theta,1})$ can be back-propagated using

[TABLE]

instead of Eq. (13), where $sign(x)$ is equal to $+1$ or $-1$ depending on the sign of $x$. We emphasize that the distances between the cumulative densities are therefore not taken into account when computing the gradients of the Wasserstein-1, although this is the case for Wasserstein-2 distances.

C.2 Logistic Regression

We now show how to simply implement our regularization model for Logistic Regression. We minimize:

[TABLE]

where $f_{\theta}(x_{i})=(1+\exp{(-\theta^{0}-\theta^{\prime}x_{i})})^{-1}$ is the logistic function and $\theta=(\theta^{1},\ldots,\theta^{p})^{\prime}$ is a vector in $\mathbb{R}^{p}$ representing the weights given to each dimension. The derivatives of the whole energy Eq. (C.2) with respect to each $\theta^{j}$, $j=0,\ldots,p$, can be directly computed using finite differences here.

We emphasize that a fundamental difference between using our Wasserstein based regularization model in Section 3.4 and here is that $p$ derivatives of the minimized energy are approximated using Logistic Regression (derivatives w.r.t. the $\theta^{j}$, $j=0,\ldots,p$), while $n$ derivatives are required when using Neural Networks with a standard gradient descent (derivatives w.r.t. the $f_{\theta}(x_{i})$, $i=0,\ldots,n$). As a cumulative histogram is computed each time the derivative of a Wasserstein-2 distance is approximated, this task can be bottleneck for common Neural-Networks applications where $n$ is large. This fully justifies the proposed batch-training regularization strategy of Section 3.4.

C.3 Automatic tuning of $\lambda$

The minimized energy Eq. (5) depends on a weight $\lambda$ which balances the influence of the regularization term $W_{2}^{2}(\mu^{n}_{\theta,0},\mu^{n}_{\theta,1})$ with respect to the data attachment term $R(\theta)$. A simple way to automatically tune $\lambda$ is the following. Compute the average derivatives of $W_{2}$ and $R$ after a predefined warming phase of several epochs, where $\lambda=0$. We denote $d_{W_{2}}$ and $d_{R}$ these values. Then tune $\lambda$ as equal to $\alpha\frac{g_{R}}{g_{W_{2}}}$, where $\alpha$ is typically in $[0.1,1]$. This makes it intuitive to tune the scale of $\lambda$.

In the disparate impact case (Section 3.4), it can be interesting to accurately adapt $\alpha$ to the machine learning problem, in order to finely tune $\lambda$ with regards to the fact that we simultaneously want fair and accurate predictions. Inspired by the hard constraints of [19] to enforce fair predictions, we update $\alpha$ based on measures of the Disparate Impact (DI), Eq. (1), and average Prediction Accuracy (Acc) at the beginning of each epoch. Remark that lowering the Wasserstein-2 distance between the predictions $f_{\theta}(x_{i})$ in groups [math] and $1$ naturally tend to make decrease $\mathbbm{1}_{f_{\theta}(x_{i},s_{i}=0)>0.5}-\mathbbm{1}_{f_{\theta}(x_{i},s_{i}=1)>0.5}$, which we empirically verified. The disparate impact therefore tends to be improved. We believe that hard constraints based on other fairness measures could also be used. Establishing a clear relation of causality between the Wasserstein-2 distance and different fairness measures is however out of the scope of this paper and hence considered as future work. Note that the same technique holds in the squared error case (Section 3.5) but the Disparate Mean Squared Error (DMSE) index, Eq. (14), is used instead of the DI.

In the experiments of Section 4.1, our hard constraints are for instance: If the prediction accuracy is too low (Acc$<0.75$), then $\alpha$ is slightly decreased to favor the predictions accuracy ($\alpha=0.9\alpha$). If the prediction accuracy is sufficiently high and the DI is too low (DI$<0.85$) then $\alpha$ is slightly increased ($\alpha=1.1\alpha$) to favor fair decisions. We empirically verified in our experiments that $\alpha$ converges to satisfactory values using this method, if the classifier is able to learn classification rules leading to sufficiently high PA. Parameter $\alpha$ converges to zero otherwise.

Appendix D Impact of $\lambda$ on the convergence

D.1 Results

In this section, we extend the results of Section 6 by discussing the convergence of the training algorithm on the training set and the test set for different values of $\lambda$. Each sub-figure specifically represents the accuracy of the prediction model mini-batch after mini-batch, where we distinguish the average accuracy obtained on observations with $S=0$ and those with $S=1$. The disparate impact is also given. We additionally compare the convergence of the training algorithm without regularization when using a mean squared error (MSE) loss, as in other experiments, and a binary cross entropy (BCE) loss, which is known to be less sensitive than MSE to outlier observations. The regularization was applied on the output predictions directly. Results are detailed in Fig. 6.

D.2 Discussion

We can first notice that the convergence properties of the training algorithm with $\lambda=0$ are very similar when using the MSE and the BCE loss. The potential outlier observations seem therefore to have no or little impact here. In both cases and both for $S=0$ and $S=1$, the observed level of accuracy is similar for the training set and the test set until about the 30th iteration. Then, the accuracy becomes clearly higher on the training set than on the test set, which means that the trained model overfits the training data. This phenomenon is particularly strong for the observations in the group $S=1$, which have a lower accuracy than in the group $S=0$ before overfitting. In all cases, the disparate impact is close to $0.25$ at convergence.

Now, when increasing $\lambda$ from $1.e-3$ to $5.e-3$ the level of accuracy becomes increasingly similar between $S=0$ and $S=1$, and the disparate impact is closer and closer to $0.4$. It can be remarked that for $\lambda=4.e-3$ and $\lambda=5.e-3$, the curves start oscillating when the training algorithm overfits the training data. The regularization therefore seems to be strong compared with the data attachment term. This phenomenon is stronger, even before overfitting, for $\lambda=6.e-3$. This explains why the results become unstable in Fig. 4-(top) for high values of $\lambda$, and confirms that a reasonable trade-off should be searched between regularization and prediction accuracy when tuning $\lambda$.