# A novel multigrid method for elliptic distributed control problems

Yunhui He Department of Computer Science, The University of British Columbia,
Vancouver, BC, V6T 1Z4, Canada<EMAIL_ADDRESS>
###### Abstract

Large linear systems of saddle-point type have arisen in a wide variety of
applications throughout computational science and engineering. The
discretizations of distributed control problems have a saddle-point structure.
The numerical solution of saddle-point problems has attracted considerable
interest in recent years. In this work, we propose a novel Braess-Sarazin
multigrid relaxation scheme for finite element discretizations of the
distributed control problems, where we use the stiffness matrix obtained from
the five-point finite difference method for the Laplacian to approximate the
inverse of the mass matrix arising in the saddle-point system. We apply local
Fourier analysis to examine the smoothing properties of the Braess-Sarazin
multigrid relaxation. From our analysis, the optimal smoothing factor for
Braess-Sarazin relaxation is derived. Numerical experiments validate our
theoretical results. The relaxation scheme considered here shows its high
efficiency and robustness with respect to the regularization parameter and
grid size.

Keywords. Multigrid methods, distributed optimal control, saddle-point
problem, Braess-Sarazin relaxation, local Fourier analysis

2000 MSC: 49M25, 49K20, 65N55, 65F10

##  1 Introduction

In this paper we consider the following two-dimensional elliptic optimal
control problem [11]: Find the state $u^{*}\in H^{1}(\Omega)$ and the control
$f^{*}\in L^{2}(\Omega)$ such that

$\mathcal{J}(u^{*},f^{*})=\displaystyle\min_{(u,f)\in H^{1}(\Omega)\times
L^{2}(\Omega)}\mathcal{J}(u,f),$ (1)

subject to the state equation

$\left\\{\begin{aligned} -\triangle u&=&f,\,\,\,\,{\rm in}\,\,\Omega,\\\
u&=&g,\,\,\,{\rm at}\,\,\,\partial\Omega,\\\ \end{aligned}\right.\\\ $ (2)

with cost functional

$\mathcal{J}(u,f)=\frac{1}{2}\|u-\hat{u}\|^{2}_{L^{2}(\Omega)}+\beta\|f\|^{2}_{L^{2}(\Omega)},$
(3)

where $\hat{u}$ is the desired state and $\beta>0$ is the weight of the cost
of the control (or simply a regularization parameter).

The discretization of the optimization problem (1) leads to a large scale
saddle-point system. Since multigrid methods offer the possibility of solving
problems with $N$ unknowns using $O(N)$ work and storage, and lead to
substantial improvements in computational efficiency over direct methods, they
have gained growing interest in the area of optimal control problems, see, for
example [1, 2, 4, 6, 12, 13, 14].

The choice of multigrid smoother is crucial to construct efficient multigrid
algorithms. Local Fourier analysis (LFA) is a useful tool to help choose
multigrid components, such as relaxation schemes and grid-transfer operators,
and can quantitatively predict the two-grid convergence factor. The main goal
of this work is to construct and analyze multigrid methods for the discrete
optimal control problems. Note that there are a few studies of LFA applied to
the optimal control problems. In [9], the smoothing factor for a one-
dimensional optimal control problem was proposed by LFA for a special
smoother. [4] reviewed some recent efforts and outlined recent developments in
the field of multigrid methods for partial differential equation optimization,
where LFA was applied to study collective Gauss–Seidel methods. The finite
difference multigrid solution of an optimal control problem associated with an
elliptic equation was considered in [3], and LFA was used to estimate the
convergence factor for the collective Gauss–Seidel relaxation.

In this work, we consider finite element monolithic multigrid solution of the
elliptic distributed control problem (1), that is, employing multigrid methods
directly to the discrete system in coupled form. In contrast to collective
smoothers, we consider a well-known Braess-Sarazin-type relaxation (BSR)
scheme, which was original designed for the Stokes equations [5]. In our
recently work [7], we have proved that the mass matrix obtained from bilinear
finite elements is a good approximation to the inverse of the scalar Laplacian
discretized by a five-point finite difference method, where the optimal
smoother factor of $\frac{1}{3}$ is derived for the Laplacian. Furthermore, we
extended this mass-based approximation to a Braess-Sarazin-type relaxation
scheme for the Stokes equations [8], where we obtained the optimal smoother
factor of $\frac{1}{3}$ for Stokes. This motivates us to consider the
stiffness matrix derived from the finite-difference discrerization of the
scalar Laplacian to approximate the mass matrix arsing from the finite element
discretizations of optimal control problems. In the resulting BSR relaxation,
we can directly apply the stiffness matrix and there is no need to invert a
matrix. LFA is applied to examine this stiffness-based BSR applied to the
optimal control problems considered here. We derive optimal smoothing factor
of $\frac{1}{3}$ for the stiffness-based BSR for the optimal control problems.
Considering practical use, an inexact version of stiffness-based Braess-
Sarazin relaxation is developed, where a 3 $V(2,2)$-cycles multigrid with
weighted Jacobi relaxation is used to approximate the solution of Schur
complement system appeared in the exact BSR process. This inexact version
preserves the optimal smoothing factor of $\frac{1}{3}$ obtained from exact
BSR. The LFA convergence estimates agree very well with results of numerical
experiments and are independent of the grid size and of the value of the
control parameter.

This paper is organized as follows. In Section 2, we propose a stiffness-based
Braess-Sarazin relaxation scheme for the elliptic distributed control problem.
In Section 3, we apply LFA to the proposed stiffness-based Braess-Sarazin
relaxation, and we obtain a highly satisfactory convergence factor. Numerical
results are reported in Section 4 to confirm our theoretical analysis.
Finally, we draw conclusions in Section 5.

##  2 Discretization and relaxation

In our work, we discretize problem (1) using bilinear quadrilateral $Q_{1}$
finite elements. The discretize formulations of (2) and (3) are given by

$\displaystyle\min_{u_{h},f_{h}}\frac{1}{2}\|u_{h}-\hat{u}\|^{2}_{L^{2}(\Omega)}+\beta\|f_{h}\|^{2}_{L^{2}(\Omega)},$
(4) $\displaystyle{\rm such\,\,that}\int\nabla u_{h}v_{h}=\int
v_{h}f_{h},\,\,\forall v_{h}\in V_{0}^{h},$ (5)

where $V_{0}^{h}$ is the $Q_{1}$ finite element space. With a Lagrange
multiplier ($\tau$) approach to (4) and (5), it typically leads to a linear
system of the form [11]

$\begin{pmatrix}2\beta M&0&-M^{T}\\\ 0&M&K^{T}\\\
-M&K&0\end{pmatrix}\begin{pmatrix}f_{h}\\\ u_{h}\\\
\tau_{h}\end{pmatrix}=b_{h},$ (6)

where $K$ is the stiffness matrix (the discrete Laplacian) and $M$ is the
discrete mass matrix. For more details, we refer to [11].

We can rewrite system (6) as a saddle-point system, that is,

$\mathcal{L}_{h}z=\begin{pmatrix}A&B^{T}\\\ B&0\\\
\end{pmatrix}\begin{pmatrix}x\\\ y\end{pmatrix}=b_{h},$ (7)

where

$A=\begin{pmatrix}2\beta M&0\\\ 0&M\end{pmatrix},\quad
B=\begin{pmatrix}-M&K\end{pmatrix}.$

Here, we consider multigrid methods for solving (7). One important process in
multgrid is the smoothing, where we update a current approximation $z_{k}$ via

$\hat{z}_{k}=z_{k}-\omega\mathcal{K}^{-1}_{h}(b_{h}-\mathcal{L}_{h}z_{k})=z_{k}-\omega\delta
z,$ (8)

where $\mathcal{K}^{-1}_{h}$ is an approximation to $\mathcal{L}_{h}$, and
$\delta z=\mathcal{K}^{-1}_{h}(b_{h}-\mathcal{L}_{h}z_{k})$. Then, the error-
propagation for the relaxation scheme (8) is

$\mathcal{S}_{h}=I-\omega\mathcal{K}^{-1}_{h}\mathcal{L}_{h}.$ (9)

The choices of $\mathcal{K}_{h}$ and $\omega$ play a crucial role in
determining the convergence speed of multigrid methods. For simplicity,
throughout the rest of this paper, we drop the subscript $h$, except when
necessary for clarity.

In this work, we consider Braess-Sarazin relaxation [5], that is,
$\mathcal{K}$ has the form of

$\mathcal{K}=\begin{pmatrix}\alpha C&B^{T}\\\ B&0\\\ \end{pmatrix},$ (10)

where $C$ is an approximation of $A$, whose inverse is easy to apply, and
$\alpha>0$.

With $\mathcal{K}$ defined in (10), we compute the update $\delta z=(\delta
x,\delta y)$ in (8) by the following two stages

$\displaystyle(BC^{-1}B^{T})\delta y$ $\displaystyle=$ $\displaystyle
BC^{-1}r_{x}-\alpha r_{y},$ (11) $\displaystyle\delta x$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha}C^{-1}(r_{x}-B^{T}\delta y),$ (12)

where $(r_{x},r_{y})=b_{h}-\mathcal{L}z_{k}$.

Note that, in (11) and (12), we can directly consider an approximation to
$C^{-1}$. In our recently work [7], we have shown that the mass matrix $M$
derived from bilinear finite elements is a good approximation to the inverse
of the scalar Laplacian discretized by the five-point finite difference
method. This inspires us to consider the discrete scalar Laplacian
approximating the inverse of the mass matrix $M$ in the diagonal part of the
coefficient matrix defined in (6). Let $A_{fd}$ be the stiffness matrix
obtained from the five-point finite difference discretization for the scalar
Laplacian. In (10), we consider

$C^{-1}=C^{-1}_{f}=\begin{pmatrix}\frac{1}{2\beta}A_{fd}&0\\\
0&A_{fd}\end{pmatrix},$ (13)

and denote the corresponding $\mathcal{K}$ as

$\mathcal{K}_{f}=\begin{pmatrix}\alpha C_{f}&B^{T}\\\ B&0\\\ \end{pmatrix}.$
(14)

To choose suitable values of the parameters $\omega$ and $\alpha$, we consider
LFA to examine the relaxation scheme with $\mathcal{K}=\mathcal{K}_{f}$ in the
following section.

##  3 Smoothing analysis

LFA [15, 16] has been widely used to analyze actual multigrid performance for
different types of problems, by investigating the spectral radius of symbol of
the underlying operator. Specifically, the LFA smoothing factor with
simplifying assumptions on boundary conditions offers a sharp convergence
estimate of multigrid methods. Thus, in this work, we employ LFA to examine
the LFA smoothing factor for BSR, and help choose relaxation parameters. Next,
we first give a brief introduction to LFA.

We consider two-dimensional infinite uniform grids $\mathbf{G}_{h}$

$\mathbf{G}_{h}=\left\\{\boldsymbol{x}:=(x_{1},x_{2})=\boldsymbol{k}h=(k_{1},k_{2})h,k\in\mathbb{Z}\right\\},$

and Fourier modes
$\varphi(\boldsymbol{\theta},\boldsymbol{x})=e^{\imath\boldsymbol{\theta}\boldsymbol{x}/h}$
on $\mathbf{G}_{h}$, where $\boldsymbol{\theta}=(\theta_{1},\theta_{2})$ and
$\imath^{2}=-1$. Since $\varphi(\boldsymbol{\theta},\boldsymbol{x})$ is
periodic in $\boldsymbol{\theta}$ with period $2\pi$, we consider
$\theta_{i}\in\big{(}-\frac{\pi}{2},\frac{3\pi}{2}\big{]}$. The coarse grid
$\mathbf{G}_{2h}$ is defined similarly.

Let $L_{h}$ be a Toeplitz operator acting on $\mathbf{G}_{h}$ as follows [15],

$\displaystyle L_{h}w_{h}(\boldsymbol{x})$ $\displaystyle=$
$\displaystyle\sum_{\boldsymbol{\kappa}\in\boldsymbol{V}}s_{\boldsymbol{\kappa}}w_{h}(\boldsymbol{x}+\boldsymbol{\kappa}h),$

with constant coefficients $s_{\boldsymbol{\kappa}}\in\mathbb{R}$, where
$w_{h}(\boldsymbol{x})$ is a function in $l^{2}(\mathbf{G}_{h})$. Here,
$\boldsymbol{V}$ is a finite index set.

######  Definition 3.1.

We call
$\widetilde{L}_{h}(\boldsymbol{\theta})=\displaystyle\sum_{\boldsymbol{\kappa}\in\boldsymbol{V}}s_{\boldsymbol{\kappa}}e^{\iota\boldsymbol{\theta}\cdot\boldsymbol{\kappa}}$
the symbol of $L_{h}$.

We consider multigrid methods for finite-element discretizations with standard
geometric grid coarsening, that is, we construct a sequence of coarse grids by
doubling the mesh size in each spatial direction. High and low frequencies for
standard coarsening are given by

$\boldsymbol{\theta}\in T^{{\rm
low}}=\left[-\frac{\pi}{2},\frac{\pi}{2}\right)^{2},\,\boldsymbol{\theta}\in
T^{{\rm
high}}=\displaystyle\left[-\frac{\pi}{2},\frac{3\pi}{2}\right)^{2}\bigg{\backslash}\left[-\frac{\pi}{2},\frac{\pi}{2}\right)^{2}.$

######  Definition 3.2.

The LFA smoothing factor for the error-propagation for the relaxation scheme
$\mathcal{S}$, see (9), is defined as

$\mu_{{\rm loc}}(\boldsymbol{p})=\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}\big{\\{}\big{|}\lambda(\widetilde{\mathcal{S}}(\boldsymbol{\theta}))\big{|}\,\,\big{\\}},$
(15)

where $\boldsymbol{p}$ is algorithmic parameters, and
$\lambda\big{(}\widetilde{\mathcal{S}}(\boldsymbol{\theta})\big{)}$ denotes
the eigenvalue of symbol $\widetilde{\mathcal{S}}(\boldsymbol{\theta})$.

Note that for the $Q_{1}$ elements considered in this work,
$\widetilde{\mathcal{S}}(\boldsymbol{\theta})$ is a $3\times 3$ matrix. Since
$\mu_{{\rm loc}}(\boldsymbol{p})$ is a function of $\boldsymbol{p}$, we can
minimize $\mu_{{\rm loc}}(\boldsymbol{p})$ to obtain a fast convergence speed.
Thus, we define the optimal smoothing factor as follows.

######  Definition 3.3.

The optimal LFA smoothing factor for the error-propagation for the relaxation
scheme $\mathcal{S}$ is defined as

$\mu_{{\rm opt}}=\min_{\boldsymbol{p}\in\Upsilon}{\mu_{{\rm
loc}}}(\boldsymbol{p}),$

where $\Upsilon$ is the set of allowable parameters.

Another important LFA factor is the LFA two-grid convergence factor. In
general, the two-grid error-propagation can be expressed as

$E_{h}=\mathcal{S}^{\nu_{2}}(I-P(\mathcal{L}_{2h})^{-1}R\mathcal{L}_{h})\mathcal{S}^{\nu_{1}},$

where $\mathcal{L}_{2h}$ is the coarse grid operator, $R$ and $P$ are the
restriction and interpolation operators, respectively. The integers $\nu_{1}$
and $\nu_{2}$ are pre- and post- smoothing steps, respectively.

######  Definition 3.4.

The LFA two-grid convergence factor for $E_{h}$ is defined as

$\rho=\max_{\boldsymbol{\theta}\in T^{{\rm
low}}}\left\\{\rho(\widetilde{E}_{h}(\boldsymbol{\theta}))\right\\},$ (16)

where $\widetilde{E}_{h}(\boldsymbol{\theta})$ is the two-grid LFA symbol of
$E_{h}$, and $\rho(\widetilde{E}_{h})$ denotes the spectral radius of matrix
$\widetilde{E}_{h}$.

Note that in our case, $\widetilde{E}_{h}(\boldsymbol{\theta})$ is a $12\times
12$ matrix. For more details on how to compute two-grid symbol, we refer to
[15]. Since it is easy to compute the LFA smoothing factor, (15), compared to
the LFA two-grid convergence factor, (16), we focus on the analysis of LFA
smoothing factor in the following.

###  3.1 Fourier representation of mass and stiffness operators

In one dimension (1D), the stiffness and mass stencils of $Q_{1}$ elements are

$K_{1}=\frac{1}{h}\begin{bmatrix}-1&2&-1\end{bmatrix}\quad\text{and}\quad
M_{1}=\frac{h}{6}\begin{bmatrix}1&4&1\end{bmatrix},$ (17)

respectively.

Then, the stiffness stencil for bilinear discretization in two dimensions (2D)
is given by

$\displaystyle K$ $\displaystyle=$ $\displaystyle K_{1}\otimes
M_{1}+M_{1}\otimes K_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\begin{bmatrix}-1&-1&-1\\\ -1&8&-1\\\
-1&-1&-1\end{bmatrix},$

and the mass stencil for bilinear discretization in 2D is given by

$\displaystyle M$ $\displaystyle=$ $\displaystyle M_{1}\otimes M_{1}$ (18)
$\displaystyle=$ $\displaystyle\frac{h^{2}}{36}\begin{bmatrix}1&4&1\\\
4&16&4\\\ 1&4&1\end{bmatrix}.$

From (17), the symbols of the stiffness and mass stencils in 1D are

$\widetilde{K}_{1}(\theta)=\frac{2}{h}(1-\cos\theta),\,\,\,\widetilde{M}_{1}(\theta)=\frac{h}{3}(2+\cos\theta).$
(19)

Based on (19), the symbol of $K$ can be written as

$\displaystyle\widetilde{K}(\theta_{1},\theta_{2})$ $\displaystyle=$
$\displaystyle\widetilde{K}_{1}(\theta_{2})\widetilde{M}_{1}(\theta_{1})+\widetilde{M}_{1}(\theta_{2})\widetilde{K}_{1}(\theta_{1})$
$\displaystyle=$
$\displaystyle\frac{2}{3}(4-\cos\theta_{1}-\cos\theta_{2}-2\cos\theta_{1}\cos\theta_{2}),$

and the symbol of $M$ is

$\displaystyle\widetilde{M}(\theta_{1},\theta_{2})$ $\displaystyle=$
$\displaystyle\widetilde{M}_{1}(\theta_{2})\widetilde{M}_{1}(\theta_{1})$
$\displaystyle=$
$\displaystyle\frac{h^{2}}{9}(4+2\cos\theta_{1}+2\cos\theta_{2}+\cos\theta_{1}\cos\theta_{2}).$

We introduce the stencil of the Laplacian $-\triangle_{h}$ discretized by the
standard five-point scheme, given by

$A_{fd}=-\triangle_{h}=\frac{1}{h^{2}}\begin{bmatrix}&-1&\\\ -1&4&-1\\\
&-1&\end{bmatrix}.$

Then, the symbol of $A_{fd}$ is

$\widetilde{A}_{fd}=\frac{4-2\cos\theta_{1}-2\cos\theta_{2}}{h^{2}}.$ (20)

For simplicity, let $a=\widetilde{M}$ and $\hat{a}=(\widetilde{A}_{fd})^{-1}$.
Using (20) gives

$\widetilde{A}_{fd}\widetilde{M}=\widetilde{M}\widetilde{A}_{fd}=\frac{a}{\hat{a}}.$
(21)

###  3.2 LFA for Braess-Sarazin relaxation scheme

Now, we examine Braess-Sarazin-type algorithms, which were originally
developed as a relaxation scheme for the Stokes equations [5]. In practice,
(11) is not solved exactly. Thus, in the following, we first consider exact
solve for (11), and then inexact one.

Exact Braess-Sarazin relaxation: In (10), we consider $C^{-1}$ defined in (13)
and $\alpha=1$. We first derive the corresponding optimal smoothing factor for
the exact Braess-Sarazin relaxation. We refer to the relaxation scheme as
stiffness-based BSR. Let $b=\widetilde{K}(\theta_{1},\theta_{2})$. The symbol
of operator $\mathcal{L}$ defined in (7) is given by

$\widetilde{\mathcal{L}}(\theta_{1},\theta_{2})=\begin{pmatrix}2\beta
a&0&-a\\\ 0&a&b\\\ -a&b&0\end{pmatrix},$

and

$\widetilde{\mathcal{K}}_{f}(\theta_{1},\theta_{2})=\begin{pmatrix}2\beta\hat{a}&0&-a\\\
0&\hat{a}&b\\\ -a&b&0\end{pmatrix}.$

Then,

$\widetilde{\mathcal{L}}-\lambda\widetilde{\mathcal{K}}_{f}=\begin{pmatrix}2\beta(a-\lambda\hat{a})&0&-a(1-\lambda)\\\
0&a-\lambda\hat{a}&b(1-\lambda)\\\ -a(1-\lambda)&b(1-\lambda)&0\end{pmatrix}.$

The determinant of
$\widetilde{\mathcal{L}}-\lambda\widetilde{\mathcal{K}}_{f}$ is

$\displaystyle|\widetilde{\mathcal{L}}-\lambda\widetilde{\mathcal{K}}_{f}|$
$\displaystyle=-2\beta(a-\lambda\hat{a})b^{2}(1-\lambda)^{2}-a^{2}(1-\lambda)^{2}(a-\lambda\hat{a})$
$\displaystyle=(1-\lambda)^{2}(a-\lambda\hat{a})(-2\beta b^{2}-a^{2}).$

It follows that the eigenvalues of
$\widetilde{\mathcal{K}}^{-1}_{f}\widetilde{\mathcal{L}}$ are $1,1$ and
$\frac{a}{\hat{a}}$.

We restate the results for solving the Laplace problem using mass matrix
approximation [7] in the following.

######  Lemma 3.1.

[7] For $\boldsymbol{\theta}\in T^{high}$,
$\frac{a}{\hat{a}}(\boldsymbol{\theta})\in[8/9,16/9]$.

######  Lemma 3.2.

[7] If we consider the mass matrix $M$ to approximate the inverse of the
scalar Laplacian, $A_{fd}$, and $\mathcal{S}=I-\omega MA_{fd}$, then the
optimal smoothing factor for $\mathcal{S}$ is

$\mu_{\rm opt}=\min_{\omega}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}\left|1-\omega\widetilde{M}\widetilde{A}_{fd}\right|=\frac{1}{3},$

where the minimum is uniquely achieved at $\omega=\omega_{\rm
opt}=\frac{3}{4}$.

Since
$\widetilde{M}\widetilde{A}_{fd}=\widetilde{A}_{fd}\widetilde{M}=\frac{a}{\hat{a}}$,
see (21), we have the following result.

######  Corollary 3.1.

If we consider the discrete Laplacian $A_{fd}$ to approximate the inverse of
the mass matrix, $M$, and $\mathcal{S}_{fd}=I-\omega A_{fd}M$, then the
optimal smoothing factor for $\mathcal{S}_{fd}$ is

$\mu_{\rm opt}=\min_{\omega}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}\left|1-\omega\widetilde{A}_{fd}\widetilde{M}\right|=\frac{1}{3},$

provided that $\omega=\omega_{\rm opt}=\frac{3}{4}$.

Since $\widetilde{\mathcal{K}}^{-1}_{f}\widetilde{\mathcal{L}}$ contains
eigenvalue $\frac{a}{\hat{a}}$, which is an eigenvalue of
$\widetilde{A}_{fd}\widetilde{M}$, Corollary 3.1 indicates that $\frac{1}{3}$
is a lower bound on the optimal smoothing factor for the stiffness-based BSR.

######  Theorem 3.1.

The optimal smoothing factor for the stiffness-based BSR for the optimal
control system (7) is

$\mu_{{\rm opt},B}=\displaystyle\min_{\omega_{B}}\max_{\boldsymbol{\theta}\in
T^{{\rm
high}}}\left|\lambda(\widetilde{\mathcal{S}}_{B}(\omega_{B},\boldsymbol{\theta}))\right|=\frac{1}{3},$

where the minimum is attained provided that $\omega_{B}=\omega_{B,\rm
opt}=\frac{3}{4}$.

###### Proof.

Since the eigenvalues of
$\widetilde{\mathcal{S}}_{B}=I-\omega_{B}\widetilde{\mathcal{K}}^{-1}_{f}\widetilde{\mathcal{L}}$
are $1-\omega_{B},1-\omega_{B}$ and $1-\omega_{B}\frac{a}{\hat{a}}$ and Lemma
3.1, the smoothing factor for $\mathcal{\widetilde{S}}_{B}$ is

$\mu_{{\rm
loc}}=\max\left\\{|1-\omega_{B}|,\left|1-\frac{16}{9}\omega_{B}\right|,\left|1-\frac{8}{9}\omega_{B}\right|\right\\},$

From Corollary 3.1, we know that

$\min_{\omega_{B}}\max\left\\{\left|1-\frac{16}{9}\omega_{B}\right|,\left|1-\frac{8}{9}\omega_{B}\right|\right\\}=\frac{1}{3},$

provided that $\omega_{B}=\omega_{B,\rm opt}=\frac{3}{4}$. Furthermore,
$|1-\omega_{B,\rm opt}|=\frac{1}{4}<\frac{1}{3}$. Thus,
$\min_{\omega_{B}}\mu_{{\rm loc}}=\frac{1}{3}$. ∎

######  Remark 3.1.

We point out that one can consider $\alpha$ as a variable in (14) rather than
$\alpha=1$. However, this does not change the optimal smoothing factor for
exact BSR obtained in Theorem 3.1, and only changes the choice of optimal
parameter $\omega_{B}$. In contrast to exact BSR, varying $\alpha$ is needed
for inexact version of BSR, which will be discussed in Section 4.

Inexact Braess-Sarazin relaxation: Exact Braess-Sarazin algorithm requires an
exact inversion of the Schur complement (11), which is very expensive. In the
literature, it has been shown that inexact solve using a few sweeps of
weighted Jacobi iteration for the Schur complement system is enough to
maintain the same convergence of exact BSR for the Stokes equations, see, for
example [8]. In (11), the Schur complement matrix is

$\displaystyle BC^{-1}B^{T}$
$\displaystyle=\begin{pmatrix}-M&K\end{pmatrix}\begin{pmatrix}\frac{1}{2\beta}A_{fd}&0\\\
0&A_{fd}\end{pmatrix}\begin{pmatrix}-M^{T}\\\ K^{T}\end{pmatrix}$
$\displaystyle=\frac{1}{2\beta}MA_{fd}M^{T}+KA_{fd}K^{T}.$

We explore inexact stiffness-based BSR, using either a few sweeps of weighted
Jacobi iteration to approximate the solution of (11) or a $V$-cycles multigrid
on the Schur complement system. We are wondering whether the inexact version
is possible to achieve the same smoothing factor of $\frac{1}{3}$. Here, we do
not seek theoretical analysis. Instead, we numerically examine the performance
of inexact BSR in Section 4.

Note that another choice for $\mathcal{K}$ is

$\mathcal{K}=\mathcal{K}_{D}=\begin{pmatrix}C_{D}&B^{T}\\\ B&0\end{pmatrix},$
(22)

where

$C_{D}=\begin{pmatrix}2\beta{\rm diag}(M)&0\\\ 0&{\rm diag}(M)\end{pmatrix}.$

From (18), it can be seen that the symbol of ${\rm diag}(M)$ is
$\frac{4h^{2}}{9}=:a_{D}$. Here, we are curious about the optimal smoothing
factor for the smoother (22). As a comparison, we derive the optimal smoothing
factor for this choice.

######  Theorem 3.2.

If we consider $\mathcal{K}_{D}$ defined in (22) and
$\mathcal{S}_{D}=I-\omega_{D}\mathcal{K}^{-1}_{D}\mathcal{L}$, them the
optimal smoothing factor for $\mathcal{S}_{D}$ is

$\mu_{{\rm opt},D}=\min_{\omega_{D}}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}\left|\lambda(\widetilde{\mathcal{S}}_{D}(\omega_{D},\boldsymbol{\theta}))\right|=\frac{5}{7},$

provided that $\omega_{D}=\frac{8}{7}$.

###### Proof.

It can easily be shown that the eigenvalues of $\widetilde{\mathcal{S}}_{D}$
are $1,1$ and $\frac{a}{a_{D}}$. Now, we consider the optimal smoothing factor
for the modes, $\frac{a}{a_{D}}$. Note that

$\frac{a}{a_{D}}=\frac{1}{4}(4+2\cos\theta_{1}+2\cos\theta_{2}+\cos\theta_{1}\cos\theta_{2}).$

For $\boldsymbol{\theta}\in T^{\rm high}$, it can easily be shown that the
extreme values of $\frac{a}{a_{D}}$ are $\frac{1}{4}$ and $\frac{3}{2}$. Thus,

$\mu_{\rm opt}=\min_{\omega_{D}}\max_{\boldsymbol{\theta}\in T^{{\rm
high}}}\left\\{|1-\omega_{D}|,\left|1-\frac{1}{4}\omega_{D}\right|,\left|1-\frac{3}{2}\omega_{D}\right|\right\\}=\frac{5}{7},$

provided that $\omega_{D}=\frac{2}{1/4+3/2}=\frac{8}{7}$. ∎

From Theorems 3.1 and 3.2, we see that the stiffness-based BSR offers a highly
satisfactory convergence speed. Thus, in our numerical tests, we only consider
the stiffness-based BSR.

##  4 Numerical experiments

In this section, we first present the LFA two-grid convergence to compare with
the optimal LFA smoothing factor derived in Subsection 3.2 for the exact
stiffness-based BSR. Second, we display actual multigrid performance to
validate our analytical result for exact BSR. To consider practical use, we
develop an inexact version of BSR, where a $V$-cycle multigrid is applied to
the Schur complement system to maintain the performance observed for exact
BSR. To demonstrate the robustness of our multigrid methods with respect to
the regularization parameters, we consider different values of $\beta$ in our
tests.

###  4.1 LFA predictions

In this subsection, we report the LFA predicted two-grid convergence factor,
$\rho_{h}$, for different values of regularization parameter $\beta$ and grid
size $h$ for exact stiffness-based BSR. For the prolongation of correction, we
consider the standard bilinear interpolation and apply the corresponding
adjoint operator for the restriction. Galerkin operator is used to define the
coarse grid operator. The number of smoothing steps $\nu$ is
$\nu=\nu_{1}+\nu_{2}$. We use $32\times 32$ evenly distributed Fourier
frequencies in the Fourier domain
$[-\frac{\pi}{2}+\tau,\frac{\pi}{2}-\tau]^{2}$ with $\tau=\frac{\pi}{64}$ to
compute LFA two-grid convergence factor $\rho_{h}$ in (16).

In Table 1, we show the LFA predicted two-grid convergence factor as a
function of $\beta$ and $h$ with $\nu=1$. We see that the LFA two-grid
convergence factor is the same as the optimal smoothing factor derived from
our theoretical analysis. Moreover, the convergence factor is robust in
regularization parameter $\beta$ and is $h$-independent.

Table 1: LFA predicted two-grid convergence factor, $\rho_{h}$, as a function of $\beta$ and $h$ with $\nu=1$ for exact stiffness-based BSR. $h$ | $\beta=10^{-2}$ | $\beta=10^{-4}$ | $\beta=10^{-6}$ | $\beta=10^{-8}$
---|---|---|---|---
$1/64$ | 0.333 | 0.333 | 0.333 | 0.333
$1/128$ | 0.333 | 0.333 | 0.333 | 0.333
$1/256$ | 0.333 | 0.333 | 0.333 | 0.333

In Table 2, we report the LFA predicted two-grid convergence factor as a
function of the number of smoothing steps $\nu$ and $h$ with $\beta=10^{-2}$.
We see that $\rho_{h}=\mu^{\nu}_{\rm opt}$ for different $h$. For other
$\beta$, we see same behaviour, but omit the results.

Table 2: LFA predicted two-grid convergence factor, $\rho_{h}$, as a function of $\nu$ and $h$ with $\beta=10^{-2}$ for exact stiffness-based BSR. $h$ | $\nu=2$ | $\nu=3$ | $\nu=4$
---|---|---|---
$1/64$ | 0.111 | 0.037 | 0.012
$1/128$ | 0.111 | 0.037 | 0.012
$1/256$ | 0.111 | 0.037 | 0.012

###  4.2 Multigrid performance

In this subsection, we report some numerical results to confirm our
theoretical results. We consider the distributed control problem on the unit
square domain with homogeneous Dirichlet boundary conditions and $\hat{u}=0$
such that the discrete solution is $z_{h}=0$. The coarsest grid is a $4\times
4$ mesh. For the stiffness matrix $A_{fd}$, we consider rediscretization
operator on the coarse meshgrid. We measure the multigrid convergence factor
by the following form [15]

$\hat{\rho}_{h}=\sqrt[n]{\frac{||d_{h}^{(n)}||_{2}}{||d_{h}^{(0)}||_{2}}},$

where $d_{h}^{(n)}=b_{h}-\mathcal{L}_{h}z_{h}^{(n)}$, and $z_{h}^{(n)}$ is the
approximation to the solution of (7) at the $n^{\rm th}$ multigrid iteration.
We display the convergence factors obtained after 100 multigrid cycles for
different values of $\nu$. The initial guess $z_{h}^{(0)}$ is chosen randomly.
We use the MATLAB code of Rees [10] to generate the linear systems.

Table 3 shows the measured $W$-cycles multigrid convergence factors for exact
stiffness-based BSR with algorithmic parameters,
$\alpha_{B}=1,\omega_{B}=\frac{3}{4}$, for different values of $\beta$ and
grid size $h$. We see that our numerical results are independent of grid size
$h$, and $\hat{\rho}_{h}$ for $\nu>1$ matches well with the LFA two-grid
convergence factor $\rho_{h}$. Although for $\nu>1$ the measured convergence
factors are slightly better than the predictions, it is not surprising since
LFA two-grid prediction is the same as the actual convergence for problems
with periodic boundary conditions, and here we consider Dirichlet boundary
conditions.

Table 3:  $W$-cycles measured multigrid convergence factors, $\hat{\rho}_{h}$, for exact stiffness-based BSR with algorithmic parameters, $\alpha_{B}=1,\omega_{B}=\frac{3}{4}$, compared with LFA two-grid predictions, $\rho_{h}$, for exact BSR with $\alpha_{B}=1,\omega_{B}=\frac{3}{4}$. $\nu$ | 1 | 2 | 3 | 4
---|---|---|---|---
$\rho_{h}$ | 0.333 | 0.111 | 0.037 | 0.012
$\beta=10^{-2}$
$\hat{\rho}_{h=1/64}$ | 0.279 | 0.093 | 0.031 | 0.010
$\hat{\rho}_{h=1/128}$ | 0.275 | 0.091 | 0.030 | 0.010
$\hat{\rho}_{h=1/256}$ | 0.272 | 0.090 | 0.030 | 0.010
$\beta=10^{-4}$
$\hat{\rho}_{h=1/64}$ | 0.268 | 0.090 | 0.030 | 0.010
$\hat{\rho}_{h=1/128}$ | 0.264 | 0.089 | 0.029 | 0.010
$\hat{\rho}_{h=1/256}$ | 0.261 | 0.088 | 0.029 | 0.010
$\beta=10^{-6}$
$\hat{\rho}_{h=1/64}$ | 0.289 | 0.094 | 0.031 | 0.012
$\hat{\rho}_{h=1/128}$ | 0.276 | 0.093 | 0.030 | 0.010
$\hat{\rho}_{h=1/256}$ | 0.272 | 0.092 | 0.030 | 0.010
$\beta=10^{-8}$
$\hat{\rho}_{h=1/64}$ | 0.302 | 0.098 | 0.032 | 0.011
$\hat{\rho}_{h=1/128}$ | 0.299 | 0.098 | 0.032 | 0.012
$\hat{\rho}_{h=1/256}$ | 0.290 | 0.096 | 0.030 | 0.012

As mentioned before, using exact BSR, we need to solve the Schur complement
system (11) exactly. In practice, an inexact solve is preferred. First, we
tried two or three sweeps of weighted Jacobi iteration and Gauss-Seidel
iteration to approximate the solution of (11). However, experiments showed
that these approaches maintain the performance observed for exact BSR only for
$\nu=1$, and there is a degradation for $\nu>1$. Moreover, these approaches
are very sensitive to algorithmic parameters. Instead, we consider solving the
Schur complement system by applying a multigrid $V(2,2)$-cycle using weighted
Jacobi relaxation with weight $\omega_{J}$. In this test, we use
$\frac{\omega_{B}}{\alpha_{B}}=\frac{3}{4}$ because
$1-\frac{\omega_{B}}{\alpha_{B}}\frac{a}{\hat{a}}$ is an eigenvalue of the
relaxation error-propagation operator for both exact and inexact version. We
found that there are many choices of
$\alpha_{B},\omega_{B}=\frac{3}{4}\alpha_{B}$ and $\omega_{J}$ resulting in
good multigrid results. In our report, we consider
$\alpha_{B}=1.5,\omega_{B}=\frac{3}{4}\alpha_{B}$ and $\omega_{J}=0.8$. We
observe that using only 2 or 3 $V(2,2)$-cycles on the approximate Schur
complement achieves convergence factors essentially matching those in Table 3
for $\nu=1,2$. We display the inexact BSR results in Table 4, where 3
$V(2,2)$-cycles are used. We see there is a slight degradation of the measured
convergence factors for inexact BSR for $\nu=3,4$ compared with the results in
Table 3. However, for $\nu=1,2$, the measured convergence factors $\hat{\rho}$
agree very well with the results in Table 3, and are $h$-independent and
robust to $\beta$. In practice, we recommand $W(1,0)$-cycle and $W(1,1)$-cycle
for inexact BSR.

Table 4: $W$-cycles measured convergence factors, $\hat{\rho}_{h}$, for inexact stiffness-based BSR with inner 3 $V(2,2)$-cycles and algorithmic parameters, $\alpha_{B}=1.5,\omega_{B}=\frac{3}{4}\alpha_{B},\omega_{J}=0.8$. $\nu$ | 1 | 2 | 3 | 4
---|---|---|---|---
$\beta=10^{-2}$
$\hat{\rho}_{h=1/64}$ | 0.280 | 0.092 | 0.046 | 0.044
$\hat{\rho}_{h=1/128}$ | 0.275 | 0.092 | 0.046 | 0.044
$\hat{\rho}_{h=1/256}$ | 0.272 | 0.090 | 0.047 | 0.045
$\beta=10^{-4}$
$\hat{\rho}_{h=1/64}$ | 0.271 | 0.089 | 0.045 | 0.043
$\hat{\rho}_{h=1/128}$ | 0.267 | 0.088 | 0.046 | 0.044
$\hat{\rho}_{h=1/256}$ | 0.264 | 0.086 | 0.047 | 0.045
$\beta=10^{-6}$
$\hat{\rho}_{h=1/64}$ | 0.283 | 0.092 | 0.039 | 0.025
$\hat{\rho}_{h=1/128}$ | 0.280 | 0.091 | 0.044 | 0.041
$\hat{\rho}_{h=1/256}$ | 0.276 | 0.090 | 0.047 | 0.045
$\beta=10^{-8}$
$\hat{\rho}_{h=1/64}$ | 0.296 | 0.096 | 0.032 | 0.010
$\hat{\rho}_{h=1/128}$ | 0.293 | 0.095 | 0.034 | 0.025
$\hat{\rho}_{h=1/256}$ | 0.289 | 0.094 | 0.043 | 0.030

Next, we consider the model problem (1) subject to the state equation [11]

$\left\\{\begin{aligned} -\triangle u&=&f,\,\,\,\,{\rm in}\,\,\Omega,\\\
u&=&\hat{u}|_{\partial\Omega},\,\,\,{\rm at}\,\,\,\partial\Omega,\\\
\end{aligned}\right.\\\ $ (23)

where $\Omega=[0,1]^{2}$, and

$\hat{u}=\left\\{\begin{aligned}
&(2x-1)^{2}(2y-1)^{2}&\quad\text{if}\quad(x,y)\in[0,1/2]^{2},\\\
&0&\quad\text{otherwise}.\\\ \end{aligned}\right.\\\ $

For this problem, we consider $W$-cycles multigrid methods with inner 3
$V(2,2)$-cycles on the Schur complement system with
$\alpha_{B}=1.5,\omega_{B}=\frac{3}{4}\alpha_{B},\omega_{J}=0.8$. In Tables 5
and 6, we report the smallest number of iterations $n$ such that
$\frac{||d_{h}^{(n)}||_{2}}{||d_{h}^{(0)}||_{2}}<10^{-10}$ for this
distributed control model problem for different values of $\beta$, with
$\nu=1$ and $\nu=2$, respectively. Again, from Tables 5 and 6, we see the
robustness of our multigrid methods with respect to the regularization
parameters, $\beta$, and grid size $h$.

Table 5: Iteration counts for the control model problem with state equation (23) using $\nu=1$ and $W$-cycles multigrid for inexact stiffness-based BSR with inner 3 $V(2,2)$-cycles and algorithmic parameters, $\alpha_{B}=1.5,\omega_{B}=\frac{3}{4}\alpha_{B},\omega_{J}=0.8$.  $h$ | $\beta=10^{-2}$ | $\beta=10^{-4}$ | $\beta=10^{-6}$ | $\beta=10^{-8}$
---|---|---|---|---
$1/64$ | 17 | 17 | 17 | 19
$1/128$ | 17 | 17 | 17 | 18
$1/256$ | 18 | 18 | 18 | 18
$1/512$ | 19 | 19 | 19 | 19
Table 6: Iteration counts for the control model problem with state equation (23) using $\nu=2$ and $W$-cycles multigrid for inexact stiffness-based BSR with inner 3 $V(2,2)$-cycles and algorithmic parameters, $\alpha_{B}=1.5,\omega_{B}=\frac{3}{4}\alpha_{B},\omega_{J}=0.8$.  $h$ | $\beta=10^{-2}$ | $\beta=10^{-4}$ | $\beta=10^{-6}$ | $\beta=10^{-8}$
---|---|---|---|---
$1/64$ | 12 | 12 | 11 | 10
$1/128$ | 12 | 12 | 12 | 11
$1/256$ | 13 | 13 | 12 | 12
$1/512$ | 13 | 13 | 13 | 13

##  5 Conclusion

In this work, we present a novel multigrid method for the elliptic distributed
control problems, where stiffness-based Braess-Sarazin multigrid relaxation
scheme is proposed to the discrete saddle-point system of the control
problems. In this relaxation scheme, the stiffness matrix obtained from the
five-point finite difference method for the Laplacian is used to approximate
the mass matrix in the saddle-point system discretized by a finite element
method. We hire LFA to minimize LFA smoothing factor for Braess-Sarazin
relaxation, leading to a highly satisfactory convergence factor. For practical
use, we develop an inexact version of Braess-Sarazin relaxation, where a 3
$V(2,2)$-cycles multigrid with weighted Jacobi relaxation is applied to the
Schur complement system. This inexact version preserves the optimal smoothing
factor observed for exact Braess-Sarazin relaxation. Numerical results show
the robustness of our multigrid methods with respect to the regularization
parameters and grid size. In future, we will extend the relaxation scheme
considered here to more complex optimal problems, such as the Stokes control
problem.

## References

  * [1] A. Borzı, Multigrid methods for parabolic distributed optimal control problems, Journal of Computational and Applied Mathematics, 157 (2003), pp. 365–382.
  * [2] A. Borzı, High-order discretization and multigrid solution of elliptic nonlinear constrained optimal control problems, Journal of Computational and Applied Mathematics, 200 (2007), pp. 67–85.
  * [3] A. Borzì, K. Kunisch, and D. Y. Kwak, Accuracy and convergence properties of the finite difference multigrid solution of an optimal control optimality system, SIAM Journal on Control and Optimization, 41 (2002), pp. 1477–1497.
  * [4] A. Borzi and V. Schulz, Multigrid methods for PDE optimization, SIAM review, 51 (2009), pp. 361–395.
  * [5] D. Braess and R. Sarazin, An efficient smoother for the Stokes problem, Appl. Numer. Math., 23 (1997), pp. 3–19.
  * [6] M. Engel and M. Griebel, A multigrid method for constrained optimal control problems, Journal of Computational and Applied Mathematics, 235 (2011), pp. 4368–4388.
  * [7] C. Greif and Y. He, A closed-form multigrid smoothing factor for an additive Vanka-type smoother applied to the Poisson equation, arXiv preprint arXiv:2111.03190, (2021).
  * [8] Y. He, New mass-based multigrid relaxation schemes for the Stokes equations, arXiv preprint arXiv:2111.04922, (2021).
  * [9] J. KEPLER, Fourier analysis and local Fourier analysis for multigrid methods, (2009).
  * [10] T. Rees, Github/poisson-control.  https://github.com/tyronerees/poisson-control, 2010.
  * [11] T. Rees, H. S. Dollar, and A. J. Wathen, Optimal solvers for PDE-constrained optimization, SIAM Journal on Scientific Computing, 32 (2010), pp. 271–298.
  * [12] J. Schöberl, R. Simon, and W. Zulehner, A robust multigrid method for elliptic optimal control problems, SIAM Journal on Numerical Analysis, 49 (2011), pp. 1482–1503.
  * [13] R. Simon and W. Zulehner, On schwarz-type smoothers for saddle point problems with applications to PDE-constrained optimization problems, Numerische Mathematik, 111 (2009), pp. 445–468.
  * [14] S. Takacs and W. Zulehner, Convergence analysis of multigrid methods with collective point smoothers for optimal control problems, Computing and Visualization in Science, 14 (2011), pp. 131–141.
  * [15] U. Trottenberg, C. W. Oosterlee, and A. Schüller, Multigrid, Academic Press, Inc., San Diego, CA, 2001.  With contributions by A. Brandt, P. Oswald and K. Stüben.
  * [16] R. Wienands and W. Joppich, Practical Fourier analysis for multigrid methods, CRC press, 2004.


Copyright for this paper by its authors. Use permitted under Creative Commons
License Attribution 4.0 International (CC BY 4.0).

Manuscript submitted to Aequitas 2024: Workshop on Fairness and Bias in AI | co-located with ECAI 2024, Santiago de Compostela, Spain

[1]

[orcid=0000-0002-4653-7685<EMAIL_ADDRESS>] [1]

[orcid=0000-0003-2788-3854<EMAIL_ADDRESS>]

[1]Corresponding author.

# Measuring and Mitigating Bias for Tabular Datasets with Multiple Protected
Attributes

Manh Khoi Duong  Heinrich Heine University, Universitätsstraße 1, 40225
Düsseldorf, Germany  Stefan Conrad

(2024)

###### Abstract

Motivated by the recital (67) of the current corrigendum of the AI Act in the
European Union, we propose and present measures and mitigation strategies for
discrimination in tabular datasets. We specifically focus on datasets that
contain multiple protected attributes, such as nationality, age, and sex. This
makes measuring and mitigating bias more challenging, as many existing methods
are designed for a single protected attribute. This paper comes with a twofold
contribution: Firstly, new discrimination measures are introduced. These
measures are categorized in our framework along with existing ones, guiding
researchers and practitioners in choosing the right measure to assess the
fairness of the underlying dataset. Secondly, a novel application of an
existing bias mitigation method, FairDo, is presented. We show that this
strategy can mitigate any type of discrimination, including intersectional
discrimination, by transforming the dataset. By conducting experiments on
real-world datasets (Adult, Bank, Compas), we demonstrate that de-biasing
datasets with multiple protected attributes is achievable. Further, the
transformed fair datasets do not compromise any of the tested machine learning
models’ performances significantly when trained on these datasets compared to
the original datasets. Discrimination was reduced by up to 83% in our
experimentation. For most experiments, the disparity between protected groups
was reduced by at least 7% and 27% on average. Generally, the findings show
that the mitigation strategy used is effective, and this study contributes to
the ongoing discussion on the implementation of the European Union’s AI Act.

###### keywords:

Machine Learning Bias Mitigation Intersectional Discrimination Fairness AI Act

##  1 Introduction

Discrimination in artificial intelligence (AI) applications is a growing
concern since the adoption of the _AI Act_ by the European Parliament on March
13, 2024 [1]. It still remains a significant challenge across numerous domains
[2, 3, 4, 5]. To prevent biased outcomes, _pre-processing_ methods are often
used to mitigate biases in datasets before training machine learning models
[6, 7, 8, 9]. These methods are even mentioned in Recital (67) at the current
corrigendum of the _AI Act_ [1]:

> _“[…] The data sets should also have the appropriate statistical properties,
> including as regards the persons or groups of persons in relation to whom
> the high-risk AI system is intended to be used, with specific attention to
> the mitigation of possible biases in the data sets […]”_

Since datasets often consist of multiple protected attributes, pre-processing
methods should be able to handle these cases. However, only a few works have
addressed it [7, 10, 11, 12] and de-biasing datasets with multiple protected
attributes is still an ongoing research topic.

Our paper mainly focuses on how to measure and mitigate discrimination in
datasets where multiple protected attributes are present. In our first
contribution, we provide a comprehensive categorization of discrimination
measuring methods. Besides introducing new measures for some of these cases,
we also categorize existing measures from the literature. Some of the listed
measures specifically address _intersectional discrimination_ and _non-binary
groups_. The second contribution deals with bias mitigation. For this, we use
an existing pre-processing framework, FairDo [9], that is _fairness-agnostic_.
The fairness-agnostic property makes it possible to define any discrimination
measure that should be minimized. By implementing the introduced measures, we
can therefore mitigate biases for multiple protected attributes. Another
advantage of FairDo is that it preserves data integrity and does not modify
the features of individuals during the optimization process, unlike other
methods [13, 3, 7].

We evaluated our methodology on popular tabular datasets with fairness
concerns, such as Adult [14], Bank [15], and Compas [16]. These datasets
contain multiple protected attributes and we mitigated the inherent biases. We
used different discrimination measures to evaluate the effectiveness of the
bias mitigation process. Because a successful mitigation process does not
guarantee that the outcomes of machine learning models are fair, we trained
machine learning models on the transformed datasets and evaluated their
predictions regarding fairness and performance.

The results of the bias mitigation process as well as the performance of the
machine learning models are promising. They indicate that achieving fairness
in datasets with multiple protected attributes is perhaps possible, and FairDo
is a proper framework for this task. Overall, our work contributes technical
solutions for stakeholders to enhance the fairness of datasets and machine
learning models, aiming for compliance with the _AI Act_ [1].

Figure 1: Stick figures can be differentiated by their color and shape. In
intersectional discrimination, attributes are intersected, which leads to new
subgroups. In non-intersectional, each attribute is treated separately, i.e.,
colors and shapes are not intersecting in this case.

##  2 Preliminaries

To handle multiple protected attributes, we define
$Z=\\{Z_{1},\ldots,Z_{p}\\}$ as a set of protected attributes. It can
represent the set of sociodemographic features such as age, gender, and
ethnicity. These factors may make individuals vulnerable to discrimination.
Each protected attribute $Z_{k}\in Z$ is formally a _discrete random variable_
that can take on values from a set $g_{k}$. In this context, we refer $g_{k}$
to groups that describe distinct social categories of a protected attribute.
For example, let $Z_{k}$ represent gender; then $g_{k}$ is a set containing
the genders male, female, and non-binary. To avoid limitations to a particular
group fairness notion, we introduce a generalized notation based on the works
of Žliobaitė [2], Duong and Conrad [9] in the following.

######  Definition 2.1 (Treatment).

Let $E_{1},E_{2}$ be events and $Z_{k}$ be a random variable that can take on
values from $g_{k}$, then we call the conditional probability

$P(E_{1}\mid E_{2},Z_{k}=i)$

_treatment_ , where $i\in g_{k}$. $E_{1}$ describes some favorable outcome,
such as getting accepted for a job, while $E_{2}$ often represents some
additional information about the individual, such as their qualifications.

######  Definition 2.2 (Fairness Criteria).

With the definition of treatment, we can define fairness criteria that demand
equal treatments for different groups. Let $P(E_{1}\mid E_{2},Z_{k}=i)$ and
$P(E_{1}\mid E_{2},Z_{k}=j)$ be two treatments, then we call the following
equation:

$P(E_{1}\mid E_{2},Z_{k}=i)=P(E_{1}\mid E_{2},Z_{k}=j)$

a _fairness criterion_ , for all $i,j\in g_{k}$.

Definition 2.2 allows us to define various group fairness criteria, including
_statistical parity_ [17], _predictive parity_ [3], _equality of opportunity_
[18], etc. They all demand some sort of equal outcome for different groups and
can be defined by configuring the events $E_{1},E_{2}$. For instance,
statistical parity [17] requires that two different groups have an equal
probability of receiving a favorable outcome $(Y=1)$.

######  Example 2.1 (Statistical Parity [17]).

To define statistical parity for the attribute $Z_{k}$ using our notation, we
set $E_{1}\coloneqq(Y=1)$ and $E_{2}\coloneqq\text{True}$. By setting $E_{2}$
to a tautology, we compare the probabilities of the event $Y=1$ across
different groups without conditioning on any additional event:

$\displaystyle P(Y=1\mid\text{True},Z_{k}=i)$
$\displaystyle=P(Y=1\mid\text{True},Z_{k}=j)$ $\displaystyle\iff P(Y=1\mid
Z_{k}=i)$ $\displaystyle=P(Y=1\mid Z_{k}=j),$

where $i,j\in g$ represent different groups.

In real-world applications, achieving equal probabilities for certain outcomes
is not always possible. Due to variations in sample sizes in the groups, it is
common to yield unequal treatments, even when they are similar. Thus, existing
literature [2] uses the absolute difference to quantify the strength of
discrimination.

######  Definition 2.3 (Disparity).

Let $P(E_{1}\mid E_{2},Z_{k}=i)$ and $P(E_{1}\mid E_{2},Z_{k}=j)$ be two
treatments, then we refer to

$\delta_{Z_{k}}(i,j,E_{1},E_{2})=|P(E_{1}\mid E_{2},Z_{k}=i)-P(E_{1}\mid
E_{2},Z_{k}=j)|$

as the _disparity_ , for all $i,j\in g_{k}$. Trivially, $\delta_{Z_{k}}$ is
commutative regarding $i,j$. In practice, it prevents reverse discrimination
due to the absolute value.

######  Definition 2.4 (Discrimination).

We use $\psi\colon\mathbb{D}\rightarrow\mathbb{R}$ to denote some
discrimination measure that quantifies the discrimination inherent in any
dataset $\mathcal{D}\in\mathbb{D}$. A dataset $\mathcal{D}$ consists of
features, protected attributes, and labels for each individual. The explicit
form of $\psi$ depends on cases introduced in Section 3.

##  3 Measuring Discrimination for Multiple Attributes

We found that numerous scenarios arise when dealing with multiple protected
attributes. We categorize these scenarios based on the number of groups,
denoted as $|g|$, and number of protected attributes, denoted as $|Z|$. By
going through all cases, we present possible approaches from the literature as
well as our own suggestions to measure discrimination.

###  3.1 Single Protected Attribute $(|Z|=1)$

In case of having only one protected attribute, i.e., $|Z|=|\\{Z_{1}\\}|=1$,
we distinguish between cases by the number of available groups $|g|$ in the
dataset. We categorize the cases by $|g|=0,1,2$, and $|g|>2$.

####  3.1.1 No Groups ($|g|=0$)

When having no groups, the measurement of discrimination is impossible if no
assumptions are being made. Discrimination can be assessed through proxy
variables [19]; however, this approach can be imprecise and may introduce new
biases. This case is equivalent to having no protected attribute, i.e.,
$|Z|=0$.

####  3.1.2 Single Group ($|g|=1$)

Similarly to the case of having no groups, discrimination cannot be measured
when having only one group. For this, we propose practices where prior
information can be incorporated:

  1. 1.

_No discrimination_ : As no difference towards any other group can be
measured, returning a discrimination score of 0 is one viable option.

$\psi(\mathcal{D})=0.$ (1)

  2. 2.

_Difference to optimal treatment_ : Another way is to return the absolute
difference of the group’s outcome to the optimal treatment. For example, group
$i$ has an 80% chance of getting the favorable treatment. Ideally, having a
100% chance of receiving the favorable treatment would represent the optimal
scenario. Therefore, the discrimination score is 20% in this case. It is given
by:

$\psi(\mathcal{D})=|P(E_{1}\mid E_{2},Z_{1}=i)-1|.$ (2)

  3. 3.

_Difference to expected treatment_ : We can use the expected treatment as a
reference point. For example, we know that a company has a 50% acceptance rate
for job applications. Now a machine learning classifier is trained to predict
whether an applicant will be accepted and the model’s predictions result in a
60% acceptance rate for group $i$. Hence, the model is positively biased
towards group $i$ by 10%. This can be formulated as:

$\psi(\mathcal{D})=|P(E_{1}\mid E_{2},Z_{1}=i)-p_{\text{expect.}}|,$ (3)

where $p_{\text{expect.}}$ is the expected treatment. It can describe the
average treatment across all groups [20] or some other prior information that
is not included in the dataset.

####  3.1.3 Binary Groups ($|g|=2$)

Without using any prior information, we can calculate the discrimination score
by taking the absolute difference between the treatments of the two groups, as
advised by Žliobaitė [2]. The discrimination measure $\psi$ is then simply
given by the disparity as mentioned in Definition (2.3).

####  3.1.4 Non-binary Groups ($|g|>2$)

While the case for binary attributes is straightforward, it becomes non-
trivial for non-binary attributes that arise naturally in real-world data. We
can fall back to $|g|=2$ by calculating the absolute difference between every
distinct group $i,j\in g$. Because the discrimination between $i$ and $j$ is
the same as between $j$ and $i$, only $\binom{|g|}{2}$ pairs need to be
compared and we use an aggregation function $\operatorname{agg}^{(1)}$ to
report the differences [2]. Lum et al. [21] refer to measures that aggregate
or summarize discrimination scores as _meta-metrics_. The aggregate can be the
_sum_ or _maximum_ function, depending on the use case. The result for a
single protected attribute $Z_{k}$ with two or more groups can be computed as
follows:

$\psi(\mathcal{D})=\underset{i,j\in
g_{k},i<j}{\operatorname{agg}^{(1)}}\delta_{Z_{k}}(i,j,E_{1},E_{2}),$ (4)

where $\delta_{Z_{k}}$ is the disparity as defined in Definition (2.3) and
$i<j$ ensures that each pair is considered only once (assuming label encoded
groups). According to Žliobaitė [2] and her personal discussions with legal
experts, she advocates using the maximum function, i.e.,

$\displaystyle\psi(\mathcal{D})$ $\displaystyle=\underset{i,j\in
g_{k},i<j}{\max}\delta_{Z_{k}}(i,j,E_{1},E_{2})$ (5)
$\displaystyle=\underset{i\in g_{k}}{\max}\,P(E_{1}\mid
E_{2},Z_{k}=i)-\underset{j\in g_{k}}{\min}\,P(E_{1}\mid E_{2},Z_{k}=j).$ (6)

Equation (5) describes the maximum discrimination obtainable between two
groups. An alternative and equivalent formulation is given in Equation (6)
[7]. The latter is computationally more efficient as it requires
$\mathcal{O}(2|g|)$ operations compared to $\mathcal{O}(|g|^{2})$ operations
for the former.

A more general approach to measure discrimination is to calculate some form of
_correlation coefficient_ between the protected attribute and the outcome. The
correlation coefficient can be calculated using Pearson’s correlation [22],
Spearman or Kendall’s rank correlation [23, 24]. The discrimination measure
can then be defined as the absolute value of the correlation coefficient:

$\psi(\mathcal{D})=|\text{Corr}(E_{1},Z_{k})|.$ (7)

This approach can be applied to any number of groups. Fairlearn provides a
pre-processing method that removes the correlation between the protected
attribute and the outcome by transforming the data [7]. However, the given
approach violates data integrity constraints as categorical attributes are
transformed into continuous values. Further, zero correlation does not imply
independency between two variables.

###  3.2 Multiple Protected Attributes $(|Z|>1)$

There are several ways to measure discrimination for multiple protected
attributes ($|Z|>1$). Based on the works of Kearns et al. [20] and Kang et al.
[12], we categorize them into two approaches: _intersectional_ and _non-
intersectional_ (see Figure 1). Intersectional approaches consider the
intersection of groups/identities. The overlapping of such identities forms
_subgroups_ [20]. Non-intersectional approaches treat each protected attribute
independently.

####  3.2.1 Intersectional Discrimination

The central idea of intersectionality is that individuals experience
overlapping forms of oppression or privilege based on the combination of
multiple social categories they belong to. In the following, we will introduce
necessary definitions to formulate intersectional discrimination, which is
based on the work of Kearns et al. [20].

######  Definition 3.1 (Subgroup [20]).

Let $Z=\\{Z_{1},\ldots,Z_{p}\\}$ be a set of discrete random variables
representing protected attributes that can take on values from corresponding
groups $g_{1},\ldots,g_{p}$. A _subgroup_ $i$ is defined as
$i=(i_{1},\ldots,i_{p})\in g_{1}\times\ldots\times g_{p}$. In other words, a
_subgroup_ encompasses multiple groups from different protected attributes.

######  Definition 3.2 (Subgroup Treatment).

Let $i$ be a subgroup as defined in Definition 3.1 and let
$Z=\\{Z_{1},\ldots,Z_{p}\\}$ be a set of discrete random variables. _Subgroup
treatment_ is then defined as:

$P(E_{1}\mid E_{2},Z_{1}=i_{1},\ldots,Z_{p}=i_{p}).$

######  Definition 3.3 (Subgroup Disparity).

Let $Z=\\{Z_{1},\ldots,Z_{p}\\}$ be a set of discrete random variables. Let
$i,j\in g_{1}\times\ldots\times g_{p}$ be two subgroups with
$i=(i_{1},\ldots,i_{p})$ and $j=(j_{1},\ldots,j_{p})$. The disparity between
two subgroups is denoted as $\hat{\delta}_{Z}$ and is given by:

$\hat{\delta}_{Z}(i,j,E_{1},E_{2})=|P(E_{1}\mid
E_{2},Z_{1}=i_{1},\ldots,Z_{p}=i_{p})-P(E_{1}\mid
E_{2},Z_{1}=j_{1},\ldots,Z_{p}=j_{p})|.$

Similarly to Equation (4), we can calculate the discrimination score for
multiple protected attributes by aggregating disparities across all subgroups.
A subgroup can be treated like a normal group. According to Definition 3.1,
there are theoretically at least $2^{p}$ subgroups, where $p$ is the number of
protected attributes. However, not all subgroups may be available in the
dataset. For unavailable subgroups, the disparity cannot be calculated as the
corresponding treatment is undefined.

Let denote the set of available subgroups as $G_{\text{avail}}\subseteq
g_{1}\times\ldots\times g_{k}$. To finally capture the discrepancies across
all available subgroup pairs, an aggregation function
$\operatorname{agg}^{(1)}$ is applied to the subgroup disparities
$\delta_{Z}$:

$\displaystyle\psi_{\text{intersect.}}(\mathcal{D})$
$\displaystyle=\underset{i,j\in
G_{\text{avail}}}{\operatorname{agg}^{(1)}}\hat{\delta}_{Z}(i,j,E_{1},E_{2}).$
(8)

Equation (8) represents the aggregated discrimination between all available
subgroups in the dataset. When using the maximum function as the aggregator,
the calculations are equivalent to Equation (5) and Equation (6). The only
difference is that the conditionals are now subgroups instead of groups:

$\displaystyle\psi_{\text{intersect.}}(\mathcal{D})$
$\displaystyle=\underset{i,j\in
G_{\text{avail}}}{\max}\hat{\delta}_{Z_{k}}(i,j,E_{1},E_{2})$ (9)
$\displaystyle=\underset{i\in G_{\text{avail}}}{\max}\,P(E_{1}\mid
E_{2},Z_{1}=i_{1},\ldots,Z_{p}=i_{p})-\underset{j\in
G_{\text{avail}}}{\min}\,P(E_{1}\mid E_{2},Z_{1}=j_{1},\ldots,Z_{p}=j_{p}).$
(10)

Kang et al. [12] also dealt with intersectional discrimination in their work
by introducing a multivariate random variable where each dimension represents
a protected attribute. Their fairness objective is to minimize the mutual
information between the outcome and the multivariate random variable. By
minimizing the mutual information, the outcome is independent of the protected
attributes, which is a desirable property for fairness [13, 25]. However, this
approach relies on expensive techniques to approximate the mutual information.
Using our notation, their formulation can be written as [12]:

$\psi_{\text{MI}}(\mathcal{D})=\text{MI}(E_{1},Z),$ (11)

where MI denotes the mutual information and $Z$ is the multivariate random
variable representing the protected attributes. We can trivially conclude that
zero mutual information implies zero intersectional discrimination [12].

Table 1: Example dataset of individuals receiving a favorable ($Y=1$) or unfavorable ($Y=0$) outcome. The dataset shows four individuals with their respective age group and sex. Individual | Age | Sex | Outcome ($Y$)
---|---|---|---
1 | Old | Male | 1
2 | Old | Female | 0
3 | Young | Male | 0
4 | Young | Female | 1

####  3.2.2 Non-intersectional Discrimination

The problem with intersectional approaches is that they can lead to an upwards
bias in the discrimination score when using meta-metrics [21]. In our case,
the number of subgroups grows exponentially with the number of protected
attributes. This leads to many subgroups where the number of samples in each
subgroup is small, resulting in larger noise in the treatment estimates [21].

Hence, we present a non-intersectional approach, which is more suitable when
dealing with a large number of subgroups or when intersectional discrimination
is not deemed important. The non-intersectional approach treats each protected
attribute independently and aggregates the discrimination scores across all
protected attributes. For this, a second aggregate function with
$\operatorname{agg}^{(2)}$ is introduced, yielding the following equation:

$\psi_{\text{multi.}}(\mathcal{D})=\underset{Z_{k}\in
Z}{\operatorname{agg}^{(2)}}\left\\{\underset{i,j\in
g_{k},i<j}{\operatorname{agg}^{(1)}}\delta_{Z_{k}}(i,j,E_{1},E_{2})\right\\}.$
(12)

The first-level aggregator $\operatorname{agg}^{(1)}$ aggregates disparities
within a protected attribute, considering unique pairs of groups $i$ and $j$
in each attribute. The second-level aggregator $\operatorname{agg}^{(2)}$ then
combines the results across all protected attributes. By applying both
operators, we obtain a discrimination measure that captures disparities
between groups across multiple attributes.

####  3.2.3 Example

Let us consider a dataset with two protected attributes, age and sex (see
Table 1). The set of protected attributes is
$Z=\\{Z_{1},Z_{2}\\}=\\{\text{Age},\text{Sex}\\}$ and the set of available
subgroups in the dataset is $G_{\text{avail}}=\\{\text{Old,
Young}\\}\times\\{\text{Male, Female}\\}$. We measure discrimination using
_statistical disparity_. For simplicity, all aggregation functions are set to
the maximum function. The _intersectional approach_ yields the following
discrimination score:

$\displaystyle\psi_{\text{intersect.}}(\mathcal{D})$
$\displaystyle=\underset{i,j\in
G_{\text{avail}}}{\max}\hat{\delta}_{Z}(i,j,(Y=1),\text{True})$ (13)
$\displaystyle=\underset{i,j\in
G_{\text{avail}}}{\max}\hat{\delta}_{\text{\\{Age,
Sex\\}}}(i,j,(Y=1),\text{True})$ (14) $\displaystyle=\underset{i\in
G_{\text{avail}}}{\max}P(Y=1\mid Z_{1}=i_{1},Z_{2}=i_{2})-\underset{j\in
G_{\text{avail}}}{\min}P(Y=1\mid Z_{1}=j_{1},Z_{2}=j_{2})$ (15)
$\displaystyle=|P(Y=1\mid\text{Age}=\text{Old},\text{Sex}=\text{Male})-P(Y=1\mid\text{Age}=\text{Young},\text{Sex}=\text{Male})|=1.$
(16)

We show in the following that the discrimination for the _non-intersectional
approach_ is different as it treats each protected attribute independently:

$\displaystyle\psi_{\text{multi.}}(\mathcal{D})$
$\displaystyle=\underset{Z_{k}\in Z}{\max}\left\\{\underset{i,j\in
g_{k},i<j}{\max}\delta_{Z_{k}}(i,j,(Y=1),\text{True})\right\\}$ (17)
$\displaystyle=\max\left\\{\delta_{\text{Age}}(\text{Old},\text{Young},(Y=1),\text{True}),\delta_{\text{Sex}}(\text{Male},\text{Female},(Y=1),\text{True})\right\\}$
(18) $\displaystyle=\max\\{|0.5-0.5|,|0.5-0.5|\\}=\max\\{0,0\\}=0.$ (19)

The non-intersectional approach yields a discrimination score of 0 because the
disparities for both protected attributes are 0. This is quite different from
the intersectional approach, which reports a discrimination score of 1. As
seen, the choice of the approach can lead to different results.

##  4 Experiments

Our experimentation follows a pipeline consisting of data pre-processing, bias
mitigation, model training, and evaluation. To mitigate bias in tabular
datasets with multiple protected attributes, we used a sampling method, FairDo
[9], that constructs fair datasets by selectively sampling data points. The
method is very flexible and only requires the user to define the
discrimination measure that should be minimized. In our case, we are
interested in yielding a fair dataset that minimizes the bias for multiple
protected attributes. The experiments revolve around two research questions:

  * •

RQ1 Is it possible to yield a fair dataset with FairDo, where bias for
multiple protected attributes is reduced?

  * •

RQ2 Are machine learning models trained on fair datasets more fair in their
predictions than those trained on original datasets?

###  4.1 Experimental Setup

##### Datasets and Pre-processing

Table 2: Overview of Datasets Dataset | Samples | Feats. |  Label |  Protected Attributes |  Description
---|---|---|---|---|---
Adult [14] | 32 561 | 21 |  Income |  Race: White, Black, Asian-Pacific-Islander, American-Indian-Eskimo, Other
Sex: Male, Female |  Indicates individuals earning over $50,000 annually
Bank [15] | 41 188 | 50 |  Term deposit
subscription |  Job: Admin, Blue-Collar, Technician, Services, Management, Retired, Entrepreneur, Self-Employed, Housemaid, Unemployed, Student, Unknown
Marital Status: Divorced, Married, Single, Unknown |  Shows whether the client has subscribed to a term deposit.
Compas [16] | 7 214 | 13 |  2-year
recidivism |  Race: African-American, Caucasian, Hispanic, Other, Asian, Native American
Sex: Male, Female
Age Category: <25, 25-45, >45 |  Displays individuals that were rearrested for a new crime within 2 years after initial arrest.

The tabular datasets used in our experiments are: Adult [14], Bank [15], and
Compas [16]. They are known for their use in fairness research and contain
multiple protected attributes. We pre-processed the datasets by applying one-
hot encoding to categorical variables and label encoding to protected
attributes. Table 2 shows important characteristics of the datasets after pre-
processing.

Each dataset was divided into training and testing sets using an 80/20 split,
respectively. We ensured that the split was stratified (if possible) based on
protected attributes to maintain representativeness across different groups in
both sets.

##### Bias Mitigation

Applying the bias mitigation method FairDo [9] to the datasets can be regarded
as a pre-processing step, too. This is because the method simply returns a
dataset that is fair with respect to the given discrimination measure. FairDo
[9] offers a variety of methods to mitigate bias and we chose the
_undersampling_ method that removes samples. The fairness objective is stated
as [9]:

$\min_{\mathcal{D}_{\text{fair}}\subseteq\mathcal{D}}\quad\psi(\mathcal{D}_{\text{fair}}),$
(20)

where $\mathcal{D}$ is the training set of Adult, Bank, or Compas, and $\psi$
is the fairness objective function. We considered $\psi_{\text{intersect.}}$
and $\psi_{\text{multi.}}$ as our objectives. The first measures
intersectional discrimination, while the latter measures discrimination for
multiple protected attributes where the intersection of groups is not
considered. Bias mitigation is only applied to the training set and the
testing set remains unchanged. FairDo internally uses genetic algorithms for
optimization. We used the same settings and operators as provided in the
package and only adjusted the population size (200) and the number of
generations (400).

##### Model Training

We utilized the scikit-learn library [26] to train various machine learning
classifiers, namely _Logistic Regression_ (LR), _Support Vector Machines_
(SVM), _Random Forest_ (RF), and _Artificial Neural Networks_ (ANNs). These
classifiers were trained on both the original and fair datasets. Classifiers
trained on the original datasets serve as a baseline for comparison. We used
the default hyperparameters given by scikit-learn for each classifier.

##### Evaluation Metrics

We evaluated the models’ predictions on fairness and performance using the
test set. For fairness, we assessed $\psi_{\text{intersect.}}$ and
$\psi_{\text{multi.}}$. For the classifiers’ performances, we report the _area
under the receiver operating characteristic curve_ (AUROC) [27], where higher
values indicate better performance. Because removing data points can
compromise the overall quality of the data, we also report the number of
subgroups before and after bias mitigation, as well as the number of samples.

##### Trials

For each dataset, the bias mitigation process was repeated 10 times. The
results were averaged over the trials to obtain a more robust evaluation.

###  4.2 Results

##### Fair Dataset Generation

Table 3: Average discrimination and number of subgroups before and after pre-processing the training sets with FairDo. |  | Disc. Before | Disc. After | Subgroups Before | Subgroups After
---|---|---|---|---|---
Dataset | Metrics |  |  |  |
Adult | $\psi_{\text{multi.}}$ | $20\%$ | $13\%$ | 10 | 10
$\psi_{\text{intersect.}}$ | $31\%$ | $16\%$ | 10 | 10
Bank | $\psi_{\text{multi.}}$ | $24\%$ | $5\%$ | 48 | 48
$\psi_{\text{intersect.}}$ | $33\%$ | $15\%$ | 48 | 46.2
Compas | $\psi_{\text{multi.}}$ | $30\%$ | $5\%$ | 34 | 34
$\psi_{\text{intersect.}}$ | $100\%$ | $17\%$ | 34 | 28.8

Table 3 shows the average discrimination before and after mitigating bias in
the training sets. On all datasets, discrimination was reduced after applying
FairDo. It can be observed that calculating the discrimination for multiple
protected attributes while considering intersectionality leads to a higher
discrimination value. Without considering group intersections, discrimination
was reduced by 7%, 19%, and 25% for Adult, Bank, and Compas, respectively.
When considering intersectionality, the discrimination was reduced by 15%,
18%, and 83%. Hence, discrimination was reduced by 28% on average across all
datasets.

We also report the number of subgroups before and after bias mitigation to
assess the impact of the undersampling method on the dataset. The removal of
subgroups can only be observed in the intersectional setting. In the Compas
dataset 5.2 out of 34 subgroups were removed on average, indicating the
largest amount of subgroups removed across all datasets. While the Bank
dataset consists of 48 subgroups, only 1.8 subgroups were removed on average.
Because the Compas dataset’s initial intersectional discrimination score is
100%, removing multiple subgroups seems inevitable to reduce discrimination.

##### Model Performance and Fairness

(a)

(a) Adult (LR)

(b) Adult (SVM)

(c) Adult (RF)

(d) Adult (ANN)

(e) Bank (LR)

(f) Bank (SVM)

(g) Bank (RF)

(h) Bank (ANN)

(i) Compas (LR)

(j) Compas (SVM)

(k) Compas (RF)

(l) Compas (ANN)

Figure 2: Results on the test set. The x-axis represents the discrimination
values (legend indicates used measure) and the y-axis represents the
classifiers’ performances. We compare the pre-processed (fair) data with the
original data. The points represent the average and the error bars show the
standard deviation of the AUROC and discrimination values over 10 trials.

Figure 2 shows the results of the classifiers’ performances on the test set.
The classifiers’ performances are displayed on the y-axis, while the
discrimination values are shown on the x-axis. We show AUROC to evaluate the
classifiers’ performances, where higher values indicate better performances.
We note that the axes are not shared across the subfigures for better
visibility.

Classifiers trained on fair datasets did not suffer a significant decline of
performance than those trained on original datasets. In all cases, only a
slight decrease of 1%-3% in performance can be observed. This indicates that
the bias mitigation process does not compromise the dataset’s fidelity and
therefore, the classifiers’ performances. Regarding discrimination, a great
reduction can be seen. The x-axis scales are much larger than the y-axis
scales, suggesting that changes in discrimination are larger than changes in
performance. For example, the RF classifier trained on the Bank dataset
(Figure 2g) shows a decrease in intersectional discrimination from 38% to 15%,
while the performance only decreases by 2%. The results suggest that FairDo
can be reliably used to mitigate bias in tabular datasets for various measures
that consider multiple protected attributes. Still, we advise users to
carefully perform similar analyses when applying the method to their datasets.

##  5 Discussion

The results of our experiments show that the presented measures detect
discrimination in datasets with multiple protected attributes differently.
When using the intersectional discrimination measure, more groups are
identified and compared to each other, which leads to a higher reported
discrimination score. While subgroups are not ignored by this measure,
exponentially more groups are considered, and measuring higher discrimination
scores by random chance becomes more likely [20, 21]. In contrast, treating
each protected attribute separately prevents this issue but may lead to an
underestimation of discrimination. The choice of measure is up to the
stakeholders and depends on the context of the dataset. It also depends on the
regulations for the application of the AI system. We generally recommend using
the intersectional discrimination measure if the interplay of protected
attributes is of interest and the number of groups is low. If a higher number
of groups is present, the number of individuals in each subgroup should be
high enough to prevent high discrimination scores that are due to noise.
Otherwise, treating each protected attribute separately is more suitable.

Using the mitigation strategy FairDo [9], the resulting datasets in the
experiments have improved statistical properties regarding fairness. Whether
intersectionality was considered or not, reducing discrimination in datasets
was possible, which led to training fairer machine learning models without
compromising their performances substantially. At the current state, the AI
Act [1] does not explicitly mention _intersectional discrimination_ nor how to
deal with multiple protected attributes generally. While recital (67) states
that datasets _“should […] have the appropriate statistical properties”_ , it
does not specify what these properties are. Hence, our work serves as an
initial guideline for what these properties could be and how to achieve them
in practice.

##  6 Conclusion

Datasets often come with multiple protected attributes, which makes measuring
and mitigating discrimination more challenging. Most existing studies only
deal with a single protected attribute, and works that consider multiple
protected attributes often focus on intersectionality. In opposition to this,
we have proposed a new measure that treats each protected attribute
separately. This is more suitable when the number of subgroups is too high and
the number of individuals in each subgroup is low. We used both measures in
our experiments and used them as objectives for the bias mitigation method
FairDo. On average, the disparity between protected groups was reduced by 28%
in the bias-mitigated datasets compared to the original datasets. Machine
learning models trained on these datasets also showed a reduction in
discrimination while maintaining their performance.

## References

  * European Commission [2024] European Commission, Artificial Intelligence Act, Corrigendum, 19 April 2024, Available online: https://www.europarl.europa.eu/doceo/document/TA-9-2024-0138-FNL-COR01_EN.pdf, 2024\. Accessed: 17 May 2024.
  * Žliobaitė [2017] I. Žliobaitė,  Measuring discrimination in algorithmic decision making,  Data Mining and Knowledge Discovery 31 (2017) 1060–1089.
  * Zafar et al. [2017] M. B. Zafar, I. Valera, M. Gomez Rodriguez, 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, International World Wide Web Conferences Steering Committee, 2017. doi:10.1145/3038912.3052660.
  * Corbett-Davies et al. [2017] S. Corbett-Davies, E. Pierson, A. Feller, S. Goel, 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, 2017, pp. 797–806.
  * Barocas et al. [2019] S. Barocas, M. Hardt, A. Narayanan, Fairness and Machine Learning, fairmlbook.org, 2019.
  * Feldman et al. [2015] M. Feldman, S. A. Friedler, J. Moeller, C. Scheidegger, S. Venkatasubramanian,  Certifying and removing disparate impact,  in: proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, 2015, pp. 259–268.
  * Bird et al. [2020] S. Bird, M. Dudík, R. Edgar, B. Horn, R. Lutz, V. Milan, M. Sameki, H. Wallach, K. Walker, Fairlearn: A toolkit for assessing and improving fairness in AI, Technical Report MSR-TR-2020-32, Microsoft, 2020. URL: https://www.microsoft.com/en-us/research/publication/fairlearn-a-toolkit-for-assessing-and-improving-fairness-in-ai/.
  * Agarwal et al. [2018] A. Agarwal, A. Beygelzimer, M. Dudík, J. Langford, H. Wallach,  A reductions approach to fair classification,  in: International Conference on Machine Learning, PMLR, 2018, pp. 60–69.
  * Duong and Conrad [2024] M. K. Duong, S. Conrad,  Towards fairness and privacy: A novel data pre-processing optimization framework for non-binary protected attributes,  in: D. Benavides-Prado, S. Erfani, P. Fournier-Viger, Y. L. Boo, Y. S. Koh (Eds.), Data Science and Machine Learning, Springer Nature Singapore, Singapore, 2024, pp. 105–120.
  * Foulds et al. [2020] J. R. Foulds, R. Islam, K. N. Keya, S. Pan, Bayesian Modeling of Intersectional Fairness: The Variance of Bias, 2020, pp. 424–432. doi:10.1137/1.9781611976236.48.
  * Celis et al. [2021] L. E. Celis, L. Huang, V. Keswani, N. K. Vishnoi,  Fair classification with noisy protected attributes: A framework with provable guarantees,  in: M. Meila, T. Zhang (Eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, PMLR, 2021, pp. 1349–1361.
  * Kang et al. [2021] J. Kang, T. Xie, X. Wu, R. Maciejewski, H. Tong,  Infofair: Information-theoretic intersectional fairness,  2022 IEEE International Conference on Big Data (Big Data) (2021) 1455–1464.
  * Zemel et al. [2013] R. Zemel, Y. Wu, K. Swersky, T. Pitassi, C. Dwork,  Learning fair representations,  in: International conference on machine learning, PMLR, 2013, pp. 325–333.
  * Kohavi [1996] R. Kohavi,  Scaling up the accuracy of naive-bayes classifiers: A decision-tree hybrid,  KDD’96, AAAI Press, 1996, p. 202–207.
  * Moro et al. [2014] S. Moro, P. Cortez, P. Rita,  A data-driven approach to predict the success of bank telemarketing,  Decision Support Systems 62 (2014) 22–31.
  * Larson et al. [2016] J. Larson, J. Angwin, S. Mattu, L. Kirchner, Machine bias, 2016. URL: https://www.propublica.org/article/machine-bias-risk-assessments-in-criminal-sentencing.
  * Calders et al. [2009] T. Calders, F. Kamiran, M. Pechenizkiy,  Building classifiers with independency constraints,  in: 2009 IEEE International Conference on Data Mining Workshops, 2009, pp. 13–18. doi:10.1109/ICDMW.2009.83.
  * Hardt et al. [2016] M. Hardt, E. Price, N. Srebro,  Equality of opportunity in supervised learning,  Advances in neural information processing systems 29 (2016).
  * Mehrabi et al. [2021] N. Mehrabi, F. Morstatter, N. Saxena, K. Lerman, A. Galstyan,  A survey on bias and fairness in machine learning,  ACM Computing Surveys (CSUR) 54 (2021) 1–35.
  * Kearns et al. [2018] M. Kearns, S. Neel, A. Roth, Z. S. Wu,  Preventing fairness gerrymandering: Auditing and learning for subgroup fairness,  in: J. Dy, A. Krause (Eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, PMLR, 2018, pp. 2564–2572.
  * Lum et al. [2022] K. Lum, Y. Zhang, A. Bower,  De-biasing “bias” measurement,  in: Proceedings of the 2022 ACM Conference on Fairness, Accountability, and Transparency, FAccT ’22, Association for Computing Machinery, New York, NY, USA, 2022, p. 379–389. doi:10.1145/3531146.3533105.
  * Pearson [1895] K. Pearson,  Notes on regression and inheritance in the case of two parents,  Proceedings of the Royal Society of London 58 (1895) 240–242.
  * Spearman [1904] C. Spearman,  The proof and measurement of association between two things,  American Journal of Psychology 15 (1904) 72–101.
  * Kendall [1938] M. G. Kendall,  A new measure of rank correlation,  Biometrika 30 (1938) 81–93.
  * Ghassami et al. [2018] A. Ghassami, S. Khodadadian, N. Kiyavash,  Fairness in supervised learning: An information theoretic approach,  in: 2018 IEEE International Symposium on Information Theory (ISIT), IEEE Press, 2018, p. 176–180. doi:10.1109/ISIT.2018.8437807.
  * Pedregosa et al. [2011] 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, E. Duchesnay,  Scikit-learn: Machine learning in Python,  Journal of Machine Learning Research 12 (2011) 2825–2830.
  * Fawcett [2006] T. Fawcett,  An introduction to ROC analysis,  Pattern Recognition Letters 27 (2006) 861–874. doi:10.1016/j.patrec.2005.10.010.