This paper considers fair probabilistic classification where the outputs of primary interest are predicted probabilities, commonly referred to as scores. We formulate the problem of transforming scores to satisfy fairness constraints that are linear in conditional means of scores while minimizing the loss in utility. The formulation can be applied either to post-process classifier outputs or to pre-process training data, thus allowing maximum freedom in selecting a classification algorithm.
We derive a closed-form expression for the optimal transformed scores and a convex optimization problem for the transformation parameters. In the population limit, the transformed score function is the fairness-constrained minimizer of cross-entropy with respect to the optimal unconstrained scores.
In the finite sample setting, we propose to approach this solution using a combination of standard probabilistic classifiers and ADMM. The transformation parameters obtained from the finite-sample procedure are shown to be asymptotically optimal. Comprehensive experiments comparing to 10 existing methods show that the proposed has advantages for score-based metrics such as Brier score and AUC while remaining competitive for binary label-based metrics such as accuracy.

§ INTRODUCTION

Recent years have seen a surge of interest in the problem of
fair classification, which is concerned with disparities in classification output or performance
when conditioned on a protected attribute such as race or gender, or ethnicity.
Many measures of fairness have been introduced [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] and fairness-enhancing interventions have been proposed to mitigate these disparities
[15]. Roughly categorized, these interventions either (i) change data used to train
a classifier (pre-processing) [16, 17, 18, 19, 20], (ii) change a classifier's output
(post-processing) [21, 22, 4, 23, 24], or (iii) directly change
a classification model to ensure fairness (in-processing) [25, 26, 5, 27, 28, 29, 30, 31, 32].

This paper differs from many of the above works in placing more emphasis on probabilistic classification in which the outputs of interest are predicted probabilities of belonging to one of the classes as opposed to binary predictions. The predicted probabilities are often referred to as scores and are desirable because they indicate confidences in predictions.
We propose an optimization formulation for  transforming scores to satisfy fairness constraints while minimizing the loss in utility. The formulation accommodates any fairness criteria that can be expressed as linear inequalities involving conditional means of scores, including variants of statistical parity (SP) [1] and equalized odds (EO) [4, 5].

We make the following theoretical and methodological contributions beyond a novel problem formulation. We derive a closed-form expression for the optimal transformed scores and a convex dual optimization problem for the Lagrange multipliers that parametrize the transformation. In the population limit, the transformed scores minimize cross-entropy with respect to the conditional distribution $p_{Y\given X}$ of the outcome $Y$ given features $X$, in other words the unconstrained optimal score, subject to the fairness constraints. In the finite sample setting, we propose a method called (FST) that uses standard probabilistic classifiers (e.g. logistic regression) to approximate $p_{Y\given X}$ and the alternating direction method of multipliers (ADMM) to solve the dual problem. The closed-form expression for the transformed scores and the low dimension of the dual problem (a small multiple of the number of protected groups) make FST computationally lightweight. Furthermore, we prove that the Lagrange multiplier parameters obtained by FST are asymptotically optimal.

FST lends itself naturally to post-processing and can also be applied in pre-processing.  As such, we envision that FST will be particularly beneficial in situations that make post- and pre-processing attractive,
as articulated by several authors
[17, 20, 29, 33, 34]:
a) when it is not possible or desirable to modify an existing classifier (only post-processing is possible); b) when freedom is desired to select the most suitable classifier for an application, whether it maximizes performance or has some other desired property such as interpretability (post- and pre-processing apply); and c) when standard training algorithms are used without the additional complexity of accounting for fairness (post- and pre-processing again). In-processing meta-algorithms [29, 32] can also support situation b) but not a) or c).

Pre-processing is further motivated by its attempt to address the problem at its source, namely historical bias and imbalance in the data, and also allows corrected data to be published as an output in its own right for use by others. Compared to existing post- and pre-processing methods, FST is considerably more flexible in handling more cases (see summary in Table <ref>).

We have conducted comprehensive experiments comparing FST to 10 existing methods, a number that compares favorably to recent meta-studies [15]. On score-based metrics such as Brier score and AUC, FST achieves better fairness-utility trade-offs and hence is indeed advantageous when scores are of interest. At the same time, it remains competitive on binary label-based metrics such as accuracy.

In summary, it is shown that enables fairness-ensuring post- and pre-processing that

* is theoretically grounded and optimal in the population limit (Sections <ref>, <ref>, and <ref>),

* is computationally lightweight (Section <ref>),

* performs favorably compared to the state-of-the-art (Section <ref> and  Supplementary Material).

§.§ Related work

Existing post-processing methods for fairness
include [21, 22, 4, 23, 35]; limitations of post-processing are studied in [24]. While these methods take predicted scores
as input, most [21, 22, 4] are designed to
produce only binary output and not scores. The method of [23]
maintains calibrated probability estimates, which is
a requirement that we do not enforce herein. Furthermore, [21, 22, 4, 23] all assume knowledge of the protected attribute at test time, and [21, 22, 35] address only SP ([21] as originally proposed)
while [4, 23] address disparities in error rates. Our approach does not have these limitations.

Pre-processing methods range from reweighing, resampling, and relabeling training data [16], to performing probability transformations on features [19], to modifying both labels and features through optimization [20] or labels and protected attributes using classification rules [17]. The above methods

only address SP or the related notion of disparate impact [19]. Learning representations that are invariant to protected attributes [18, 36, 37, 38, 39] can also be seen as pre-processing,
and recent adversarial approaches [40, 41, 33] permit control of EO as well as SP. Representation learning however does not preserve the original data domain and its semantics, while adversarial algorithms
can be computationally challenging.

Several works [29, 32, 42, 43, 44] have technical similarities to the approach herein but
focus on binary outputs, with $0$-$1$ risk [32, 29] or cost-sensitive risk [42, 43] as the objective function, and/or lead to in-processing algorithms. The closest of these is [32], which also solves a fairness-constrained classification problem via the dual problem. However, [32] along with [29] propose in-processing algorithms that solve multiple instances of a subproblem whereas we solve only one instance. [32] also addresses a larger class of fairness measures that are linear-fractional in the classifier output.
Similar to us, [42, 43] also characterize optimal fair classifiers in the population limit in which probability distributions are known; however, they do not propose algorithms for computing the Lagrange multipliers or thresholds that parametrize
the solution.

§ PROBLEM FORMULATION

We represent one or more protected attributes such as gender and race by a random variable $A$ and an outcome variable by $Y$. We make the common assumption that $Y \in \{0,1\}$ is binary-valued.
It is assumed that $A$ takes a finite number of values in a set $\cA$, corresponding to protected groups.  Let $X$ denote features used to predict $Y$ in a supervised classification setting.  We consider two scenarios in which $X$ either includes or does not include $A$, like in other works in fair classification
(e.g. [16, 29, 31]).  While it is recognized that
the former scenario can achieve better trade-offs between utility and
fairness, the latter is needed in applications where disparate treatment laws and regulations forbid the explicit use of $A$.  To develop our approach in this section and Section <ref>, we work in the population limit and make use of probability distributions involving $A$, $X$, $Y$.  Section <ref> discusses how these distributions are approximated using a training sample.

As stated earlier, we focus more heavily on probabilistic classification in which the output of interest is the predicted probability of being in the positive class $Y=1$ rather than a binary prediction.  The optimal probabilistic classifier is the conditional probability $r(x) \equiv p_{Y\given X}(1\given x)$, which we refer to as the
original score.  Bayes-optimal binary classifiers can be derived from $r(x)$ by thresholding, specifically at level $c \in [0,1]$ if $c$ and $1-c$ are the relative costs of false positive and false negative errors.  Score functions will thus play the central role in our development.

We propose a mathematical formulation and method called  (FST) that can be applied to both post-processing and pre-processing.  In both cases, the goal is to transform $r(x)$ into a transformed score $r'(x)$ that satisfies fairness conditions while minimizing the loss in optimality compared to $r(x)$.  We elaborate on the utility and fairness measures considered in Sections <ref> and <ref>.

The application of FST to post-processing is straightforward: $r'(x)$ is used directly as the classification output and
can be thresholded to provide a binary prediction.

In the pre-processing case, we additionally define a transformed outcome variable $Y' \in \{0,1\}$ and let $r'(x) = p_{Y'\given X}(1\given x)$ be the conditional probability associated with it.
The overall procedure consists of two steps, performed in general by two different parties: 1) The data owner transforms the outcome variable from $Y$ to $Y'$; 2) The modeler trains a classifier with $Y'$ as target variable and $X$ as input, without regard for fairness.  The transformed score $r'(x)$ plays two roles in this procedure.  The first is to specify the (randomized) mapping from $X$ to $Y'$ in step 1).  As will be seen, the mapping depends only indirectly on $Y$ through the original score $r(x)$.  The second role stems from the main challenge faced by pre-processing methods, namely that the predominant fairness metrics depend on the output of the classifier trained in step 2) but this classifier is not under direct control of the pre-processing in step 1).
In recognition of this challenge, we make the following assumption, also discussed in [33, 34]:

The classifier trained by the modeler approximates
the transformed score $r'(x)$ if it is a probabilistic classifier or a thresholded version of $r'(x)$ if it is a binary classifier.

This assumption is satisfied for modelers who are “doing their job” in learning to predict $Y'$ from $X$ since the optimal classifier in this case is $r'(x)$ or a function thereof.  Given the assumption, we will use $r'(x)$ as a surrogate for the actual classifier output.  The assumption is not satisfied if the modeler is not competent or, worse, malicious in trying to discriminate against certain protected groups.

§.§ Utility measure

We propose to measure the loss in optimality, i.e. utility, between the transformed score $r'(x)$ and original score $r(x)$ using the following cross-entropy:

\begin{equation}\label{eqn:utility}
\mathbb{E}\bigl[-\log p_{Y'\given X}(Y\given X)\bigr] = \mathbb{E} \bigl[ -r(X) \log r'(X) - (1-r(X)) \log(1 - r'(X)) \bigr],
\end{equation}

where the right-hand side results from expanding the expectation over $Y$ conditioned on $X$, and $p_{Y'\given X}$ is used only as notational shorthand in the post-processing case since $Y'$ is not generated. For simplicity, we shall also use the following notation for cross-entropy:

\begin{equation}\label{eqn:crossentropy}
H_b(p,q )\defined -p\log q-(1-p)\log(1-q).
\end{equation}

The utility measure in (<ref>) is equivalent to $\EE{H_b\left(r(X),r'(X)\right)}.$

One way to arrive at (<ref>)
is to assume that
$r'(x)$, which is the classifier output in the post-processing case and a surrogate thereof in the pre-processing case, is evaluated against the original outcomes $y_1,\dots,y_n$ in a training set using the cross-entropy a.k.a. log loss.  This yields the empirical version of the left-hand side of (<ref>),

\[
-\frac{1}{n} \sum_{i=1}^n \log p_{Y'\given X}(y_i\given x_i).
\]

The use of log loss is well-motivated by the desire for $r'(x)$ to be close to the true conditional probability $r(x)$.

An equivalent way to motivate (<ref>) in the pre-processing context is to measure the utility lost in transformation by
the Kullback-Leibler (KL) divergence between the original and transformed joint
distributions $p_{X,Y}$, $p_{X,Y'}$:

\begin{equation}
\label{eqn:utilityKL}
D_{\mathrm{KL}}\bigl(p_{X,Y} \KLsep p_{X,Y'}\bigr) = \mathbb{E}_{p_{X,Y}} \left[ \log\frac{p_{X,Y}}{p_{X,Y'}} \right]= \mathbb{E}_{p_{X,Y}}[\log p_{Y\given X}] - \mathbb{E}_{p_{X,Y}}[\log p_{Y'\given X}].
\end{equation}

The first term depends on the data distribution but not $r'(x)$ and the second term is exactly (<ref>).

Starting from a different premise, [44] proposed a similar mathematical formulation in which the arguments of the KL divergence are reversed from those in (<ref>), i.e. the given distribution is the second argument while the distribution to be determined is the first.
The form of the solution in [44] is therefore different from the one presented herein.  The order of arguments in (<ref>) is justified by the connection to log loss in classification discussed above.

§.§ Fairness measures

We consider fairness criteria expressible
as linear inequalities involving
conditional means of scores,

\begin{equation}\label{eqn:fairnessLinIneq}
\sum_{j=1}^J b_{lj} \mathbb{E}\bigl[ r'(X) \given \cE_{lj} \bigr] \leq c_l, \quad l = 1,\dots,L,
\end{equation}

where $\{b_{lj}\}$ and $\{c_l\}$ are real-valued coefficients
and the conditioning events $\cE_{lj}$ are defined in terms of   $(A, X, Y)$ but do not depend on $r'$.
Special cases of (<ref>) correspond to the well-studied notions of statistical parity (SP) and equalized odds (EO).

More precisely, we focus on the following variant of SP:

\begin{equation}\label{eqn:MSP}
-\epsilon \leq \mathbb{E}[r'(X)\given A=a] - \mathbb{E}[r'(X)] \leq \epsilon \qquad \forall a \in \cA,
\end{equation}

which we refer to as mean score parity (MSP) following [45]. Condition (<ref>) corresponds to approximate mean independence of random variable $R' = r'(X)$ with respect to $A$.

Similar notions can also be put in the form of (<ref>), for example bounds on the ratio

\[
1-\epsilon \leq \frac{\mathbb{E}[r'(X)\given A=a]}{\mathbb{E}[r'(X)]} \leq 1+\epsilon,
\]

referred to as disparate impact in [19], as well as conditional statistical parity [3, 43].

For EO, we add the condition $Y=y$ to the conditioning events in (<ref>), resulting in
\begin{equation}\label{eqn:MEO}
-\epsilon \leq \mathbb{E}[r'(X)\given A=a, Y=y] - \mathbb{E}[r'(X) \given Y=y] \leq \epsilon \qquad \forall a \in \cA, \; y \in \{0,1\}.
\end{equation}

For $y = 0$ (respectively $y = 1$), $\mathbb{E}[r'(X) \given Y=y]$ is the false (true) positive rate (FPR, TPR) generalized for a probabilistic classifier, and $\mathbb{E}[r'(X) \given A=a, Y=y]$ is the corresponding group-specific rate.  Following [23], we refer to (<ref>) for $y=0$ or $y=1$ alone as approximate equality in generalized FPRs or TPRs, and to (<ref>) for $y=0$ and $y=1$ together as generalized EO (GEO).

For later use in proving Proposition <ref>, we specify the exact correspondences between (<ref>), (<ref>) and (<ref>). The MSP constraint (<ref>) can be obtained from (<ref>) by setting $J = 2$, $l = (a,\pm)$ for $a\in\cA$ where $+$ corresponds to the $\leq \epsilon$ constraint and $-$ to the $\geq -\epsilon$ constraint, $L = 2\abs{\cA}$, $\cE_{(a,\pm),1} = \{A = a\}$, $\cE_{(a,\pm),2} = \Omega$ (the entire sample space), $c_l = \epsilon$, and $b_{(a,\pm),j} = \mp (-1)^j$.
For the GEO constraint (<ref>), set $J = 2$, $l = (a,y,\pm)$ for $a\in\cA$, $y\in\{0,1\}$ and the same $\pm$ correspondences, $L = 4\abs{\cA}$, $\cE_{(a,y,\pm),1} = \{A = a, Y = y\}$, $\cE_{(a,y,\pm),2} = \{Y = y\}$, $c_l = \epsilon$, and $b_{(a,y,\pm),j} = \mp (-1)^j$.

The fairness measures (<ref>) in our formulation are defined in terms of probabilistic scores.  Parallel notions defined for binary predictions, i.e. by replacing $r'(X)$ with a thresholded version $\ones(r'(X) > t)$, are more common in the literature.  For example, the counterpart to (<ref>) is (non-generalized) EO while the counterpart to (<ref>) is called thresholded score parity in [45].  While our formulation does not optimize for these binary prediction measures, we nevertheless use them for evaluation in Section <ref>.

The form of (<ref>) is inspired by but is less general than the linear conditional moment
constraints in [29], which replace $r'(X)$ in (<ref>) by
an arbitrary bounded function $g_j(A, X, Y, r'(X))$.
We have restricted ourselves to (<ref>) so that a closed-form optimal solution can be derived
in Section <ref>.  We note however that in both of
the examples in [29] and many
fairness measures, $g_j(A, X, Y, r'(X)) = r'(X)$ and the additional generality is not required.

§.§ Optimization problem

The transformed score $r'(x)$ is obtained by minimizing the cross-entropy in (<ref>) (equivalently maximizing its negative) subject to fairness constraints (<ref>):

\begin{equation}\label{eqn:primalGeneral}
%    \begin{split}
\max_{r'} \;\; -\EE{ H_b\left(r(X),r'(X)\right) } \quad
\text{s.t.} \quad \sum_{j=1}^J b_{lj} \mathbb{E}\bigl[ r'(X) \given \cE_{lj} \bigr] \leq c_l, \quad l = 1,\dots,L.
%    \end{split}
\end{equation}

The next section characterizes the optimal solution to this  problem.

On the sufficiency of pre-processing scores.

In the pre-processing case, the proposed optimization (<ref>) transforms only scores and uses them to generate a weighted dataset, as described further in Section <ref>. Can a better trade-off between utility and fairness  be achieved by also pre-processing features $X$, i.e., mapping each pair $(x,r(x))$ into a new $(x',r'(x))$? Note that pre-processing both scores/labels and input features is suggested in [17, 19, 20]. When utility and fairness are measured according to the objective and constraints in (<ref>), the answer is negative, since both the objective and the constraints only depend on the marginals of $r'(X)$ on events given in terms of $A$ and $Y$. Thus, for the metrics considered here, pre-processing the scores is sufficient.

§ CHARACTERIZATION OF OPTIMAL FAIRNESS-CONSTRAINED SCORE

We derive a closed-form expression for the optimal solution to problem (<ref>) using the method of Lagrange multipliers.  We then state the dual optimization problem that determines the Lagrange multipliers.  These results are specialized to the cases of MSP (<ref>) and GEO (<ref>).

Define Lagrange multipliers $\lambda_l \geq 0$, $l = 1,\dots,L$ for the constraints in (<ref>), and let $\lambda \defined (\lambda_1,\dots,\lambda_L)$.  Then the Lagrangian function is given by

\begin{equation}\label{eqn:Lagrangian1}
L(r', \lambda) = -\mathbb{E} \bigl[ H_b\left(r(X),r'(X)\right) \bigr] - \sum_{l=1}^L \sum_{j=1}^J \lambda_l b_{lj} \mathbb{E}\bigl[ r'(X) \given \cE_{lj} \bigr] + \sum_{l=1}^L c_l \lambda_l.
\end{equation}

The dual optimization problem corresponding to (<ref>) is

\[
\min_{\lambda \geq 0} \max_{r'} L(r', \lambda).
\]

Note that $L(r',\lambda)$ is a strictly concave function of $r'$ and the fairness constraints in (<ref>) are affine functions of $r'$. Consequently, as long as the constraints in (<ref>) are feasible, the optimal transformed score $r^*$ can be found by (i) maximizing $L(r',\lambda)$ with respect to $r'$, resulting in an optimal solution $r^*$ that is a function of $\lambda$, and then (ii) minimizing $L(r^*,\lambda)$ with respect to $\lambda$ <cit.>. Substituting the optimal $\lambda^*$ into the solution for $r^*$ found in the first step then yields the optimal transformed score. Note that this procedure would not necessarily be correct if a linear objective function were considered (e.g., 0-1 loss in [32]) due to lack of strict concavity.  The next proposition states the general form of the solution to the inner maximization (i) above. Its proof is in Appendix <ref>.

Let $L(r',\lambda)$ be as given in (<ref>). Then for fixed $\lambda$,

$r^*(\lambda) = \arg \max_{r'} L(r',\lambda)$

is given by

\begin{equation}\label{eqn:r*}
r^*\bigl(\mu(x); r(x)\bigr) = \begin{cases}
\dfrac{1 + \mu(x) - \sqrt{(1 + \mu(x))^2 - 4 r(x) \mu(x)}}{2\mu(x)}, & \mu(x) \neq 0\\
r(x), & \mu(x) = 0,
\end{cases}
\end{equation}

\begin{equation}\label{eqn:mu}
\mu(x) \defined \sum_{l=1}^L \sum_{j=1}^J \lambda_l b_{lj} \frac{\Pr(\cE_{lj} \given X=x)}{\Pr(\cE_{lj})}.
\end{equation}

We can interpret the optimal primal solution (<ref>) as a prescription for score transformation controlled by
$\mu(x)$, which is in turn a linear function of $\lambda$. When $\mu(x) = 0$, the score is unchanged from the original $r(x)$, and as $\mu(x)$ increases or decreases away from zero, the score $r^*(\mu(x); r(x))$ decreases or increases smoothly from $r(x)$ (as can be seen by plotting the function).  It can also be shown that the transformed score
$r^*$ has a rank-preserving property
for fixed $\mu$ in the sense that if $r_1 < r_2$ then $r^*(\mu; r_1) < r^*(\mu; r_2)$.  This can be shown by obtaining the partial derivative of $r^*$ with respect to $r(x)$ and noting that it is non-negative, i.e. $r^*$ is an increasing function of $r(x)$ for fixed $\mu(x)$.

It is shown in Appendix <ref> that the result of substituting the optimal primal solution (<ref>) into the first two terms of the Lagrangian (<ref>) is the expectation of the function

\begin{equation}
\label{eqn:g}
g\bigl(\mu(x); r(x)\bigr) \defined -H_b\Bigl(r(x),r^*\left(\mu(x);r(x) \right) \Bigr)  - \mu(x) r^*\bigl(\mu(x); r(x)\bigr).
\end{equation}

The dual problem is therefore

\begin{equation}\label{eqn:dualGeneral}
\begin{split}
\min_{\lambda} \quad &\mathbb{E}\left[ g\bigl(\mu(X); r(X)\bigr) \right] + \sum_{l=1}^L c_l \lambda_l \\
\text{s.t.} \quad &\mu(X) = \sum_{l=1}^L \sum_{j=1}^J \lambda_l b_{lj} \frac{\Pr(\cE_{lj} \given X)}{\Pr(\cE_{lj})}, \qquad \lambda \geq 0.
\end{split}
\end{equation}

The solution to the above minimization provides the values of $\lambda^*$ for the optimal transformed score (<ref>). Like all Lagrangian duals, (<ref>) is a convex optimization (although it is no longer apparent from (<ref>) that this is the case).  Furthermore, (<ref>) is typically low-dimensional in cases where the number of dual variables $L$ is a small multiple of the number of protected groups $\abs{\cA}$.

We now specialize and simplify (<ref>) to MSP (<ref>) and GEO (<ref>) fairness constraints. The following proposition follows from the correspondences between (<ref>), (<ref>) and (<ref>) identified in Section <ref> and is proved in Appendix <ref>.

Under the MSP constraint (<ref>), the dual optimization (<ref>) reduces to

\begin{equation}\label{eqn:dualMSP}
\begin{split}
\min_{\lambda} \quad &\mathbb{E}\left[ g\bigl(\mu(X); r(X)\bigr) \right] + \epsilon \norm{\lambda}_1\\
\text{s.t.} \quad &\mu(X) = \sum_{a\in\cA} \lambda_a \left(\frac{p_{A\given X}(a\given X)}{p_A(a)} - 1\right).
\end{split}
\end{equation}

For the GEO constraint (<ref>), (<ref>) reduces to

\begin{equation}\label{eqn:dualMEO}
\begin{split}
\min_{\lambda} \quad &\mathbb{E}\left[ g\bigl(\mu(X); r(X)\bigr) \right] + \epsilon \norm{\lambda}_1, \\
\text{s.t.} \quad &\mu(X) = \sum_{y\in\{0,1\}} \frac{p_{Y\given X}(y\given X)}{p_{Y}(y)} \sum_{a\in\cA} \lambda_{a,y} \left( \frac{p_{A\given X,Y}(a\given X,y)}{p_{A\given Y}(a\given y)} - 1\right).
\end{split}
\end{equation}

In (<ref>), (<ref>), there is no longer a non-negativity constraint on $\lambda$ but instead an $\ell_1$ norm, and the problem dimension is only $|\cA|$
in (<ref>) and $2|\cA|$ in (<ref>).

Moreover, both dual formulations are well-suited for decomposition using the alternating direction method of multipliers (ADMM), as discussed further in Section <ref>.

In the case where the features $X$ include the protected attribute $A$, we have $p_{A\given X}(a\given X) = p_{A\given X,Y}(a\given X,y) = \ones(a = A)$, where $A$ is the component of $X$ that is given. The constraints in (<ref>) and (<ref>) then simplify to

\begin{equation}\label{eqn:muMSPAtrue}
\mu(X) = \frac{\lambda_A}{p_A(A)} - \sum_{a\in\cA} \lambda_a,
\end{equation}

\begin{equation}\label{eqn:muMEOAtrue}
\mu(X) = \sum_{y\in\{0,1\}} \frac{p_{Y\given X}(y\given X)}{p_{Y}(y)} \left( \frac{\lambda_{A,y}}{p_{A\given Y}(A\given y)} - \sum_{a\in\cA} \lambda_{a,y} \right)
\end{equation}

respectively. Interestingly, the only difference between the cases of including or excluding $A$ is that in the latter, the constraints in (<ref>), (<ref>) indicate that $A$ should be inferred from the available features $X$ and possibly $Y$, whereas in the former, $A$ can be used directly.

§ PROPOSED FAIRSCORETRANSFORMER PROCEDURE

We now consider the finite sample setting in which the probability distributions of $A,X,Y$ are not known and we have instead a training set $\mathcal{D}_n\defined \{(a_i, x_i, y_i), i=1,\dots,n\}$.  This section presents the proposed
(FST) procedure that approximates
the optimal fairness-constrained score in Section <ref>. We focus on the cases of MSP and GEO. The procedure consists of the following steps:

* Estimate the original score and other probabilities required to define the dual problem (<ref>) or (<ref>).

* Solve the dual problem to obtain dual variables $\lambda^*$ (the “fit” step).

* Transform scores using (<ref>) and (<ref>) (“transform” step).

* For pre-processing, modify the training data.

* For binary-valued predictions, binarize scores.

The following subsections elaborate on steps 1, 2, and 4. Step 5 is done simply by selecting a threshold $t \in [0,1]$ to maximize accuracy on the training set.

§.§ Estimation of original score and other probabilities

In some post-processing applications, original scores $r(x)$
may already be estimated by an existing base classifier.  If no suitable base classifier exists, any probabilistic classification algorithm may be used to estimate $r(x)$.  We experiment with logistic regression and gradient boosting machines in Section <ref>.  We naturally recommend selecting a model and any hyperparameter values to maximize performance in this regard, i.e. to yield accurate and calibrated probabilities.  This can be done through cross-validation on the training set using an appropriate metric such as Brier score [47].

In the case where $A$ is one of the features in $X$, the other probabilities required are $p_A(a)$ for MSP (<ref>) and $p_Y(y)$, $p_{A\given Y}(a\given y)$ for GEO (<ref>) ($p_{Y\given X}(y\given x)$ is already given by $r(x)$ and $p_{A\given X}$, $p_{A\given X,Y}$ are delta functions).  Since $Y$ is binary and $\abs{\cA}$ is typically small, it suffices to use the empirical estimates of these probabilities.  If $A$ is not included in $X$, then it is also necessary to estimate it using $p_{A\given X}(a\given X)$ for MSP (<ref>) and $p_{A\given X,Y}(a\given X,y)$ for GEO (<ref>).  Again, any probabilistic classification algorithm can be used, provided that it can handle more than two classes if $\abs{\cA} > 2$.

We highlight that FST
translates the effort of ensuring fair classification into training well-calibrated models for predicting $Y$ and, if necessary, $A$. This echoes the plug-in approach advocated by [42].

§.§ ADMM for optimizing dual variables

In the finite sample case, we solve an empirical version of the dual problem in Proposition <ref>. We write $\mu(x) = \lambda^T \bf(x)$, where $\mathbf{f}:\mathcal{X}\to \mathbb{R}^d$ is defined by the expression for $\mu(x)$ in (<ref>) or (<ref>), and $d$ is the dimension of $\lambda$. Let $\hat{r}(x)$ denote the estimate of $r(x)$ from Section <ref>, and $\hat{\bf}(x)$ be an empirical version of $\bf(x)$ in which all probabilities (e.g. $p_A(a)$ for MSP (<ref>)) are replaced by their estimates. With these definitions, both optimizations in Proposition <ref> have the general form

\begin{equation}
\label{eqn:dualSpecF}
\min_{\lambda\in\mathbb{R}^d} \quad \frac{1}{n}\sum_{i=1}^n g\bigl(\mu(x_i); \hat{r}(x_i)\bigr)  + \epsilon \norm{\lambda}_1 \quad
\text{s.t.} \quad \mu(x_i) = \lambda^T\hat{\mathbf{f}}(x_i), \quad i=1,\dots,n,
\end{equation}

where the expectation in the objective has also been approximated by the average over the training dataset. To simplify notation in this subsection, we suppress the hat symbols on $\hat{r}(x)$ and $\hat{\bf}(x)$.

Formulation (<ref>) is well-suited for ADMM because the objective function is separable between $\mu(x)$ and $\lambda$, which are linearly related through the constraint. We present one ADMM decomposition here and alternatives in Appendix <ref>. Under the first decomposition, application of the scaled ADMM algorithm <cit.> to (<ref>) yields the following three steps in each iteration $k = 0,1,\dots$:

\begin{align}
\mu^{(k+1)}(x_i) &= \arg\min_{\mu} \; \frac{1}{n} g\bigl(\mu; r(x_i)\bigr) + \frac{\rho}{2}\left(\mu- (\lambda^{(k)})^T\mathbf{f}(x_i)+c^{(k)}(x_i) \right)^2 \quad \forall i=1,\dots,n \label{eq:iter-mu}\\
\lambda^{(k+1)} &= \arg\min_{\lambda}~ \epsilon\|\lambda \|_1+\frac{\rho}{2}\sum_{i=1}^n\left(\mu^{(k+1)}(x_i)- \lambda^T\mathbf{f}(x_i) +c^{(k)}(x_i)\right)^2 \label{eq:iter-lambda}\\
c^{(k+1)}(x_i) &= c^{(k)}(x_i)+\mu^{(k+1)}(x_i)- \left(\lambda^{(k+1)}\right)^T\mathbf{f}(x_i) \quad \forall i=1,\dots,n.\label{eq:iter-c}
\end{align}

The first update (<ref>)  can be computed in parallel for each sample $x_i$ in the dataset. Given an $x_i$, finding $\mu(x_i)$ is a single-parameter optimization where the objective possesses   closed-form expressions for its derivatives. For simplicity of notation, let $r_i \defined r(x_i)$, $a_i\defined (\lambda^{(k)})^T\mathbf{f}(x_i)+c(x_i)$,   and
$$\mathsf{obj}(\mu)\defined \frac{1}{n}g\bigl(\mu; r_i\bigr) + \frac{{\rho}}{2}\left(\mu- a_i \right)^2.$$
The first two derivatives of $\mathsf{obj}(\mu)$ are

\begin{align*}
\frac{\partial \mathsf{obj}(\mu)}{\partial \mu} = -r^*(\mu;r_i)+{\rho}(\mu-a_i),\qquad
\frac{\partial^2 \mathsf{obj}(\mu)}{\partial \mu^2} &
\begin{cases}
-\frac{1}{2\mu^2}\left(1-\frac{1+\mu(1-2r_i)}{\sqrt{(1+\mu)^2-4r_i\mu}} \right)+{\rho},& \mu\neq 0, \\
\end{cases}%\label{eq:step1}
\end{align*}

The above expressions are useful whether (<ref>) is solved by a generic optimization solver or a root-finding method (e.g., Newton's method).

The second update (<ref>) reduces to an $\ell_1$-penalized quadratic minimization over (at most) $2|A|$ variables. Specifically,

\begin{align}
\lambda^{(k+1)} &= \arg\min_{\lambda}~ \epsilon\|\lambda \|_1+\lambda^T\mathbf{v} +\lambda^T\mathbf{F}\lambda, \label{eqn:step2}
\end{align}

\begin{equation}
\mathbf{v} \defined -\rho\sum_{i=1}^n\mathbf{f}(x_i)\left(\mu^{(k+1)}(x_i)+c^{(k)}(x_i)\right),\qquad \mathbf{F}\defined \frac{\rho}{2}\sum_{i=1}^n\mathbf{f}(x_i)\mathbf{f}(x_i)^T.
\end{equation}

The values of $\mathbf{v}$ and $\mathbf{F}$ above can be pre-computed prior to solving (<ref>). In fact, $\mathbf{F}$ can be computed once at the start of the iterations. The ensuing minimization only involves $|\cA|$ variables under the MSP constraint (<ref>), and $2|\cA|$ variables under the GEO constraint (<ref>).

From (<ref>)–(<ref>), it is seen that the computational complexity of each ADMM iteration scales linearly with
$n$. We have fixed the ADMM penalty parameter $\rho = 1$ and have not attempted to tune it for faster convergence.

§.§ Additional steps for pre-processing

In the pre-processing scenario, the transformed score $r'(x)$ is used to generate samples of a transformed outcome
$Y'$.  Since $r'(x) = p_{Y'\given X}(1\given x)$ is a probabilistic mapping, we propose generating a weighted dataset $\cD' = \{(x_i, y'_i, w_i)\}$ with weights $w_i$ that reflect the conditional distribution $p_{Y'\given X}$.  Specifically, $\cD' = \cD'_0 \cup \cD'_1$ with $\cD'_0 = \{(x_i, 0, 1-r'(x_i)), i=1,\dots,n\}$ and $\cD'_1 = \{(x_i, 1, r'(x_i)), i=1,\dots,n\}$ so that $\cD'$ is twice the size of the original dataset.  The data owner passes the transformed dataset
$\cD'$ to the modeler, who uses it to train a classifier for $Y'$ given $X$
without fairness constraints.

§ CONSISTENCY OF FAIRSCORETRANSFORMER

In this section, we present a result toward guaranteeing the consistency of the FST procedure of Section <ref>, again focusing on the cases of MSP and GEO. We consider in particular step 2 of the procedure.
Consistency in this case means that

optimal solutions to the empirical dual problem (<ref>) become asymptotically optimal for the population dual problem (<ref>) or (<ref>) as the sample size $n \to \infty$ and the estimates $\hat{r}(x)$, $\hat{\bf}(x)$ converge to their respective true quantities.

We first discuss the assumptions that are made before formally stating and proving the result.

§.§ Assumptions

To simplify the proofs, we assume in this section that $A$ is available at test time, as stated below for easy reference:

The protected attributes $A$ are known at test time.

The proofs can be extended to the case in which $A$ is not known by also assuming a consistent estimator of the conditional probability $p_{A\given X}$ in the MSP case or $p_{A\given X,Y}$ in the GEO case and accounting for the error of this estimator.

We use the following assumption to show that it is sufficient to consider a bounded feasible set for the dual problem.

The fairness constraint parameter $\epsilon > 0$.

To ensure consistency of FST, we naturally assume that $\hat{r}(X)$ is a consistent estimator of the original score $r(X)$. The precise sense in which $\hat{r}(X)$ is consistent is stated below.
We view $\hat{r}(X)$ and $r(X)$ as random variables induced by $X$ and define $D_{\mathrm{TV}}(\hat{r}(X), r(X))$ to be the total variation distance between them.

The estimate $\hat{r}(X)$ converges to the original score $r(X)$ in distribution. More specifically,

\[
D_{\mathrm{TV}}\bigl(\hat{r}(X), r(X)\bigr) \overset{p}{\to} 0.
\]

Lastly, we state assumptions regarding the probabilities $p_A(a)$ (MSP case) and $p_{A,Y}(a,y)$ (GEO case) together with their estimates.

* For the MSP case, $p_A(a)$ and its estimate $\hat{p}_A(a)$ are bounded away from zero, i.e. $p_A(a) \geq \delta$ and $\hat{p}_A(a) \geq \delta$ for all $a \in \cA$ and some $\delta > 0$. For the GEO case, $p_{A,Y}(a,y) \geq \delta$ and $\hat{p}_{A,Y}(a,y) \geq \delta$ for all $a \in \cA$, $y \in \{0,1\}$, and some $\delta > 0$.

* For the MSP case, $\abs{\hat{p}_A(a) - p_A(a)} \overset{p}{\to} 0$ for all $a \in \cA$. For the GEO case, $\abs{\hat{p}_{A,Y}(a,y) - p_{A,Y}(a,y)} \overset{p}{\to} 0$ for all $a \in \cA$ and $y \in \{0,1\}$.

Assumption <ref>.1 is very reasonable in that if a protected group is to be considered, it should represent a constant fraction of the population (and have non-negligible probabilities of being in classes $0$ and $1$). The boundedness of the estimated probabilities can be ensured by truncating them, i.e. setting $\hat{p}_A(a) \leftarrow \max\{\hat{p}_A(a), \delta\}$. Assumption <ref>.2 is hardly an assumption since for finite $\abs{\cA}$, it is satisfied by the empirical probabilities. Furthermore if the true probability is lower bounded by $\delta$, convergence is not affected by truncating the estimated probability.

§.§ Asymptotic dual optimality

We now formalize the statement
regarding asymptotic optimality of empirical dual solutions. Let $J(\lambda)$ and $\hat{J}(\lambda)$ denote the objective functions in the population dual (<ref>), (<ref>) and empirical dual (<ref>) respectively.

Let $\hat{\lambda} \in \argmin \hat{J}(\lambda)$ and $\lambda^* \in \argmin J(\lambda)$ be optimal solutions to the empirical dual problem (<ref>) and population dual problem (<ref>), (<ref>) respectively. Under Assumptions <ref>, <ref>, <ref>, and <ref> and as $n \to \infty$, we have
$J(\hat{\lambda}) \to J(\lambda^*)$, i.e. $\hat{\lambda}$ is asymptotically optimal for the population dual.

We prove the theorem by showing that $\hat{J}(\lambda)$ converges uniformly in probability to $J(\lambda)$. As a consequence, for any deviation $\epsilon > 0$ and any $\delta > 0$, with probability $1-\delta$ we have

\[
J(\hat{\lambda}) \leq \hat{J}(\hat{\lambda}) + \epsilon \leq \hat{J}(\lambda^*) + \epsilon \leq J(\lambda^*) + 2\epsilon,
\]

where the second inequality is by definition of $\hat{\lambda}$.

Toward proving uniform convergence, we first establish that it suffices to solve the dual problem over a closed and bounded (and hence compact) feasible set. The same argument applies to both the population and empirical duals and we will illustrate with the former. Indeed, it always suffices to restrict to a sub-level set defined by the objective value of an initial solution. We take $\lambda = 0$ as the initial solution and define

\begin{equation}\label{eqn:Lambda0}
\Lambda_0 = \{\lambda: J(\lambda) \leq J(0)\}.
\end{equation}

The $\ell_1$ norm of elements in $\Lambda_0$ can be bounded explicitly as proved in Appendix <ref>.

Under Assumption <ref>,

\[
\norm{\lambda}_1 \leq \frac{\log 2}{\epsilon} \quad \forall \, \lambda \in \Lambda_0.
\]

We then consider the supremum over $\Lambda_0$ of the absolute difference $\abs{\hat{J}(\lambda) - J(\lambda)}$ as the quantity of interest for uniform convergence. We use the triangle inequality and separate suprema to decompose this into three terms:

\begin{align}\label{eqn:JunifConv}
\sup_{\lambda\in\Lambda_0} \abs[\big]{\hat{J}(\lambda) - J(\lambda)} \leq& \sup_{\lambda\in\Lambda_0} \abs*{\frac{1}{n} \sum_{i=1}^n g\bigl(\lambda^T \hat{\bf}(a_i, \hat{r}(x_i)); \hat{r}(x_i)\bigr) - \EE{g\bigl(\lambda^T \hat{\bf}(A, \hat{r}(X)); \hat{r}(X)\bigr)}}\nonumber\\
&{} + \sup_{\lambda\in\Lambda_0} \abs*{\EE{g\bigl(\lambda^T \hat{\bf}(A, \hat{r}(X)); \hat{r}(X)\bigr)} - \EE{g\bigl(\lambda^T \bf(A, \hat{r}(X)); \hat{r}(X)\bigr)}}\nonumber\\
&{} + \sup_{\lambda\in\Lambda_0} \abs*{\EE{g\bigl(\lambda^T \bf(A, \hat{r}(X)); \hat{r}(X)\bigr)} - \EE{g\bigl(\lambda^T \bf(A, r(X)); r(X)\bigr)}}.
\end{align}

Under Assumption <ref>, $\mu(X) = \lambda^T \bf(X)$ is given by (<ref>) or (<ref>) and the notation $\bf(A, r(X))$ above makes clear that $\bf$ depends on $X$ through $A$ and $r(X)$ in this case. Recall also from Section <ref> that in $\hat{\bf}$, the probabilities that parametrize $\bf$ ($p_A$ for MSP, $p_{A,Y}$ for GEO) are replaced by their estimates. The first right-hand side quantity in (<ref>) is the difference between the empirical average and expectation of the same quantity. The second difference is due to having $\hat{\bf}$ instead of $\bf$, and the third is due to having $\hat{r}(X)$ instead of $r(X)$.

The following three lemmas, proven in Appendix <ref>, show that the three right-hand side terms in (<ref>) all converge uniformly to zero, thus completing the proof of the theorem.

Under Assumptions <ref>, <ref>, and <ref>.1,

\[
\lim_{n\to\infty} \sup_{\lambda\in\Lambda_0} \abs*{\frac{1}{n} \sum_{i=1}^n g\bigl(\lambda^T \hat{\bf}(a_i, \hat{r}(x_i)); \hat{r}(x_i)\bigr) - \EE{g\bigl(\lambda^T \hat{\bf}(A, \hat{r}(X)); \hat{r}(X)\bigr)}} \overset{p}{\to} 0.
\]

Under Assumptions <ref>, <ref>, and <ref>,

\[
\sup_{\lambda\in\Lambda_0} \abs*{\EE{g\bigl(\lambda^T \hat{\bf}(A, \hat{r}(X)); \hat{r}(X)\bigr)} - \EE{g\bigl(\lambda^T \bf(A, \hat{r}(X)); \hat{r}(X)\bigr)}} \overset{p}{\to} 0.
\]

Under Assumptions <ref>, <ref>, <ref>, and <ref>.1,

\[
\sup_{\lambda\in\Lambda_0} \abs*{\EE{g\bigl(\lambda^T \bf(A, \hat{r}(X)); \hat{r}(X)\bigr)} - \EE{g\bigl(\lambda^T \bf(A, r(X)); r(X)\bigr)}} \overset{p}{\to} 0.
\]

§ EMPIRICAL EVALUATION

This section discusses experimental evaluation of the proposed FST methods for MSP and GEO constraints and both pre- and post-processing.

Datasets Four datasets were used,
the first three of which are standard in the fairness literature: 1) adult income, 2) ProPublica's COMPAS recidivism, 3) German credit risk, 4) Medical Expenditure Panel Survey (MEPS).
Specifically, we used versions pre-processed by
an open-source library for algorithmic
fairness [49].
To facilitate comparison with other methods, we used binary-valued protected attributes and consider gender and race for both adult and COMPAS, age and gender for German, and race for MEPS.  Each dataset was randomly split $10$ times into training ($75\%$) and test ($25\%$) sets and all methods were subject to the same splits.

Methods compared
Since FST is intended for post- and pre-processing,
comparisons to other post- and pre-processing methods are most natural as they accommodate situations a)–c) in Section <ref>. For post-processing, we have chosen [4] (HPS) and the reject option method of [21], both as implemented in [49], as well as the Wass-1 Post-Process $\hat{p}_S$ method (WPP) of [35]. For pre-processing, the massaging and reweighing methods of [16] and the optimization method of [20] were chosen.  Among in-processing methods, meta-algorithms that work with essentially any base classifier can handle situation b). The reductions method [29] (`red') was selected from this class. We also compared to in-processing methods specific to certain types of classifiers, which do not allow for any of a)–c): fairness constraints (FC) [27], disparate mistreatment (DM) [5], and fair empirical risk minimization (FERM) [31]. Last but not least, availability of code was an important criterion.

Capabilities of methods in comparison. $\star$ refers to an extension implemented in [49].

method     pre     in     post     SP     EO     no $A$ at test time     scores     approx fairness     any classifier

massage, reweigh [16]

OPP [20]

HPS [4]

reject [21]                     $\star$

WPP [35]

FC [27]

DM [5]

FERM [31]

reductions [29]

proposed FST

The methods in the previous paragraph have various limitations, summarized by Table <ref>, that affect the design of the experiments.
First, the post-processing methods [4, 21, 35] (specifically the WPP variant for [35]) require knowledge of the protected attribute $A$ at test time. Accordingly, the experiments presented in this section include $A$ in the features $X$ to make it available to all methods; experiments without $A$
at test time (excluding [4, 21, 35]) are presented in Appendix <ref>.  We also encountered computational problems with
[20, 5] and thus perform separate comparisons
with FST on reduced feature sets, also reported in Appendix <ref>.

Three versions of FST were evaluated: post-processing (FSTpost), pre-processing (FSTpre), and a second post-processing version (FSTbatch) that assumes that test instances can be processed in a batch rather than one by one.  In this case, the fitting of the dual variables that parametrize FST (Section <ref>)
can actually be done on test data since it does not depend on labels $y_i$ (and uses only predicted probabilities for $A$ if $A$ is unavailable at test time).

Base classifiers We used $\ell_1$-regularized logistic regression (LR) and gradient boosted classification trees (GBM) from scikit-learn [50] as base classifiers. These are used in different ways depending on the method:
Post-processing methods operate on the scores produced by the base classifier, pre-processing methods train the base classifier after modifying the training data, and the reductions method repeatedly calls the base classification algorithm with different instance-specific costs.
In Appendix <ref>, we used linear
SVMs (with Platt scaling [51] to output probabilities) to compare with FERM [31]. We found it impractical to train nonlinear SVMs on the larger datasets for reductions and FERM since reductions needs to do so repeatedly and FERM uses a slower specialized algorithm.  For a similar reason, $5$-fold cross-validation to select parameters for LR (regularization parameter $C$ from $[10^{-4}, 10^4]$) and GBM (minimum number of samples per leaf from
$\{5,10,15,20,30\}$) was done only once per training set. All other parameters were set to the scikit-learn
defaults. The base classifier was then instantiated with the best parameter value for use by all methods.

Metrics Classification performance and fairness were evaluated using both score-based metrics (Brier score, AUC, differences in mean scores (MSP) and GEO) and binary label-based metrics (accuracy, differences in mean binary predictions (SP) and non-generalized EO).
We account for the fact that the reductions method [29] returns a randomized classifier, i.e. a probability distribution over a set of classifiers. For the binary label-based metrics, we used the methods provided with the code[<https://github.com/microsoft/fairlearn>] for reductions. The score-based metrics were computed by evaluating the metric for each classifier in the distribution and then averaging weighted by their probabilities.

Trade-offs between fairness and classification performance for a subset of dataset-protected attribute combinations and the case in which the features include the protected attributes (see Appendix <ref> for full set). Pareto efficient points are connected by line segments. Horizontal and vertical bars represent standard errors in the means.

Results Figure <ref> shows the trade-offs between classification performance and fairness obtained in a subset of the experiments we conducted.
The full set with other dataset-protected attribute combinations, protected attributes excluded from features, and reduced feature sets is in Appendix <ref>. Each dataset occupies two rows with the first showing score-based measures (Brier score vs. MSP or GEO differences, AUC is in Appendix <ref>) and the second showing binary label-based measures (accuracy vs. SP or EO differences). The columns correspond to combinations of base classifier (LR, GBM) and fairness measure targeted (SP, EO). Markers indicate mean values over the $10$ splits, error bars indicate standard errors in the means, and Pareto-optimal points have been connected with line segments to ease visualization.

Considering first the score-based plots (odd rows), FSTpost and FSTbatch achieve trade-offs that are at least as good as all other methods, with the slight exception of the GBM case on MEPS. In all cases, the advantage of FST lies in extending the Pareto frontiers farther to the left, attaining smaller MSP or GEO differences; this is especially apparent for GEO. FSTpre sometimes performs less well, e.g. with GBM on adult and MEPS, likely due to the loss incurred in approximating the transformed score $r'(x)$ with the output of a classifier fit to the pre-processed data.

Turning to the binary label-based plots (even rows), the trade-offs for FSTpost and FSTbatch generally coincide with or are close to the trade-offs of the best method, and are even sometimes the best, despite not optimizing for binary metrics beyond tuning the binarization threshold for accuracy. Again FSTpre with GBM is worse on adult, but FSTpre with LR is the best performer on COMPAS. The main disadvantage of FST is that its trade-off curves may not extend as far to the left as other methods, in particular on adult. This is the converse of its advantage for score-based metrics.

Among the existing methods, reductions is the strongest and also the most versatile, handling all cases that FST does. However, it is an in-processing method and far more computationally expensive, requiring an average of nearly $30$ calls to the base classification algorithm compared to one for FSTpost, FSTbatch and two for FSTpre. Reductions also returns a randomized classifier, which may not be desirable in some applications. The post-processing methods [21, 4] are not designed to output scores and hence are omitted from the score-based plots. Reject option [21] performs close to the best except on MEPS and at small unfairness values. HPS is limited to EO, does not have a parameter to vary the trade-off, and is less competitive.
WPP and the pre-processing methods of [16],
massaging and reweighing, likewise do not have a trade-off parameter and are limited to SP. As also observed in [29], massaging is often dominated by other methods while reweighing lies on the Pareto frontier but with substantial disparity. WPP results in low disparity but its Brier score or accuracy is sometimes less competitive. FC applies only to the LR-SP column and could not substantially reduce unfairness, possibly due to the larger feature dimension.

Limitations We caution the reader against some of the limitations of FST. First, the method inherently depends on well-calibrated classifiers  that approximate $p_{Y\given X}$ and, if necessary, $p_{A|X}$ or $p_{A|X,Y}$.
A poorly calibrated model
(e.g., due lack of samples) may lead to transformed scores that do not achieve the target fairness criteria. Second, thresholding the transformed scores may have an adverse impact on fairness guarantees, as seen throughout Figures <ref>–<ref>.
Finally, like most pre- and post-processing methods, the score transformation found by the  FST is vulnerable to distribution shifts between training and deployment.

§ CONCLUSION

We introduced a fairness-ensuring score transformation method called $\mathsf{FairScoreTransformer}$. This method can handle fairness criteria given in terms of linear inequalities involving conditional means of scores (cf. (<ref>)) and minimizes the cross-entropy between the original and the transformed scores (<ref>). FST is based on the  optimization problem (<ref>) which, in turn, has a dual form  that can be efficiently solved using the ADMM algorithm (Section <ref>). We provided explicit instantiations of the method for MSP and GEO fairness constraints in Prop. <ref>. Moreover, via a comprehensive set of experiments (cf. Section <ref>), we numerically demonstrated that FST is either as competitive or outperforms several existing fairness intervention mechanisms over a range of constraints and datasets. Future directions include characterizing convergence rates for the ADMM iterations and adapting to fairness criteria that are not based on conditional expectation of scores (e.g., calibration across groups [23]) as well as to non-binary outcomes $Y$.

[1]
D. Pedreschi, S. Ruggieri, and F. Turini, “A study of top-k measures for
discrimination discovery,” in Proceedings of the 27th Annual ACM
Symposium on Applied Computing.1em plus 0.5em minus 0.4emACM,
2012, pp. 126–131.

[2]
C. Dwork, M. Hardt, T. Pitassi, O. Reingold, and R. Zemel, “Fairness through
awareness,” in Proceedings of the 3rd Innovations in Theoretical
Computer Science Conference.1em plus 0.5em minus 0.4emACM,
2012, pp. 214–226.

[3]
F. Kamiran, I. Žliobaitė, and T. Calders, “Quantifying explainable
discrimination and removing illegal discrimination in automated decision
making,” Knowledge and Information Systems, vol. 35, no. 3, pp.
613–644, Jun. 2013. [Online]. Available:

[4]
M. Hardt, E. Price, and N. Srebro, “Equality of opportunity in supervised
learning,” in Advances in Neural Information Processing Systems
(NeurIPS), Dec. 2016, pp. 3323–3331. [Online]. Available:

[5]
M. B. Zafar, I. Valera, M. Gomez Rodriguez, and K. P. Gummadi, “Fairness
beyond disparate treatment & disparate impact: Learning classification
without disparate mistreatment,” in Proceedings of the 26th
International Conference on World Wide Web (WWW).1em plus 0.5em minus
0.4emInternational World Wide Web Conferences Steering Committee,
2017, pp. 1171–1180.

[6]
A. Chouldechova, “Fair prediction with disparate impact: A study of bias in
recidivism prediction instruments,” Big Data, vol. 5, no. 2, pp.
153–163, 2017.

[7]
J. Kleinberg, S. Mullainathan, and M. Raghavan, “Inherent trade-offs in the
fair determination of risk scores,” in Proceedings of the 8th
Innovations in Theoretical Computer Science Conference (ITCS), 2017, pp.

[8]
N. Kilbertus, M. Rojas-Carulla, G. Parascandolo, M. Hardt, D. Janzing, and
B. Schölkopf, “Avoiding discrimination through causal reasoning,” in
Advances in Neural Information Processing Systems (NeurIPS), Dec.
2017, pp. 656–666. [Online]. Available:

[9]
M. J. Kusner, J. Loftus, C. Russell, and R. Silva, “Counterfactual fairness,”
in Advances in Neural Information Processing Systems (NeurIPS), Dec.
2017, pp. 4066–4076. [Online]. Available:

[10]
M. B. Zafar, I. Valera, M. Rodriguez, K. Gummadi, and A. Weller, “From parity
to preference-based notions of fairness in classification,” in
Advances in Neural Information Processing Systems (NeurIPS), Dec.
2017, pp. 229–239. [Online]. Available:

[11]
R. Nabi and I. Shpitser, “Fair inference on outcomes,” in AAAI
Conference on Artificial Intelligence (AAAI), Feb. 2018, pp. 1931–1940.
[Online]. Available:

[12]
M. Kearns, S. Neel, A. Roth, and Z. S. Wu, “Preventing fairness
gerrymandering: Auditing and learning for subgroup fairness,” in
International Conference on Machine Learning, 2018, pp. 2569–2577.

[13]
H. Heidari, C. Ferrari, K. Gummadi, and A. Krause, “Fairness behind a veil of
ignorance: A welfare analysis for automated decision making,” in
Advances in Neural Information Processing Systems (NeurIPS), Dec.
2018, pp. 1265–1276. [Online]. Available:

[14]
S. Chiappa, “Path-specific counterfactual fairness,” in AAAI Conference
on Artificial Intelligence (AAAI), Jan. 2019.

[15]
S. A. Friedler, C. Scheidegger, S. Venkatasubramanian, S. Choudhary, E. P.
Hamilton, and D. Roth, “A comparative study of fairness-enhancing
interventions in machine learning,” in Proceedings of the Conference
on Fairness, Accountability, and Transparency.1em plus 0.5em minus
0.4emACM, 2019, pp. 329–338.

[16]
F. Kamiran and T. Calders, “Data preprocessing techniques for classification
without discrimination,” Knowledge and Information Systems, vol. 33,
no. 1, pp. 1–33, Oct. 2012. [Online]. Available:

[17]
S. Hajian and J. Domingo-Ferrer, “A methodology for direct and indirect
discrimination prevention in data mining,” IEEE Transactions on
Knowledge and Data Engineering, vol. 25, no. 7, pp. 1445–1459, July 2013.

[18]
R. Zemel, Y. Wu, K. Swersky, T. Pitassi, and C. Dwork, “Learning fair
representations,” in Proceedings of the 30th International Conference
on Machine Learning (ICML), Atlanta, Georgia, USA, Jun. 2013, pp. 325–333.
[Online]. Available: <http://proceedings.mlr.press/v28/zemel13.html>

[19]
M. Feldman, S. A. Friedler, J. Moeller, C. Scheidegger, and
S. Venkatasubramanian, “Certifying and removing disparate impact,” in
Proceedings of the 21st ACM SIGKDD International Conference on
Knowledge Discovery and Data Mining (KDD), Aug. 2015, pp. 259–268.
[Online]. Available: <http://doi.acm.org/10.1145/2783258.2783311>

[20]
F. Calmon, D. Wei, B. Vinzamuri, K. Natesan Ramamurthy, and K. R. Varshney,
“Optimized pre-processing for discrimination prevention,” in Advances
in Neural Information Processing Systems (NeurIPS), Dec. 2017, pp.
3992–4001. [Online]. Available:

[21]
F. Kamiran, A. Karim, and X. Zhang, “Decision theory for
discrimination-aware classification,” in IEEE 12th International
Conference on Data Mining (ICDM), Dec 2012, pp. 924–929.

[22]
B. Fish, J. Kun, and Ádám D. Lelkes, “A confidence-based approach for
balancing fairness and accuracy,” in Proceedings of the 2016 SIAM
International Conference on Data Mining (SDM), May 2016, pp. 144–152.

[23]
G. Pleiss, M. Raghavan, F. Wu, J. Kleinberg, and K. Q. Weinberger, “On
fairness and calibration,” in Advances in Neural Information
Processing Systems (NeurIPS), Dec. 2017, pp. 5680–5689. [Online].

[24]
B. Woodworth, S. Gunasekar, M. I. Ohannessian, and N. Srebro, “Learning
non-discriminatory predictors,” in Proceedings of the 2017 Conference
on Learning Theory (COLT), Jul. 2017, pp. 1920–1953. [Online]. Available:

[25]
T. Calders and S. Verwer, “Three naive Bayes approaches for
discrimination-free classification,” Data Mining and Knowledge
Discovery, vol. 21, no. 2, pp. 277–292, Sep. 2010. [Online]. Available:

[26]
T. Kamishima, S. Akaho, H. Asoh, and J. Sakuma, “Fairness-aware classifier
with prejudice remover regularizer,” in Joint European Conference on
Machine Learning and Knowledge Discovery in Databases (ECML PKDD), Sep.
2012, pp. 35–50.

[27]
M. B. Zafar, I. Valera, M. G. Rogriguez, and K. P. Gummadi, “Fairness
Constraints: Mechanisms for Fair Classification,” in Proceedings of
the 20th International Conference on Artificial Intelligence and Statistics
(AISTATS), Apr. 2017, pp. 962–970. [Online]. Available:

[28]
C. Dwork, N. Immorlica, A. T. Kalai, and M. Leiserson, “Decoupled classifiers
for group-fair and efficient machine learning,” in Proceedings of the
1st Conference on Fairness, Accountability and Transparency (FAT*), Feb.
2018, pp. 119–133. [Online]. Available:

[29]
A. Agarwal, A. Beygelzimer, M. Dudík, J. Langford, and H. Wallach, “A
reductions approach to fair classification,” in Proceedings of the
35th International Conference on Machine Learning (ICML), Jul. 2018, pp.
60–69. [Online]. Available:

[30]
E. Krasanakis, E. Spyromitros-Xioufis, S. Papadopoulos, and Y. Kompatsiaris,
“Adaptive sensitive reweighting to mitigate bias in fairness-aware
classification,” in Proceedings of the 2018 World Wide Web Conference
(WWW), 2018, pp. 853–862.

[31]
M. Donini, L. Oneto, S. Ben-David, J. S. Shawe-Taylor, and M. Pontil,
“Empirical risk minimization under fairness constraints,” in Advances
in Neural Information Processing Systems (NeurIPS), Dec. 2018, pp.
2791–2801. [Online]. Available:

[32]
L. E. Celis, L. Huang, V. Keswani, and N. K. Vishnoi, “Classification with
fairness constraints: A meta-algorithm with provable guarantees,” in
Proceedings of the Conference on Fairness, Accountability, and
Transparency (FAT*), Jan. 2019, pp. 319–328. [Online]. Available:

[33]
D. Madras, E. Creager, T. Pitassi, and R. Zemel, “Learning adversarially fair
and transferable representations,” in Proceedings of the 35th
International Conference on Machine Learning (ICML), Jul. 2018, pp.
3384–3393. [Online]. Available:

[34]
B. Salimi, L. Rodriguez, B. Howe, and D. Suciu, “Capuchin: Causal database
repair for algorithmic fairness,” Feb. 2019, arXiv e-print

[35]
R. Jiang, , A. Pacchiano, T. Stepleton, H. Jiang, and S. Chiappa, “Wasserstein
fair classification,” in Conference on Uncertainty in Artificial
Intelligence (UAI), Jul. 2019, pp. 1–10.

[36]
C. Louizos, K. Swersky, Y. Li, M. Welling, and R. Zemel, “The variational fair
encoder,” in Proceedings of the International Conference on Learning
Representations (ICLR), May 2016, pp. 1–11.

[37]
H. Edwards and A. Storkey, “Censoring representations with an adversary,” in
Proceedings of the International Conference on Learning Representations
(ICLR), May 2016, pp. 1–14.

[38]
Q. Xie, Z. Dai, Y. Du, E. Hovy, and G. Neubig, “Controllable invariance
through adversarial feature learning,” in Advances in Neural
Information Processing Systems (NeurIPS), Dec. 2017, pp. 585–596. [Online].

[39]
D. Xu, S. Yuan, L. Zhang, and X. Wu, “FairGAN: Fairness-aware
generative adversarial networks,” in IEEE International Conference on
Big Data (Big Data), Dec. 2018, pp. 570–575.

[40]
A. Beutel, J. Chen, Z. Zhao, and E. H. Chi, “Data decisions and theoretical
implications when adversarially learning fair representations,” in
Workshop on Fairness, Accountability, and Transparency in Machine
Learning (FATML), Aug. 2017, pp. 1–5.

[41]
B. H. Zhang, B. Lemoine, and M. Mitchell, “Mitigating unwanted biases with
adversarial learning,” in Proceedings of the AAAI/ACM Conference on
Artifical Intelligence, Ethics, and Society (AIES), Feb. 2018, pp. 335–340.
[Online]. Available: <http://doi.acm.org/10.1145/3278721.3278779>

[42]
A. K. Menon and R. C. Williamson, “The cost of fairness in binary
classification,” in Proceedings of the 1st Conference on Fairness,
Accountability and Transparency (FAT*), Feb. 2018, pp. 107–118. [Online].
Available: <http://proceedings.mlr.press/v81/menon18a.html>

[43]
S. Corbett-Davies, E. Pierson, A. Feller, S. Goel, and A. Huq, “Algorithmic
decision making and the cost of fairness,” in Proceedings of the 23rd
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining
(KDD), Aug. 2017, pp. 797–806. [Online]. Available:

[44]
H. Jiang and O. Nachum, “Identifying and correcting label bias in machine
learning,” Jan. 2019, arXiv e-print <https://arxiv.org/abs/1901.04966>.

[45]
A. Coston, K. N. Ramamurthy, D. Wei, K. R. Varshney, S. Speakman, Z. Mustahsan,
and S. Chakraborty, “Fair transfer learning with missing protected
attributes,” in Proceedings of the AAAI/ACM Conference on Artifical
Intelligence, Ethics, and Society (AIES), Jan. 2019, pp. 1–8.

[46]
S. Boyd and L. Vandenberghe, Convex Optimization.1em plus 0.5em
minus 0.4emCambridge, UK: Cambridge University Press, 2004.

[47]
J. Hernández-Orallo, P. Flach, and C. Ferri, “A unified view of
performance metrics: Translating threshold choice into expected
classification loss,” J. Mach. Learn. Res., vol. 13, pp. 2813–2869,
Oct. 2012.

[48]
S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein et al.,
“Distributed optimization and statistical learning via the alternating
direction method of multipliers,” Foundations and
Trends® in Machine learning, vol. 3, no. 1, pp. 1–122,

[49]
R. K. E. Bellamy, K. Dey, M. Hind, S. C. Hoffman, S. Houde, K. Kannan,
P. Lohia, J. Martino, S. Mehta, A. Mojsilovic, S. Nagar, K. N. Ramamurthy,
J. Richards, D. Saha, P. Sattigeri, M. Singh, K. R. Varshney, and Y. Zhang,
“AI Fairness 360: An extensible toolkit for detecting, understanding, and
mitigating unwanted algorithmic bias,” Oct. 2018,
<https://arxiv.org/abs/1810.01943>. [Online]. Available:

[50]
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel,
M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos,
D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn:
Machine learning in Python,” Journal of Machine Learning Research,
vol. 12, pp. 2825–2830, 2011.

[51]
J. C. Platt, “Probabilistic outputs for support vector machines and
comparisons to regularized likelihood methods,” in Advances in Large
Margin Classifiers, 1999, pp. 61–74.

[52]
D. Bertsimas and J. N. Tsitsiklis, Introduction to Linear Optimization,
2nd ed.1em plus 0.5em minus 0.4emBelmont, MA, USA: Athena
Scientific, 1997.

[53]
R. I. Jennrich, “Asymptotic properties of non-linear least squares
estimators,” Annals of Mathematical Statistics, vol. 40, no. 2, pp.
633–643, 1969.

[54]
D. P. Bertsekas, Nonlinear Programming, 2nd ed.1em plus 0.5em
minus 0.4emBelmont, MA, USA: Athena Scientific, 1999.


* Ali-Haïmoud (2018) Ali-Haïmoud, Y.  Correlation Function of High-Threshold Regions and Application to the Initial Small-Scale Clustering of Primordial Black Holes.  Phys. Rev. Lett. 2018, 121, 081304, [arXiv:astro-ph.CO/1805.05912].  doi:black10.1103/PhysRevLett.121.081304.
  * Petač et al. (2022) Petač, M.; Lavalle, J.; Jedamzik, K.  Microlensing constraints on clustered primordial black holes.  arXiv e-prints 2022, p. arXiv:2201.02521, [arXiv:astro-ph.CO/2201.02521].
  * Suyama and Yokoyama (2019) Suyama, T.; Yokoyama, S.  Clustering of primordial black holes with non-Gaussian initial fluctuations.  Progress of Theoretical and Experimental Physics 2019, 2019, 103E02, [arXiv:astro-ph.CO/1906.04958].  doi:black10.1093/ptep/ptz105.
  * Young and Byrnes (2020) Young, S.; Byrnes, C.T.  Initial clustering and the primordial black hole merger rate.  Journal of Cosmology and Astroparticle Physics 2020, 2020, 004, [arXiv:astro-ph.CO/1910.06077].  doi:black10.1088/1475-7516/2020/03/004.
  * Atal et al. (2020) Atal, V.; Sanglas, A.; Triantafyllou, N.  LIGO/Virgo black holes and dark matter: the effect of spatial clustering.  Journal of Cosmology and Astroparticle Physics 2020, 2020, 036, [arXiv:astro-ph.CO/2007.07212].  doi:black10.1088/1475-7516/2020/11/036.
  * Abbott and et al. (2021) Abbott, R.; et al..  Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run.  Phys. Rev. D 2021, 104, 022004.  doi:black10.1103/PhysRevD.104.022004.
  * Gorton and Green (2022) Gorton, M.; Green, A.M.  Effect of clustering on primordial black hole microlensing constraints.  arXiv e-prints 2022, p. arXiv:2203.04209, [arXiv:astro-ph.CO/2203.04209].
  * Roncadelli et al. (2009) Roncadelli, M.; Treves, A.; Turolla, R.  Primordial black holes are again on the limelight.  arXiv e-prints 2009, p. arXiv:0901.1093, [arXiv:astro-ph.CO/0901.1093].
  * Capela et al. (2013) Capela, F.; Pshirkov, M.; Tinyakov, P.  Constraints on primordial black holes as dark matter candidates from star formation.  Phys. Rev. D 2013, 87, 023507, [arXiv:astro-ph.CO/1209.6021].  doi:black10.1103/PhysRevD.87.023507.
  * Capela et al. (2014) Capela, F.; Pshirkov, M.; Tinyakov, P.  Adiabatic contraction revisited: Implications for primordial black holes.  Phys. Rev. D 2014, 90, 083507, [arXiv:astro-ph.CO/1403.7098].  doi:black10.1103/PhysRevD.90.083507.
  * Montero-Camacho et al. (2019) Montero-Camacho, P.; Fang, X.; Vasquez, G.; Silva, M.; Hirata, C.M.  Revisiting constraints on asteroid-mass primordial black holes as dark matter candidates.  Journal of Cosmology and Astro-Particle Physics 2019, 2019, 031, [arXiv:astro-ph.CO/1906.05950].  doi:black10.1088/1475-7516/2019/08/031.
  * Oncins et al. (2022) Oncins, M.; Miralda-Escudé, J.; Gutiérrez, J.L.; Gil-Pons, P.  Primordial black holes capture by stars and induced collapse to low-mass stellar black holes.  arXiv e-prints 2022, p. arXiv:2205.13003, [arXiv:astro-ph.GA/2205.13003].
  * Pani and Loeb (2014) Pani, P.; Loeb, A.  Tidal capture of a primordial black hole by a neutron star: implications for constraints on dark matter.  Journal of Cosmology and Astroparticle Physics 2014, 2014, 026, [arXiv:astro-ph.CO/1401.3025].  doi:black10.1088/1475-7516/2014/06/026.
  * Capela et al. (2014) Capela, F.; Pshirkov, M.; Tinyakov, P.  A comment on “Exclusion of the remaining mass window for primordial black holes …”, arXiv:1401.3025.  arXiv e-prints 2014, p. arXiv:1402.4671, [arXiv:astro-ph.CO/1402.4671].
  * Defillon et al. (2014) Defillon, G.; Granet, E.; Tinyakov, P.; Tytgat, M.H.G.  Tidal capture of primordial black holes by neutron stars.  Phys. Rev. D 2014, 90, 103522, [arXiv:gr-qc/1409.0469].  doi:black10.1103/PhysRevD.90.103522.
  * Hütsi et al. (2021) Hütsi, G.; Raidal, M.; Vaskonen, V.; Veermäe, H.  Two populations of LIGO-Virgo black holes.  Journal of Cosmology and Astroparticle Physics 2021, 2021, 068, [arXiv:astro-ph.CO/2012.02786].  doi:black10.1088/1475-7516/2021/03/068.
  * De Luca et al. (2021) De Luca, V.; Franciolini, G.; Pani, P.; Riotto, A.  Bayesian evidence for both astrophysical and primordial black holes: mapping the GWTC-2 catalog to third-generation detectors.  Journal of Cosmology and Astro-Particle Physics 2021, 2021, 003, [arXiv:astro-ph.CO/2102.03809].  doi:black10.1088/1475-7516/2021/05/003.
  * Franciolini et al. (2022) Franciolini, G.; Baibhav, V.; De Luca, V.; Ng, K.K.Y.; Wong, K.W.K.; Berti, E.; Pani, P.; Riotto, A.; Vitale, S.  Searching for a subpopulation of primordial black holes in LIGO-Virgo gravitational-wave data.  Phys. Rev. D 2022, 105, 083526, [arXiv:gr-qc/2105.03349].  doi:black10.1103/PhysRevD.105.083526.
  * Eroshenko (2021) Eroshenko, Y.  Mergers of primordial black holes in extreme clusters and the H0 tension.  Physics of the Dark Universe 2021, 32, 100833, [arXiv:astro-ph.CO/2105.03704].  doi:black10.1016/j.dark.2021.100833.
  * Kusenko et al. (2020) Kusenko, A.; Sasaki, M.; Sugiyama, S.; Takada, M.; Takhistov, V.; Vitagliano, E.  Exploring Primordial Black Holes from the Multiverse with Optical Telescopes.  Phys. Rev. Lett. 2020, 125, 181304, [arXiv:astro-ph.CO/2001.09160].  doi:black10.1103/PhysRevLett.125.181304.
  * Hawkins (2022) Hawkins, M.R.S.  New evidence for a cosmological distribution of stellar mass primordial black holes.  MNRAS 2022, 512, 5706–5714, [arXiv:astro-ph.CO/2204.09143].  doi:black10.1093/mnras/stac863.
  * Churazov et al. (2005) Churazov, E.; Sunyaev, R.; Sazonov, S.; Revnivtsev, M.; Varshalovich, D.  Positron annihilation spectrum from the Galactic Centre region observed by SPI/INTEGRAL.  MNRAS 2005, 357, 1377–1386, [arXiv:astro-ph/astro-ph/0411351].  doi:black10.1111/j.1365-2966.2005.08757.x.
  * Kierans et al. (2020) Kierans, C.A.; Boggs, S.E.; Zoglauer, A.; Lowell, A.W.; Sleator, C.; Beechert, J.; Brandt, T.J.; Jean, P.; Lazar, H.; Roberts, J.; Siegert, T.; Tomsick, J.A.; Ballmoos, P.v.  Detection of the 511 keV Galactic Positron Annihilation Line with COSI.  ApJ 2020, 895, 44, [arXiv:astro-ph.HE/1912.00110].  doi:black10.3847/1538-4357/ab89a9.
  * Keith and Hooper (2021) Keith, C.; Hooper, D.  511 keV excess and primordial black holes.  Phys. Rev. D 2021, 104, 063033, [arXiv:astro-ph.CO/2103.08611].  doi:black10.1103/PhysRevD.104.063033.
  * Mukhopadhyay et al. (2021) Mukhopadhyay, U.; Majumdar, D.; Paul, A.  Discriminating and Constraining the Synchrotron and Inverse Compton Radiations from Primordial Black Hole and Dark Matter at the Galactic Centre Region.  arXiv e-prints 2021, p. arXiv:2109.14955, [arXiv:astro-ph.CO/2109.14955].
  * Arzoumanian and et al. (2020) Arzoumanian, Z.; et al..  The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background.  ApJ 2020, 905, L34, [arXiv:astro-ph.HE/2009.04496].  doi:black10.3847/2041-8213/abd401.
  * De Luca et al. (2021) De Luca, V.; Franciolini, G.; Riotto, A.  NANOGrav Data Hints at Primordial Black Holes as Dark Matter.  Phys. Rev. Lett. 2021, 126, 041303, [arXiv:astro-ph.CO/2009.08268].  doi:black10.1103/PhysRevLett.126.041303.
  * Domènech and Pi (2020) Domènech, G.; Pi, S.  NANOGrav Hints on Planet-Mass Primordial Black Holes.  arXiv e-prints 2020, p. arXiv:2010.03976, [arXiv:astro-ph.CO/2010.03976].
  * Kohri and Terada (2021) Kohri, K.; Terada, T.  Solar-mass primordial black holes explain NANOGrav hint of gravitational waves.  Physics Letters B 2021, 813, 136040, [arXiv:astro-ph.CO/2009.11853].  doi:black10.1016/j.physletb.2020.136040.
  * Sugiyama et al. (2021) Sugiyama, S.; Takhistov, V.; Vitagliano, E.; Kusenko, A.; Sasaki, M.; Takada, M.  Testing stochastic gravitational wave signals from primordial black holes with optical telescopes.  Physics Letters B 2021, 814, 136097, [arXiv:astro-ph.CO/2010.02189].  doi:black10.1016/j.physletb.2021.136097.
  * Green and Kavanagh (2021) Green, A.M.; Kavanagh, B.J.  Primordial black holes as a dark matter candidate.  Journal of Physics G Nuclear Physics 2021, 48, 043001, [arXiv:astro-ph.CO/2007.10722].  doi:black10.1088/1361-6471/abc534.
  * Bird et al. (2022) Bird, S.; Albert, A.; Dawson, W.; Ali-Haimoud, Y.; Coogan, A.; Drlica-Wagner, A.; Feng, Q.; Inman, D.; Inomata, K.; Kovetz, E.; Kusenko, A.; Lehmann, B.V.; Munoz, J.B.; Singh, R.; Takhistov, V.; Tsai, Y.D.  Snowmass2021 Cosmic Frontier White Paper:Primordial Black Hole Dark Matter.  arXiv e-prints 2022, p. arXiv:2203.08967, [arXiv:hep-ph/2203.08967].
  * LIGO Scientific Collaboration and Harms (2017) LIGO Scientific Collaboration.; Harms, J.  Exploring the sensitivity of next generation gravitational wave detectors.  Classical and Quantum Gravity 2017, 34, 044001, [arXiv:astro-ph.IM/1607.08697].  doi:black10.1088/1361-6382/aa51f4.
  * Reitze and et al. (2019) Reitze, D.; et al..  Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO.  Bulletin of the American Astronomical Society, 2019, Vol. 51, p. 35, [arXiv:astro-ph.IM/1907.04833].
  * Maggiore and et al. (2020) Maggiore, M.; et al..  Science case for the Einstein telescope.  Journal of Cosmology and Astro-Particle Physics 2020, 2020, 050, [arXiv:astro-ph.CO/1912.02622].  doi:black10.1088/1475-7516/2020/03/050.
  * Ng et al. (2022) Ng, K.K.Y.; Franciolini, G.; Berti, E.; Pani, P.; Riotto, A.; Vitale, S.  Constraining high-redshift stellar-mass primordial black holes with next-generation ground-based gravitational-wave detectors.  arXiv e-prints 2022, p. arXiv:2204.11864, [arXiv:astro-ph.CO/2204.11864].
  * De Luca et al. (2021) De Luca, V.; Franciolini, G.; Pani, P.; Riotto, A.  The minimum testable abundance of primordial black holes at future gravitational-wave detectors.  Journal of Cosmology and Astro-Particle Physics 2021, 2021, 039, [arXiv:astro-ph.CO/2106.13769].  doi:black10.1088/1475-7516/2021/11/039.
  * Amaro-Seoane and et al. (2017) Amaro-Seoane, P.; et al..  Laser Interferometer Space Antenna.  arXiv e-prints 2017, p. arXiv:1702.00786, [arXiv:astro-ph.IM/1702.00786].
  * Vaskonen and Veermäe (2021) Vaskonen, V.; Veermäe, H.  Did NANOGrav See a Signal from Primordial Black Hole Formation?  Phys. Rev. Lett. 2021, 126, 051303, [arXiv:astro-ph.CO/2009.07832].  doi:black10.1103/PhysRevLett.126.051303.
  * NANOGrav Collaboration (2018) NANOGrav Collaboration.  Science with the Next-Generation VLA and Pulsar Timing Arrays.  arXiv e-prints 2018, p. arXiv:1810.06594, [arXiv:astro-ph.IM/1810.06594].
  * Singh and et al (2014) Singh, K.P.; et al.  ASTROSAT mission.  Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray; Takahashi, T.; den Herder, J.W.A.; Bautz, M., Eds., 2014, Vol. 9144, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 91441S.  doi:black10.1117/12.2062667.
  * Barcons et al. (2015) Barcons, X.; Nandra, K.; Barret, D.; den Herder, J.W.; Fabian, A.C.; Piro, L.; Watson, M.G.; the Athena Team.  Athena: the X-ray observatory to study the hot and energetic Universe.  Journal of Physics Conference Series, 2015, Vol. 610, Journal of Physics Conference Series, p. 012008.  doi:black10.1088/1742-6596/610/1/012008.
  * The Lynx Team (2018) The Lynx Team.  The Lynx Mission Concept Study Interim Report.  arXiv e-prints 2018, p. arXiv:1809.09642, [arXiv:astro-ph.IM/1809.09642].
  * Zhang and et al (2016) Zhang, S.N.; et al.  eXTP: Enhanced X-ray Timing and Polarization mission.  Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray; den Herder, J.W.A.; Takahashi, T.; Bautz, M., Eds., 2016, Vol. 9905, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 99051Q, [arXiv:astro-ph.IM/1607.08823].  doi:black10.1117/12.2232034.