Fair Regression: Quantitative Definitions and Reduction-based Algorithms

Introduction

As machine learning touches increasingly critical aspects of our life, including education, healthcare, criminal justice and lending, there is a growing focus to ensure that the algorithms treat various subpopulations fairly (see, e.g., Barocas & Selbst, 2016; Podesta et al., 2014; Corbett-Davies & Goel, 2018; and references therein). These questions have been particularly extensively researched in the context of classification, where several quantitative measures of fairness have been proposed (Berk et al., 2017; Chouldechova, 2017; Hardt et al., 2016; Kleinberg et al., 2017), leading to a variety of algorithms that aim to satisfy them (see, e.g., Corbett-Davies & Goel, 2018, for an overview of the literature).

These classifier-based formulations appear to fit the settings where the decision space is discrete and small, such as accept/reject decisions in hiring, school admissions, or lending. However, in practice, the decision makers work with tools that estimate a continuous quantity, such as success on the job, GPA in the first year of college, or risk of default on a loan. Predictions of these quantities are treated as scores, which are used by human decision makers, perhaps in the context of a partly automated workflow, to reach final decisions (see, e.g., Waters & Miikkulainen, 2014; US Federal Reserve, 2007; Northpointe, 2010; Lowenkamp et al., 2012). While, in principle, a fair classification tool could be used to recommend the yes/no decision directly, such tools are often resisted by practitioners, because they limit their autonomy, whereas ranking or scoring tools do not have this drawback (Veale et al., 2018). In such situations, it is desirable to work with real-valued scores that satisfy some notion of fairness. Yet, despite ample motivation and use cases, the prior work on designing fair continuous predictors is quite limited in its scope compared with the generality of methods for fair classification (e.g., Hardt et al., 2016; Agarwal et al., 2018).

This paper seeks to diminish this gap by developing efficient algorithms for a substantially broader set of regression tasks and model classes than done before, in many cases providing the first method with theoretical performance guarantees.

We consider the problem of predicting a real-valued target, where the prediction quality is measured by any Lipschitz-continuous loss function. Each example contains a protected attribute, such as race or gender, with respect to which we seek to guarantee fairness. We study two definitions of fairness from previous literature: statistical parity (SP), which asks that the prediction be statistically independent of the protected attribute, and bounded group loss (BGL), which asks that the prediction error restricted to any protected group stay below some pre-determined level. We define fair regression as the task of minimizing the expected loss of our real-valued predictions, subject to either of these fairness constraints. By choosing the appropriate loss, we obtain a wide range of standard prediction tasks including least-squares, logistic, Poisson, and quantile regression (with labels and predictions restricted to a bounded set to obtain Lipschitz continuity). While we seek to solve the regression tasks under fairness constraints, our schemes only require access to standard risk minimization algorithms such as standard classification or least-squares regression.

Several prior works also seek predictors that exhibit some form of independence from the protected attribute similar to statistical parity. Calders et al. (2013), Johnson et al. (2016) and Komiyama et al. (2018) consider a more limited form of independence, expressed via a small number of moment constraints, such as lack of correlation, and design specific algorithms for linear least squares. Berk et al. (2017) study notions of individual and group fairness specialized to linear regression. Pérez-Suay et al. (2017) seek zero correlation in a reproducing kernel Hilbert space (RKHS), which can capture statistical independence, but it only yields predictors in the same RKHS and the loss is limited to least squares. Kamishima et al. (2012) and Fukuchi et al. (2013) seek to fit a probabilistic model that satisfies statistical independence, but they do not present efficient algorithms or statistical guarantees. In contrast, we consider full statistical independence, arbitrary model classes and Lipschitz losses, and our algorithms are efficient and come with statistical guarantees.

Our second fairness definition, bounded group loss, fits into the general framework of Alabi et al. (2018), whose goal is to minimize a general function of group-wise prediction losses, but their algorithm is less efficient (albeit still polynomial), and they do not provide statistical guarantees.

We design a separate algorithm for each of the two fairness definitions. For BGL, our insight is that the problem of loss minimization subject to a loss bound in each subpopulation can be algorithmically reduced to a weighted loss minimization problem for which standard approaches exist. For SP, the main obstacle is that the number of constraints is uncountable. Here, the main insight that allows us to design and analyze the algorithm is that if we discretize the real-valued prediction space, then the task of fair regression can be reduced to cost-sensitive classification under certain constraints. We build on the recent work of Agarwal et al. (2018), and use the special structure of our discretization scheme to develop several algorithms reducing to standard classification or regression problems without fairness constraints. We provide theoretical results to bound the computational cost, generalization error and fairness violation of the returned predictor for both of our fairness measures with arbitrary Lipschitz-continuous loss functions and with arbitrary regression-function classes of bounded complexity, again building on the analysis of Agarwal et al. (2018). Prior works in the regression setting lack such guarantees.

Empirically, we evaluate our method on several standard datasets, on the tasks of least-squares and logistic regression under statistical parity, with linear and tree-ensemble learners, and compare it with the unconstrained baselines as well as the technique of Johnson et al. (2016). Our method uncovers fairness–accuracy frontiers and provides the first systematic scheme for enforcing fairness in a significantly broader class of learning problems than prior work.

Usage guidelines. We envision the use of our algorithms in uncovering fairness–accuracy frontiers in a variety of applications. Any substantial tradeoffs along the frontier need to be analyzed. They might point to data issues requiring non-algorithmic interventions, such as gathering of additional (less biased) data or introduction of new features (Chen et al., 2018). As with other algorithmic fairness tools, in order to successfully use our algorithms in practice, it is essential to consider the societal context of the application (Selbst et al., 2018). In some contexts, the best fairness intervention might be to avoid a technological intervention altogether.

Problem Formulation

We consider a general prediction setting where the training examples consist of triples $(X,A,Y)$, where $X\in\mathcal{X}$ is a feature vector, $A\in\mathcal{A}$ is a protected attribute and $Y\in\mathcal{Y}\subseteq$ is the label. Throughout, we focus on the protected attribute taking a small number of discrete values, i.e., $\mathcal{A}$ is finite, but $\mathcal{X}$ is allowed to be continuous and high-dimensional. We make no specific assumptions about whether the protected attribute is included in the feature vector $X$ or not; also the set of labels $\mathcal{Y}$ can be discrete (but embedded in $$) or continuous. Given a set of predictors$\mathcal{F}$containing functions$f:\mathcal{X}\to$, our goal is to find$f\in\mathcal{F}$which is accurate in predicting$Y$given$X$while satisfying some fairness condition such as statistical parity or bounded group loss (formally defined below). Note that the functions$f$do not explicitly depend on$A$unless it is included in$X$.

We consider two quantitative definitions of fairness appearing in prior work on fair classification and regression.

The first definition, called statistical (or demographic) parity, says that the prediction should be independent of the protected attribute. In classification, it corresponds to the practice of affirmative action (see, e.g., Holzer & Neumark, 2006, and references therein) and it is also invoked to address disparate impact under the US Equal Employment Opportunity Commission’s “four-fifths rule,” which requires that the “selection rate for any race, sex, or ethnic group [must be at least] four-fifths (4/5) (or eighty percent) of the rate for the group with the highest rate.”See the Uniform Guidelines on Employment Selection Procedures, 29 C.F.R. §1607.4(D) (2015).

The characterization through the properties of the CDF of $f(X)$ is particularly useful when $f(X)$ can take any real values in $$, because it allows us to design efficient algorithms. It also makes it obvious that if$f$satisfies SP, then any classifier induced by thresholding$f$ will also satisfy SP.

Our second fairness definition, called bounded group loss, formalizes the requirement that the predictor’s loss remain below some acceptable level for each protected group. In settings such as speech or face recognition, this corresponds to the requirement that all groups receive good service (cf. Buolamwini & Gebru, 2018). In other settings, such as lending and hiring, it aims to prevent situations when the predictor has a high error on some of the groups (cf. Section 3.3 of Corbett-Davies & Goel, 2018).

Hence, fair regression with BGL minimizes the overall loss, while controlling the worst loss on any protected group. By Lagrangian duality, this is equivalent to minimizing the worst loss on any group while maintaining good overall loss (referred to as max-min fairness). Unlike overall accuracy equality in classification (Dieterich et al., 2016), which requires the losses on all groups to be equal, BGL does not force an artificial decrease in performance on every group just to match the hardest-to-predict group. BGL can be used as a diagnostic for the potential shortcomings of a chosen featurization or dataset. If it is not possible to achieve a loss below $\zeta$ on some group, then to achieve fairness we need to collect more data for that group, or develop more informative features for individuals in that group.

2 Fair Regression

Statistical parity. Similar to prior works on fair classification (Agarwal et al., 2018), it is frequently desirable to have a tunable knob for navigating the fairness-accuracy tradeoff, such as $\zeta$ in the definition of bounded group loss. To allow such a tradeoff in SP, we consider slack parameters $\varepsilon_{a}$ for each attribute and define the fair regression task under SP as

Bounded group loss. In this case, the constrained optimization formulation follows directly from the definition. For the sake of flexibility, we allow specifying a different bound $\zeta_{a}$ for each attribute value, leading to the formulation

Randomized predictors. Similar to fair classification, in order to achieve better fairness–accuracy tradeoffs, we consider randomized predictors which first pick $f$ according to some distribution $Q$ and then predict according to $f$. We first introduce additional notation for the objective and constraints appearing in (1) and (2):

For a randomized predictor represented by a distribution $Q$, we have $\textup{loss}(Q)=\sum_{f}Q(f)\textup{loss}(f)$, $\gamma_{a,z}(Q)=\sum_{f}Q(f)\gamma_{a,z}(f)$, and $\gamma^{\textup{BGL}}_{a}(Q)=\sum_{f}Q(f)\gamma^{\textup{BGL}}_{a}(f)$.

Supervised Learning Oracles

(3) Cost-sensitive classification (CS). Our third type of oracle optimizes over classifiers $h:\mathcal{X}^{\prime}\to\{0,1\}$ from some class $\mathcal{H}$. As input, we are given a dataset $\{(X^{\prime}_{i},C_{i})\}_{i=1}^{n}$, where $X^{\prime}_{i}$ is a feature vector and $C_{i}$ indicates the difference between the cost (i.e., the loss) of predicting 1 versus 0; positive $C_{i}$ means that 0 is favored, negative $C_{i}$ means that 1 is favored. The goal is to find a classifier $h\in\mathcal{H}$, which minimizes the empirical cost relative to the cost of predicting all zeros: $\sum_{i=1}^{n}C_{i}h(X^{\prime}_{i})$.

CS reduces to weighted binary classification on the data $\{(W_{i},X^{\prime}_{i},Y_{i})\}_{i=1}^{n}$ with $Y_{i}=\mathbf{1}\{C_{i}\leq 0\}$ and $W_{i}=\left\lvert C_{i}\right\rvert$, where we minimize $\sum_{i=1}^{n}W_{i}\mathbf{1}\{h(X^{\prime}_{i})\neq Y_{i}\}$. Weighted classification oracles exist for many classifier families $\mathcal{H}$.

Fair Regression under Statistical Parity

We next show how to solve the fair regression problem (3) using a CS oracle. We begin by recasting the problem (3) as a constrained (and cost-sensitive) classification problem, which we then solve via the reduction approach of Agarwal et al. (2018), by repeatedly invoking the CS oracle.

We proceed in two steps. First we discretize our prediction space and show that a loss function in the discretized space approximates our original loss well, owing to its Lipschitz continuity. We then show how the fair regression problem in this discretized space can be turned into a constrained classification problem, which we solve via reduction.

where $\bar{z}=\lceil z\rceil_{\alpha}$ is the value of $z$ rounded up to the nearest integer multiple of $\alpha$. This allows us to replace the uncountable set of constraints indexed by $z\in$ with the finite set indexed by $z\in\mathcal{Z}$. Thus, denoting \underaccent{\bar}{\calF}=\{\underaccent{\bar}{f}:\>f\in\mathcal{F}\}, we have argued that the solution of (3), can be approximated by

Let $Q^{\star}$ be any feasible point of (3) and $\underline{\smash{Q}\!\!}\,\,^{\star}$ be the solution of (8). Then $\textup{loss}(\underline{\smash{Q}\!\!}\,\,^{\star})\leq\textup{loss}(Q^{\star})+\alpha$ and $\left\lvert\gamma_{a,z}(\underline{\smash{Q}\!\!}\,\,^{\star})\right\rvert\leq\varepsilon_{a}$ for all $a\in\mathcal{A}$, $z\in$.

2 Reduction to Constrained Classification

We next show that (8) can be rewritten as a constrained classification problem for the family of classifiers $\mathcal{H}=\{h_{f}:\>f\in\mathcal{F}\}$ defined in Section 3.

Given any distribution $D$ over $(X,A,Y)$ and any $f\in\mathcal{F}$, the cost and constraints satisfy \textup{cost}(h_{f})=\textup{loss}_{\alpha}(\underaccent{\bar}{f})+c_{0}, where $c_{0}$ is independent of $f$, and \gamma_{a,z}(h_{f})=\gamma_{a,z}(\underaccent{\bar}{f}) for all $a\in\mathcal{A}$, $z\in\mathcal{Z}$.

By linearity of expectation, the lemma implies analogous equalities also for distributions over $f$. Thus, in problem (8), we can replace the optimization over \underline{\smash{Q}\!\!}\,\,\in\Delta(\underaccent{\bar}{\calF}) with $Q\in\Delta(\mathcal{H})$. Notice that while we started from discretized regressors in problem (8), Lemma 1 allows us to work with the full classifier family $\{h_{f}:\>f\in\mathcal{F}\}$, which is important as we typically only have computational oracles for non-discretized classes $\mathcal{H}$ and $\mathcal{F}$. We next show how to solve an empirical version of this classification problem.

3 Algorithm and Generalization Bounds

We are interested in the following empirical optimization problem, which is, according to Lemma 1, an empirical approximation of the original problem (8):

The slacks $\widehat{\varepsilon}_{a}$ are slightly larger than $\varepsilon_{a}$ to compensate for finite-sample errors in measuring constraint violations (more on that below). This problem is a special case of that studied by Agarwal et al. (2018) with a key difference. Since the distribution of $Z$ is known, we can take expectation according to $Z$ rather than a sample, which leads to substantially better estimates of constraint violations. Thus, our objective uses a product of an empirical distribution over $(X,A,Y)$ with the uniform distribution over $Z$ rather than an i.i.d. sample as assumed by Agarwal et al.. However, their algorithm and generalization bounds still apply (as we show in our proofs).

It solves the saddle-point problem $\min_{Q}\max_{\boldsymbol{\lambda}}L(Q,\boldsymbol{\lambda})$ over $Q\in\Delta(\mathcal{H})$ and $\boldsymbol{\lambda}\geq\mathbf{0}$, $\lVert\boldsymbol{\lambda}\rVert_{1}\leq B$, by treating it as a two-player zero-sum game (see Algorithm 1 for details).

We bound the suboptimality and fairness of the returned solution largely following their analysis. Let $R_{n}(\mathcal{H})$ denote the Rademacher complexity of $\mathcal{H}$ (see Eq. 17 in Appendix C). To state the bounds, recall an assumption from their paper on the setting of the empirical slacks $\widehat{\varepsilon}_{a}$:

There exist $C,C^{\prime}>0$ and $\beta\leq 1/2$ such that $R_{n}(\mathcal{H})\leq Cn^{-\beta}$ and $\widehat{\varepsilon}_{a}=\varepsilon_{a}+C^{\prime}n_{a}^{-\beta}$, where $n_{a}$ is the number of samples with $A=a$.

Under this assumption, we obtain the following guarantees.The notation $\widetilde{O}(\cdot)$ suppresses polynomial dependence on $\ln n$, $\ln\,\lvert\mathcal{A}\rvert$, and $\ln(1/\delta)$.

Let Assumption 1 hold for $C^{\prime}\geq 2C+2+\sqrt{2\ln(4\lvert\mathcal{A}\rvert N/\delta)}$, where $\delta>0$. Let $Q^{\star}$ be any feasible distribution for the fair regression problem (3). Then Algorithm 1 with $\nu\propto n^{-\beta}$, $B\propto n^{\beta}$, and $N\propto n^{\beta}$ terminates in $O\left(n^{4\beta}\ln(n^{\beta}\lvert\mathcal{A}\rvert)\right)$ iterations and returns $\widehat{Q}$, which, when viewed as a distribution over \underaccent{\bar}{\calF}, satisfies with probability at least 1-$\delta$,

4 Efficient Implementation of Algorithm 1

It is not too difficult to show that each iteration of Algorithm 1 can be implemented in time $O(n\log n+\lvert\mathcal{A}\rvert N)$ plus the complexity of two calls to $\textsc{Best}_{h}$, on which we focus here, while deferring the remaining details to Appendix F.

This corresponds to a CS problem with $nN$ instances $\{(X^{\prime}_{i,z},C_{i,z})\}_{i\leq n,\,z\in\mathcal{Z}}$ defined as

The sum $\sum_{a}N\lambda_{a,Z}$ in the definition of $c_{\boldsymbol{\lambda}}$ can be precalculated once for each value of $Z$ in the overall time $O(N\lvert\mathcal{A}\rvert)$. After that the construction of the dataset takes time $O(nN)$.

Based on Assumption 1 and Theorem 2, we expect $N\propto n^{\beta}$, so this reduction to cost-sensitive classification takes time $\Omega(n^{1+\beta})$ and creates a dataset of size $\Omega(n^{1+\beta})$. This is substantially larger than the original problem of size $n$, given the typical value $\beta\approx 1/2$. We next describe two alternatives that run faster and only create datasets of size $n$.

Reduction to least-squares regression. The main overhead in the CS reduction above comes from the summation over $z\in\mathcal{Z}$, implicit in the expectation over $Z$ in Eq. (12). In order to eliminate this overhead, suppose we have access to a function $g_{\boldsymbol{\lambda}}$ such that

Fair Regression with Bounded Group Loss

The approach still follows the scheme of Agarwal et al. (2018), but thanks to the matched loss function between the objective and the constraints, fair regression can be reduced directly to regression, without the need for discretization. We first replace the problem (4) by its empirical version

We give a detailed pseudocode for our approach in Algorithm 2 in Appendix D, and describe the main differences from Algorithm 1 here. As before, the algorithm alternates between exponentiated gradient updates on $\boldsymbol{\lambda}$ and best responses for $Q$ to compute an approximate saddle point:

The saddle point always exists. However, unlike the fair regression problem under SP, the fair regression problem under BGL, i.e., Eq. (13), might be infeasible. Therefore, Algorithm 2 explicitly checks whether the distribution $\widehat{Q}$ that it finds satisfies constraints of Eq. (13).

The other main difference between Algorithms 1 and 2 is in the computation of the best response $f$ to any given $\boldsymbol{\lambda}$, which requires solving the problem

Denoting by $n_{a}$ the number of samples with $A_{i}=a$, this minimization can be written as

which can be solved using one call to the weighted risk-minimization oracle.

We finish this section with the optimality and fairness guarantees for $\smash{\widehat{Q}}$ returned by Algorithm 2. We assume that $\smash{\widehat{\zeta}_{a}}$ are set according to the Rademacher complexity of $\mathcal{F}$:

There exist $C,C^{\prime}>0$ and $\omega\leq 1/2$ such that $R_{n}(\mathcal{F})\leq Cn^{-\omega}$ and $\smash{\widehat{\zeta}_{a}}=\zeta_{a}+C^{\prime}n_{a}^{-\omega}$, where $n_{a}$ is the number of samples where $A=a$.

Let Assumption 2 hold for $C^{\prime}\geq 4C+2+\sqrt{2\ln(4\lvert\mathcal{A}\rvert/\delta)}$, where $\delta>0$. Then Algorithm 2 with $\nu\propto n^{-\omega}$ and $B\propto n^{\omega}$ terminates in $O(n^{4\omega}\ln\,\lvert\mathcal{A}\rvert)$ iterations and returns $\smash{\widehat{Q}}$ such that, with probability at least 1-$\delta$, one of the following holds:

$\widehat{Q}\neq\textit{null}$ and, for any $Q^{\star}$ feasible in problem (4), ${}\quad\textup{loss}(\widehat{Q})\leq\textup{loss}(Q^{\star})+\widetilde{O}(n^{-\omega})$ {}\quad\gamma^{\textup{BGL}}_{a}(\widehat{Q})\leq\zeta_{a}+\widetilde{O}(n_{a}^{-\omega})\quad\text{for alla\in\mathcal{A}.}

$\widehat{Q}=\textit{null}$ and problem (4) is infeasible.

Experiments

We evaluate our method on the tasks of least-squares regression and logistic regression under statistical parity. We use the following three datasets:

Adult: The adult income dataset (Lichman, 2013) has 48,842 examples. The task is to predict the probability that an individual makes more than $50k per year via logistic loss minimization, with gender as the protected attribute.

Law school: Law School Admissions Council’s National Longitudinal Bar Passage Study (Wightman, 1998) has 20,649 examples. The task is to predict a student’s GPA (normalized to $$) via square loss minimization, with race as the protected attribute (white versus non-white).

Communities & crime: The dataset contains socio-economic, law enforcement, and crime data about communities in the US (Redmond & Baveja, 2002) with 1,994 examples. The task is to predict the number of violent crimes per 100,000 population (normalized to $$) via square loss minimization, with race as the protected attribute (whether the majority population of the community is white).

For the two larger datasets (adult and law school), we also created smaller (subsampled) versions by picking random 2,000 points. Thus we ended up with a total of five datasets, and split each into 50% for training and 50% for testing.

We ran Algorithm 1 on each training set over a range of constraint slack values $\hat{\varepsilon}$, with a fixed discretization grid of size 40: $\mathcal{Z}=\{1/40,2/40,\ldots,1\}$. Among the solutions for different $\hat{\varepsilon}$, we selected the ones on the Pareto front based on their training losses and SP disparity $\max_{a,z}\{\hat{\gamma}_{a,z}\}$. We then evaluated the selected predictors on the test set, and show the resulting Pareto front in Figure 1.

We ran our algorithm with the three types of reductions from Section 4.4: reductions to cost-sensitive (CS) oracles, least-squares (LS) oracles, and logistic-loss minimization (LR) oracles. Our CS oracle sought the linear model minimizing weighted hinge-loss (as a surrogate for weighted classification error). Because of unfavorable scaling of the cost-sensitive problem sizes (see Section 4.4), we only ran the CS oracle on the three small datasets. We considered two variants of LS and LR oracles: linear learners from scikit-learn (Pedregosa et al., 2011), and tree ensembles from XGBoost (Chen & Guestrin, 2016). Tree ensembles heavily overfitted smaller datasets, so we only show their performance on two larger datasets. We only used LR oracles when the target loss was logistic, whereas we used LS oracles across all datasets.

In addition to our algorithm, we also evaluated regression without any fairness constraints, and two baselines from the fair classification and fair regression literature.

On the three datasets where the task was least-squares regression, we evaluated the full substantive equality of opportunity (SEO) estimate of Johnson et al. (2016). It can be obtained in a closed form by solving for the linear model that minimizes least-squares error while having zero correlation with the protected attribute. In contrast, our method seeks to limit not just correlation, but statistical dependence.

On the two datasets where the task was logistic regression, we ran the fair classification (FC) reduction of Agarwal et al. (2018) with the same LR oracles as in our algorithm. For this choice of oracles, the classifiers returned by FC are implemented by logistic models and return real-valued scores, which we evaluated. We ran FC across a range of trade-offs between classification accuracy and statistical parity (in the classification sense) and show the resulting Pareto front. Note that FC only enforces statistical parity (SP) when the scores are thresholded at zero, whereas our method enforces SP across all thresholds.

In Figure 1, we see that all of our reductions are able to significantly reduce disparity, without strongly impacting the overall loss. On communities & crime, there is a more substantial accuracy–fairness tradeoff, which can be used as a starting point to diagnose the data quality for the two racial subgroups. Our methods dominate SEO in least-squares tasks, but are slightly worse than FC in logistic regression. The difference is statistically significant only on adult, where it points to the limitations of our LS and LR reduction heuristics. However, for the most part, LR and LS reductions achieve tradeoffs on par with the CS reduction, and are substantially faster to run (see Appendix G). The results on adult and adult subsampled suggest that reducing to a matching loss is preferable over reducing to another loss.

In summary, we have shown that our scheme efficiently handles a range of losses and regressor classes and, where possible, diminishes disparity while maintaining the overall accuracy. The emergence of FC as a strong baseline for logistic regression suggests that our regression-based reduction heuristics can be further improved, which we leave open for future research.

Acknowledgements

ZSW is supported in part by a Google Faculty Research Award, a J.P. Morgan Faculty Award, and a Facebook Research Award. Part of this work was completed while ZSW was at Microsoft Research-New York City.

References

Appendix A Proof of Lemma 1

Now plugging in u=\underaccent{\bar}{f}(x) and using the fact that for $z\in\mathcal{Z}$, we have \mathbf{1}\{\underaccent{\bar}{f}(x)\geq z\}=\mathbf{1}\{f(x)\geq z\}=h_{f}(x,z), we obtain

For $z\in\mathcal{Z}$, we can rewrite \gamma_{a,z}(\underaccent{\bar}{f}) as

Appendix B Iteration Complexity of Algorithm 1

Algorithm 1 terminates in at most $\frac{16B^{2}\log(2\lvert\mathcal{A}\rvert N+1)}{\nu^{2}}$ iterations. Furthermore, if $Q$ is any feasible point of (11) then the solution $\widehat{Q}$ returned by Algorithm 1 satisfies:

This result is essentially a corollary of Theorem 1, and Lemmas 2 and 3 of Agarwal et al. (2018). Specifically, we note that the constraints appearing in our problem (11) can be cast in their general framework along the lines of their Example 1, with a total of $2\lvert\mathcal{A}\rvert N$ constraints. Following their Example 3, we obtain that the maximal constraint violation $\rho$, needed in their Theorem 1, is at most 2. We further observe that the violation of the i.i.d. structure by explicit averaging over $z$ values does not impact their optimization analysis in any way. Therefore, their Theorem 1 with $\rho=2$ implies that our Algorithm 1 finds a $\nu$-approximate saddle point of the Lagrangian in at most $\frac{16B^{2}\log(2\lvert\mathcal{A}\rvert N+1)}{\nu^{2}}$ iterations as desired.

To bound $\widehat{\textup{cost}}(\widehat{Q})$ and \bigl{\lvert}\widehat{\gamma}_{a,z}(\widehat{Q})\bigr{\rvert} we appeal to their Lemmas 2 and 3. First note that their approach applies to the objective of our problem (11) as long as the costs $c(y,z)$ are in $$(see their footnote 4). However, in our case, these can be in$$ (see Eq. 9). This does not affect their Theorem 1 and Lemma 2, but their Lemma 3 now holds with the right-hand side equal to $\frac{2+2\nu}{B}$ instead of $\frac{1+2\nu}{B}$. Their Lemma 2 immediately yields the bound $\widehat{\textup{cost}}(\widehat{Q})\leq\widehat{\textup{cost}}(Q)+2\nu$, whereas the modified Lemma 3 yields the bound \bigl{\lvert}\widehat{\gamma}_{a,z}(\widehat{Q})\bigr{\rvert}\leq\widehat{\varepsilon}_{a}+\frac{2+2\nu}{B} for all $a,z$, finishing the proof. ∎

Appendix C Proof of Theorem 2

The Rademacher complexity of a class $\mathcal{G}$ can be used to obtain uniform bounds of any Lipschitz continuous transformations of $g\in\mathcal{G}$ as follows:

Let $D$ be a distribution over a pair of random variables $(S,U)$ taking values in $\mathcal{S}\times\mathcal{U}$. Let $\mathcal{G}$ be a class of functions $g:\mathcal{U}\to$, and let $\varphi:\mathcal{S}\times\to$ be a contraction in its second argument, i.e., for all $s\in\mathcal{S}$ and all $t,t^{\prime}\in$, $\left\lvert\varphi(s,t)-\varphi(s,t^{\prime})\right\rvert\leq\left\lvert t-t^{\prime}\right\rvert$. Then with probability at least $1-\delta$, for all $g\in\mathcal{G}$,

where the expectation is with respect to $D$ and the empirical expectation is based on $n$ i.i.d. draws from $D$. If $\varphi$ is also linear in its second argument then a tighter bound holds, with $4R_{n}(\mathcal{G})$ replaced by $2R_{n}(\mathcal{G})$.

Let $\Phi\coloneqq\{\varphi_{g}\}_{g\in\mathcal{G}}$ be the class of functions $\varphi_{g}:(s,u)\mapsto\varphi(s,g(u))$. By Theorem 3.2 of Boucheron et al. (2005), we then have with probability at least $1-\delta$, for all $g$,

We will next bound $R_{n}(\Phi)$ in terms of $R_{n}(\mathcal{G})$. For a fixed tuple $(s_{1},u_{1}),\dotsc,(s_{n},u_{n})$, we have

where the first inequality follows from Theorem 12(5) of Bartlett & Mendelson (2002) and the last inequality follows from the contraction principle of Ledoux & Talagrand (1991), specifically their Theorem 4.12. Dividing by $n$ and taking a supremum over $(s_{1},u_{1}),\dotsc,(s_{n},u_{n})$ yields the bound

Together with the bound (18), this proves the lemma for an arbitrary contraction $\varphi$. If $\varphi$ is linear in its second argument, we get a tighter bound by invoking Theorem 4.4 of Ledoux & Talagrand (1991) instead of their Theorem 4.12. ∎

Our proof largely follows the proof of Theorems 2 and 3 of Agarwal et al. (2018). We first use Lemma 2 to show that by solving the empirical problem (11), we also obtain an approximate solution of the corresponding population problem:

The theorem will then follow by invoking the equivalence between problem (19) and the discretized fair regression (8), and adding up various approximation errors.

To bound the deviations in the cost, we need to be a bit careful, because the definition of $\widehat{\textup{cost}}$ mixes the empirical expectation over the data with the averaging over $z\in\mathcal{Z}$. For the analysis, we therefore define

Since c(\underaccent{\bar}{Y},z)\in, we can invoke Lemma 2 with S=c(\underaccent{\bar}{Y},z), $U=(X,z)$, $\mathcal{G}=\mathcal{H}$, and $\varphi(s,t)=st$ to obtain that with probability at least $1-\delta/2$ for all $z\in\mathcal{Z}$ and all $h\in\mathcal{H}$

where the last equality follows by Assumption 1 and the setting $N\propto n^{\beta}$. Taking an average over $z\in\mathcal{Z}$ and a convex combination according to any $Q\in\Delta(\mathcal{H})$, we obtain by Jensen’s inequality that with probability at least $1-\delta/2$ for all $Q\in\Delta(\mathcal{H})$

To bound the deviations in the constraints, we invoke Lemma 2 with $S=1$, $U=(X,z)$, $\mathcal{G}=\mathcal{H}$, and $\varphi(s,t)=st$, but apply it to the data distribution conditioned on $A=a$. We thus obtain that with with probability at least $1-\delta/2$ for all $a\in\mathcal{A}$, $z\in\mathcal{Z}$, and $h\in\mathcal{H}$

By Jensen’s inequality this also means that with probability at least $1-\delta/2$ for all $a\in\mathcal{A}$, $z\in\mathcal{Z}$, and $Q\in\Delta(\mathcal{H})$

In the remainder of the analysis, we assume that Eqs. (20) and (21) both hold, which occurs with probability at least $1-\delta$ by the union bound.

Putting it all together.

Given the settings of $\nu$, $B$ and $N$, we obtain by Theorem 4 that Algorithm 1 terminates in O\bigl{(}n^{4\beta}\ln(n^{\beta}\lvert\mathcal{A}\rvert)\bigr{)} iterations, as desired, and returns a distribution $\widehat{Q}$ which compares favorably with any feasible point $Q$ of the empirical problem (11), meaning that for any such $Q$, we have

Now bounding $\widehat{\textup{cost}}(\widehat{Q})$ and $\widehat{\textup{cost}}(Q)$ in Eq. (22) via the uniform convergence bound (20), we obtain

Above, we assumed that $Q$ was a feasible point of the empirical problem (11). However, assuming that Eq. (21) holds, any feasible solution of the population problem (19) is also feasible in the empirical problem (11) thanks to our setting of $C^{\prime}$. Thus, Eqs. (24) and (25) show that the solution $\widehat{Q}$ is approximately feasible and approximately optimal in the population problem (19). It remains to relate $\widehat{Q}$ to the original fair regression problem (3).

First, by Lemma 1 and Eqs. (24) and (25), we can interpret $\widehat{Q}$ as a distribution over the set of discretized regressors \underaccent{\bar}{\calF} and obtain that for all \underline{\smash{Q}\!\!}\,\,\in\Delta(\underaccent{\bar}{\calF}) that are feasible in the discretized fair regression problem (8), we have

where in Eq. (27) we have expanded the domain $z\in\mathcal{Z}$ to $z\in$ thanks to Eq. (7). Finally, by substituting the solution $\underline{\smash{Q}\!\!}\,\,^{*}$ of problem (8) as $\underline{\smash{Q}\!\!}\,\,$ in Eq. (26) and applying Theorem 1, we obtain that for any $Q^{*}\in\Delta(\mathcal{F})$ that is feasible in the discretized fair regression problem (8), we have

where the last equality follows by our setting $\alpha=1/N=O(n^{-\beta})$. The theorem now follows from Eqs. (28) and (27).

Appendix D Algorithm for Fair Regression under Bounded Group Loss

In this section we provide a detailed pseudocode of our algorithm for fair regression under BGL, described at a high level in Section 5.

Appendix E Proof of Theorem 3

The analysis of Algorithm 2 proceeds similarly to the analysis of Algorithm 1.

Similarly to our analysis of Algorithm 1 in Theorem 4, we can appeal to Theorem 1, and Lemmas 2 and 3 of Agarwal et al. (2018). While our objective and constraints are for the distributions $Q$ over $$-valued predictors$f\in\mathcal{F}$, whereas their analysis is for distributions over$\{0,1\}$-valued classifiers, we can still directly use their Theorem 1, and Lemmas 2 and 3, because they only rely on the bilinear structure of the Lagrangian with respect to$Q$and$\boldsymbol{\lambda}$and the boundedness of the objective and constraints, which all hold in our setting. The maximal constraint violation, needed in their Theorem 1, is$\rho\leq 1$. Therefore, their Theorem 1 implies that Algorithm 2 terminates in at most$4B^{2}\ln(\lvert\mathcal{A}\rvert+1)/\nu^{2}=O(n^{4\omega}\ln\lvert\mathcal{A}\rvert)$iterations and finds a$\nu$-approximate saddle point$\widehat{Q}$of problem (14), albeit sometimes it ends up returning null instead of$\widehat{Q}$. To prove the theorem we consider two cases.

Given the settings of $\nu$ and $B$, and using Lemmas 2 and 3 of Agarwal et al. (2018), we obtain that the $\nu$-approximate saddle point $\widehat{Q}$ of the empirical problem (14) satisfies

for any distribution $Q$ feasible in the empirical problem (13). Eq. (30) implies that in this case the algorithm returns $\widehat{Q}\neq\textit{null}$. It remains to argue that statements similar to (29) and (30) hold for true expectations rather than just empirical expectations.

By Assumption 2 and our setting of $C^{\prime}$, this implies that with probability at least $1-\delta$, for all $Q\in\Delta(\mathcal{F})$,

We continue the analysis assuming that the uniform convergence bounds (31) and (32) both hold. Instantiating the bound (31) for $\widehat{\textup{loss}}(Q)$ and $\widehat{\textup{loss}}(\widehat{Q})$ in Eq. (29) yields, for any $Q$ feasible in the empirical problem (13),

where the last inequality follows by our setting of $B$, $\nu$, and $C^{\prime}$ as well as the bound $\widehat{\zeta}_{a}\leq\zeta_{a}+C^{\prime}n_{a}^{-\omega}$ from Assumption 2.

Above, we assumed that $Q$ was a feasible point of the empirical problem (13). However, assuming that Eq. (32) holds, any feasible solution of the population problem (4) is also feasible in the empirical problem (13) thanks to our setting of $C^{\prime}$. Thus, Eqs. (33) and (34) hold for any $Q^{*}$ feasible in (4), proving the theorem in this case.

In this case, the $\nu$-approximate saddle point $\widehat{Q}$ that the algorithm finds may still satisfy

in which case the algorithm returns $\widehat{Q}$ and the theorem holds vacuously since here is no feasible point $Q^{*}$. If the found approximate saddle point does not satisfy Eq. (35), then the algorithm returns null and the theorem also holds.

Appendix F Details for Efficient Implementation of Algorithm 1

The algorithm checks the suboptimality of this solution and terminates once the convergence $\nu$ is reached (Steps 1–1).

The updates of $\boldsymbol{\theta}_{t}^{+}$ and $\boldsymbol{\theta}_{t}^{-}$ (Step 1) run in time $O(\lvert\mathcal{A}\rvert N)$ assuming that $\boldsymbol{\gamma}(h_{t})$ has already been calculated. Similarly, the transformation of $\boldsymbol{\theta}_{t}^{+}$ and $\boldsymbol{\theta}_{t}^{-}$ to $\boldsymbol{\lambda}_{t}$ (Steps 1–1) runs in time $O(\lvert\mathcal{A}\rvert N)$, because the denominator is the same across all coordinates and so it needs to be computed only once. We next show that the remaining operations, except for the two $\textsc{Best}_{h}$ calls, run in time $O(n\log n+\lvert\mathcal{A}\rvert N)$.

Computation of $\widehat{\textup{cost}}(h_{f})$ and $\widehat{\boldsymbol{\gamma}}(h_{f})$. This computation is implicit in the calculation of Lagrangian in Steps 1–1, and also in the updates of $\boldsymbol{\theta}_{t}^{+}$ and $\boldsymbol{\theta}_{t}^{-}$ in Step 1. By Eq. (15),

so it can be calculated in time $O(n)$ in a single pass over examples. To calculate $\widehat{\boldsymbol{\gamma}}(h_{f})$, we keep the training examples partitioned into $\lvert\mathcal{A}\rvert$ disjoint sets according to their protected attribute $A$. Let $n_{a}$ be the number of examples with $A=a$. We sort these examples according to $f(X)$ in time $O(n_{a}\log n_{a})$. Now going through these examples from the largest $f(X)$ value to the smallest allows us to calculate the conditional expectations

Computation of $L(h_{f},\boldsymbol{\lambda})$ and $L(\widehat{Q}_{t},\boldsymbol{\lambda})$. Lagrangian is evaluated in Steps 1–1 to determine whether the algorithm has converged. Note that $L(Q,\boldsymbol{\lambda})$ depends on $Q$ only through $\widehat{\textup{cost}}(Q)$ and $\widehat{\boldsymbol{\gamma}}(Q)$. If we have already computed these, $L(Q,\boldsymbol{\lambda})$ can be evaluated in time $O(\lvert\mathcal{A}\rvert N)$. To calculate $L(h_{f},\boldsymbol{\lambda})$ for an arbitrary $h_{f}$, we first need to calculate $\widehat{\textup{cost}}(h_{f})$ and $\widehat{\boldsymbol{\gamma}}(h_{f})$, so the overall running time is $O(n\log n+\lvert\mathcal{A}\rvert N)$. To calculate $L(\widehat{Q}_{t},\boldsymbol{\lambda})$, note that $\widehat{\textup{cost}}(\widehat{Q}_{t})=\frac{1}{t}\sum_{t^{\prime}=1}^{t}\widehat{\textup{cost}}(h_{t^{\prime}})$, so we can obtain $\widehat{\textup{cost}}(\widehat{Q}_{t})$ from $\widehat{\textup{cost}}(\widehat{Q}_{t-1})$ at the cost of evaluation of $\widehat{\textup{cost}}(h_{t})$ and similarly for $\widehat{\boldsymbol{\gamma}}(\widehat{Q}_{t})$. Therefore, the first evaluation of the form $L(\widehat{Q}_{t},\boldsymbol{\lambda})$ takes time $O(n\log n+\lvert\mathcal{A}\rvert N)$ and each consequent evaluation takes time $O(\lvert\mathcal{A}\rvert N)$.

Computation of $\textsc{Best}_{\boldsymbol{\lambda}}(\widehat{Q}_{t})$. The best response of the $\boldsymbol{\lambda}$-player is used in Step 1 to determine the suboptimality of the current solution. Given an arbitrary $Q$, $\textsc{Best}_{\boldsymbol{\lambda}}(Q)$ returns $\boldsymbol{\lambda}$ maximizing $L(Q,\boldsymbol{\lambda})$ over $\boldsymbol{\lambda}\geq\mathbf{0},\lVert\boldsymbol{\lambda}\rVert_{1}\leq B$. By first-order optimality, the optimizing $\boldsymbol{\lambda}$ is either $\mathbf{0}$ or puts all of its mass on the most violated constraint among $\widehat{\gamma}_{a,z}(Q)\leq\widehat{\varepsilon}_{a}$, $\widehat{\gamma}_{a,z}(Q)\geq-\widehat{\varepsilon}_{a}$. In particular, let $\mathbf{e}_{a,z}^{+}$ and $\mathbf{e}_{a,z}^{-}$ denote the basis vectors corresponding to coordinates $\lambda_{a,z}^{+}$ and $\lambda_{a,z}^{-}$. The call to $\textsc{Best}_{\boldsymbol{\lambda}}(Q)$ first calculates

Thus, $\textsc{Best}_{\boldsymbol{\lambda}}(\widehat{Q}_{t})$ can be calculated in time $O(\lvert\mathcal{A}\rvert N)$ as long as we have $\widehat{\boldsymbol{\gamma}}(\widehat{Q}_{t})$, whose computation we have already accounted for within the computation of the Lagrangian of the form $L(\widehat{Q}_{t},\boldsymbol{\lambda})$.

By definition of $\widehat{\gamma}_{a,z}(h)$, we have

In particular, Eqs. (10) and (37) imply that the minimization of $L(h,\boldsymbol{\lambda})$ is indeed equivalent to minimizing the right-hand side of Eq. (12) as claimed.

F.2 Details for Reduction to Least-squares Regression

The construction of the least-squares regression data set begins with Eq. (12), which states that for $h_{f}\in\mathcal{H}$

Plugging this back into Eq. (38), we see that the minimization of Eq. (38) over $h\in\mathcal{H}$ is equivalent to the minimization of the empirical loss under $g_{\boldsymbol{\lambda}}$ among $f\in\mathcal{F}$:

Solving this problem directly seems to require access to a generic optimization oracle. We instead use a heuristic, where we first pick

and then seek to solve the least-squares regression problem

Hence, we create two weighted examples for each $(X_{i},A_{i},Y_{i})$ triple in our dataset, and solve

Appendix G Additional Experimental Results

In this section, we include further details on our experimental evaluation.

In Figure 2 we include the training performances of our algorithm and the baseline methods, including the SEO method and the unconstrained regressors. Our method generally dominated or closely matched the baseline methods. The SEO method provided solutions that were not Pareto optimal on the law school data set.

Implementation of the cost-sensitive oracle.

Given an instance of cost-sensitive classification problem, CS oracle optimizes the equivalent weighted binary classification problem on the data $\{W_{i},X_{i}^{\prime},Y_{i}\}_{i=1}^{n}$ with each $X_{i}^{\prime}=(x_{i},z_{i})$ (see Section 3 for the transformation). The oracle aims to solve

where every function $h$ in the class $\mathcal{H}$ is parameterized by a vector $\beta$ and defined as $h(x,z)=\mathbf{1}\left\{\langle\beta,x\rangle\geq z\right\}$ for any input $(x,z)$. Instead of optimizing over the objective in (39), we will consider the following minimization problem with hinge loss. It will be convenient to consider the labels $Y_{i}$ take values $\{\pm 1\}$ and each predictor $h$ predicts in $\{\pm 1\}$. Then the optimization becomes:

In our experiments, this optimization problem was solved with the Gurobi Optimizer (Gurobi Optimization, 2018).

Runtime comparison.

We performed a comparison on the running time of a single call of the three supervised learning oracles. On a subsampled law school data set with 1,000 examples, we ran the oracles to solve an instance of the $\textsc{Best}_{h}$ problem, optimizing over either the linear models or tree ensemble models. The details are listed in Table 1. We also compare the number of oracle calls for different specified values of fairness slackness.