Individual Fairness Guarantees for Neural Networks

1 Introduction

Reservations have been raised about the application of neural networks (NN) in contexts where fairness is of concern [Barocas and Selbst, 2016]. Because of inherent biases present in real-world data, if unchecked, these models have been found to discriminate against individuals on the basis of sensitive features, such as race or sex [Bolukbasi et al., 2016; Angwin et al., 2016]. Recently, the topic has come under the spotlight, with technologies being increasingly challenged for bias [Kirk et al., 2021], leading to the introduction of a range of definitions and techniques for capturing the multifaceted properties of fairness.

Fairness approaches are broadly categorised into: group fairness [Hardt et al., 2016], which inspects the model over data demographics; and individual fairness (IF) [Dwork et al., 2012], which considers the behaviour over each individual. Despite its wider adoption, group fairness is only concerned with statistical properties of the model so that a situation may arise where predictions of a group-fair model can be perceived as unfair by a particular individual. In contrast, IF is a worst-case measure with guarantees over every possible individual in the input space. However, while techniques exist for group fairness of NNs [Albarghouthi et al., 2017; Bastani et al., 2019], research on IF has thus far been limited to designing training procedures that favour fairness [Yurochkin et al., 2020; Yeom and Fredrikson, 2021; McNamara et al., 2017] and verification over specific individuals [Ruoss et al., 2020]. To the best of our knowledge, there is currently no work targeted at global certification of IF for NNs.

∗Corresponding Author

We develop an anytime algorithm with provable bounds for the certification of IF on NNs. We build on the ϵ-δ-IF formalisation employed by [John et al., 2020]. That is, given ϵ, δ ≥ 0 and a distance metric dfair that captures the similarity between individuals, we ask that, for every pair of points x′ and x′′ in the input space with dfair(x′ , x′′ ) ≤ ϵ, the NN's output does not differ by more than δ. Although related to it, IF certification on NNs poses a different problem than adversarial robustness [Tjeng et al., 2018], as both x′ and x′′ are here problem variables, spanning the whole space. Hence, local approximation techniques developed in the adversarial literature cannot be employed in the context of IF.

Nevertheless, we show how this global, non-linear requirement can be encoded in Mixed-Integer Linear Programming (MILP) form, by deriving a set of global upper and lower piecewise-linear (PWL) bounds over each activation function in the NN over the whole input space, and performing linear encoding of the (generally non-linear) similarity metric dfair(x′ , x′′ ). The formulation of our optimisation as a MILP allows us to compute an anytime, worst-case bound on IF, which can thus be computed using standard solvers from the global optimisation literature [Dantzig, 2016]. Furthermore, we demonstrate how our approach can be embedded into the NN training so as to optimise for individual fairness at training time. We do this by performing gradient descent on a weighted loss that also accounts for the maximum δ-variation in dfair-neighborhoods for each training point, similarly to what is done in adversarial learning [Goodfellow et al., 2015; Wicker et al., 2021].

We apply our method on four benchmarks widely employed in the fairness literature, namely, the Adult, German, Credit and Crime datasets 1 , and an array of similarity metrics learnt from data that include ℓ∞, Mahalanobis, and NN embeddings. We empirically demonstrate how our method is able to provide the first, non-trivial IF certificates for NNs commonly employed for tasks from the IF literature, and even larger NNs comprising up to thousands of neurons. Furthermore, we find that our MILP-based fair training approach consistently outperforms, in terms of IF guarantees, NNs trained with a competitive state-of-the-art technique by orders of magnitude, albeit at an increased computational cost.

The paper makes the following main contributions: 2

• We design a MILP-based, anytime verification approach for the certification of IF as a global property on NNs.
• We demonstrate how our technique can be used to modify the loss function of a NN to take into account certification of IF at training time.
• On four datasets, and an array of metrics, we show how our techniques obtain non-trivial IF certificates and train NNs that are significantly fairer than state-of-the-art.

Related Work. A number of works have considered IF by employing techniques from adversarial robustness. [Yeom and Fredrikson, 2021] rely on randomized smoothing to find the highest stable per-feature difference in a model. Their method, however, provides only (weak) guarantees on model statistics. [Yurochkin et al., 2020] present a method for IF training that builds on projected gradient descent and optimal transport. While the method is found to decrease model bias to state-of-the-art results, no formal guarantees are obtained. [Ruoss et al., 2020] adapted the MILP formulation for adversarial robustness to handle fair metric embeddings. However, rather than tackling the IF problem globally as introduced by [Dwork et al., 2012], the method only works iteratively on a finite set of data, hence leaving open the possibility of unfairness in the model. In contrast, the MILP encoding we obtain through PWL bounding of activations and similarity metrics allows us to provide guarantees over any possible pair of individuals. [Urban et al., 2020] employ static analysis to certify causal fairness. While this method yields global guarantees, it cannot be straightforwardly employed for IF, and it is not anytime, making exhaustive analysis impractical. [John et al., 2020] present a method for the computation of IF, though limited to linear and kernel models. MILP and linear relaxation have been employed to certify NNs in local adversarial settings [Ehlers, 2017; Tjeng et al., 2018; Wicker et al., 2020]. However, local approximations cannot be employed for the global IF problem. While [Katz et al., 2017; Leino et al., 2021] consider global robustness, their methods are restricted to ℓp metrics. Furthermore, they require the knowledge of a Lipschitz constant or are limited to ReLU.

2 Individual Fairness

We focus on regression and binary classification with NNs with real-valued inputs and one-hot encoded categorical features. 3 Such frameworks are often used in automated decision-making, e.g. for loan applications [Hardt et al., 2016]. Formally, given a compact input set X ⊆ Rn and an output set Y ⊆ R, we consider an L layer fully-connected NN f w : X → Y , parameterised by a vector of weights w ∈ Rnw trained on D = {(xi, yi), i ∈ {1, ..., nd}}. For an input x ∈ X, i = 1, . . . , L and j = 1, . . . , ni, the NN is defined as:

1 http://archive.ics.uci.edu/ml

2 Proofs and additional details can be found in Appendix of an extended version of the paper available at http://www.fun2model. org/bibitem.php?key=BPW+22.

where ζ (0) j = xj. Here, ni is the number of units in the ith layer, W (i) jk and b(i) j are its weights and biases, σ(i) is the activation function, ϕ(i) is the pre-activation and ζ (i) the activation. The NN output is the result of these computations, f w (x) := ζ (L) . In regression, f w (x) is the prediction, while for classification it represents the class probability. In this paper we focus on fully-connected NNs as widely employed in the IF literature[Yurochkin et al., 2020; Urban et al., 2020; Ruoss et al., 2020]. However, we should stress that our framework, being based on MILP, can be easily extended to convolutional, max-pool and batch-norm layers or res-nets by using embedding techniques from the adversarial robustness literature (see e.g. [Boopathy et al., 2019].

Individual Fairness. Given a NN f w , IF [Dwork et al., 2012] enforces the property that similar individuals are similarly treated. Similarity is defined according to a taskdependent pseudometric, dfair : X × X ↦→ R≥0, provided by a domain expert (e.g., a Mahalanobis distance correlating each feature to the sensitive one), whereas similarity of treatment is expressed via the absolute difference on the NN output f w (x). We adopt the ϵ-δ-IF formulation of [John et al., 2020] for the formalisation of input-output IF similarity.

Definition 1 (ϵ-δ-IF [John et al., 2020]). Consider ϵ ≥ 0 and δ ≥ 0. We say that f w is ϵ-δ-individually fair w.r.t. dfair iff

Here, ϵ measures similarity between individuals and δ is the difference in outcomes (class probability for classification). We emphasise that individual fairness is a global notion, as the condition in Definition 1 must hold for all pairs of points in X. We remark that the ϵ-δ-IF formulation of [John et al., 2020] (which is more general than IF formulation typically used in the literature [Yurochkin et al., 2020; Ruoss et al., 2020]) is a slight variation on the Lipschitz property introduced by [Dwork et al., 2012]. While introducing greater flexibility thanks to its parametric form, it makes an IF parametric analysis necessary at test time. In Section 4 we analyse how ϵ-δ-IF of NNs is affected by variations of ϵ and δ. A crucial component of IF is the similarity dfair. The intuition is that sensitive features, or their sensitive combination, should not influence the NN output. While a number of metrics has been discussed in the literature [Ilvento, 2020], we focus on the following representative set of metrics which can be automatically learnt from data [John et al., 2020;

3 Multi-class can be tackled with component-wise analyses.

Ruoss et al., 2020; Mukherjee et al., 2020; Yurochkin et al., 2020]. Details on metric learning is given in Appendix B.

Weighted ℓp. In this case dfair(x′ , x′′ ) is defined as a weighted version of an ℓp metric, i.e. dfair(x′ , x′′ ) = p√∑ni=1 θi|x′i -x′′ i |p . Intuitively, we set the weights θi related to sensitive features to zero, so that two individuals are considered similar if they only differ with respect to those. The weights θi for the remaining features can be tuned according to their degree of correlation to the sensitive features.

Mahalanobis. In this case we have dfair(x′ , x′′ ) = √(x′ -x′′ )T S(x′ -x′′ ), for a given positive semi-definite (SPD) matrix S. The Mahalanobis distance generalises the ℓ2 metric by taking into account the intra-correlation of features to capture latent dependencies w.r.t. the sensitive features.

Feature Embedding. The metric is computed on an embedding, so that dfair(x′ , x′′ ) = ˆ d(φ(x′ ), φ(x′′ )), where ˆ d is either the Mahalanobis or the weighted ℓp metric, and φ is a feature embedding map. These allow for greater modelling flexibility, at the cost of reduced interpretability.

2.1 ProblemFormulation

We aim at certifying ϵ-δ-IF for NNs. To this end we formalise two problems: computing certificates and training for IF.

Problem 1 (Fairness Certification). Given a trained NN f w , a similarity dfair and a distance threshold ϵ ≥ 0, compute

Problem 1 provides a formulation in terms of optimisation, seeking to compute the maximum output change δmax for any pair of input points whose dfair distance is no more than ϵ. One can then compare δmax with any threshold δ: if δmax ≤ δ holds then the model f w has been certified to be ϵ-δ-IF.

While Problem 1 is concerned with an already trained NN, the methods we develop can also be employed to encourage IF at training time. Similarly to the approaches for adversarial learning [Goodfellow et al., 2015], we modify the training loss L(f w (x), y) to balance between the model fit and IF.

Problem 2 (Fairness Training). Consider an NN f w , a training set D, a similarity metric dfair and a distance threshold ϵ ≥ 0. Let λ ∈ [0, 1] be a constant. Define the IF-fair loss as

where x∗i = argmaxx∈X s.t. dfair(xi,x)≤ϵ |f w (xi) - f w (x)|. The ϵ-IF training problem is defined as finding wfair s.t.:

In Problem 2 we seek to train a NN that not only is accurate, but whose predictions are also fair according to Definition 1. Parameter λ balances between accuracy and IF. In particular, for λ = 1 we recover the standard training that does not account for IF, while for λ = 0 we only consider IF.

3 AMILP Approach For Individual Fairness

Certification of individual fairness on a NN thus requires us to solve the following global, non-convex optimisation problem:

We develop a Mixed-Integer Linear Programming (MILP) over-approximation (i.e., providing a sound bound) to this problem. We notice that there are two sources of non-linearity here, one induced by the NN (Equation (2)), which we refer to as the model constraint, and the other by the fairness metric (Equation (3)), which we call fairness constraint. In the following, we show how these can be modularly bounded by piecewise-linear functions. In Section 3.3 we bring the results together to derive a MILP formulation for ϵ-δ-IF.

3.1 Model Constraint

We develop a scheme based on piecewise-linear (PWL) upper and lower bounding for over-approximating all commonly used non-linear activation functions. An illustration of the PWL bound is given in Figure 1. Let ϕ(i)L j and ϕ(i)U j ∈ R be lower and upper bounds on the pre-activation ϕ(i) j . 4 We proceed by building a discretisation grid over the ϕ(i) j values on M grid points: ϕgrid = [ϕ(i) j,0 , . . . , ϕ(i) j,M ], with ϕ(i) j,0 := ϕ(i)L j and ϕ(i) j,M := ϕ(i)U j , such that, in each partition interval [ϕ(i) j,l , ϕ(i) j,l+1 ], we have that σ(i) is either convex or concave. We then compute linear lower and upper bound functions for σ(i) in each [ϕ(i) j,l , ϕ(i) j,l+1 ] as follows. If σ(i) is convex (resp. concave) in [ϕ(i) j,l , ϕ(i) j,l+1 ], then an upper (resp. lower) linear bound is given by the segment connecting the two extremum points of the interval, and a lower (resp. upper) linear bound is given by the tangent through the mid-point of the interval. We then compute the values of each linear bound in each of its grid points, and select the minimum of the lower bounds and the maximum of the upper bound values, which we store in two vectors ζ PWL,(i),U j = [ζ PWL,(i),U j,0 , . . . , ζ PWL,(i),U j,M ] and ζ PWL,(i),L j = [ζ PWL,(i),L j,0 , . . . , ζ PWL,(i),L j,M ]. The following lemma is a consequence of this construction.

4 Computed by bound propagation over X [Ehlers, 2017].

that is, ζ PWL,(i),L j and ζ PWL,(i),U j define continuous PWL lower and upper bounds for ϕ in [ϕ(i)L j , ϕ(i)U j ].

Lemma 3.1 guarantees that we can bound the non-linear activation functions using PWL functions. Crucially, PWL functions can then be encoded into the MILP constraints.

Proposition 1. Let y(i) j,l for l = 1, . . . ,M, be binary variables, and η(i) j,l ∈ [0, 1] be continuous ones. Consider ϕ(i) j ∈ [ϕ(i)L j , ϕ(i)U j ] then it follows that ζ (i) j = σ(i) (ϕ(i) j ) implies:

A proof can be found in Appendix A. Proposition 1 ensures that the global behaviour of each NN neuron can be over-approximated by 5 linear constraints using 2M auxiliary variables. Employing Proposition 1 we can encode the model constraint of Equation (2) into the MILP form in a sound way.

The over-approximation error does not depend on the MILP formulation (which is exact), but on the PWL bounding, and is hence controllable through the selection of the number of grid points M, and becomes exact in the limit. Notice that in the particular case of ReLU activation functions the over-approximation is exact for any M > 0.

Proposition 2. Assume σ(i) to be continuously differentiable everywhere in [ϕ(i)L j , ϕ(i)U j ], except possibly in a finite set. Then PWL lower and upper bounding functions of Lemma 3.1 converge uniformly to σ(i) as M goes to infinity.

Furthermore, define ∆M = (ϕ(i)U j - ϕ(i)L j )/M, then for finite values of M the error on the lower (resp. upper) bounding in convex (resp. concave) regions of σ(i) for ϕ ∈ [ϕ(i) j,l , ϕ(i) j,l+1 ] is given by:

and upper (resp. lower) in concave (resp. convex) regions:

A proof of Proposition 2 is given in Appendix A, alongside an experimental analysis of the convergence rate.

We remark that the PWL bound can be used over all commonly employed activation functions σ. The only assumption made is that σ has a finite number of inflection points over any compact interval of R. For convergence (Prop. 2) we require continuous differentiability almost everywhere, which is satisfied by commonly used activations.

3.2 Fairness Constraint

The encoding of the fairness constraint within the MILP formulation depends on the specific form of the metric dfair.

Weighted ℓp Metric: The weighted ℓp metric can be tackled by employing rectangular approximation regions. While this is straightforward for the ℓ∞ metric, for the remaining cases interval abstraction can be used [Dantzig, 2016].

Mahalanobis Metric: We first compute an orthogonal decomposition of S as in U T SU = Λ, where U is the eigenvector matrix of S and Λ is a diagonal matrix with S eigenvalues as entries. Consider the rotated variables z′ = U T x′ and z′′ = U T x′′ , then we have that Equation (3) can be re-written as (z′ -z′′ )T Λ(z′ -z′′ ) ≤ ϵ2 . By simple algebra we thus have that, for each i, (z′i -z′′ i )2 ≤ ϵ2 Λii . By transforming back to the original variables, we obtain that Equation (3) can be overapproximated by: - ϵ √diag(Λ) ≤ U T x′ -U T x′′ ≤ ϵ √diag(Λ) . Feature Embedding Metric We tackle the case in which φ used in the metric definition, i.e. dfair(x′ , x′′ ) = ˆ d(φ(x′ ), φ(x′′ )), is a NN embedding. This is straightforward as φ can be encoded into MILP as for the model constraint.

3.3 Overall Formulation

We now formulate the MILP encoding for the overapproximation δ∗ ≥ δmax of ϵ-δ-IF. For Equation (2), we proceed by deriving a set of approximating constraints for the variables x′ and x′′ by using the techniques described in Section 3.1. We denote the corresponding variables as ϕ′(i) j , ζ ′(i) j and ϕ′′(i) j , ζ ′′(i) j , respectively. The NN final output on x′ and on x′′ will then respectively be ζ ′(L) and ζ ′′(L) , so that δ = ζ ′(L) - ζ ′′(L) . Finally, we over-approximate Equation (3) as described in Section 3.2. In the case of Mahalanobis distance, we thus obtain:

Though similar, the above MILP is significantly different from those used for adversarial robustness (see e.g. [Tjeng et al., 2018]). First, rather than looking for perturbations around a fixed a point, here we have both x′ and x′′ as variables. Furthermore, rather than being local, the MILP problem for ϵ-δ-IF is global, over the whole input space X. As such, local approximations of non-linearities cannot be used, as the bounding needs to be valid simultaneously over the whole input space. Finally, while in adversarial robustness one can ignore the last sigmoid layer, for IF, because of the two optimisation variables, one cannot simply map from the last preactivation value to the class probability, so that even for ReLU NNs one needs to employ bounding of non-piecewise activations for the final sigmoid.

By combining the results from this section, we have:

Theorem 1. Consider ϵ ≥ 0, a similarity dfair and a NN f w . Let x′∗ and x′′ ∗ be the optimal points for the optimisation problem in Equation (4). Define δ∗ = |f w (x′∗ ) - f w (x′′ ∗ )|. Then f w is ϵ-δ-individually fair w.r.t. dfair for any δ ≥ δ∗.

Theorem 1, whose proof can be found in Appendix A, states that a solution of the MILP problem provides us with a sound estimation of individual fairness of an NN. Crucially, it can be shown that branch-and-bound techniques for the solution of MILP problems converge in finite time to the optimal solution [Del Pia and Weismantel, 2012], while furthermore providing us with upper and lower bounds for the optimal value at each iteration step. Therefore, we have:

Corollary 1. Let δL k and δU k lower and upper bounds computed by a MILP solver at step k > 0. Then we have that: δL k ≤ δ∗ ≤ δU k . Furthermore, given a precision, τ, there exist a finite k∗ such that δU k∗ - δL k∗ ≤ τ.

That is, our method is sound and anytime, as at each iteration step in the MILP solving we can retrieve a lower and an upper bound on δ∗, which can thus be used to provide provable guarantees while converging to δ∗ in finite time.

Complexity Analysis. The encoding of the model constraint can be done in O(LMnmax), where nmax is the maximum width of f w , L is the number of layers, and M is the number of grid points used for the PWL bound. The computational complexity of the fairness constraints depends on the similarity metric employed. While for ℓ∞ no processing needs to be done, the computational complexity is O(n3 ) for the Mahalanobis distance and again O(LMnmax) for the feature embedding metric. Each iteration of the MILP solver entails the solution of a linear programming problem and is hence O((MnmaxL)3 ). Finite time convergence of the MILP solver to δ∗ with precision τ is exponential in the number of problem variables, in τ and ϵ.

3.4 Fairness Training for Neural Networks

The ϵ-δ-IF MILP formulation introduced in Section 3 can be adapted for the solution of Problem 2. The key step is the computation of x∗i in the second component of the modified loss introduced in Problem 2, which is used to introduce fairness directly into the loss of the neural network. This computation can be done by observing that, for every training point xi drawn from D, the computation of x∗i = argmaxx∈X s.t. dx(xi,x)≤ϵ |f w (xi) - f w (x)| is a particular case of the formulation described in Section 3, where, instead of having two variable input points, only one input point is a problem variable while the other is given and drawn from the training dataset D. Therefore, x∗i can be computed by solving the MILP problem, where we fix a set of the problem variables to xi, and can be subsequently used to obtain the value of the modified loss function. Note that these constraints are not cumulative, since they are built for each minibatch, and discarded after optimization is solved to update the weights.

Algorithm 1 Fair Training with MILP.

Input: NN architecture: f w , Dataset: D, Learning rate: α, Iterations: nepoch, Batch Size: nbatch, Similarity metric: dfair, Maximum similarity: ϵ, Fairness Loss Weighting: λ.

Output: wfair: weight values balancing between accuracy and fairness.

B1: wfair ← InitWeights(f w ) 2: for t = 1, . . . , nepoch do 3: for b = 1, . . . , ⌈|D|/nbatch⌉ do 4: {X, Y } ← {xi, yi}nbatch i=0 ∼ D #Sample Batch 5: Yclean ← f w (X) #Standard forward pass 6: [ϕ′ , ζ ′ ,ϕ′′ , ζ ′′ ] ← InitMILP(f w , dfair, ϵ) # Section 3 7: XMILP ← ∅ 8: for i = 0, . . . nbatch do 9: ϕ′i, ζ ′ i ← F ixV arConst(xi) #Fix constraints 10: x∗i ←MILP(xi,ϕ′i, ζ ′ i) # Solve 'local' MILP prob. 11: XMILP ← XMILP ⋃{x∗i } 12: end for 13: YMILP ← f w (XMILP) #MILP inputs forward pass 14: l ← Lfair(Yclean, Y, YMILP, λ) #Fair Loss 15: wfair ←wfair -α∇wl #Optimizer step (here, SGD) 16: end for 17: end for 18: return wfair #Weights optimized for fairness & accuracy

We summarise our fairness training method in Algorithm 1. For each batch in each of the nepoch training epochs, we perform a forward pass of the NN to obtain the output, Yclean (line 5). We then formulate the MILP problem as in Section 3 (line 6), and initialise an empty set variable to collect the solutions to the various sub-problems (line 7). Then, for each training point xi in the mini-batch, we fix the MILP constraints to the variables associated with xi (line 9), solve the resulting MILP for x∗i , and place x∗i in the set that collects the solutions, i.e. XMILP. Finally, we compute the NN predictions on XMILP (line 13); the result is used to compute the modified loss function (line 14) and the weights are updated by taking a step of gradient descent. The resulting set of weights wfair balances the empirical accuracy and fairness around the training points.

The choice of λ affects the relative importance of standard training w.r.t. the fairness constraint: λ = 1 is equivalent to standard training, while λ = 0 only optimises for fairness. In our experiments we keep λ = 1 for half of the training epochs, and then change it to λ = 0.5.

4 Experiments

In this section, we empirically validate the effectiveness of our MILP formulation for computing ϵ-δ-IF guarantees as well as for fairness training of NNs. We perform our experiments on four UCI datasets: the Adult dataset (predicting income), the Credit dataset (predicting payment defaults), the German dataset (predicting credit risk) and the Crime dataset (predicting violent crime). In each case, features encoding information regarding gender or race are considered sensitive. In the certification experiments we employ a precision τ for the MILP solvers of 10-5 and a time cutoff of 180 seconds. We compare our training approach with two different learningmethods: Fairness-Through-Unawareness (FTU), in which the sensitive features are simply removed, and SenSR [Yurochkin et al., 2020]. Exploration of the cutoff, group fairness, certification of additional NNs, scalability of the methods and additional details on the experimental settings are given in Appendix C and D. 5

Fairness Certification. We analyse the suitability of our method in providing non-trivial certificates on ϵ-δ-IF with respect to the similarity threshold ϵ (which we vary from 0.01 to 0.25), the similarity metric dfair , the width of the NN (from 8 to 64), and its number of layers (from 1 to 4). These reflect the characteristics of NNs and metrics used in the IF literature [Yurochkin et al., 2020; Ruoss et al., 2020; Urban et al., 2020]; for experiments on larger architectures, demonstrating the scalability of our approach, see Appendix D.3. For each dataset we train the NNs by employing the FTU approach. The results for these analyses are plotted in Figure 2 for the Adult and the Crime datasets (results for Credit and German datasets can be found in Appendix D.1). Each heat map depicts the variation of δ∗ as a function of ϵ and the NN architecture. The top row in the figure was computed by considering the Mahalanobis similarity metric; the bottom row was computed for a weighted ℓ∞ metric (with coefficients chosen as in [John et al., 2020]) and results for the feature embedding metrics are given in Appendix D.2. As one might expect, we observe that, across all the datasets and architectures, increasing ϵ correlates with an increase in the values for δ∗, as higher values of ϵ allow for greater feature changes. Interestingly, δ∗ tends to decrease (i.e., the NN becomes more fair) as we increase the number of NN layers. This is the opposite to what is observed for the adversarial robustness, where increased capacity generally implies more fragile models [Madry et al., 2018]. In fact, as those NNs are trained via FTU, the main sensitive features are not accessible to the NN. A possible explanation is that, as the number of layers increases, the NN's dependency on the specific value of each feature diminishes, and the output becomes dependent on their nonlinear combination. The result suggests that overparametrised NNs could be more adept at solving IF tasks though this would come with a loss of model interpretability, and exploration would be needed to assess under which condition this holds. Finally, we observe that our analysis confirms how FTU training is generally insufficient in providing fairness on the model behaviour for ϵ-δ-IF. For each model, individuals that are dissimilar by ϵ ≥ 0.25 can already yield a δ∗ > 0.5, meaning they would get assigned to different classes if one was using the standard classification threshold.

5 An implementation of the method and of the experiments can be found at https://github.com/eliasbenussi/nn-cert-individual-fairness.

Fairness Training. We investigate the behaviour of our fairness training algorithm for improving ϵ-δ-IF of NNs. We compare our method with FTU and SenSR [Yurochkin et al., 2020]. For ease of comparison, in the rest of this section we measure fairness with dfair equal to the Mahalanobis similarity metric, with ϵ = 0.2, for which SenSR was developed. The results for this analysis are given in Figure 3, where each point in the scatter plot represents the values obtained for a given NN architecture. We train architectures with up to 2 hidden layers and 64 units, in order to be comparable to those trained by [Yurochkin et al., 2020]. As expected, we observe that FTU performs the worst in terms of certified fairness, as simple omission of the sensitive features is unable to obfuscate latent dependencies between the sensitive and non-sensitive features. As previously reported in the literature, SenSR significantly improves on FTU by accounting for features latent dependencies. However, on all four datasets, our MILP-based training methodology consistently improves IF by orders of magnitude across all the architectures when compared to SenSR. In particular, for the architectures with more than one hidden layer, on average, MILP outperforms FTU by a factor of 78598 and SenSR by 27739. Intuitively, while SenSR and our approach have a similar formulation, the former is based on gradient optimisation so that no guarantees are provided in the worst case for the training loss. In contrast, by relying on MILP, our method optimises the worst-case behaviour of the NN at each step, which further encourages training of individually fair models. The cost of the markedly improved guarantees is, of course, a higher computational costs. In fact, the training of the models in Figure 3 with MILP had an average training time of about 3 hours. While the increased cost is significant, we highlight that this is a cost that is only paid once and may be justified in sensitive applications by the necessity of fairness at deployment time. We furthermore notice that, while our implementation is sequential, parallel per-batch solution of the MILP problems during training would markedly reduce the computational time and leave for future work the parallelisation and tensorisation of the techniques. Interestingly, we find that balanced accuracy also slightly improved with SenSR and MILP training in the tasks considered here, possibly as a result of the bias in the class labels w.r.t. sensitive features. Finally, in Figure 4 we further analyse the certified δ∗-profile w.r.t. to the input similarity ϵ, varying the value of ϵ used in for the ceritification of ϵ-δ-IF. In the experiment, both SenSR and MILP are trained with ϵ = 0.2, which means that our method, based on formal IF certificates, is guaranteed to outperform SenSR up until ϵ = 0.2 (as in fact is the case). Beyond 0.2, no such statement can be made, and it is still theoretically possible for SenSR to outperform MILP in particular circumstances. Empirically, however, MILP-based training still largely outperforms SenSR in terms of certified fairness obtained.

5 Conclusion

We introduced an anytime MILP-based method for the certification and training of ϵ-δ-IF in NNs, based on PWL bounding and MILP encoding of non-linearities and similarity metrics. In an experimental evaluation comprising four datasets, a selection of widely employed NN architectures and three types of similarity metrics, we found that our method is able to provide the first non-trivial certificates for ϵ-δ-IF in NNs and yields NNs which are orders of magnitude more fair than those obtained by a competitive techniques.

Acknowledgements

This project was funded by the ERC European Union's Horizon 2020 research and innovation programme (FUN2MODEL, grant agreement No. 834115).