Fairness for Robust Log Loss Classification

TL;DR

This paper introduces a new fairness-aware classification method based on distributional robustness, formulated as a minimax game, which offers theoretical and practical improvements over existing fairness constraints.

Contribution

The authors derive a novel classifier from first principles of distributional robustness that integrates fairness into a worst-case logarithmic loss minimization framework.

Findings

Convexity and asymptotic convergence of the proposed method.

Improved fairness and accuracy on benchmark datasets.

Produces a parametric exponential family distribution resembling truncated logistic regression.

Abstract

Developing classification methods with high accuracy that also avoid unfair treatment of different groups has become increasingly important for data-driven decision making in social applications. Many existing methods enforce fairness constraints on a selected classifier (e.g., logistic regression) by directly forming constrained optimizations. We instead re-derive a new classifier from the first principles of distributional robustness that incorporates fairness criteria into a worst-case logarithmic loss minimization. This construction takes the form of a minimax game and produces a parametric exponential family conditional distribution that resembles truncated logistic regression. We present the theoretical benefits of our approach in terms of its convexity and asymptotic convergence. We then demonstrate the practical advantages of our approach on three benchmark fairness datasets.

Figures

Equations

$$\displaystyle\Gamma:\Big{\{}\mathbb{P}\mid$$

$$\displaystyle=\tfrac{1}{p_{\gamma_{0}}}\mathbb{E}_{\begin{subarray}{c}\widetilde{P}({\bf x},a,y)\\ \mathbb{P}(\widehat{y}|{\bf x},a,y)\end{subarray}}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{0}(A,Y))]\Big{\}},$$

$${\Gamma }_{\text{dp}}⟺{\gamma }_j(A,Y)=\mathbb{I}(A=j);$$

$${\Gamma }_{\text{e.opp}}⟺{\gamma }_j(A,Y)=\mathbb{I}(A=j\wedge Y=1);$$

$$\displaystyle\Gamma_{\text{e.odd}}\iff\gamma_{j}(A,Y)=\left[\begin{array}[]{c}\mathbb{I}(A=j\wedge Y=1)\\ \mathbb{I}(A=j\wedge Y=0)\end{array}\right].$$

$$\underset{\mathbb{P}(\hat{y}|\mathbf{x})\in \Delta }{\min }\underset{\mathbb{Q}(\hat{y}|\mathbf{x})\in \Delta \cap \Xi }{\max }-\sum _{\mathbf{x},\hat{y}}\tilde{P}(\mathbf{x})\mathbb{Q}(\hat{y}|\mathbf{x})\log \mathbb{P}(\hat{y}|\mathbf{x})$$

$$=\underset{\mathbb{P}(\hat{y}|\mathbf{x})\in \Xi }{\max }-\sum _{\mathbf{x},\hat{y}}\tilde{P}(\mathbf{x})\mathbb{P}(\hat{y}|\mathbf{x})\log \mathbb{P}(\hat{y}|\mathbf{x})=\underset{\mathbb{P}(\hat{y}|\mathbf{x})\in \Xi }{\max }H(\hat{Y}|\mathbf{X}),$$

$$\displaystyle\Xi:\Big{\{}\mathbb{Q}\mid\mathbb{E}_{\widetilde{P}({\bf x});\mathbb{Q}(\widehat{y}|{\bf x})}[\phi({\bf X},\widehat{Y})]=\mathbb{E}_{\widetilde{P}({\bf x},y)}\left[\phi({\bf X},Y)\right]\Big{\}},$$

$$\underset{\mathbb{P}\in \Delta \cap \Gamma }{\min }\underset{\mathbb{Q}\in \Delta \cap \Xi }{\max }$$

$$\underset{\theta }{\min }$$

$$+{\theta }^{\top }\left({\mathbb{E}}_{{\mathbb{Q}}_{\theta ,\lambda }(\hat{y}|\mathbf{x},a,y)}[\varphi (\mathbf{x},\hat{Y})]-\varphi (\mathbf{x},y)\right)$$

$$\displaystyle+\lambda\big{(}\tfrac{1}{p_{\gamma_{1}}}\mathbb{E}_{\mathbb{P}_{\theta,\lambda}(\widehat{y}|{\bf x},a,y)}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{1}(A,Y))]$$

$$\displaystyle-\tfrac{1}{p_{\gamma_{0}}}\mathbb{E}_{\mathbb{P}_{\theta,\lambda}(\widehat{y}|{\bf x},a,y)}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{0}(A,Y))]\big{)}\bigg{\}},$$

$${\mathbb{P}}_{\theta ,\lambda }(\hat{y}=1|\mathbf{x},a,y)=$$

$$\displaystyle\begin{cases}\min\big{\{}e^{\theta^{\top}\phi({\bf x},1)}/Z_{\theta}({\bf x}),\frac{p_{\gamma_{1}}}{\lambda}\big{\}}&\text{if }\gamma_{1}(a,y)\land\lambda>0\\ \max\big{\{}e^{\theta^{\top}\phi({\bf x},1)}/Z_{\theta}({\bf x}),1\!-\!\frac{p_{\gamma_{0}}}{\lambda}\big{\}}&\text{if }\gamma_{0}(a,y)\land\lambda>0\\ \max\big{\{}e^{\theta^{\top}\phi({\bf x},1)}/Z_{\theta}({\bf x}),1\!+\!\frac{p_{\gamma_{1}}}{\lambda}\big{\}}&\text{if }\gamma_{1}(a,y)\land\lambda<0\\ \min\big{\{}e^{\theta^{\top}\phi({\bf x},1)}/Z_{\theta}({\bf x}),-\frac{p_{\gamma_{0}}}{\lambda}\big{\}}&\text{if }\gamma_{0}(a,y)\land\lambda<0\\ e^{\theta^{\top}\phi({\bf x},1)}/Z_{\theta}({\bf x})&\text{otherwise},\end{cases}$$

$${\mathbb{Q}}_{\theta ,\lambda }(\hat{y}=1|\mathbf{x},a,y)={\mathbb{P}}_{\theta ,\lambda }(\hat{y}=1|\mathbf{x},a,y)\times$$

$$\displaystyle\qquad\begin{cases}\big{(}1+\frac{\lambda}{p_{\gamma_{1}}}\mathbb{P}_{\theta,\lambda}(\widehat{y}=0|{\bf x},a,y)\big{)}&\text{if }\gamma_{1}(a,y)\\ \big{(}1-\frac{\lambda}{p_{\gamma_{0}}}\mathbb{P}_{\theta,\lambda}(\widehat{y}=0|{\bf x},a,y)\big{)}&\text{if }\gamma_{0}(a,y)\\ 1&\text{otherwise}.\end{cases}$$

$${\mathbb{Q}}_{\theta ,\lambda }(\hat{y}=1|\mathbf{x},a,y)=\rho (1+\frac{\lambda }{p_{{\gamma }_1}}(1-\rho ))=(1+\frac{\lambda }{p_{{\gamma }_1}})\rho -\frac{\lambda }{p_{{\gamma }_1}}{\rho }^2,$$

$$\underset{\theta }{\min }\frac{1}{n}{\sum }_{(\mathbf{x},a,y)\in \mathcal{D}}{\ell }_{\theta ,{\lambda }_{\theta }^{\ast }}(\mathbf{x},a,y),\text{ where: }$$

$${\ell }_{\theta ,{\lambda }^{\ast }}(\mathbf{x},a,y)=-{\theta }^{\top }\varphi (\mathbf{x},y)+$$

$$\{\begin{array}{ll}-\log (\frac{p_{{\gamma }_1}}{{\lambda }_{\theta }^{\ast }})+{\theta }^{\top }(\varphi (\mathbf{x},1))& \text{if }{\gamma }_1(a,y)\wedge T(\mathbf{x},\theta )\wedge {\lambda }_{\theta }^{\ast }>0\\ -\log (\frac{p_{{\gamma }_0}}{{\lambda }_{\theta }^{\ast }})+{\theta }^{\top }(\varphi (\mathbf{x},0))& \text{if }{\gamma }_0(a,y)\wedge T(\mathbf{x},\theta )\wedge {\lambda }_{\theta }^{\ast }>0\\ -\log (-\frac{p_{{\gamma }_1}}{{\lambda }_{\theta }^{\ast }})+{\theta }^{\top }(\varphi (\mathbf{x},0))& \text{if }{\gamma }_1(a,y)\wedge T(\mathbf{x},\theta )\wedge {\lambda }_{\theta }^{\ast }<0\\ -\log (-\frac{p_{{\gamma }_0}}{{\lambda }_{\theta }^{\ast }})+{\theta }^{\top }(\varphi (\mathbf{x},1))& \text{if }{\gamma }_0(a,y)\wedge T(\mathbf{x},\theta )\wedge {\lambda }_{\theta }^{\ast }<0\\ \log Z_{\theta }(\mathbf{x})& \text{otherwise}.\end{array}$$

$$\mathbb{P}(\hat{y}|\mathbf{x},a)$$

$$+\mathbb{P}(\hat{y}|\mathbf{x},a,y=0)P(y=0|\mathbf{x},a).$$

$$\mathbb{Q}(\hat{y}|\mathbf{x},a)$$

$$+\mathbb{Q}(\hat{y}|\mathbf{x},a,y=0)\mathbb{Q}(\hat{y}=0|\mathbf{x},a).$$

$$\mathbb{Q}(\hat{y}=1|\mathbf{x},a)$$

$$+\mathbb{Q}(\hat{y}=1|\mathbf{x},a,y=0)),$$

$$\underset{\mathbb{P}\in \Delta \cap \Gamma }{\min }\underset{\mathbb{Q}\in \Delta \cap \Xi }{\max }{\mathbb{E}}_{\begin{array}{c}\tilde{P}(\mathbf{x},a,y)\\ \mathbb{Q}(\hat{y}|\mathbf{x},a,y)\end{array}}\left[-\log \mathbb{P}(\hat{Y}|\mathbf{X},A,Y)\right]$$

$$\stackrel{(a)}{=}$$

$$\stackrel{(b)}{=}$$

$$\stackrel{(c)}{=}$$

$$\stackrel{(d)}{=}$$

$$\stackrel{(e)}{=}$$

$$\displaystyle\qquad\qquad+\lambda\Big{(}\tfrac{1}{p_{\gamma_{1}}}\mathbb{E}_{\begin{subarray}{c}\widetilde{P}({\bf x},a,y)\\ \mathbb{P}(\widehat{y}|{\bf x},a,y)\end{subarray}}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{1}(A,Y))]-\tfrac{1}{p_{\gamma_{0}}}\mathbb{E}_{\begin{subarray}{c}\widetilde{P}({\bf x},a,y)\\ \mathbb{P}(\widehat{y}|{\bf x},a,y)\end{subarray}}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{0}(A,Y))]\Big{)}$$

$$\stackrel{(f)}{=}$$

$$\displaystyle\qquad\qquad+\lambda\Big{(}\tfrac{1}{p_{\gamma_{1}}}\mathbb{E}_{\begin{subarray}{c}\widetilde{P}({\bf x},a,y)\\ \mathbb{P}(\widehat{y}|{\bf x},a,y)\end{subarray}}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{1}(A,Y))]-\tfrac{1}{p_{\gamma_{0}}}\mathbb{E}_{\begin{subarray}{c}\widetilde{P}({\bf x},a,y)\\ \mathbb{P}(\widehat{y}|{\bf x},a,y)\end{subarray}}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{0}(A,Y))]\Big{)}$$

$$\stackrel{(g)}{=}$$

$$\displaystyle\qquad\qquad+\lambda\Big{(}\tfrac{1}{p_{\gamma_{1}}}\mathbb{E}_{\mathbb{P}(\widehat{y}|{\bf x},a,y)}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{1}(A,Y))]-\tfrac{1}{p_{\gamma_{0}}}\mathbb{E}_{\mathbb{P}(\widehat{y}|{\bf x},a,y)}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{0}(A,Y))]\Big{)}\Bigg{]}$$

$$\underset{\mathbb{P}\in \Delta }{\min }\underset{\mathbb{Q}\in \Delta }{\max }$$

$$\displaystyle\qquad\qquad+\lambda\Big{(}\tfrac{1}{p_{\gamma_{1}}}\mathbb{E}_{\mathbb{P}(\widehat{y}|{\bf x},a,y)}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{1}(A,Y))]-\tfrac{1}{p_{\gamma_{0}}}\mathbb{E}_{\mathbb{P}(\widehat{y}|{\bf x},a,y)}[\mathbb{I}(\widehat{Y}\!=\!1\wedge\gamma_{0}(A,Y))]\Big{)}\Bigg{]}$$

$$\underset{\mathbb{P}}{\min }\underset{\mathbb{Q}}{\max }\underset{Z_{\mathbb{Q}},\beta \ge 0}{\min }\underset{Z_{\mathbb{P}}}{\max }$$

$$+$$

Full text

Introduction

Though maximizing accuracy has been the principal objective for classification tasks, competing priorities are also often of key concern in practice. Fairness properties that guarantee equivalent treatment to different groups in various ways are a prime example. These may be desirable—or even legally required—when making admissions decisions for universities (?; ?), employment and promotion decisions for organizations (?), medical decisions for hospitals and insurers (?; ?), sentencing guidelines within the judicial system (?; ?), loan decisions for the financial industry (?; ?) and in many other applications. Group fairness criteria generally partition the population based on a protected attribute into groups and mandate equal treatment of members across groups based on some defined statistical measures. We focus on three prevalent group fairness measures in this paper: demographic parity (?), equalized odds, and equalized opportunity (?).

Techniques for constructing predictors with group fairness properties can be categorized into pre-, post-, and in-processing methods. Pre-processing methods use reweighting and relabeling (?; ?) or other transformations of input data (?; ?; ?; ?; ?; ?) to remove unfair dependencies with protected attributes. Post-processing methods adjust the class labels (or label distributions) provided from black box classifiers to satisfy desired fairness criteria (?; ?; ?). In-processing methods integrate fairness criteria into the optimization procedure of the classifier with constraints/penalties (?; ?; ?; ?; ?; ?; ?; ?; ?; ?), meta-algorithms (?; ?), reduction-based methods (?), or generative-adversarial training (?; ?; ?; ?; ?).

Unlike many existing methods that directly form a constrained optimization from base classifiers, we take a step back and re-derive prediction from the underlying formulation of logistic regression. Working from the first principles of distributionally robust estimation (?; ?; ?), we incorporate fairness constraints into the formulation of the predictor. We pose predictor selection as a minimax game between a predictor that is fair on a training sample and a worst-case approximator of the training data labels that maintains some statistical properties of the training sample. Like post-processing methods, our approach reshapes its predictions for each group to satisfy fairness requirements. However, our approach is inherently an in-process method that jointly optimizes this fairness transformation and linear feature-based parameters for an exponential family distribution that can be viewed as truncated logistic regression. Our method assumes group membership attributes are given at training and testing time, which matches many real-world applications. We leave the extension of our approach to settings with inferred group attributes as future work.

Our method reduces to a convex optimization problem with a unique solution for resolving unfairness between groups that asymptotically minimizes the KL divergence from the true distribution. In contrast, many existing methods are susceptible to the local optima of non-convex optimization or to the approximation error from relaxations (?; ?; ?; ?; ?; ?; ?), do not have unique solutions for ensuring fairness (?), or produce mixtures of predictors (?) rather than a single coherent predictor. For fairness criteria that include the true label (e.g., equalized opportunity, equalized odds), we introduce a method for making predictions from label-conditioned distributions and establish desirable asymptotic properties. We demonstrate the practical advantages of our approach compared to existing fair classification methods on benchmark data-driven decision tasks.

Background

Measures of fairness for decision making

Several useful measures have been proposed to quantitatively assess fairness in decision making. Though our approach can be applied to a wider range of fairness constraints, we focus on three prominent ones: Demographic Parity (?), Equality of Opportunity (?) and Equality of Odds (?). These are defined for binary decision settings with examples drawn from a population distribution, $({\bf X},A,Y)\sim P$, with $\widetilde{P}({\bf x},a,y)$ denoting this empirical sample distribution, $\{{\bf x}_{i},a_{i},y_{i}\}_{i=1:n}$. Here, $y=1$ is the “advantaged” class for positive decisions. Each example also possesses a protected attribute $a\in\{0,1\}$ that defines membership in one of two groups. The general decision task is to construct a probabilistic prediction, $\mathbb{P}(\widehat{y}|{\bf x},a)$ over the decision variable $\widehat{y}\in\{0,1\}$ given input ${\bf x}\in\boldsymbol{\mathcal{X}}$ and training data $\widetilde{P}({\bf x},a,y)$. We similarly notate an adversarial conditional distribution that approximates the labels as $\mathbb{Q}$. $\mathbb{P}$ and $\mathbb{Q}$ are the key objects being optimized in our formulation.

Fairness requires treating the different groups equivalently in various ways. Unfortunately, the naïve approach of excluding the protected attribute from the decision function, e.g., restricting to $\mathbb{P}(\widehat{y}|{{\bf x}})$, does not guarantee fairness because the protected attribute $a$ may still be inferred from ${\bf x}$ (?). Instead of imposing constraints on the predictor’s inputs, definitions of fairness require statistical properties on its decisions to hold.

Definition 1**.**

A classifier satisfies Demographic Parity (D.P.) if the output variable $\widehat{Y}$ is statistically independent of the protected attribute $A$: $P(\widehat{Y}\!=\!1|A\!=\!a)=P(\widehat{Y}\!=\!1),\hskip 5.0pt\forall a\in\{0,1\}$.

Definition 2**.**

*A classifier satisfies Equalized Odds (E.Odd.) if the output variable $\widehat{Y}$ is conditionally independent of the protected attribute $A$ given the true label $Y$: $P(\widehat{Y}\!=\!1|A\!=\!a,Y\!=\!y)=P(\widehat{Y}\!=\!1|Y\!=\!y),\hskip 5.0pt\forall y,a\in\{0,1\}.$ *

Definition 3**.**

*A classifier satisfies Equalized Opportunity (E.Opp.) if the output variable $\widehat{Y}$ and protected attribute $A$ are conditionally independent given $Y=1$: $P(\widehat{Y}\!=\!1|A\!=\!a,Y\!=\!1)=\;P(\widehat{Y}\!=\!1|Y\!=\!1),\hskip 5.0pt\forall a\in\{0,1\}$. *

The sets of decision functions $\mathbb{P}$ satisfying these fairness constraints are convex and can be defined using linear constraints (?). The general form for these constraints is:

[TABLE]

where $\gamma_{1}$ and $\gamma_{0}$ denote some combination of group membership and ground-truth class for each example, while $p_{\gamma_{1}}$ and $p_{\gamma_{0}}$ denote the empirical frequencies of $\gamma_{1}$ and $\gamma_{0}$: $p_{\gamma_{i}}=\mathbb{E}_{\widetilde{P}(a,y)}[\gamma_{i}(A,Y)]$. We specify $\gamma_{1}$ and $\gamma_{0}$ in (1) for fairness constraints (Definitions 1, 2, and 3) as:

[TABLE]

Robust log-loss minimization, maximum entropy, and logistic regression

The logarithmic loss, $-\sum_{{\bf x},y}P({\bf x},y)\log\mathbb{P}(y|{\bf x})$, is an information-theoretic measure of the expected amount of “surprise” (in bits for $\log_{2}$) that the predictor, $\mathbb{P}(y|{\bf x})$, experiences when encountering labels $y$ distributed according to $P({\bf x},y)$. Robust minimization of the logarithmic loss serves a fundamental role in constructing exponential probability distributions (e.g., Gaussian, Laplacian, Beta, Gamma, Bernoulli (?)) and predictors (?). For conditional probabilities, it is equivalent to maximizing the conditional entropy (?):

[TABLE]

after simplifications based on the fact that the saddle point solution is $\mathbb{P}=\mathbb{Q}$. When the loss maximizer $\mathbb{Q}$ is constrained to match the statistics of training data (specified using vector-valued feature function $\phi$),

[TABLE]

the robust log loss minimizer/maximum entropy predictor (Eq. (7)) is the logistic regression model, $P(y|{\bf x})\propto e^{\theta^{\mathrm{T}}\phi({\bf x},y)}$, with $\theta$ estimated by maximizing data likelihood (?). While this distribution technically needs to only be defined at input values in which training data exists (i.e., $\widetilde{P}({\bf x})>0$), we employ an inductive assumption that generalizes the form of the distribution to other inputs.

This formulation has been leveraged to provide robust predictions under covariate shift (i.e., difference in training and testing distributions) (?) and for constructing consistent predictors for multiclass classifications (?) and graphical models (?). Our approach similarly extends this fundamental formulation by imposing fairness constraints on $\mathbb{P}$. However, since the fairness constraints and statistic-matching constraints are often not fully compatible (i.e., $\Gamma\not\subseteq\Xi$), the saddle point solution is no longer simple (i.e., $\mathbb{P}\neq\mathbb{Q}$).

Formulation and Algorithms

Given fairness requirements for a predictor (Eq. (1)) and partial knowledge of the population distribution provided by a training sample (Eq. (8)), how should a fair predictor be constructed? Like all inductive reasoning, good performance on a known training sample does not ensure good performance on the unknown population distribution. We take a robust estimation perspective by seeking the best solution for the worst-case population distribution under these constraints.

Robust and fair log loss minimization

We formulate the robust fair predictor’s construction as a minimax game between the predictor and a worst-case approximator of the population distribution. We assume the availability of a set of training samples, $\{({\bf x}_{i},a_{i},y_{i})\}_{i=1:n}$, which we equivalently denote by probability distribution $\widetilde{P}({\bf x},a,y)$.

Definition 4**.**

The Fair Robust Log-Loss Predictor, $\mathbb{P}$, minimizes the worst-case log loss—as chosen by approximator $\mathbb{Q}$ constrained to reflect training statistics (denoted by set $\Xi$ of Eq. (8))—while providing empirical fairness guarantees111$\Delta$ is the set of conditional probability simplexes (i.e., $\mathbb{P}(y|{\bf x},a)\geq 0,\sum_{y^{\prime}}\mathbb{P}(y^{\prime}|{\bf x},a)=1,\forall{\bf x},y,a$). (denoted by set $\Gamma$ of Eq. (1)):

[TABLE]

Though conditioning the decision variable $\widehat{Y}$ on the true label $Y$ would appear to introduce a trivial solution ($\widehat{Y}=Y$), instead, $Y$ only influences $\widehat{Y}$ based on fairness properties due to the robust predictor’s construction. Note that if the fairness constraints do not relate $Y$ and $\widehat{Y}$, the resulting distribution is conditionally independent (i.e., $\mathbb{P}(\widehat{Y}|{\bf X},A,Y=0)=\mathbb{P}(\widehat{Y}|{\bf X},A,Y=1)$), and when all fairness constraints are removed, this formulation reduces to the familiar logistic regression model (?). Conveniently, this saddle point problem is convex-concave in $\mathbb{P}$ and $\mathbb{Q}$ with additional convex constraints ($\Gamma$ and $\Xi$) on each distribution.

Parametric Distribution Form

By leveraging strong minimax duality in the “log-loss game” (?; ?) and strong Lagrangian duality (?), we derive the parametric form of our predictor. ​​222The proofs of Theorem 1 and other theorems in the paper are available in the supplementary material.

Theorem 1**.**

The Fair Robust Log-Loss Predictor (Definition 4) has equivalent dual formulation:

[TABLE]

*with Lagrange multipliers $\theta$ and $\lambda$ for moment matching and fairness constraints, respectively, and $n$ samples in the dataset. The parametric distribution of $\mathbb{P}$ is: *

[TABLE]

where $Z_{\theta}({\bf x})=e^{\theta^{\top}\phi({\bf x},1)}+e^{\theta^{\top}\phi({\bf x},0)}$ is the normalization constant. The parametric distribution of $\mathbb{Q}$ is defined using the following relationship with $\mathbb{P}$:

[TABLE]

Note that the predictor’s distribution is a member of the exponential family that is similar to standard binary logistic regression, but with the option to truncate the probability based on the value of $\lambda$. The truncation of $\mathbb{P}_{\theta,\lambda}(\widehat{y}=1|{\bf x},a,y)$ is from above when $0\!<\!{p_{\gamma_{1}}}/{\lambda}\!<\!1$ and $\gamma_{1}(a,y)\!=\!1$, and from below when $-1\!<\!{p_{\gamma_{1}}}/{\lambda}\!<\!0$ and $\gamma_{1}(a,y)\!=\!1$. The approximator’s distribution is computed from the predictor’s distribution using the quadratic function in Eq. (12), e.g., in the case where $\gamma_{1}(a,y)\!=\!1$:

[TABLE]

where $\rho\triangleq\mathbb{P}_{\theta,\lambda}(\widehat{y}\!=\!1|{\bf x},a,y)$.

Figure 1 illustrates the relationship between $\mathbb{P}_{\theta,\lambda}(\widehat{y}\!=\!1|{\bf x},a,y)$ and $\mathbb{Q}_{\theta,\lambda}(\widehat{y}\!=\!1|{\bf x},a,y)$ for decisions influencing the fairness of group one (i.e., $\gamma_{1}(a,y)=1$). When ${\lambda}/{p_{\gamma_{1}}}\!=\!0$, the approximator’s probability is equal to the predictor’s probability as shown in the plot as a straight line. Positive values of $\lambda$ curve the function upward (e.g., ${\lambda}/{p_{\gamma_{1}}}\!=\!1$) as shown in the plot. For larger $\lambda$ (e.g., ${\lambda}/{p_{\gamma_{1}}}\!=\!2$), some of the valid predictor probabilities ($0<\mathbb{P}<1$) map to invalid approximator probabilities (i.e., $\mathbb{Q}\geq 1$) according to the quadratic function. In this case (e.g., ${\lambda}/{p_{\gamma_{1}}}\!=\!2$ and $\mathbb{P}_{\theta,\lambda}(\widehat{y}\!=\!1|{\bf x},a,y)>0.5$), the predictor’s probability is truncated to ${p_{\gamma_{1}}}/{\lambda}\!=\!0.5$ according to Eq. (11). Similarly, for negative $\lambda$, the curve is shifted downward and the predictor’s probability is truncated when the quadratic function mapping results in a negative value of $\mathbb{Q}$. When $\gamma_{0}(a,y)=1$, the reverse shifting is observed, i.e., shifting downward when $\lambda>0$ and shifting upward when $\lambda<0$.

We contrast our reshaping function of the decision distribution (Figure 1) with the post-processing method of ? (?) shown in Figure 3. Here, we use $\mathbb{Q}(\widehat{Y}=1|{\bf x},a)$ to represent the estimating distributions (the approximator’s distribution in our method, and the standard logistic regression in ? (?)) and the post-processed predictions as $\mathbb{P}(\widehat{Y}=1|{\bf x},a)$. Both shift the positive prediction rates of each group to provide fairness. However, our approach provides a monotonic and parametric transformation, avoiding the criticisms that ? (?)’s modification (flipping some decisions) is partially random, creating an unrestricted hypothesis class (?). Additionally, since our parametric reshaping function is learned within an in-processing method, it avoids the noted suboptimalities that have been established for certain population distributions when employing post-processing alone (?).

Enforcing fairness constraints

The inner maximization in Eq. (10) finds the optimal $\lambda$ that enforces the fairness constraint. From the perspective of the parametric distribution of $\mathbb{P}$, this is equivalent to finding threshold points (e.g., ${p_{\gamma_{1}}}/{\lambda}$ and $1-{p_{\gamma_{0}}}/{\lambda}$) in the $\min$ and $\max$ function of Eq. (11) such that the expectation of the truncated exponential probabilities of $\mathbb{P}$ in group $\gamma_{1}$ match the one in group $\gamma_{0}$. Given the value of $\theta$, we find the optimum $\lambda^{*}$ directly by finding the threshold points. We first compute the exponential probabilities $P_{e}(\widehat{y}=1|{\bf x},a,y)=\exp(\theta^{\top}\phi({\bf x},1))/{Z_{\theta}({\bf x})}$ for each examples in $\gamma_{1}$ and $\gamma_{0}$. Let $E_{1}$ and $E_{0}$ be the sets that contain $P_{e}$ for group $\gamma_{1}$ and $\gamma_{0}$ respectively. Finding $\lambda^{*}$ given the sets $E_{1}$ and $E_{0}$ requires sorting the probabilities for each set, and then iteratively finding the threshold points for both sets simultaneously. We refer to the supplementary material for the detailed algorithm.

Learning

Our learning process seeks parameters $\theta,\lambda$ for our distributions ($\mathbb{P}_{\theta,\lambda}$ and $\mathbb{Q}_{\theta,\lambda}$) that match the statistics of the approximator’s distribution with training data ($\theta$) and provide fairness ($\lambda$), as illustrated in Eq. (10). Using our algorithm from the previous subsection to directly compute the best $\lambda$ given arbitrary values of $\theta$, denoted $\lambda^{*}_{\theta}$, the optimization of Eq. (10) reduces to a simpler optimization solely over $\theta$, as described in Theorem 2.

Theorem 2**.**

Given the optimum value of $\lambda^{*}_{\theta}$ for $\theta$, the dual formulation in Eq. (10) reduces to:

[TABLE]

Here, $T({\bf x},\theta)\triangleq 1$ if the exponential probability is truncated (for example when $e^{\theta^{\top}\phi({\bf x},1)}/{Z_{\theta}({\bf x})}>{p_{\gamma_{1}}}/{\lambda^{*}_{\theta}}$, $\gamma_{1}(a,y)=1$, and $\lambda^{*}_{\theta}>0$), and is [math] otherwise.

We present an important optimization property for our objective function in the following theorem.

Theorem 3**.**

The objective function in Theorem 2 (Eq. (13)) is convex with respect to $\theta$.

To improve the generalizability of our parametric model, we employ a standard L2 regularization technique that is common for logistic regression models: $\theta^{*}=\textstyle\mathop{\mathrm{argmin}}_{\theta}\textstyle\sum_{({\bf x},a,y)\in\mathcal{D}}\ell_{\theta,\lambda^{*}_{\theta}}({\bf x},a,y)+\tfrac{C}{2}\|\theta\|_{2}^{2},$ where $C$ is the regularization constant. We employ a standard batch gradient descent optimization algorithm (e.g., L-BFGS) to obtain a solution for $\theta^{*}$.​444We refer the reader to the supplementary material for details. We also compute the corresponding solution for the inner optimization, $\lambda^{*}_{\theta^{*}}$, and then construct the optimal predictor and approximator’s parametric distributions based on the values of $\theta^{*}$ and $\lambda^{*}_{\theta^{*}}$.

Inference

In the inference step, we apply the optimal parametric predictor distribution $\mathbb{P}_{\theta^{*},\lambda^{*}_{\theta^{*}}}$ to new example inputs $({\bf x},a)$ in the testing set. Given the value of $\theta^{*}$ and $\lambda^{*}_{\theta^{*}}$, we calculate the predictor’s distribution for our new data point using Eq. (11). Note that the predictor’s parametric distribution also depends on the group membership of the example. For fairness constraints not based on the actual label $Y$, e.g., D.P., this parametric distribution can be directly applied to make predictions. However, for fairness constraints that depend on the true label, e.g., E.Opp. and E.Odd., we introduce a prediction procedure that estimates the true label using the approximator’s parametric distribution.

For fairness constraints that depend on the true label, our algorithm outputs the predictor and approximator’s parametric distributions conditioned on the value of true label, i.e., $\mathbb{P}(\widehat{y}|{\bf x},a,y)$ and $\mathbb{Q}(\widehat{y}|{\bf x},a,y)$. Our goal is to produce the conditional probability of $\widehat{y}$ that does not depend on the true label, i.e., $\mathbb{P}(\widehat{y}|{\bf x},a)$. We construct the following procedure to estimate this probability. Based on the marginal probability rule, $\mathbb{P}(\widehat{y}|{\bf x},a)$ can be expressed as:

[TABLE]

However, since we do not have access to $P(y|{\bf x},a)$, we cannot directly apply this expression. Instead, we approximate $P(y|{\bf x},a)$ with the approximator’s distribution $\mathbb{Q}(\widehat{y}|{\bf x},a)$. Using the similar marginal probability rule, we express the estimate as:

[TABLE]

By rearranging the terms above, we calculate the estimate as:

[TABLE]

which is directly computed from the approximator’s parametric distribution produced by our model using Eq. (12). Finally, to obtain the predictor’s conditional probability estimate ($\mathbb{P}(\widehat{y}|{\bf x},a)$), we replace $P(y|{\bf x},a)$ in Eq. (14) with $\mathbb{Q}(\widehat{y}|{\bf x},a)$ calculated from Eq. (16).

Asymptotic convergence property

The ideal behavior of an algorithm is an important consideration in its design. Asymptotic convergence properties consider a learning algorithm when it is provided with access to the population distribution $P({\bf x},a,y)$ and a fully expressive feature representation. We show in Theorem 4 that in the limit, our method finds a predictor distribution that has a desirable characteristic in terms of the Kullback-Leibler (KL) divergence from the true distribution.

Theorem 4**.**

Given the population distribution $P({\bf x},a,y)$ and a fully expressive feature representation, our formulation (Def. 4) finds the fair predictor with the minimal KL-divergence from $P({\bf x},a,y)$.

We next show in Theorem 5 that for the case where the fairness constraint depends on the true label (e.g., E.Opp. and E.Odd.), our prediction procedure outputs a predictor distribution with the same desired characteristic, after being marginalized over the true label.

Theorem 5**.**

For fairness constraints that depend on the true label, our inference procedure in Eq. (14) produces the marginal predicting distribution $\mathbb{P}$ of the fair predictor distribution with the closest KL-divergence to $P({\bf x},a,y)$ in the limit.

Experiments

Illustrative behavior on synthetic data

We illustrate the key differences between our model and logistic regression with demographic parity requirements on 2D synthetic data in Figure 3. The predictive distribution includes different truncated probabilities for each group: raising the minimum probability for group $A=1$ and lowering the maximum probability for group $A=0$. This permits a decision boundary that differs significantly from the logistic regression decision boundary and better realizes the desired fairness guarantees. In contrast, post-processing methods using logistic regression as the base classifier (?) are constrained to reshape the given unfair logistic regression predictions without shifting the decision boundary orientation, often leading to suboptimality (?).

Datasets

We evaluate our proposed algorithm on three benchmark fairness datasets:

1. (1)

The UCI Adult (?) dataset includes 45,222 samples with an income greater than $50k considered to be a favorable binary outcome. We choose gender as the protected attribute, leaving 11 other features for each example. 2. (2)

The ProPublica’s COMPAS recidivism dataset (?) contains 6,167 samples, and the task is to predict the recidivism of an individual based on criminal history, with the binary protected attribute being race (white and non-white) and an additional nine features. 3. (3)

The dataset from the Law School Admissions Council’s National Longitudinal Bar Passage Study (?) has 20,649 examples. Here, the favorable outcome for the individual is passing the bar exam, with race (restricted to white and black only) as the protected attribute, and 13 other features.

Comparison methods

We compare our method (Fair Log-loss) against various baseline/fair learning algorithms that are primarily based on logistic regression as the base classifier:

1. (1)

Unconstrained logistic regression is a standard logistic regression model that ignores all fairness requirements. 2. (2)

The cost sensitive reduction approach by ? (?) reduces fair classification to learning a randomized hypothesis over a sequence of cost-sensitive classifiers. We use the sample-weighted implementation of Logistic Regression in scikit-learn as the base classifier, to compare the effect of the reduction approach. We evaluate the performance of the model by varying the constraint bounds across the set $\epsilon\in\{.001,.01,.1\}$. 3. (3)

The constraint-based learning method555https://github.com/mbilalzafar/fair-classification of (?; ?) uses a covariance proxy measure to achieve equalized odds (under the name disparate mistreatment) (?), and improve the disparate impact ratio (?), which we use as a baseline method to evaluate demographic parity violation. They cast the resulting non-convex optimization as a disciplined convex-concave program in training time. We use the logistic regression as the base classifier. 4. (4)

For demographic parity, we compare with the reweighting method (reweighting) of ? (?), which learns weights for each combination of class label and protected attribute and then uses these weights to resample from the original training data which yields a new dataset with no statistical dependence between class label and protected attribute. The new balanced dataset is then used for training a classifier. We use IBM AIF360 toolkit to run this method. 5. (5)

For equalized odds, we also compare with the post-processing method of ? (?) which transforms the classifier’s output by solving a linear program that finds a prediction minimizing misclassification errors and satisfying the equalized odds constraint from the set of probability formed by the convex hull of the original classifier’s probabilities and the extreme point of probability values (i.e., zero and one).

Evaluation measures and setup

Data-driven fair decision methods seek to minimize both prediction error rates and measures of unfairness. We consider the misclassification rate (i.e., the 0-1 loss, $\mathbb{E}[\widehat{Y}\neq Y]$) on a withheld test sample to measure prediction error. To quantify the unfairness of each method, we measure the degree of fairness violation for demographic parity (D.P.) as: \big{|}\mathbb{E}[\mathbb{I}(\widehat{Y}=1)|A=1]-\mathbb{E}[\mathbb{I}(\widehat{Y}=1)|A=0]\big{|}, and the sum of fairness violations for each class to measure the total violation for equalized odds (E.Odd.) as: \sum_{y\in\{0,1\}}\big{(}\big{|}\mathbb{E}[\mathbb{I}(\widehat{Y}=1)|A=1,Y=y]-\mathbb{E}[\mathbb{I}(\widehat{Y}=1)|A=0,Y=y]\big{|}\big{)}, to obtain a level comparison across different methods. We follow the methodology of ? (?) to give all methods access to the protected attribute both at training and testing time by including the protected attribute in the feature vector. We perform all of our experiments using 20 random splits of each dataset into a training set (70% of examples) and a testing set (30%). We record the averages over these twenty random splits and the standard deviation. We cross validate our model on a separate validation set using the best logloss to select an L2 penalty from ($\{.001,.005,.01,.05,.1,.2,.3,.4,.5\}$).

Experimental Results

Figure 4 provides the evaluation results (test error and fairness violation) of each method for demographic parity and equalized odds on test data from each of the three datasets Fairness can be vacuously achieved by an agnostic predictor that always outputs labels according to independent (biased) coin flips. Thus, the appropriate question to ask when considering these results is: “how much additional test error is incurred compared to the baseline of the unfair logistic regression model for how much of an increase in fairness?”

For demographic parity on the Adult dataset, our Fair Log-loss approach outperforms all baseline methods on average for both test error rate and for fairness violation, and on COMPAS dataset it achieves the lowest ratio of increased fairness over increased error. Additionally, the increase in test error over the unfair unconstrained logistic regression model is small. For demographic parity on the Law dataset, the relationship between methods is not as clear, but our Fair Log-loss approach still resides in the Pareto optimal set, i.e., there are no other methods that are significantly better than our result on both criteria. For equalized odds, Fair Log-loss provides the lowest ratios of increased fairness over increased error rate for the Adult and COMPAS datasets, and competitive performance on the Law dataset. The post-processing method provides comparable or better fairness at the cost of significantly higher error rates. This shows that the approximation in our prediction procedure does not significantly impact the performance of our method. In terms of the running time, our method is an order of magnitude faster than comparable methods (e.g., the train and test running time on one random split of the Adult dataset takes approximately 5 seconds by our algorithm, 80 seconds for the constraint-based method (?), and 100 seconds for the reduction-based method (?)).

Conclusions and Future Work

We have developed a novel approach for providing fair data-driven decision making in this work by deriving a new classifier from the first principles of distributionally robust estimation (?; ?; ?). We formulated a learning objective that imposes fairness requirements on the predictor and views uncertainty about the population distribution pessimistically while maintaining a semblance of the training data characteristics through feature-matching constraints. This resulted in a parametric exponential family conditional distribution that resemble a truncated logistic regression model.

In future work, we plan to investigate the setting in which group membership attributes are not available at testing time. Extending our approach using a plug-in estimator of $P(a|{\bf x})$ in the fairness constraints introduces this estimator in the parametric form of the model. Understanding the impact of error from this estimator on predictor fairness in both theory and practice is an important direction of research.

Acknowledgments

This work was supported, in part, by the National Science Foundation under Grant No. 1652530 and by the Future of Life Institute (futureoflife.org) FLI-RFP-AI1 program.

Appendix A Proofs

Proof of Theorem 1

Proof of Theorem 1.

Following the strong minimax duality analysis in "log-loss game" (?; ?), and strong Lagrange duality in convex optimization (?), we perform the following transformations:

[TABLE]

The transformation steps above are described as follows:

1. (a)

We flip the min and max order using strong minimax duality in "log-loss game" (?; ?). 2. (b)

We introduce the Lagrange dual variable $\theta$ to directly incorporate the moment matching constraints over $\mathbb{Q}$ into the objective function. 3. (c)

The objective is concave on $\mathbb{Q}$ for all $\theta$, while $\Xi$ is a convex set. Given a feasible solution on the relative interior of $\Xi$, strong Lagrange duality holds (?) and thus we can flip the optimization order of $\mathbb{Q}$ and $\theta$. 4. (d)

We flip the inner min and max over $\mathbb{P}$ and $\mathbb{Q}$ using the minimax duality, as in (a). 5. (e)

We introduce the Lagrange dual variable $\lambda$ to directly incorporate the fairness constraints over $\mathbb{P}$ into the objective function. 6. (f)

Similar to (c), we use strong Lagrange duality theorem to flip the optimization order of $\lambda$ and $\mathbb{P}$. 7. (g)

We group the expectation with respect to the empirical training data.

We now focus on the inner minimax formulation over $\mathbb{P}$ and $\mathbb{Q}$, given the value of $\theta$ and $\lambda$, i.e.:

[TABLE]

We aim to find the analytical solutions for $\mathbb{P}$ and $\mathbb{Q}$ in the equation above. First, we write the Lagrangian by incorporating the probability simplex constraints into the objective, i.e.:

[TABLE]

[TABLE]

We now take the derivative of the Lagrangian with respect to $\mathbb{P}(\widehat{y}|{\bf x},a,y)$:

[TABLE]

By setting Eq. (28) to zero, we rewrite $\mathbb{P}$ in terms of $\mathbb{Q}$:

[TABLE]

Using Eq. (29) we rewrite Eq. (26) as:

[TABLE]

Replacing $\mathbb{P}$ in Lagrangian we get:

[TABLE]

We now calculate the derivative with respect to $\mathbb{Q}$.

[TABLE]

Setting Eq. (32) to 0 yields:

[TABLE]

We analytically solve the normalization constraint for $\mathbb{P}$, i.e., $\sum_{\widehat{y}\in\mathcal{Y}}\mathbb{P}(\widehat{y}|{\bf x},a,y)=1$

[TABLE]

which yields following parametric form of the predictor distribution:

[TABLE]

Notice the similarity to standard logistic regression. Where in contrast, here the probability for each class is adjusted with terms $\frac{\beta({\bf x},a,y,\widehat{y})}{\widetilde{P}({\bf x},a,y)}$ to satisfy the fairness constraints.

From Eq. (29) we get the relation of $\mathbb{Q}$ and $\mathbb{P}$. Solving the normalization constraint for $\mathbb{Q}(\widehat{y}|{\bf x},a,y)$ yields:

[TABLE]

Thus, we can rewrite $\mathbb{Q}$ as:

[TABLE]

We consider the binary classification $\widehat{y},y=\{0,1\}$, and expand $\mathbb{Q}$ as:

[TABLE]

The above equation shows that the approximator’s distribution $\mathbb{Q}$ is a quadratic function of predictor $\mathbb{P}$, for example in the case where $\gamma_{1}(a,y)=1$:

[TABLE]

where $\rho=\mathbb{P}_{\theta,\lambda}(\widehat{y}\!=\!1|{\bf x},a,y)$. For the region where the function goes above 1 (or below 1 depending on sign of $F_{\lambda}$), the predictor’s probability must be truncated in terms of fairness function such that $\mathbb{Q}=1$ (or zero). We derive these cases in the following by considering that the complementary slackness ensures non-negativity of $\mathbb{Q}$.

The complementary slackness from the KKT condition requires:

[TABLE]

Suppose that $\mathbb{Q}(\widehat{y}|{\bf x},a,y)=0$, then:

[TABLE]

Since $F_{\lambda}(a,y,0)=0$, then the equation above reduces to:

[TABLE]

Observe that the above equation can only hold if $\gamma_{1}(a,y)=1$, or $\gamma_{0}(a,y)=1$. For the other cases, complementary slackness requires that $\beta({\bf x},a,y,\widehat{y})=0$ and $\mathbb{Q}(\widehat{y}|{\bf x},a,y)=\mathbb{P}(\widehat{y}|{\bf x},a,y)=\frac{e^{\theta^{\top}\phi({\bf x},\widehat{y})}}{Z_{\theta}({\bf x})}$.

Thus, we have the following cases:

[TABLE]

Therefore, if $\lambda\geq 0$, we have the following parametric form for the predictor distribution:

[TABLE]

and if $\lambda<0$:

[TABLE]

Note that if $\lambda=0$, all of the cases collapse to a single case $\mathbb{P}(\widehat{y}=1|{\bf x},a,y)=\frac{e^{\theta\phi({\bf x},1)}}{Z_{\theta}({\bf x})}$.

∎

Proof of Theorem 2

Proof of Theorem 2.

Given the optimum $\lambda^{*}_{\theta}$ for each $\theta$, Eq. (10) reduces to:

[TABLE]

Plugging the parametric distribution forms of $\mathbb{P}$ and $\mathbb{Q}$, for $\lambda^{*}_{\theta}>0$, we get:

[TABLE]

and for $\lambda^{*}_{\theta}<0$, we get:

[TABLE]

and for $\lambda^{*}_{\theta}=0$, we get:

[TABLE]

∎

Proof of Theorem 3

Proof of Theorem 3.

It is easy to see that given $\theta$ and $\lambda_{\theta}$, the function $\ell_{\theta,\lambda_{\theta}}({\bf x},a,y)$ is convex for each sample. It is a linear function with respect to $\theta$ for the truncated cases. For the “otherwise” case, we know that $\log Z_{\theta}({\bf x})$ is convex. Hence, the full objective given $\theta$ and $\lambda_{\theta}$ is convex (non-negative weighted sum of convex functions is convex).

Now, the objective function given $\theta$ can be written as:

[TABLE]

Since for each $\lambda$, the function $\sum_{({\bf x},a,y)\in\mathcal{D}}\ell_{\theta,\lambda_{\theta}}({\bf x},a,y)$ is convex, the objective in Eq. (59) is also convex with respect to $\theta$ (pointwise supremum of convex functions is convex). ∎

Proof of Theorem 4

Proof of Theorem 4.

A fully expressive feature representation constrains the approximator’s distribution in our primal formulation Eq. (9) to match the population distribution. Then, the optimization simplifies to:

[TABLE]

This means that the optimal solution of our method when learning from the population distribution with a fully expressive feature representation is the fair predictive distribution that has the minimum KL-divergence from the population distribution. ∎

Proof of Theorem 5

Proof of Theorem 5.

For fairness constraints that depend on the true label (e.g., E.Opp. and E.Odd.), as described in §Inference, we compute $\mathbb{P}^{*}(\widehat{y}|{\bf x},a)$ using Eq. (14) with the input of $\mathbb{P}^{*}(\widehat{y}|{\bf x},a,y)$ and the approximator’s distribution to approximate the true distribution. Based on the proof of Theorem 4, we know that, in the limit, the approximator’s distribution matches with the true distribution $P({\bf x},a,y)$. Hence, our prediction becomes the standard marginal probability rule (it is no longer an approximation), i.e.:

[TABLE]

Therefore, our predictor is the marginal predictor distribution computed from the fair predictor’s distribution with the closest KL-divergence from the true distribution, marginalized over the true label. ∎

Appendix B Optimization Details

Algorithm for finding the optimal fairness parameters

The inner maximization in Eq. (10) finds the optimal $\lambda$ that enforces the fairness constraint. From the perspective of the parametric distribution of $\mathbb{P}$, this is equivalent to finding threshold points (e.g., ${p_{\gamma_{1}}}/{\lambda}$ and $1-{p_{\gamma_{0}}}/{\lambda}$) in the $\min$ and $\max$ function of Eq. (11) such that the expectation of the truncated exponential probabilities of $\mathbb{P}$ in group $\gamma_{1}$ match the one in group $\gamma_{0}$. Given the value of $\theta$, we find the optimum $\lambda^{*}$ directly by finding the threshold. We first compute the exponential probabilities $P_{e}(\widehat{y}=1|{\bf x},a,y)=\exp(\theta^{\top}\phi({\bf x},1))/{Z_{\theta}({\bf x})}$ for each examples in $\gamma_{1}$ and $\gamma_{0}$. Let $E_{1}$ and $E_{0}$ be the sets that contains $P_{e}$ for group $\gamma_{1}$ and $\gamma_{0}$ respectively, and let $\bar{e}_{1}$ and $\bar{e}_{0}$ be the average of $E_{1}$ and $E_{0}$ respectively.

Finding $\lambda^{*}$ given the sets $E_{1}$ and $E_{0}$ requires sorting the probabilities for each set, and then iteratively finding the threshold points for both sets ($t_{1}$ and $t_{0}$ respectively) simultaneously as described in Algorithm 1. Without loss of generality666For the case when $\bar{e}_{1}<\bar{e}_{0}$, we flip the group membership and then $\lambda^{*}$ is the negative of the solution produced by the algorithm. When $\bar{e}_{1}=\bar{e}_{0}$, the exponential probabilities are already fair, and we set $\lambda^{*}=0$.​, the algorithm assumes that $\bar{e}_{1}>\bar{e}_{0}$. It find the threshold points by traversing the sorted list of points in $E_{1}$ and $E_{0}$ until it find the thresholds that ensure the equality of the average truncated probabilities in both groups.

The runtime of Algorithm 1 is dominated by sorting, i.e., $O(n\log n)$ time. However, if we perform subgradient based optimization on $\theta$, the value of the current $\theta$ in each iteration does not change much, and neither do the exponential probabilities in each group. By maintaining the sorted list in each iteration as the basis index, the next iteration will have the probabilities in a nearly sorted order. Therefore, we can improve the sorting cost requirement by running sorting algorithms that work best on a nearly sorted list (e.g. insertion sort, Timsort, or $P^{3}$-sort) with run times approaching $O(n)$ as the list is close to being fully sorted (?).

Subgradient-based Optimization

Our optimization objective is:

[TABLE]

Taking the gradient of the objective with respect to $\theta$, for $\lambda^{*}_{\theta}>0$, we get:

[TABLE]

and for $\lambda^{*}_{\theta}<0$, we get:

[TABLE]

and for $\lambda^{*}_{\theta}=0$, we get:

[TABLE]

This can be simplified as:

[TABLE]

or using our $T({\bf x},\theta)$ notation as:

[TABLE]