# Pole structure of the electronic self-energy with coexistence of Charge
order and Superconductivity

M. Grandadam  Institut de Physique Théorique, Université Paris-Saclay, CEA,
CNRS, F-91191 Gif-sur-Yvette, France  C. Pépin  Institut de Physique
Théorique, Université Paris-Saclay, CEA, CNRS, F-91191 Gif-sur-Yvette, France

###### Abstract

We compare the pole structure of the electronic Green’s function obtained by
Cluster Dynamical Mean Field Theory to the results from the fractionalized
Pair Density Wave idea. In the superconducting phase, we can consider the
system in a state with coexistence of Superconducting and Charge order.
Writing the Green’s function in a way analogous to the previously proposed
“hidden-fermions” model from Ref.[Sakai  _et al._ , 2016a] leads to a similar
pole structure for the self-energy. The fractionalization of the Pair Density
Wave order also describes the pseudogap phase as a superposition of
superconducting and charge order fluctuations. Considering a phenomenological
lifetime for the particle-particle and particle-hole pairs leads to an
electronic spectral function that matches the numerical results.

##  I Introduction

Since their discovery, cuprate superconductors have been the focus of many
investigations to understand the nature of the relation between the high-
temperature superconducting phase and the pseudogap that appears when they are
doped away from half-fillingAlloul  _et al._ (1989); Warren  _et al._ (1989).
They have been the subject of an important amount of work on the theoretical,
experimental and numerical side. In the latter case, a lot of efforts have
been invested into the development of the Dynamical Mean Field Theory that
allows for solving the Hubbard model by mapping it onto an impurity
problemGeorges  _et al._ (1996). The technique has been extended to include
multiple sites in the Cluster Dynamical Mean Field Theory (CDMFT)Kotliar  _et
al._ (2001); Maier  _et al._ (2005) and gives precise results on the energy
dependence of the electronic Green’s function at specific points of the
Brillouin zoneGull  _et al._ (2010); Sordi  _et al._ (2010, 2012). Recent
results from CDMFT on 2x2 clusters have shown a peculiar link between the
structure of the electronic self-energy in the pseudogap (PG) and in the
superconducting (SC) phaseSakai  _et al._ (2016a, b, 2018) that deviate from
the standard BCS theory. Multiple propositions have been made to explain this
with an analytical model such as a hidden-fermion modelSakai  _et al._ (2016a)
or a $SU\left(2\right)$ gauge theoryScheurer  _et al._ (2018).
On the experimental side, there is growing evidence that the charge ordering
observed above the superconducting temperatureDoiron-Leyraud  _et al._ (2007);
Sebastian  _et al._ (2012); Chang  _et al._ (2012); Blanco-Canosa  _et al._
(2013); Blackburn  _et al._ (2013); Ghiringhelli  _et al._ (2012); Gerber  _et
al._ (2015); Chang  _et al._ (2016a); Wu  _et al._ (2011, 2015); Julien (2015)
plays an important role in the pseudogap physics. In fact, recent Raman
experimentLoret  _et al._ (2019) suggests that the Charge Density Wave (CDW)
gap is of the same order as the superconducting gap. This charge order is also
known to coexist with the superconductivity at low temperature and long-range
CDW is observed when superconductivity is destroyed by a magnetic fieldGerber
_et al._ (2015); Chang  _et al._ (2016a); LeBoeuf  _et al._ (2013). A recent
proposal linked these two orders by the fractionalization of modulated
particle-particle pairs, _i.e._ of a Pair Density Wave (PDW) orderChakraborty
_et al._ (2019); Grandadam  _et al._ (2020a), leading to a constraint between
the resulting CDW and SC amplitudes. The resulting superposition of
fluctuations in both channels was used to explainGrandadam  _et al._ (2020b)
the momentum, energy and temperature dependence of the electronic spectral
function observed in Bi2201 by angle-resolved photo-emission spectroscopyHe
_et al._ (2011) (ARPES).
We are here interested in the possibility that the fine structure in the self-
energy obtained by CDMFT calculationHaule and Kotliar (2007); Sakai  _et al._
(2016a); Gull  _et al._ (2013) on the Hubbard model could be explained by the
same scenario. In the superconducting phase, the electronic Green’s function
is given by considering a coexistence of SC and CDW order. We write the
electronic Green’s function in the superconducting phase in a way analogous to
the hidden-fermion proposition that was used to describe CDMFT resultsSakai
_et al._ (2016a, b, 2018) by coupling electrons to “hidden” fermionic
excitations on top of pairing. We identify the hidden-fermion with a
modulating order and show that it leads to a pole structure in the self-energy
close to the one obtained via numerical calculation below $T_{c}$. The
fractionalized PDW scenario gives a description of the pseudogap that relies
on the superposition of SC and CDW fluctuations in the antinodal region and is
different from the hidden-fermion model mentioned previously. We show that the
fractionalized PDW still reproduces the electronic Green’s function in the
pseudogap phase.

##  II Electronic spectral function from CDMFT

Cluster extension of CDMFT with 2x2 sites gives detailed energy-dependence of
the electronic Green’s function in the Hubbard model at low doping. Due to the
small system size, this is however limited to some specific points of the
Brillouin zone. One of the most important results concerns the anti-nodal
point at $k=\left(0,\pi\right)$ where the pseudogap is observed above the
superconducting state. CDMFT allows us to decompose the deviation from the
non-interacting Green’s function into two parts : the normal and anomalous
self-energy. The final form of the electronic Green’s function is given by :

$\displaystyle G\left(k,\omega\right)$
$\displaystyle=\left[\omega-\epsilon_{k}-\Sigma_{N}\left(k,\omega\right)-W\left(k,\omega\right)\right]^{-1},$
$\displaystyle W\left(k,\omega\right)$
$\displaystyle=\frac{\Sigma_{AN}\left(k,\omega\right)^{2}}{\omega+\epsilon_{k}-\Sigma_{N}\left(k,-\omega\right)^{*}}$
(1)

The anomalous part of the self-energy is directly related to the
superconducting order and vanishes outside of the ordered phase but the normal
part presents features that lead to the pseudogap in the anti-nodal region at
higher temperature. Another quantity of interest is the electronic spectral
function
$A\left(k,\omega\right)=\frac{-1}{\pi}\text{Im}\left(G\left(k,\omega\right)\right)$
which relate directly to ARPES experiment for occupied states ($\omega<0$).
Interestingly there are strong links between the structure in the normal self-
energy above and below $T_{c}$ as well as connection with the anomalous self-
energy. These seems to indicate that the pseudogap physics is strongly related
to the superconducting state.

Figure 1: (a) Normal self-energy (red line) and electronic spectral function
(green dotted lines) obtained by CDMFT for momenta on the line
$k=\left(0,\pi\right)-\left(\pi/2,\pi/2\right)$. The self-energy shows two
isolated poles in the anti-nodal region at $\omega=\pm\omega_{0}$ but with a
clear asymmetry in weight. In the nodal region the peak at positive energy
splits. The pole at negative energy is damped is difficult to follow.(b)
Anomalous self-energy obtained by CDMFT for the same momenta. It shows pairs
of poles at the same position as in the normal part but with an anti-symmetric
weight. The same splitting occurs when going from the anti-nodal to the nodal
region. The weight close to $k=\left(\pi/2,\pi/2\right)$ is close to zero,
indicating that the coupling to the hidden-fermion is specific to the anti-
nodal region.

The energy dependence of the normal self-energy (red lines) and of the
electronic spectral function (green dotted line) for momenta on a line from
$k=\left(0,\pi\right)$ to $k=\left(\pi/2,\pi/2\right)$ is shown in Fig.1(a) in
the superconducting phase.The distinctive feature of the self-energies in the
anti-nodal region (at $k=\left(0,\pi\right)$) is the two symmetric peaks at
low energy. The two poles of $\Sigma_{N}$ have a marked weight asymmetry with
the pole at $\omega>0$ having a larger spectral weight than the one at
negative energies. The spectral function plotted at the same momentum (green
dotted line in Fig.1(a)) presents 3 poles, one at negative energy which is
well defined and 2 at positive energies which are broader. In contrast, the
two poles of $\Sigma_{AN}$ shown in Fig.1(b) for the same momenta have the
same position as the poles of $\Sigma_{N}$ but opposite weight. When going
toward the nodal region we observe that the pole at positive energy split into
two poles of approximately equal weight. This is also the case for the pole at
negative energy even though the vanishing weight close to
$k=\left(\pi/2,\pi/2\right)$ makes it harder to pinpoint the position of the
two peaks.
Another crucial information given by CDMFT is the cancellation that occurs
between the poles of $\Sigma_{N}$ and $W$. Indeed it has been observed that,
besides having poles at the same energy, the residue of both the normal and
the anomalous the self-energy at these poles is such that they cancel in the
expression of the single-particle Green’s function. This has consequences for
the electronic spectral function as the poles of the self-energy should
correspond to zeros of the spectral function, effectively leading to the
splitting of the non-interacting band.
These results were interpreted in terms of ”hidden-fermions” which couple to
the original electronic degree of freedom and are themself susceptible to
pairingSakai  _et al._ (2016a, b, 2018). It is then possible to write a mean-
field Hamiltonian for this model and extracts analytical expression for the
normal and anomalous self-energy. As we will show in the next section, this
construction guarantees the cancellation between $\Sigma_{N}$ and $W$
mentioned earlier. It also displays a structure for both parts of the self-
energy in agreement with the numerical results with proper choice for the
hidden-fermion. We will present the model and its main features in the next
section and then show how to incorporate charge order in the hidden-fermion
formalism and compare the results to CDMFT in both the superconducting and the
pseudogap state.

##  III Hidden-fermion Model

We start with the same hidden-fermion model from Ref.Sakai  _et al._ , 2016a.
This model describe electrons ($c^{\dagger}_{k,\sigma}$) on a square lattice
coupled to other fermionic excitation ($f^{\alpha\dagger}_{k,\sigma}$) :

$\displaystyle H=$
$\displaystyle\sum_{k,\sigma}\epsilon_{k}c^{\dagger}_{k,\sigma}c_{k,\sigma}+\sum_{k,\sigma}\left(\sigma\Delta_{k}c^{\dagger}_{k,\sigma}c^{\dagger}_{-k,-\sigma}+h.c\right)$
$\displaystyle+\sum_{\alpha,k,\sigma}\epsilon^{f,\alpha}_{k}f^{\alpha\dagger}_{k,\sigma}f^{\alpha}_{k,\sigma}+\sum_{\alpha,k,\sigma}\left(\sigma\Delta^{f,\alpha}_{k}f^{\alpha\dagger}_{k,\sigma}f^{\alpha\dagger}_{-k,-\sigma}+h.c\right)$
$\displaystyle+\sum_{\alpha,k,\sigma}\left(V^{\alpha}_{k}c^{\dagger}_{k,\sigma}f^{\alpha}_{k,\sigma}+h.c\right)$
(2)

The electronic spectral function can then be obtain and the coupling to the
different $f^{\alpha}_{k,\sigma}$ excitation leads to a normal and anomalous
self-energy given by :

$\displaystyle\Sigma_{N}\left(k,\omega\right)=\sum_{\alpha}\frac{{V_{k}^{\alpha}}^{2}\left(\omega+\epsilon_{k}^{f,\alpha}\right)}{\omega^{2}-{\epsilon_{k}^{f,\alpha}}^{2}-{\Delta_{k}^{f,\alpha}}^{2}}$
$\displaystyle\Sigma_{AN}\left(k,\omega\right)=\Delta_{k}+\sum_{\alpha}\frac{{-V_{k}^{\alpha}}^{2}\Delta_{k}^{f,\alpha}}{\omega^{2}-{\epsilon_{k}^{f,\alpha}}^{2}-{\Delta_{k}^{f,\alpha}}^{2}}$
(3)

As can be seen from Eq.(3) the normal and anomalous part of the self-energy
share a similar pole structure which is primarily governed by the choice of
$\Delta^{\alpha}_{k}$ and $\epsilon^{\alpha}_{k}$. The value of the coupling
between the $c^{\dagger}_{k,\sigma}$ and the $f^{\alpha\dagger}_{k,\sigma}$
fermions only enters in the spectral weight at the different poles as can be
seen by computing the residue of the self-energy :

$\displaystyle\text{Res}\left(\Sigma_{N},\pm\omega^{\alpha}\right)=\frac{{V_{k}^{\alpha}}^{2}}{2}\left(1\pm\frac{\epsilon_{k}^{f,\alpha}}{\sqrt{{\epsilon_{k}^{f,\alpha}}^{2}+{\Delta_{k}^{f,\alpha}}^{2}}}\right)$
$\displaystyle\text{Res}\left(\Sigma_{AN},\pm\omega^{\alpha}\right)=\mp\frac{{V_{k}^{\alpha}}^{2}\Delta_{k}^{\alpha}}{2\omega^{\alpha}}$
(4)

We can see that the self-energies have poles due to the hidden-fermions at the
energies
$\omega^{\alpha}=\pm\sqrt{{\epsilon^{f,\alpha}_{k}}^{2}+{\Delta_{k}^{f,\alpha}}^{2}}$.
There is however an important distinction between the normal part and the
anomalous part of the self-energy in the weight associated with each of these
poles. In fact, the anomalous self-energy has equal and opposite weight at
each pole while it depends on the value of $\epsilon^{f,\alpha}_{k}$ for the
normal self-energy. In particular the sign of $\epsilon^{f,\alpha}_{k}$ will
determine which of $\omega^{\alpha}$ or $-\omega^{\alpha}$ has a higher
weight. It is also possible to check the cancellation mentioned previously :

$\displaystyle\text{Res}\left(W,\pm\omega^{\alpha}\right)$
$\displaystyle=\text{Res}\left(\frac{\Sigma_{AN}\left(k,\omega\right)^{2}}{\omega+\epsilon_{k}-\Sigma_{N}\left(k,-\omega\right)^{*}},\pm\omega^{\alpha}\right)$
$\displaystyle=-\frac{{V_{k}^{\alpha}}^{2}}{2}\left(1\pm\frac{\epsilon_{k}^{f,\alpha}}{\sqrt{{\epsilon_{k}^{f,\alpha}}^{2}+{\Delta_{k}^{f,\alpha}}^{2}}}\right)$
$\displaystyle=-\text{Res}\left(\Sigma_{N},\pm\omega^{\alpha}\right)$ (5)

The original proposition was that only one of these hidden fermionic
excitations is relevant to explain the structure of the self-energy at low
energy and that identifying
$\epsilon^{f}_{k}=z_{k}\epsilon_{k+\bm{\pi}}-\mu_{f}$, where $z_{k}$ is the
renormalization factor due to the strong interactions and
$\bm{\pi}=\left(\pi,\pi\right)$, leads to a pole structure analogous to the
CDMFT results.

Figure 2: Fermi surface for free electrons with dispersion given by
$\epsilon_{k}=-2t\left(\cos\left(k_{x}\right)+\cos\left(k_{y}\right)\right)-4t^{\prime}\cos\left(k_{x}\right)\cos\left(k_{y}\right)-\mu$
where $t=1$, $t^{\prime}=-0.2$ and $\mu=-0.56$. The black arrow indicate the
charge order wave-vector taken to link different part of the Fermi surface
while the red line shows the path between $k=\left(0,\pi\right)$ and
$k=\left(\pi/2,\pi/2\right)$ where the self energy is computed.

We show here that another choice of hidden fermion can lead to satisfactory
comparison with the results obtained by CDMFT but have a different
interpretation. Our choice of hidden fermions is based on the experimental
observation of charge order in cuprates (by X-RayChang  _et al._ (2012);
Blanco-Canosa  _et al._ (2013); Campi  _et al._ (2015); Blackburn  _et al._
(2013); Ghiringhelli  _et al._ (2012); Gerber  _et al._ (2015); Chang  _et
al._ (2016b) or NMRWu  _et al._ (2011, 2013, 2015) measurements) that persists
in the superconducting phase. This coexistence of CDW and SC has been shown to
successfully reproduce the single-particle electronic Green’s function
observed in Bi2201 by ARPESGrandadam  _et al._ (2020b) and the softening
observed in the phonon spectrum at $T_{c}$Sarkar  _et al._ (2020). Our
analysis is thus based on the previous hidden-fermion model Eq.(2) but we
consider that two fermionic excitation will contribute and identify them with
$f^{\alpha\dagger}_{k,\sigma}=c^{\dagger}_{k+Q,\sigma}$ and
$f^{\beta\dagger}_{k,\sigma}=c^{\dagger}_{k-Q,\sigma}$.

We focus on an axial wave-vector that links different parts of the Fermi
surface, _i.e._ $Q=\left(\pm Q_{0},0\right)$ as shown in Fig.2. We also take
into account strong correlations by adding a renormalization factor
$z\left(k,\omega\right)$ that will modify the dispersion of the hidden fermion
so that $\epsilon^{f,\alpha}_{k}=z\epsilon_{k+Q}$ and
$\epsilon^{f,\beta}_{k}=z\epsilon_{k-Q}$ where we took
$z\left(k,\omega\right)$ to be a constant. The bare dispersion is given by
taking nearest and second neighbours hopping
$\epsilon_{k}=-2t\left(\cos\left(k_{x}\right)+\cos\left(k_{y}\right)\right)-4t^{\prime}\cos\left(k_{x}\right)\cos\left(k_{y}\right)-\mu$.
We took $t=1$ as our energy scale, $t^{\prime}=-0.2$ and $\mu$ is chosen to
fix the electron density to $n=0.95$ using the Luttinger sum rule in the non-
interacting case. Lastly, as we are interested in the superconducting phase,
we take the SC gap with a d-wave form factor
$\Delta_{k}=\frac{\Delta_{0}}{2}\left(\cos\left(k_{x}\right)-\cos\left(k_{y}\right)\right)$
and the other gaps are given by $\Delta^{f,\alpha}_{k}=\Delta_{k+Q}$ and
$\Delta^{f,\beta}_{k}=\Delta_{k-Q}$.

##  IV Results

###  IV.1 Normal self-energy in the superconducting state

To relate the different features of the normal self-energy to the properties
of the hidden-fermions we need to study the structure of the normal part of
the self-energy given by Eq.(3). The positions of the poles are given by
$\omega^{\alpha}_{\pm}=\pm\sqrt{{\epsilon_{k}^{f,\alpha}}^{2}+{\Delta_{k}^{f,\alpha}}^{2}}$
and thus come by pair with symmetric position in energy. The weight associated
with each pole is obtained by computing the residue at each of them and is
given in Eq.(4). The asymmetry in weight is thus controlled by the sign of
$\epsilon^{f,\alpha}_{k}$ and the value of $\Delta^{f,\alpha}_{k}$. As noted
in Ref.[Sakai  _et al._ , 2016a], the observation that the pole at positive
energy has a larger weight than the pole at negative energy in the anti-nodal
region impose the fact that $\epsilon^{f,\alpha}_{\left(0,\pi\right)}$ has to
be positive. This constraint is naturally fulfilled in our model when we
choose a CDW wave-vector larger than the anti-nodal Fermi wave-vector which
also fits experimental observation. This is however no longer the case for the
fermion identified with $c_{k-Q,\sigma}$ when we approach the nodal region
which explains the discrepancy between our model and the CDMFT results near
$k=\left(\frac{\pi}{2},\frac{\pi}{2}\right)$.
The results for the original hidden-fermion ideaSakai  _et al._ (2016a) that
identifies $f^{\alpha\dagger}_{k,\sigma}=c^{\dagger}_{k+\bm{\pi},\sigma}$ are
shown in Fig.3(a) and match multiple characteristics of the numerical results.
Notably, the asymmetry in the spectral weight in the anti-nodal region is well
reproduced and the loss of weight for the peak at negative energy is found.
Note that, following Eq.(4), this happens when $\epsilon^{\alpha}_{k}$ is
positive. There is thus a crucial role of the added $\mu_{f}$ term as the
hidden-fermion dispersion $z_{k}\epsilon_{k+\pi}$ alone is negative for
$k=\left(0,\pi\right)$. The combination of the renormalization induced by
$z_{k}$ and the positive shift of $\mu_{f}$ gives an effective dispersion that
is small relative to $\Delta_{k}$ and that does not change sign. However,
because there is only one hidden fermion it is not possible to capture the
splitting in the nodal region for $\omega>0$ and the weight of the positive
peak is maximal at $k=\left(\pi/2,\pi/2\right)$ when the d-wave
superconducting gap vanishes, in contrast to the numerical results.

Figure 3: a Normal self-energy in the hidden-fermion model of Ref.Sakai  _et
al._ , 2016a. The asymmetry in the anti-nodal region is well reproduced. As
there is only one hidden-fermion coupled to electrons it is not possible to
recover the splitting in the nodal region. The weight of the positive energy
pole is also maximum for $k=\left(\pi/2,\pi/2\right)$ in contrast to the
numerical results shown in Fig.2. (b) Normal self-energy obtained while
considering hidden-fermions due to CDW order with a wave-vector along the $x$
axis. Due to the symmetry between $\epsilon_{k+Q}=\epsilon_{k-Q}$ at
$k=\left(0,\pi\right)$ we only have two visible poles with an asymmetry
consistent with the CDMFT results. This symmetry is lost when going to the
nodal region and we observe a splitting of the poles. The change of sign of
$\epsilon_{k-Q}$ close to $k=\left(\pi/2,\pi/2\right)$ explains the
significant weight for the pole at negative energy. We choose here
$Q_{0}=0.27\pi$, $z=0.22$, $\Delta_{0}=0.55$, $V_{0}=0.7$ and a numerical
broadening factor $i\eta=0.03i$.

We now turn to the fractionalized PDW hypothesis which gives the same form for
the electronic Green’s function, in the superconducting phase, to the hidden-
fermion model with $f^{\alpha\dagger}_{k,\sigma}=c^{\dagger}_{k+Q,\sigma}$ and
$f^{\beta\dagger}_{k,\sigma}=c^{\dagger}_{k-Q,\sigma}$. Our results for the
two hidden fermions coming from the charge order are shown in Fig.3(b). At
$k=\left(\pi,\pi\right)$ we can only see a pair of symmetric peaks due to the
fact that $\epsilon_{k+Q}=\epsilon_{k-Q}$. Moreover, the charge order wave-
vector being larger than the anti-nodal Fermi momentum means that we have
$\epsilon_{k\pm Q}>0$. Using the result of Eq.(4), this leads to the same
weight asymmetry observed in the CDMFT study with the pole at
$+\omega^{\alpha}$ having a higher weight than the pole at $-\omega^{\alpha}$.
When going towards the nodal region the previous symmetry between
$\epsilon_{k+Q}$ and $\epsilon_{k-Q}$ is lost and we observe a splitting of
the poles at $\omega>0$ analogous to the numerical results. It is important to
point out that the behaviours of the spectral weight for the two peaks are
different. In fact, the pole due to the coupling to $c^{\dagger}_{k+Q,\sigma}$
will have a higher weight at $-\omega^{\alpha}$ in the nodal region. This is
because $\epsilon_{k-Q}$ will change sign and become negative. Because the
pairing gap for the hidden-fermion is taken as
$\Delta_{k}^{f,\alpha}=\Delta_{k\pm Q}$, there is no cancellation of the
weight at $k=\left(\pi/2,\pi/2\right)$ which result in a non-vanishing weight
in the nodal region for the pole at $\omega<0$. Note that we took here a
coupling between the electrons and the hidden-fermions to be independent of
momentum, _i.e_ $V_{k}^{f,\alpha}=V_{0}$. This does not impact the position of
the poles as shown by Eq.(3) but only the weight at each pole (see Eq.(4)).
Furthermore, results from CDMFT away from the antinodal point
$k=\left(0,\pi\right)$ are extrapolated from the available momentum points.
Thus, the discrepancy between the spectral weight in the nodal region can be
due to multiple factors.

###  IV.2 Pseudogap phase

We showed that both the antiferromagnetic hidden-fermion and the
fractionalized PDW models give similar results in the superconducting phase.
The two models differ strongly in the pseudogap phase where the hidden-fermion
model restores the charge conservation symmetry but keeps the coupling between
electrons and the hidden fermions unchanged. In contrast, the
fractionalization of a Pair Density Wave leads to the superposition of SC and
CDW fluctuations that have a strong impact on the electronic spectral
functions.

####  IV.2.1 Pseudogap in the hidden-fermion model

In the hidden-fermion model, the pseudogap phase is obtained by setting the
superconducting order parameters $\Delta_{k}$ and $\Delta_{k}^{f,\alpha}$ to
zero. The remaining part of the self-energy is due to the hybridization to the
hidden-fermions through $V^{\alpha}$. The self-energy thus take the form:

$\displaystyle\Sigma_{N}\left(k,\omega\right)=\sum_{\alpha}\frac{{V^{\alpha}_{k}}^{2}}{\omega-{\epsilon^{f,\alpha}}_{k}}$
$\displaystyle\Sigma_{AN}\left(k,\omega\right)=0$ (6)

Figure 4: (a) Electronic spectral function (full line) and self-energy (dotted
line) in the pseudogap (red) and superconducting phase (blue) in the hidden-
fermion model. Above $T_{c}$, the self-energy is reduced to a single pole at
$\omega=\epsilon_{k\pm Q}>0$ and the resulting electronic spectral function
shows 2 poles. There is no spectral weight at $\omega=0$ but there will be
gapless states away from $k_{x}=0$ in contrast to experimental observations.
(b) Electronic spectral function (full line) and self-energy (dotted line) in
the pseudogap phase from a fractionalized PDW at $k=\left(0,\pi\right)$. Green
lines show effects of CDW fluctuations ($\Gamma_{1}$) while the red lines show
the effect of superconducting fluctuations ($\Gamma_{2}$). Due to the
superposition of SC and CDW in the fractionalized PDW idea, the electronic
spectral function retains the 4 poles structure seen in the superconducting
phase. Finite lifetime for the particle-particle or particle-hole pairs have a
much stronger dampening effect in the two poles closest to $\omega=0$ which
could explain the observed two peaks structure in CDMFT or the single occupied
band in ARPES above $T_{c}$. We choose here $Q_{0}=0.27\pi$, $z=0.22$,
$|\Psi_{0}|=V_{0}=0.15$, $\Gamma_{1}=\Gamma_{2}=0.1$ and a numerical
broadening factor $i\eta=0.03i$.

The structure of the self-energy is then reduced to peaks located at the
energy of the hidden-fermions as shown in Fig.4(a) (red dotted line). In our
case, the self-energy in the anti-nodal region will display a single peak at
$\omega>0$ because $\epsilon_{\bm{\pi}+Q}=\epsilon_{\bm{\pi}-Q}>0$. This leads
to the electronic spectral function having two poles with a vanishing spectral
weight at $\omega=0$ for $k=\left(0,\pi\right)$ (red line in Fig.4(a)). Note
that a non-zero value of hybridization between $c_{k,\sigma}$ and
$c_{k+Q,\sigma}$ should be interpreted as a long-range CDW at the mean-field
level. As such, we should observe a Fermi surface reconstruction with gapless
states in the pseudogap phase. This long-range order has only been observed by
applicatying strong magnetic fieldsGerber  _et al._ (2015); Chang  _et al._
(2016a); LeBoeuf  _et al._ (2013) in the superconducting phase. Furthermore,
the short-range charge order above $T_{c}$ is only observed below a transition
temperature $T_{CO}$ which is lower than the pseudogap temperature $T^{*}$.
The band structure observed by ARPES does not however show significant change
for temperatures in the pseudogap region. It is thus inconsistent to consider
that the pseudogap phase is solely due to the charge order.

####  IV.2.2 Fractionalized Pair Density Wave

We start by giving details on the Pair Density Wave fractionalization as
described in Ref.[Chakraborty  _et al._ , 2019; Grandadam  _et al._ , 2020a].
The fractionalization of the PDW operator is based on the relation between the
$\eta$-mode, which is a modulated particle-particle pair and the SC
($\hat{\Delta}_{ij}=J\hat{d}_{ij}\sum_{\sigma}\sigma
c_{i,\sigma}c_{j,-\sigma}$) and the CDW
($\hat{\chi}_{ij}=J\hat{d}_{ij}\sum_{\sigma}c_{i,\sigma}^{\dagger}c_{j,\sigma}e^{i\mathbf{Q.\left(\mathbf{r}_{i}+\mathbf{r}_{j}\right)/2}}$)
operators where $J$ is the typical interaction strength giving rise to the SC
and CDW orders. In particular we can write :

$\displaystyle\hat{\eta}$
$\displaystyle=[\hat{\Delta}_{ij},\hat{\chi}_{ij}^{\dagger}],$
$\displaystyle\hat{\eta}^{\dagger}=[\hat{\chi}_{ij},\hat{\Delta}_{ij}^{\dagger}],$
(7)

where $\bm{Q}$ is the modulation wave-vector of the PDW. In a way analogous to
the electron’s fractionalization in strong coupling theoryBaskaran and
Anderson (1988); Lee and Nagaosa (1992), the $\eta$-operators are then
invariant with the following gauge structure

$\displaystyle\hat{\Delta}_{ij}$ $\displaystyle\rightarrow
e^{i\theta}\hat{\Delta}_{ij},$ $\displaystyle\hat{\chi}_{ij}\rightarrow
e^{i\theta}\hat{\chi}_{ij}.$ (8)

The effective field theory for the fluctuating PDW thus involved an emergent
$U(1)$ gauge field whose fluctuation produces a constraint between the two
fields:

$\displaystyle|\hat{\Delta}_{ij}|^{2}+|\hat{\chi}_{ij}|^{2}\equiv\left|\Psi_{ij}\right|^{2}=\text{const},$
(9)

where $\Psi_{ij}=(\hat{\Delta}_{ij},\hat{\chi}_{ij})^{t}$. The energy scale
associated with Eq.(9) is typically the scale at which the fractionalization
occurs.

Earlier descriptions of the pseudogap phase as a fractionalized PDW led to
good agreement with experimental ARPES results for the electronic spectral
functionGrandadam  _et al._ (2020b). The main characteristic of the pseudogap
in this theory is to be an equal superposition of SC and CDW fluctuations in
the anti-nodal regions. This is achieved at the mean-field level by
considering the previous model but constraining the amplitude for the SC and
CDW amplitude to be equal while adding finite lifetime $\Gamma_{1}$ and
$\Gamma_{2}$ for the particle-hole and for the particle-particle pairs
respectively. Because none of the parameters is put to zero, the self-energy
and the spectral function have the same pole structure in both the pseudogap
and the superconducting phases. In fact, the electronic spectral function at
$k=\left(0,\pi\right)$ given in Fig.4(b) with both damping rate $\Gamma_{1}$
and $\Gamma_{2}$ being turned on alternatively clearly shows 4 poles indicated
by arrows. The main observation is that there is a significant reduction in
the spectral weight of the two poles at lower energy for both positive and
negative energy. This is the same phenomenon that was used to describe the
observation of a “flat band” below $T_{c}$ by ARPES in Bi2201Grandadam  _et
al._ (2020b). Similarly to this previous study, this description of the
pseudogap effectively leads to Fermi arcs and no long-range order is expected.
We show the resulting electronic spectral function at $\omega=0$ in the
antinodal region in Fig.5(b). We can see that part of the Fermi surface is
washed out even in the absence of fluctuations ($\Gamma_{1}=\Gamma_{2}=0$).
Despite having a similar structure for the self-energy in the pseudogap phase,
there is an important difference between the structure of the Green’s function
in our approach and in previous studies based on fractionalization such as the
one using a $SU(2)$ theoryScheurer  _et al._ (2018) of fluctuating
antiferromagnetism. In our case, the line of zeros for the real part of the
Green’s function (black dotted line in Fig.5) is very close to the non-
interacting Fermi surface. This is in strong contrast to other cases where the
Luttinger surface created by interaction with the antiferromagnetic
fluctuations intersect the Fermi surface to create small pockets in the nodal
regionYang  _et al._ (2006).
It is important to note that in order to obtain a good fit with the numerical
result obtained by CDMFT, it is necessary here to take a value of $V^{\alpha}$
(in the superconducting phase) and $\Psi$ (in the pseudogap phase)
significantly lower than the value of the pairing amplitude $\Delta$ in the
superconducting phase. Indeed, an agreement between the experimental ARPES
results and the fractionalized PDW scenario was obtained with the SC and CDW
order being very close, a fact supported by Raman experiment in mercury-based
cuprateLoret  _et al._ (2019). This would indicate that only a part of the
particle-particle pairs come from the fractionalization while the other part
can be described by a standard superconducting condensate. This can also be
compared to the two behaviours observed in Raman experiment for the nodal and
antinodal part of the superconducting gapLe Tacon  _et al._ (2006). In fact,
the nodal gap follows the critical temperature with doping and decreases away
from optimal doping while the antinodal gap decreases linearly with doping,
following the pseudogap temperature. This imbalance could be here a result of
the finite size of the cluster (2x2 sites) used in the CDMFT calculation which
is smaller than the wavelength of the charge modulation, thus hindering the
formation of modulated orders.

Figure 5: (a) Electronic spectral function at $\omega=0$ for the
fractionalized PDW model. The white dotted line indicates the non-interacting
Fermi surface while the black dotted line indicates the surface of zeros of
the Green’s function. We see that the Fermi surface is washed out in the
antinodal region. We choose here $Q_{0}=0.27\pi$, $z=0.22$, $|\Psi_{0}|=0.15$
and a numerical broadening factor $i\eta=0.03i$.

##  V Discussion

The results obtained from CDMFT calculations give valuable information on the
energy dependence of the electronic Green’s function in the doped Hubbard
model. This is often limited to specific points in the Brillouin zone which
depend on the cluster size. From this, it seems difficult to account for
potential modulating orders with a wavelength larger than the cluster size.
There is however a large number of experimental evidenceDoiron-Leyraud  _et
al._ (2007); Sebastian  _et al._ (2012); Chang  _et al._ (2012); Blanco-Canosa
_et al._ (2013); Blackburn  _et al._ (2013); Ghiringhelli  _et al._ (2012);
Gerber  _et al._ (2015); Chang  _et al._ (2016a); Wu  _et al._ (2011, 2015);
Julien (2015); Loret  _et al._ (2019) for these orders to exist and compete in
the pseudogap phase of cuprates superconductors. We argued in this paper that
despite the numerical limitation, there are signatures of the charge order in
the electronic self-energy calculated in CDMFT.

Much like in ARPES experiments, the observations are indirect consequences of
the charge order and there is no Fermi surface reconstruction as the charge
order remain short-ranged and fluctuating. In the case of ARPES in Bi2201, the
main indicators for a modulating order were the back-bending of the band at
the Brillouin zone’s edge and the specific way the pseudogap was closing close
to the nodal region (for more details see Ref.[Lee, 2014; He  _et al._ ,
2011]). The limited momentum resolution of CDMFT does not allow for such
distinction but the singular pole structure of the self-energy in the anti-
nodal region could itself contain evidence for this modulating order. It is
important to note the similarities between the original hidden-fermion from
Ref.[Sakai  _et al._ , 2016a] and our proposition which contain the charge
order. In particular both of them give the same result for
$k\left(0,\pi\right)$. There are however distinctions to be made when going
close to the nodal region where the two models differ. The main argument for
the CDW scenario is that it allows for the double pole structure observed for
$\omega>0$ (see Fig.3) which is observed numerically. This is to be taken with
caution as both the numerical results and the idea of CDW order reaches their
limit in the nodal region.

A more convincing argument may come from the electronic spectral function in
the pseudogap phase. The original idea for the hidden-fermion is to describe
the pseudogap by the absence of pairing. At the mean-field level, this would
lead to a CDW long-range order with the associated Fermi surface
reconstruction. The spectral function at the anti-nodal point does show the 2
poles observed by CDMFT with a vanishing weight at $\omega=0$ but there are
still gapless states away from $k_{x}=0$. In contrast, the fractionalized PDW
scenario is characterized by the equal superposition of SC and CDW
fluctuations in the pseudogap. As such the spectral function has the same
structure as in the superconducting state and is fully gapped at the Brillouin
zone’s edge, leaving Fermi arcs in the nodal region if we include a momentum
dependence for the pseudogap orderGrandadam  _et al._ (2020a, b). Taking into
account fluctuations in the SC and CDW orders leads to a drastic reduction in
the spectral weight of specific poles, giving an effective spectral function
with two apparent peaks. The same argument was given to explain the
observation of a “flat band” in Bi2201 below $T_{c}$.

There are still many differences between the results from CDMFT and the
experimental spectral function measured by ARPES. It is however interesting to
find common ground to link these two essential techniques to investigate the
pseudogap phase of cuprates. The capacity of CDMFT to explore the positive
energies which are inaccessible to ARPES and the wide range of doping
available to the experimental are complementary to elucidate the pseudogap
phase of cuprates.

##  VI ACKNOWLEDGMENTS

We thank M. Civelli, A. Banerjee, S. Sarkar and Y. Sidis for valuable
discussions. This work has received financial support from the ERC, under
grant agreement AdG-694651-CHAMPAGNE.

## References

  * Sakai  _et al._ (2016a) S. Sakai, M. Civelli, and M. Imada, Phys. Rev. Lett. 116, 057003 (2016a).
  * Alloul  _et al._ (1989) H. Alloul, T. Ohno, and P. Mendels, Phys. Rev. Lett. 63, 1700 (1989).
  * Warren  _et al._ (1989) W. W. Warren, R. E. Walstedt, G. F. Brennert, R. J. Cava, R. Tycko, R. F. Bell, and G. Dabbagh, Phys. Rev. Lett. 62, 1193 (1989).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Reviews of Modern Physics 68, 13 (1996).
  * Kotliar  _et al._ (2001) G. Kotliar, S. Y. Savrasov, G. Palsson, and G. Biroli, Phys. Rev. Lett. 87, 186401 (2001).
  * Maier  _et al._ (2005) T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
  * Gull  _et al._ (2010) E. Gull, M. Ferrero, O. Parcollet, A. Georges, and A. J. Millis, Phys. Rev. B 82, 155101 (2010).
  * Sordi  _et al._ (2010) G. Sordi, K. Haule, and A. Tremblay, Phys. Rev. Lett. 104, 226402 (2010).
  * Sordi  _et al._ (2012) G. Sordi, P. Sémon, K. Haule, and A. Tremblay, Phys. Rev. Lett. 108, 216401 (2012).
  * Sakai  _et al._ (2016b) S. Sakai, M. Civelli, and M. Imada, Physical Review B 94, 115130 (2016b).
  * Sakai  _et al._ (2018) S. Sakai, M. Civelli, and M. Imada, Physical Review B 98, 195109 (2018).
  * Scheurer  _et al._ (2018) M. S. Scheurer, S. Chatterjee, W. Wu, M. Ferrero, A. Georges, and S. Sachdev, Proceedings of the National Academy of Sciences 115, E3665 (2018), https://www.pnas.org/content/115/16/E3665.full.pdf .
  * Doiron-Leyraud  _et al._ (2007) N. Doiron-Leyraud, C. Proust, D. LeBoeuf, J. Levallois, J.-B. Bonnemaison, R. Liang, D. A. Bonn, W. N. Hardy, and L. Taillefer, Nature 447, 565 (2007).
  * Sebastian  _et al._ (2012) S. E. Sebastian, N. Harrison, R. Liang, D. A. Bonn, W. N. Hardy, C. H. Mielke, and G. G. Lonzarich, Phys. Rev. Lett. 108, 196403 (2012).
  * Chang  _et al._ (2012) J. Chang, E. Blackburn, A. T. Holmes, N. B. Christensen, J. Larsen, J. Mesot, R. Liang, D. A. Bonn, W. N. Hardy, A. Watenphul, M. v. Zimmermann, E. M. Forgan, and S. M. Hayden, Nat. Phys. 8, 871 (2012).
  * Blanco-Canosa  _et al._ (2013) S. Blanco-Canosa, A. Frano, T. Loew, Y. Lu, J. Porras, G. Ghiringhelli, M. Minola, C. Mazzoli, L. Braicovich, E. Schierle, E. Weschke, M. Le Tacon, and B. Keimer, Phys. Rev. Lett. 110, 187001 (2013).
  * Blackburn  _et al._ (2013) E. Blackburn, J. Chang, M. Hücker, A. T. Holmes, N. B. Christensen, R. Liang, D. A. Bonn, W. N. Hardy, U. Rütt, O. Gutowski, M. v. Zimmermann, E. M. Forgan, and S. M. Hayden, Phys. Rev. Lett. 110, 137004 (2013).
  * Ghiringhelli  _et al._ (2012) G. Ghiringhelli, M. Le Tacon, M. Minola, S. Blanco-Canosa, C. Mazzoli, N. B. Brookes, G. M. De Luca, A. Frano, D. G. Hawthorn, F. He, T. Loew, M. M. Sala, D. C. Peets, M. Salluzzo, E. Schierle, R. Sutarto, G. A. Sawatzky, E. Weschke, B. Keimer, and L. Braicovich, Science 337, 821 (2012).
  * Gerber  _et al._ (2015) S. Gerber, H. Jang, H. Nojiri, S. Matsuzawa, H. Yasumura, D. A. Bonn, R. Liang, W. N. Hardy, Z. Islam, A. Mehta, S. Song, M. Sikorski, D. Stefanescu, Y. Feng, S. A. Kivelson, T. P. Devereaux, Z.-X. Shen, C. C. Kao, W. S. Lee, D. Zhu, and J. S. Lee, Science 350, 949 (2015).
  * Chang  _et al._ (2016a) J. Chang, E. Blackburn, O. Ivashko, A. T. Holmes, N. B. Christensen, M. Huecker, R. Liang, D. A. Bonn, W. N. Hardy, U. Ruett, M. V. Zimmermann, E. M. Forgan, and S. M. Hayden, Nat. Commun. 7 (2016a), 10.1038/ncomms11494.
  * Wu  _et al._ (2011) T. Wu, H. Mayaffre, S. Krämer, M. Horvatic, C. Berthier, W. N. Hardy, R. Liang, D. A. Bonn, and M.-H. Julien, Nature 477, 191 (2011).
  * Wu  _et al._ (2015) T. Wu, H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, W. N. Hardy, R. Liang, D. A. Bonn, and M.-H. Julien, Nature Communications 6, 6438 (2015).
  * Julien (2015) M.-H. Julien, Science 350, 914 (2015).
  * Loret  _et al._ (2019) B. Loret, N. Auvray, Y. Gallais, M. Cazayous, A. Forget, D. Colson, M.-H. Julien, I. Paul, M. Civelli, and A. Sacuto, Nature Physics , 1 (2019).
  * LeBoeuf  _et al._ (2013) D. LeBoeuf, S. Kramer, W. N. Hardy, R. Liang, D. A. Bonn, and C. Proust, Nat. Phys. 9, 79 (2013).
  * Chakraborty  _et al._ (2019) D. Chakraborty, M. Grandadam, M. H. Hamidian, J. C. S. Davis, Y. Sidis, and C. Pépin, Phys. Rev. B 100, 224511 (2019).
  * Grandadam  _et al._ (2020a) M. Grandadam, D. Chakraborty, and C. Pépin, J. Supercond. Nov. Magn. 33, 2361 (2020a).
  * Grandadam  _et al._ (2020b) M. Grandadam, D. Chakraborty, X. Montiel, and C. Pépin, Phys. Rev. B 102, 121104 (2020b).
  * He  _et al._ (2011) R.-H. He, M. Hashimoto, H. Karapetyan, J. D. Koralek, J. P. Hinton, J. P. Testaud, V. Nathan, Y. Yoshida, H. Yao, K. Tanaka, W. Meevasana, R. G. Moore, D. H. Lu, S.-K. Mo, M. Ishikado, H. Eisaki, Z. Hussain, T. P. Devereaux, S. A. Kivelson, J. Orenstein, A. Kapitulnik, and Z.-X. Shen, Science 331, 1579 (2011).
  * Haule and Kotliar (2007) K. Haule and G. Kotliar, Phys. Rev. B 76, 104509 (2007).
  * Gull  _et al._ (2013) E. Gull, O. Parcollet, and A. J. Millis, Phys. Rev. Lett. 110, 216405 (2013).
  * Campi  _et al._ (2015) G. Campi, A. Bianconi, N. Poccia, G. Bianconi, L. Barba, G. Arrighetti, D. Innocenti, J. Karpinski, N. D. Zhigadlo, S. M. Kazakov, M. Burghammer, M. V. Zimmermann, M. Sprung, and A. Ricci, Nature (London) 525, 359 (2015), arXiv:1509.05002 [cond-mat.supr-con] .
  * Chang  _et al._ (2016b) J. Chang, E. Blackburn, O. Ivashko, A. T. Holmes, N. B. Christensen, M. Hucker, R. Liang, D. A. Bonn, W. N. Hardy, U. Rutt, M. v. Zimmermann, E. M. Forgan, and H. S. M., Nat. Commun. 7, 11494 (2016b).
  * Wu  _et al._ (2013) T. Wu, H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, P. L. Kuhns, A. P. Reyes, R. Liang, W. N. Hardy, D. A. Bonn, and M.-H. Julien, Nat. Commun. 4, 2113 (2013).
  * Sarkar  _et al._ (2020) S. Sarkar, M. Grandadam, and C. Pépin, (2020), arXiv:2009.02975 [cond-mat.supr-con] .
  * Baskaran and Anderson (1988) G. Baskaran and P. W. Anderson, Phys. Rev. B 37, 580 (1988).
  * Lee and Nagaosa (1992) P. A. Lee and N. Nagaosa, Phys. Rev. B 46, 5621 (1992).
  * Yang  _et al._ (2006) K.-Y. Yang, T. M. Rice, and F.-C. Zhang, Phys. Rev. B 73, 174501 (2006).
  * Le Tacon  _et al._ (2006) M. Le Tacon, A. Sacuto, A. Georges, G. Kotliar, Y. Gallais, D. Colson, and A. Forget, Nat. Phys. 2, 537 (2006).
  * Lee (2014) P. A. Lee, Phys. Rev. X 4, 031017 (2014).


# Learning Optimal Fair Decision Trees

Nathanael Jo  USC Center for AI in SocietyPalo AltoUSA<EMAIL_ADDRESS>,
Sina Aghaei<EMAIL_ADDRESS>USC Center for AI in Society3715 McClintock Ave
GER 240Los AngelesCAUSA90089 ,  Jack Benson<EMAIL_ADDRESS>University of
Southern California3715 McClintock Ave GER 240Los AngelesCAUSA90089 ,  Andrés
Gómez<EMAIL_ADDRESS>University of Southern California3715 McClintock Ave
GER 240Los AngelesCAUSA90089 and  Phebe Vayanos<EMAIL_ADDRESS>USC
Center for AI in Society3715 McClintock Ave GER 240Los AngelesCAUSA90089

###### Abstract.

The increasing use of machine learning in high-stakes domains – where people’s
livelihoods are impacted – creates an urgent need for interpretable and fair
algorithms. In these settings it is also critical for such algorithms to be
accurate. With these needs in mind, we propose a mixed integer optimization
(MIO) framework for learning optimal classification trees of fixed depth that
can be conveniently augmented with arbitrary domain specific fairness
constraints. We benchmark our method against state-of-the-art approaches for
building fair trees on popular datasets; given a fixed discrimination
threshold, our approach improves out-of-sample (OOS) accuracy by 2 percentage
points on average and obtains a higher OOS accuracy on 100% of the
experiments. We also incorporate various algorithmic fairness notions into our
method, showcasing its versatile modeling power that allows decision makers to
fine-tune the trade-off between accuracy and fairness.

fairness, decision trees, mixed-integer optimization

††copyright: rightsretained††conference: EAAMO ’22: ACM Conference on Equity
and Access in Algorithms, Mechanisms, and Optimization; October 6–9, 2022;
Washington DC, USA††booktitle: EAAMO ’22: ACM Conference on Equity and Access
in Algorithms, Mechanisms, and Optimization. October 6–9, 2022. Washington DC,
USA

##  1\. Introduction

There is growing interest in using machine learning (ML) to make decisions in
high-stakes domains. For instance, ML algorithms are now commonly used to
determine a criminal’s risk of recidivism in the United States (Angwin et al.,
2016). There is also a growing literature in designing algorithms to determine
the best course of action for homeless individuals (Azizi et al., 2018), to
diagnose and treat of various illnesses (Fatima et al., 2017), and many more.
In these contexts, it is necessary for such models to be both accurate (in
order to minimize erroneous predictions that negatively affect stakeholders)
and interpetable (so that decisions are transparent and hence accountable). We
therefore focus our attention to the problem of learning optimal
classification trees. A classification tree takes the form of a binary tree.
In each branching node of the tree, a binary test is performed on a specific
feature. Two branches emanate from each branching node, with each branch
representing the outcome of the test. If a datapoint passes (resp. fails) the
test, it is directed to the left (resp. right) branch. A predicted label is
assigned to all leaf nodes. Thus, each path from root to leaf represents a
classification rule that assigns a unique label to all datapoints that reach
that leaf. Classification trees are among the most interpretable of models
(Rudin, 2019), and optimal trees – rather than ones that are built using
heuristics – maximize predictive accuracy. They belong to a broader set of
models known as decision trees (the other type of decision tree being
regression trees, which apply to datasets with real valued labels). However,
for convenience, we will use the terms “classification trees” and “decision
trees” interchangeably throughout this paper.

A crucial consideration in learning optimal decision trees for use in high-
impact situations is fairness. After all, an algorithm that affects people’s
well-being should be aware of the particular historical and/or social contexts
that surround the learning problem. However, what constitutes as fair may vary
widely depending on one’s goals, domains of interest, etc. This discussion
surrounding “algorithmic fairness” consists of a rich literature that bridges
philosophy, economics, and computer science. Consider, for example, a naïve
approach where a classifier learns from data where the sensitive attributes
are omitted (such as race, sex, and others). Also known as “fairness through
unawareness”, this approach often leads to discrimination because some other
non-sensitive attribute(s) may be correlated to the withheld features (e.g.,
race is inextricably linked to ZIP code and income) (Pedreshi et al., 2008).

Another approach – and one that this paper adopts – is the notion of “group”
statistical fairness. A classifier satisfying group fairness is now aware of
the sensitive features in the data but will only produce decisions that
enforce parity between segments of the population. This is in contrast to
“individual” fairness, which requires that individuals with similar
characteristics be classified similarly (Dwork et al., 2012). Both
perspectives have their advantages and drawbacks, but in this work we will
focus on group fairness primarily because of its simplicity and
interpreability. In the following sections, we will discuss the use of many
notions of group fairness that can sufficiently satisfy more nuanced and fine-
grained fairness considerations.

###  1.1. Problem Statement

We now formalize the problem we study. Let
$\mathcal{D}:=\\{(x^{i},y^{i})\\}_{i\in\mathcal{I}}$ be training data indexed
in the set $\mathcal{I}\in\\{1,\dots,I\\}$. Each datapoint $i\in\mathcal{I}$
consists of a vector $x^{i}$ with $F$ binary features (i.e.,
$x^{i}\in\\{0,1\\}^{F}$), and a class $y^{i}\in\mathcal{K}$, where
$\mathcal{K}$ is the finite set of possible classes. Throughout this paper, we
will consider the case of binary classes, i.e., $\mathcal{K}:=\\{0,1\\}$ where
$y=1$ will be referred to as the positive class. Note that the literature in
fairness for multi-class learning is limited, see e.g. (Denis et al., 2021).
The goal is to learn, over all possible trees of maximum depth $d$, the tree
that maps $x^{i}$ to $y^{i}$ and maximizes out-of-sample accuracy, using in-
sample performance as a proxy. Further, suppose the population can be divided
into different sensitive groups (whereby discrimination exists or is a concern
within the classification problem, e.g. race, gender). Let $\mathcal{P}$
denote levels of a sensitive attribute, and let each datapoint $i$ have a
value $p^{i}\in\mathcal{P}$. The features in $\mathcal{P}$ may or may not be
included in the vector $x^{i}$ (Chen et al., 2019; Gajane and Pechenizkiy,
2017).

###  1.2. Common Notions of Group Fairness

In this section, we formally define five common notions of group fairness in
the machine learning literature that we will use in this work. Let $X$, $Y$,
and $P$ be random variables taken from an unknown distribution $\mathbb{P}$,
and let $\hat{Y}$ be the classifier’s prediction.

#### Statistical Parity

A classifier satisfies statistical parity if the probability of receiving a
positive class is equal across all protected groups (Dwork et al., 2012). For
instance, consider a company that wants to make its automated hiring algorithm
more fair; they may require that males and females be hired at equal rates.
Formally, statistical parity means

$\mathbb{P}[\hat{Y}=1|P=p]=\mathbb{P}[\hat{Y}=1|P=p^{\prime}]\quad\forall
p,p^{\prime}\in\mathcal{P}.$

#### Conditional Statistical Parity

Similarly, a classifier satisfies conditional statistical parity if the
probability of receiving a positive class is equal across all protected
groups, conditional on some legitimate feature(s) indicative of risk (Corbett-
Davies et al., 2017). For example, the authors of (Corbett-Davies et al.,
2017) consider the use of algorithms to help decide whether or not defendants
should be detained while awaiting trial. They stipulate that the detention
rate should be equal between Black and White defendants who have the same
number of prior convictions. This may be considered as a fairer notion to
statistical parity because it takes into account the distribution of risk
factors within each sensitive group. Under statistical parity, it is possible
for a classifier to detain more Black individuals with few prior convictions
than White individuals with few prior convictions as long as the detention
rate between the two groups are equal. However, one may argue that this
decision rule is not fair as it disproportionately detains Black people who do
not deserve to be detained based on their criminal history. Letting $L$
represent the random variable taken from a set $\mathcal{L}$ of legitimate
features, conditional statistical parity is satisfied if

$\mathbb{P}[\hat{Y}=1|P=p,L=\ell]=\mathbb{P}[\hat{Y}=1|P=p^{\prime},L=\ell]\quad\forall
p,p^{\prime}\in\mathcal{P},\ell\in\mathcal{L}.$

#### Predictive Equality

A classifier satisfies predictive equality if all protected groups have the
same false positive rates (FPR) (Chouldechova, 2017). For instance, consider
the use of machine learning in predicting good or bad credit scores, with
gender as the protected group. Predictive equality requires that the rate in
which we predict good credit when someone in reality has bad credit is equal
among males and females. Formally, this means

$\mathbb{P}[\hat{Y}=1|P=p,Y=0]=\mathbb{P}[\hat{Y}=1|P=p^{\prime},Y=0]\quad\forall
p,p^{\prime}\in\mathcal{P}.$

#### Equal Opportunity

A classifier satisfies equal opportunity if all protected groups have the same
true positive rate (TPR) (Hardt et al., 2016). In the credit prediction
example mentioned above, we may require that males and females are correctly
classified with good credit at the same rate. The formal definition of equal
opportunity is as follows:

$\mathbb{P}[\hat{Y}=1|P=p,Y=1]=\mathbb{P}[\hat{Y}=1|P=p^{\prime},Y=1]\quad\forall
p,p^{\prime}\in\mathcal{P}.$

Ultimately the decision to prioritize TPR or FPR depends on the learning
problem. If we only enforce equal opportunity in the credit prediction
example, then a classifier may erroneously assign males good credit at a much
higher rate than females; while the TPR may be equal between males and
females, the overall decision rule still favors males. The reverse argument
can also be made when only enforcing predictive equality.

#### Equalized Odds

Equalized odds combines predictive equality and equal opportunity such that
both FPR and TPR must be similar across all protected groups (Hardt et al.,
2016). For example, in the college admissions process, one may require that
the rate of admission for those who are qualified and not qualified are the
same among Black and White applicants. Equalized odds is, of course, a
stronger condition than predictive equality and equal opportunity
individually, and bypasses the need to choose which fairness metric to
prioritize. However, it may be computationally harder to solve, if there even
exists a feasible solution at all. Formally, equalized odds is satisfied if

$\mathbb{P}[\hat{Y}=1|P=p,Y=y]=\mathbb{P}[\hat{Y}=1|P=p^{\prime},Y=y]\quad\forall
p,p^{\prime}\in\mathcal{P},y\in\\{0,1\\}.$

###  1.3. Related Works

In this section we briefly review related works in the machine learning
literature.

Discrimination Prevention in Machine Learning. The literature on fairness in
machine learning is extensive – we point the interested reader to (Mehrabi et
al., 2021) and (Žliobaitė, 2017) for a survey of this topic. In general, there
are three approaches that existing works employ in order to promote fairness.
The first is pre-processing, which entails eliminating discrimination within
the training data, see (Zemel et al., 2013; Feldman et al., 2015; Hajian and
Domingo-Ferrer, 2012; Kamiran et al., 2013), among others. The second is post-
processing, which takes the model’s output and alters its decisions in ways
that promote some fairness metric, see (Kamiran et al., 2010; Hardt et al.,
2016; Pleiss et al., 2017). Finally, one may also employ in-processing by
modifying existing models so that fairness is integrated within its learning
goal; for instance, approaches have utilized regularization (Berk et al.,
2017; Kamishima et al., 2012; Aghaei et al., 2019) or different heuristics
(Kamiran et al., 2010; Calders and Verwer, 2010). Our work falls under the
category of in-processing discrimination prevention techniques.

Fair Decision Trees. Within the stream of literature in discrimination
prevention, our work most closely relates to approaches for learning fair
decision trees. For instance, Ranzato et al. (Ranzato et al., 2021) adapt a
genetic algorithm for learning robust decision trees in order to prioritize
individual fairness. In the online setting, Zhang et al. (Zhang and Ntoutsi,
2019) propose a fairness-aware Hoeffding tree that builds decision trees over
streams of data. Grari et al. (Grari et al., 2019) take an adversarial
approach to training fair gradient-boosted trees that promote statistical
parity. Kanamori and Arimura (Kanamori and Arimura, 2019) propose a post-
processing step that uses MIO to edit the branching thresholds of the tree’s
internal nodes to satisfy some fairness constraint; in the paper they use
statistical parity and equalized odds. In a similar manner, Zhang et al.
(Zhang et al., 2020) flip the outcomes of different paths of a learned tree in
order to improve fairness. Closely related to our work is that of Aghaei et
al. (Aghaei et al., 2019), who propose a mixed-integer optimization framework
to build optimal and fair decision trees that prevent disparate impact and/or
disparate treatment. Finally, Kamiran et al. (Kamiran et al., 2010) propose a
two-pronged approach: the first is to incorporate sensitivity gain (IGS) with
information gain (IGC) as a heuristic to build a fairer decision tree, while
the second is to relabel the predictions in the tree’s leaf nodes to further
promote statistical parity.

Optimal Decision Trees. Mixed-integer optimization (MIO) has recently gained
traction as a framework to solve various machine learning problems.
Particularly related to this paper are approaches that use MIO to learn
optimal decision trees in order to improve on the traditional
classification/regression tree (CART) algorithms (Breiman et al., 1984), which
rely on a heuristic. Many works introduce novel formulations for learning
optimal classification trees (Aghaei et al., 2021, 2019; Verwer and Zhang,
2019; Bertsimas and Dunn, 2017). Elmachtoub et al. (Elmachtoub et al., 2020)
use MIO to learn decision trees that minimize a loss function derived from the
“predict-then-optimize” framework. There exist several extensions to learning
decision trees using MIO as well. For instance, Mišić (Mišić, 2020) uses MIO
to solve tree ensemble models to optimality; Jo et al. (Jo et al., 2021) to
learn optimal prescriptive trees from observational data; and Justin et al.
(Justin et al., 2021) and Bertsimas et al. (Bertsimas et al., 2019) to learn
optimal robust decision trees. In this paper, we build on the MIO method
introduced by Aghaei et al. (Aghaei et al., 2021) by conducting extensive
experiments when we add various fairness constraints to the model.

###  1.4. Contributions

In this work, we present an MIO formulation to learn optimal decision trees
initially proposed in (Aghaei et al., 2021), and showcase its flexible
modeling power in considering various notions of fairness. We conduct a
comprehensive evaluation of its performance on several popular datasets,
incorporating various fairness constraints where applicable. We also compare
our method to one of the most popular (heuristic-based) algorithms for
learning fair trees proposed in (Kamiran et al., 2010) and show that our
optimal tree does, in fact, lead to significant performance improvements.
Given a fixed discrimination threshold, our approach improves out-of-sample
accuracy by 2.3 percentage points on average and obtains a higher out-of-
sample accuracy on 88.9% of our experiments. The code to learn fair trees with
arbitrary constraints is made publicly available.

##  2\. Formulation

In this section, we present the MIO formulation to learn optimal decision
trees proposed in (Aghaei et al., 2021). The original paper introduced the
formulation without an emphasis on fairness, but in this work we will be using
the formulation as a building block to which we add various fairness
constraints. From hereon, we will refer to the combination of the MIO
formulation and the fairness constraints as FairOCT.

###  2.1. From Decision Tree to Flow Graph

Before presenting the formulation for FairOCT, we introduce the modeling
framework we use for decision trees. The key idea is to convert a decision
tree into a directed, acyclic graph where all arcs “flow” from the tree’s root
down to its leaves. Formally, we start with a perfect binary tree of depth
$d$, whose nodes are labeled $1$ through $(2^{d+1}-1)$ in the order in which
they appear in a breadth-first search. Let
$\mathcal{B}:=\\{1,\dots,2^{d}-1\\}$ denote the set of branching nodes and
$\mathcal{T}:=\\{2^{d},\dots,2^{d+1}-1\\}$ the set of terminal nodes. We then
convert all arcs in the tree to point from the parent node to its child node
(see Figure 1, left). From the binary tree, we connect a source $s$ to the
root and all nodes $n\in\mathcal{B}\cup\mathcal{T}$ to sinks $t_{k}$, one for
every class $k\in\mathcal{K}$ (see Figure 1, right). We will again assume that
we have binary classes, $\mathcal{K}=\\{0,1\\}$, although this formulation
generalizes to arbitrary finite classes. All arcs have a capacity of 1, so
each datapoint is weighted equally and flows from the source $s$ to one sink
$t_{k}$, where $k$ is the class that the decision tree assigns to that
datapoint. From hereon, we will refer to this structure as a “flow graph”.

Figure 1. A decision tree of depth 2 (left) and its associated flow graph
(right). In this example, the flow graph contains 2 sinks to reflect 2
classes.

###  2.2. MIO Formulation

With the flow graph at hand, we present the MIO formulation for learning
optimal classification trees. We encode the branching structure of the tree
with variables $b_{nf}\in\\{0,1\\}$, for all
$n\in\mathcal{B},f\in\mathcal{F}$, which indicate if feature $f$ is selected
for branching at node $n$. We also encode the prediction scheme of the tree
using variables $w_{nk}\in\\{0,1\\}$, for all
$n\in\mathcal{B}\cup\mathcal{T},k\in\mathcal{K}$, which equal 1 if and only if
(iff) node $n$ assigns class $k$ to all datapoints that land on that node.
Accordingly, we let $p_{n}\in\\{0,1\\}$ indicate if node $n$ is a prediction
node, in which case it must assign a class to all datapoints that land on that
node and no further branching is allowed. We also define “flow variables” $z$
to capture the flow of datapoints, where $z^{i}_{a(n),n}\in\\{0,1\\}$ equals 1
iff datapoint $i$ flows from the direct ancestor of node n, $a(n)$, to $n$.

The formulation is as follows:

(1a) maximize
$\displaystyle\displaystyle\sum_{i\in\mathcal{I}}\sum_{n\in\mathcal{B}\cup\mathcal{T}}z^{i}_{n,t_{y^{i}}}$
(1b) subject to
$\displaystyle\displaystyle\sum_{f\in\mathcal{F}}b_{nf}+p_{n}+\sum_{m\in\mathcal{A}(n)}p_{m}=1$
$\displaystyle\forall n\in\mathcal{B}$ (1c) $\displaystyle
p_{n}+\sum_{m\in\mathcal{A}(n)}p_{m}=1$ $\displaystyle\forall n\in\mathcal{T}$
(1d) $\displaystyle\displaystyle
z^{i}_{a(n),n}=z^{i}_{n,\ell(n)}+z^{i}_{n,r(n)}+\sum_{k\in\mathcal{K}}z^{i}_{n,t_{k}}$
$\displaystyle\forall n\in\mathcal{B},i\in\mathcal{I}$ (1e)
$\displaystyle\displaystyle
z^{i}_{a(n),n}=\sum_{k\in\mathcal{K}}z^{i}_{n,t_{k}}$ $\displaystyle\forall
i\in\mathcal{I},n\in\mathcal{T}$ (1f) $\displaystyle\displaystyle
z^{i}_{s,1}\leq 1$ $\displaystyle\forall i\in\mathcal{I}$ (1g)
$\displaystyle\displaystyle
z^{i}_{n,\ell(n)}\leq\sum_{f\in\mathcal{F}:x_{f}^{i}=0}b_{nf}$
$\displaystyle\forall n\in\mathcal{B},i\in\mathcal{I}$ (1h)
$\displaystyle\displaystyle
z^{i}_{n,r(n)}\leq\sum_{f\in\mathcal{F}:x_{f}^{i}=1}b_{nf}$
$\displaystyle\forall n\in\mathcal{B},i\in\mathcal{I}$ (1i)
$\displaystyle\displaystyle z^{i}_{n,t_{k}}\leq w^{n}_{k}$
$\displaystyle\forall
i\in\mathcal{I},n\in\mathcal{B}\cup\mathcal{T},k\in\mathcal{K}$ (1j)
$\displaystyle\displaystyle\sum_{k\in\mathcal{K}}w^{n}_{k}=p_{n}$
$\displaystyle\forall n\in\mathcal{B}\cup\mathcal{T}$ (1k)
$\displaystyle\displaystyle w^{n}_{k}\in\\{0,1\\}$ $\displaystyle\forall
n\in\mathcal{B}\cup\mathcal{T},k\in\mathcal{K}$ (1l)
$\displaystyle\displaystyle b_{nf}\in\\{0,1\\}$ $\displaystyle\forall
n\in\mathcal{B},f\in\mathcal{F}$ (1m) $\displaystyle\displaystyle
p_{n}\in\\{0,1\\}$ $\displaystyle\forall n\in\mathcal{B}\cup\mathcal{T}$ (1n)
$\displaystyle\displaystyle z^{i}_{a(n),n},z^{i}_{n,t_{k}}\in\\{0,1\\}$
$\displaystyle\forall
n\in\mathcal{B}\cup\mathcal{T},i\in\mathcal{I},k\in\mathcal{K},$

where $\ell(n)$ (resp. $r(n)$) is the left (resp. right) descendant of $n$ and
$\mathcal{A}(n)$ is the set of all ancestors of node
$n\in\mathcal{B}\cup\mathcal{T}$. The objective (1a) maximizes the number of
correctly classified datapoints. Note that we can add a regularization term in
the objective to control for overfitting, see (Aghaei et al., 2021).
Constraints (1b) ensure that each node is either a branching node, a
prediction node, or neither of the two because one of its ancestors is already
a prediction node; the last option means that the node is pruned out.
Constraints (1c) similarly ensure that each terminal node either makes a
prediction or has an ancestor that is a prediction node. Constraints (1d) and
(1e) enable flow conservation, whereby every datapoint that flows into node
$n$ must exit to its left (or right) descendant, or flow directly to a sink
$t_{k}$ that corresponds to a class $k$. Constraints (1f) ensure that the flow
value of each datapoint entering source $s$ is at most 1. Constraints (1g) and
(1h) enforce datapoints to flow to a node’s left (resp. right) child if their
branching feature is 0 (resp. 1). Constraints (1i) require that a datapoint
must be directed to the sink corresponding to the class that the prediction
node assigns. Finally, constraints (1j) in conjunction with (1k) ensure that
if we make a prediction at node $n$, then exactly one class is associated with
its prediction. In sum, each datapoint flows into the graph and lands on
exactly one of the sinks via a path from source to sink depending on its
feature vector.

###  2.3. Fairness Constraints

One main advantage of formulation (1) is its flexible modeling power. In the
following, we showcase the variety of constraints that can be added to the
formulation to learn trees that are fair, and how we can tune the trade-off
between fairness and accuracy of the learned trees to satisfy the definitions
of fairness we introduced in Section 1.2. We previously defined the various
notions of fairness using strict equalities; in practice, however, we may
relax this condition by introducing a bias $\delta$, where $\delta$ is the
maximum disparity allowed between groups.

#### Statistical Parity.

Statistical parity is satisfied up to a bias of $\delta$ when we add the
following constraint to (1):

(2)
$\begin{split}\displaystyle\left|\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}:p^{i}=p}z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p\\}|}-\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}:p^{i}=p^{\prime}}z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p^{\prime}\\}|}\right|\leq\delta\quad\forall
p,p^{\prime}\in\mathcal{P}:p\neq p^{\prime},\end{split}$

where the left-hand side of the inequality is the absolute difference between
the proportion of positive classes assigned by the classification tree to
groups $p$ and $p^{\prime}$, respectively, in the training data.

#### Conditional Statistical Parity.

In order to describe the CSP constraint, we will let $\ell^{i}$ denote the
value of datapoint $i$’s legitimate factor(s). To ensure that the learned tree
satisfies conditional statistical parity up to a bias $\delta$ and for all
$\ell\in\mathcal{L}$, we may augment (1) with the constraint

(3)
$\begin{split}\displaystyle\left|\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}}\mathbb{I}(p^{i}=p\land\ell^{i}=\ell)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p\land\ell^{i}=\ell\\}|}-\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}}\mathbb{I}(p^{i}=p^{\prime}\land\ell^{i}=\ell)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p^{\prime}\land\ell^{i}=\ell\\}|}\right|\leq\delta\quad\forall
p,p^{\prime}\in\mathcal{P}:p\neq p^{\prime},\;\ell\in\mathcal{L},\end{split}$

where the left-hand side of the inequality is the absolute difference between
the proportions of positive classes assigned by the classification tree to
groups $p$ and $p^{\prime}$ respectively, in the training data, and whose
legitimate feature(s) equal(s) $\ell$.

#### Predictive Equality.

The following constraint may be added to (1) in order to satisfy predictive
equality up to a bias $\delta$,

(4)
$\begin{split}\left|\frac{\displaystyle\sum_{i\in\mathcal{I}}\sum_{n\in\mathcal{B}\cup\mathcal{T}}\mathbb{I}(p^{i}=p\land
y^{i}=0)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p\land
y^{i}=0\\}|}-\frac{\displaystyle\sum_{i\in\mathcal{I}}\sum_{n\in\mathcal{B}\cup\mathcal{T}}\mathbb{I}(p^{i}=p^{\prime}\land
y^{i}=0)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p^{\prime}\land
y^{i}=0\\}|}\right|\leq\delta\quad\forall p,p^{\prime}\in\mathcal{P}:p\neq
p^{\prime},\end{split}$

where the left-hand side of the inequality is the absolute difference between
the proportions of false positive assignments made by the classification tree
in groups $p$ and $p^{\prime}$, respectively, in the training data.

#### Equal Opportunity.

Similarly, the following constraint may be added to (1) in order to satisfy
equal opportunity up to a bias $\delta$,

(5)
$\begin{split}\left|\frac{\displaystyle\sum_{i\in\mathcal{I}}\sum_{n\in\mathcal{B}\cup\mathcal{T}}\mathbb{I}(p^{i}=p\land
y^{i}=1)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p\land
y^{i}=1\\}|}-\frac{\displaystyle\sum_{i\in\mathcal{I}}\sum_{n\in\mathcal{B}\cup\mathcal{T}}\mathbb{I}(p^{i}=p^{\prime}\land
y^{i}=1)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p^{\prime}\land
y^{i}=1\\}|}\right|\leq\delta\quad\forall p,p^{\prime}\in\mathcal{P}:p\neq
p^{\prime}.\end{split}$

The above constraint is similar to (4), with the exception that we take the
absolute difference between the proportions of true positive assignments made
by the tree in the training data.

#### Equalized Odds.

We may augment (1) to satisfy equalized odds up to a bias $\delta$ using the
constraint

(6)
$\begin{split}\displaystyle\left|\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}}\mathbb{I}(p^{i}=p\land
y^{i}=k)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p\land
y^{i}=k\\}|}-\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}}\mathbb{I}(p^{i}=p^{\prime}\land
y^{i}=k)z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p^{\prime}\land
y^{i}=k\\}|}\right|\leq\delta\quad\forall
k\in\mathcal{K},p,p^{\prime}\in\mathcal{P}:p\neq p^{\prime}.\end{split}$

The above constraint is equivalent to enforcing (4) and (5) simultaneously,
with $\mathcal{K}=\\{0,1\\}$.

Note that any of the above constraints can be easily linearized by decomposing
them into two. In the general case where we have constraint
$|f(x)|\leq\delta$, we can reformulate it into $f(x)\leq\delta$ and
$-f(x)\leq\delta$. If $f$ is affine (as is the case here), then these
constraints are linear in the decision variables of the problem and the
resulting problem can be solved with powerful off-the-shelf mixed-integer
_linear_ optimization solvers such as Gurobi111See
https://www.gurobi.com/products/gurobi-optimizer/. We also emphasize that the
above constraints are illustrative examples of the approach’s flexibility, and
that they may be amended or combined with other fairness considerations. For
example, instead of restricting the absolute value of the difference between
groups to be bounded by $\delta$, we may impose that the minority group should
be better off by a margin $\delta$, see Section 3.3.

###  2.4. Greedy Implementation

Since problem (1) is NP-hard to solve, any algorithm that solves it to
optimality, including using MIO technology, requires (in the worst case) a
runtime that grows exponentially with the depth of the decision tree of the
problem. In practice, in several high-stakes domains (including most of the
example discussed in Section 2.3), interpretable (and thus shallow) decision
trees are needed, and ample time is available for training; in some cases,
runtimes in the order of hours are adequate. As a consequence, large runtimes
are often not a concern, if it results in better or fairer decisions.

Nonetheless, in some cases (e.g., if the decision tree is meant to be used as
part of a black-box ensemble method), solving (1) to optimality is
impractical. For such situations, we propose to use MIO algorithms as a
heuristic, as follows:

  1. (1)

Solve Problem (1) for an initial depth $k$

  2. (2)

Solve Problem (1) for depth $k+1$, fixing the variables $b$ to the values
obtained from the depth $k$ problem. Increase $k$ and repeat until a desired
depth $d$ is reached.

In Section 3 we show that this greedy MIO-based approach outperforms other
heuristic methods proposed in the literature.

##  3\. Experiments

We now evaluate the empirical performance of the model outlined in Section 2.
We first compare its performance to a suite of methods for learning
discrimination-aware decision trees proposed by (Kamiran et al., 2010); from
hereon, we will refer to the paper’s general model as DADT, although it
comprises of a variety of heuristic models. Critically, DADT only allows for
two sensitive groups and only considers statistical parity as its fairness
metric. We therefore conduct further experiments on the proposed model where
more than two sensitive groups exist and other fairness definitions are
considered.

Note that methods other than DADT in the literature either: incorporate
different fairness metrics, which would not yield a fair comparison (e.g.,
Aghaei et al. (Aghaei et al., 2019)); consider individual fairness as opposed
to group fairness (e.g., Ranzato et al. (Ranzato et al., 2021)); learn trees
in an online setting (e.g., Zhang et al. (Zhang and Ntoutsi, 2019)); or are
not interpretable (e.g., Grari et al. (Grari et al., 2019)).

###  3.1. Datasets

COMPAS. The Correctional Offender Management Profiling for Alternative
Sanctions (COMPAS) is a popular algorithm used by judges to predict a
criminal’s risk of recidivism after 2 years. Although this algorithm was meant
to eliminate human bias, ProPublica published a seminal article analyzing that
COMPAS is nonetheless biased in favor of White criminals (Angwin et al.,
2016). The data has since become a popular baseline for fairness models to
improve upon. There are around 6,000 datapoints, and we let race be the
sensitive attribute. The original race attribute has 6 levels: White, Black,
Asian, Hispanic, Native American, and “Other”. However, Asian and Native
American criminals are quite rare in the data, which means we cannot obtain
good estimates of fairness metrics for these two groups. We therefore group
them under the “Other” category to avoid this estimation issue.

Adult. The UCI Adult dataset contains around 45,000 entries from a 1994 Census
database (Dua and Graff, 2017). The goal is to predict whether or not
someone’s income exceeds $50,000 per year, and we treat sex as the sensitive
attribute in accordance with (Agarwal et al., 2018), (Grari et al., 2019), and
(Kamishima et al., 2012), among others – females being the marginalized group.

German. The German dataset comprises of 20 attributes and classifies people as
having good or bad credit. We use age as the sensitive attribute, and we
follow the works of (Kamiran et al., 2010) and (Feldman et al., 2015) to split
the population into people 25 and younger and older than 25, with the former
as the marginalized group since younger people are often assigned worse credit
under the basis of age.

We use 75% of the data for training and 25% for testing on all the above
datasets, split five times with different seeds to account for random
variability. The Adult dataset is both large and high-dimensional, so we
additionally train FlowOCT on 2,700 datapoints due to computational
limitations; we will further discuss this choice in Section 3.6.

###  3.2. Description of DADT

As briefly mentioned in Section 1.3, DADT consists of a two-pronged approach.
The first is an in-processing step that uses sensitivity gain (IGS) and
information gain (IGC) as heuristics to build a fair decision tree. If the
tree still has a discrimination level that is greater than some bias $\delta$,
a post-processing step can be used that relabels the predictions of the tree’s
leaf nodes until the discrimination reaches $\delta$. Note that the original
paper (Kamiran et al., 2010) proposes to use a combination of these steps
where appropriate; for example, one may opt to use different heuristic
combinations (e.g. IGC+IGS, IGC-IGS) or not use IGS entirely (which is
equivalent to the CART (Breiman et al., 1984) algorithm) and just execute the
relabeling step. We compare our method with respect to all three heuristics,
which we will refer to as IGC+IGS_Relab, IGC-IGS_Relab, and IGC_Relab,
respectively.

DADT has two modeling constraints: it only considers statistical parity as its
fairness metric, and it assumes only two sensitive groups exist. In further
contrast to FairOCT, DADT does not take the absolute difference in the
probabilities of receiving a positive outcome between sensitive groups,
instead positing that we subtract the probability of the marginalized group
from the dominant group. A “fair” result under this definition may include a
scenario in which the marginalized group receives better outcomes at a much
higher rate than the dominant group, which constraint (2) does not allow.

###  3.3. Methodology

#### Comparison to DADT (Fixed Depth)

In order to compare FairOCT fairly with DADT, we collapse the protected
variables in our datasets to two groups where applicable and only consider
statistical parity as our fairness measure. We also adapt our statistical
parity constraint (2) to fit the fairness definition as implemented by DADT:

$\begin{split}\displaystyle\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}:p^{i}=p_{d}}z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p_{d}\\}|}-\frac{\displaystyle\sum_{n\in\mathcal{B}\cup\mathcal{T}}\sum_{i\in\mathcal{I}:p^{i}=p_{m}}z^{i}_{n,t_{1}}}{|\\{i\in\mathcal{I}:p^{i}=p_{m}\\}|}\leq\delta,\end{split}$

where $p_{d}$ and $p_{m}$ refer to the dominant and marginalized groups,
respectively. The only dataset that requires further processing is COMPAS,
which has multiple levels for race. We will let the sensitive groups be White
versus non-White criminals, White being the dominant group. For both DADT and
FairOCT, we grow trees of depths $d\in\\{2,3\\}$ in order to produce
interpretable results.

#### Comparison to DADT (Tuning Depth)

We also compare DADT with the greedy implementation of Formulation 1 with
initial depth $k=2$ (see Section 2.4) on deeper trees. Recall that the
training set is 75% of the data; only for this set of experiments, we further
split this set so that 50% of the data is used to train and the other 25% is
used to cross-validate. We learn trees of depths $d\in\\{3,4,5,6\\}$ for both
DADT and FairOCT; for each discrimination level and for each sample, we choose
the tree with the highest average validation accuracy. Then, we average the
performance across all five samples.

#### Additional Experiments

In addition to our comparison with DADT, we conduct further experiments
incorporating other fairness definitions as they fit in the learning problem.
We also impose absolute differences as outlined in Section 2.3 and relax the
assumption that only two sensitive groups exist. For both the COMPAS and Adult
datasets, we augment problem (1) with constraints (2), (4), (5), and (6) in
separate experiments to incorporate statistical parity, predictive equality,
equal opportunity, and equalized odds, respectively. We also consider
conditional statistical parity – constraint (3) – on COMPAS, conditioned on
the number of prior crimes. Similar to our comparison with DADT, we learn
trees of maximum depth $d=2$.

For both sets of experiments, we vary the fairness bound
$\delta\in[0.01,0.55]$, increasing in $0.01$ increments. We choose an upper
bound of 0.55 because the results above this threshold stayed the same
consistently, but for completeness we also include the case where $\delta=1$
(which is equivalent to omitting the constraint entirely). All experiments
have a time limit of 3 hours.

###  3.4. Results: Comparison to DADT (Fixed Depth)

Figure 2. Comparison of the accuracy and discrimination of FairOCT and DADT
for trees of depth 2 (top) and 3 (bottom) when the discrimination bounds
($\delta$) are varied on three datasets: Adult (left column), COMPAS (middle
column), and German (right column) – averaged over 5 random train-test splits.
For each dataset, we display both in-sample (top row) and out-of-sample
(bottom row) performances. Figure 3. Comparison of a resulting decision tree
produced by FairOCT ($d=2$) and by IGC+IGS_Relab ($d=3$). Both trees result in
the same out-of-sample accuracy and discrimination.

Figure 2 shows the accuracy and discrimination for both FairOCT and DADT when
we vary the discrimination threshold $\delta$, on trees of depths 2 and 3.
Overall, for any given discrimination level, FairOCT consistently has better
in-sample and out-of-sample accuracies. This is expected since FairOCT finds
an optimal solution whereas DADT relies on a heuristic. Given a fixed
discrimination threshold, FairOCT improves out-of-sample (OOS) accuracy by
1.6, 2.8, 1.7 percentage points on average over IGC_Relab, IGC+IGS_Relab, and
IGC-IGS_Relab, respectively. Our method also obtains a higher OOS accuracy in
100% of the experiments compared to all the DADT methods. To further
illustrate our performance improvements, we show a sample decision tree from
FairOCT (depth 2) and IGC+IGS_Relab (depth 3) on the Adult dataset, both of
which result in the same out-of-sample accuracy and discrimination
performances (see Figure 3). FairOCT produces a tree of lower depth that is
far more interpretable than IGC+IGS_Relab, whose branching decisions might
seem circular and repetitive.

Since the Adult dataset is both large and high-dimensional, FairOCT struggled
to reach an optimal solution on the full training set within the time limit,
which resulted in suboptimal performances throughout. Therefore, the results
shown in Figure 2 are from another set of experiments where FairOCT learns
from 2,700 datapoints but is evaluated on the same test set to allow for a
fair comparison with DADT. The optimal solution in fact consistently
outperforms DADT (see Figure 2, left column). This means that, in spite of the
difficulties in solving MIO problems, opting for the optimal solution (even on
a smaller training set) may yield better results than relying on a heuristic.

One distinct advantage FairOCT has over DADT is the ability to fine-tune the
accuracy-discrimination tradeoff. While both methods were trained on the same
fairness bounds, FairOCT in general produced more distinct results, which
gives one more freedom to choose a policy that best suits their needs. This is
because our fairness requirement is ingrained within the optimization problem,
so a slight change in $\delta$ may change both the branching and labelling
decisions. DADT, in contrast, does not produce many distinct solutions because
the fairness requirement only comes in the relabelling step; there is no way
of controlling how fair the branching step will be, leading to more
homogeneous trees. This fact is most evident in IGC+IGS_Relab on the Adult
dataset, refer to the lone datapoint to the left of all the Adult graphs in
Figure 2. This result arises because the heuristic step reduces discrimination
so much that the relabelling step is never executed. While the performance
itself is decent, this static result is not desirable if we want to tune the
trade-off between discrimination and accuracy.

Figure 4. Distribution of solve times for FairOCT on the COMPAS, German, and
Adult datasets when trained on depths 2 and 3. Due to computational
limitations, we display solve times for experiments trained on 2,700
datapoints of Adult. The solve times for DADT are consistently below 1 second
for both COMPAS and German, and around 2 seconds for the full Adult training
set. Figure 5. Comparison of the accuracy and discrimination of FairOCT and
DADT when depth is tuned via cross-validation. Each datapoint is the best
performing tree in terms of validation accuracy for a fixed discrimination
threshold. The results of three datasets are displayed: Adult (left column),
COMPAS (middle column), and German (right column) performances. For each
dataset, we display both in-sample (top row) and out-of-sample (bottom row)
performances. Figure 6. Accuracy and discrimination of FairOCT on the COMPAS
dataset with varying fairness metrics (top, from left to right – statistical
parity, conditional statistical parity, predictive equality; bottom, from left
to right – equal opportunity, and equalized odds), averaged over 5 random
samples. Each datapoint in the graph corresponds to a different fairness bound
$\delta$ that is imposed in the learning problem. All figures show averages
over trees of maximum depth $d=2$.

Naturally, as we have observed with the full Adult training set, obtaining
provably optimal solutions comes at a computational cost. FairOCT has average
solve times of 1.1 hours for COMPAS (see Figure 4, left), 1.5 hours for German
(Figure 4, middle), and 2.5 hours for Adult (Figure 4, right) averaged over
treed of depths 2 and 3. On the other hand, DADT consistently finds a tree in
under 1 second for COMPAS and German, and around 3 seconds for Adult. DADT may
be better suited for settings where continuous recalibration of decisions is
necessary (e.g., in an active learning setting with streams of data). However,
apart from the ability fine-tune accuracy and discrimination as discussed
previously, FairOCT is also useful in high-impact domains because it maximizes
decision accuracy. For instance, if we desire full statistical parity between
White and non-White criminals in COMPAS, a difference in accuracy of 2.3%
results in a significant number of criminals who are correctly classified as
high or low risk to recidivate, which impacts their probation, sentence, and
potentially other aspects of their lives.

###  3.5. Results: Comparison to DADT (Tuning Depth)

Figure 5 shows the accuracy and discrimination for both FairOCT and DADT when
we fine-tune the depth via cross-validation. FairOCT performs noticeably
better out-of-sample (OOS) for both the Adult and German; given a fixed
discrimination threshold, FairOCT improves OOS accuracy by 1.3, 1.7, and 2.7
percentage points on average compared to IGC_Relab, IGC+IGS_Relab, IGC-
IGS_Relab, respectively. The performance improvements are even more pronounced
in low discrimination thresholds (see Figure 5, bottom left), which is often
more desirable in practice. On the COMPAS dataset, FairOCT’s performance seems
to be slightly better but the difference is by no means pronounced. Note that
we are using a heuristic adaptation of FairOCT in these sets of experiments,
so the performance improvements we are seeing are promising nonetheless.

###  3.6. Results: Additional Fairness Measures

Figure 7. Accuracy and discrimination of FairOCT trained on a subset of the
Adult dataset with varying fairness metrics top, from left to right –
statistical parity, predictive equality; bottom, from left to right – equal
opportunity, and equalized odds), averaged over 5 random samples. Each
datapoint in the graph corresponds to a different fairness bound $\delta$ that
is imposed in the learning problem. All figures show averages over trees of
maximum depth $d=2$.

We plot the accuracy and discrimination of FairOCT on the COMPAS and Adult
datasets in Figures 6 and 7, respectively. In general, the in-sample
trendlines are almost monotonically increasing, which is expected since there
is often a penalty in accuracy when we we want to promote more equality. The
in-sample trendlines also perform consistently better than the out-of-sample
lines for a given discrimination level. Critically, however, this difference
is quite minimal, which indicates that the trees generalize well and are not
overfitted to the training data. When we consider equalized odds on COMPAS,
for instance, the optimal trees average to around 55% in-sample accuracy at
perfect in-sample equality; the out-of-sample performance is about the same.
The trade-off in accuracy and discrimination is quite stark considering the
optimal trees at 25% in-sample discrimination has around 67% accuracy.
However, this range gives decision makers the freedom to fine-tune the
accuracy and discrimination trade-off to fit their needs. We observe the same
trend for all other fairness measures except for conditional statistical
parity, which only decreases by 2 percentage points in accuracy at full parity
compared to 25% discrimination. In the Adult dataset, we observe a less
drastic trade-off, but this is because the optimal tree with no fairness bound
is quite fair to begin with; for instance, its worst discrimination for
statistical parity is 11%, and for equal opportunity it is 6%. Nonetheless,
there are sufficient points in the trendline for decision makers to similarly
make judgments on the accuracy-discrimination trade-off.

##  4\. Summary

We presented an MIO formulation for learning optimal classification trees that
can be modeled to consider a variety of algorithmic fairness notions. We
compared our method’s performance to a state-of-the-art algorithm for building
fair trees on popular datasets, improving out-of-sample accuracy by 2.3
percentage points on average for a fixed discrimination threshold making it
more appropriate to use in high stakes settings. Our method enables decision-
makers to incorporate any fairness definition as they see fit, as well as to
better fine-tune the trade-off between accuracy and fairness.

###### Acknowledgements.

P. Vayanos and S. Aghaei gratefully acknowledge support from the Hilton C.
Foundation, the Homeless Policy Research Institute, and the Home for Good
foundation under the “C.E.S. Triage Tool Research & Refinement” grant. P.
Vayanos is also grateful for the support of the National Science Foundation
under CAREER award number 2046230. A. Gómez is funded in part by the National
Science Foundation under grant 2006762.

## References

  * (1)
  * Agarwal et al. (2018) Alekh Agarwal, Alina Beygelzimer, Miroslav Dudík, John Langford, and Hanna Wallach. 2018\.  A reductions approach to fair classification. In _International Conference on Machine Learning_. PMLR, ICML, Stockholm, Sweden, 60–69.
  * Aghaei et al. (2019) Sina Aghaei, Mohammad Javad Azizi, and Phebe Vayanos. 2019.  Learning optimal and fair decision trees for non-discriminative decision-making. In _33rd AAAI Conference on Artificial Intelligence_ , Vol. 33. AAAI, Honolulu, USA, 1418–1426.
  * Aghaei et al. (2021) Sina Aghaei, Andrés Gómez, and Phebe Vayanos. 2021\.  Strong Optimal Classification Trees.
  * Angwin et al. (2016) Julia Angwin, Jeff Larson, Surya Mattu, and Lauren Kirchner. 2016.  Machine bias.  https://www.propublica.org/article/machine-bias-risk-assessments-in-criminal-sentencing
  * Azizi et al. (2018) Mohammad Javad Azizi, Phebe Vayanos, Bryan Wilder, Eric Rice, and Milind Tambe. 2018\.  Designing fair, efficient, and interpretable policies for prioritizing homeless youth for housing resources.  _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ 10848 LNCS (2018), 35–51.  https://doi.org/10.1007/978-3-319-93031-2_3
  * Berk et al. (2017) Richard Berk, Hoda Heidari, Shahin Jabbari, Matthew Joseph, Michael Kearns, Jamie Morgenstern, Seth Neel, and Aaron Roth. 2017\.  A convex framework for fair regression.
  * Bertsimas and Dunn (2017) Dimitris Bertsimas and Jack Dunn. 2017.  Optimal classification trees.  _Machine Learning_ 106, 7 (2017), 1039–1082.
  * Bertsimas et al. (2019) Dimitris Bertsimas, Jack Dunn, Colin Pawlowski, and Ying Daisy Zhuo. 2019.  Robust classification.  _INFORMS Journal on Optimization_ 1, 1 (2019), 2–34.
  * Breiman et al. (1984) Leo Breiman, Jerome H Friedman, Richard A Olshen, and Charles J Stone. 1984.  _Classification and regression trees_.  Routledge, Boca Raton, USA.
  * Calders and Verwer (2010) Toon Calders and Sicco Verwer. 2010.  Three naive bayes approaches for discrimination-free classification.  _Data mining and knowledge discovery_ 21, 2 (2010), 277–292.
  * Chen et al. (2019) Jiahao Chen, Nathan Kallus, Xiaojie Mao, Geoffry Svacha, and Madeleine Udell. 2019.  Fairness under unawareness: Assessing disparity when protected class is unobserved. In _Proceedings of the 2019 Conference on Fairness, Accountability, and Transparency_. Association for Computing Machinery, Atlanta, USA, 339–348.
  * Chouldechova (2017) Alexandra Chouldechova. 2017\.  Fair prediction with disparate impact: A study of bias in recidivism prediction instruments.  _Big data_ 5, 2 (2017), 153–163.
  * Corbett-Davies et al. (2017) Sam Corbett-Davies, Emma Pierson, Avi Feller, Sharad Goel, and Aziz Huq. 2017.  Algorithmic decision making and the cost of fairness. In _Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , Vol. Part F129685. ACM, Halifax, Canada, 797–806.  https://doi.org/10.1145/3097983.3098095
  * Denis et al. (2021) Christophe Denis, Romuald Elie, Mohamed Hebiri, and François Hu. 2021.  Fairness guarantee in multi-class classification.
  * Dua and Graff (2017) Dheeru Dua and Casey Graff. 2017.  UCI Machine Learning Repository.  http://archive.ics.uci.edu/ml
  * Dwork et al. (2012) Cynthia Dwork, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Richard Zemel. 2012.  Fairness Through Awareness. In _Proceedings of the 3rd Innovations in Theoretical Computer Science Conference_ (Cambridge, Massachusetts) _(ITCS ’12)_. ACM, New York, NY, USA, 214–226.  https://doi.org/10.1145/2090236.2090255
  * Elmachtoub et al. (2020) Adam Elmachtoub, Jason Cheuk Nam Liang, and Ryan McNellis. 2020\.  Decision trees for decision-making under the predict-then-optimize framework. In _37th International Conference on Machine Learning_. PMLR, ICML, Vienna, Austria, 2858–2867.
  * Fatima et al. (2017) Meherwar Fatima, Maruf Pasha, et al. 2017\.  Survey of machine learning algorithms for disease diagnostic.  _Journal of Intelligent Learning Systems and Applications_ 9, 01 (2017), 1.
  * Feldman et al. (2015) Michael Feldman, Sorelle A Friedler, John Moeller, Carlos Scheidegger, and Suresh Venkatasubramanian. 2015\.  Certifying and removing disparate impact. In _Proceedings of the 21st ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. Association for Computing Machinery, Sydney, Australia, 259–268.
  * Gajane and Pechenizkiy (2017) Pratik Gajane and Mykola Pechenizkiy. 2017.  On formalizing fairness in prediction with machine learning.
  * Grari et al. (2019) Vincent Grari, Boris Ruf, Sylvain Lamprier, and Marcin Detyniecki. 2019.  Fair adversarial gradient tree boosting. In _2019 IEEE International Conference on Data Mining (ICDM)_. IEEE, IEEE, Beijing, China, 1060–1065.
  * Hajian and Domingo-Ferrer (2012) Sara Hajian and Josep Domingo-Ferrer. 2012.  A methodology for direct and indirect discrimination prevention in data mining.  _IEEE transactions on knowledge and data engineering_ 25, 7 (2012), 1445–1459.
  * Hardt et al. (2016) Moritz Hardt, Eric Price, and Nathan Srebro. 2016.  Equality of opportunity in supervised learning. In _Advances in Neural Information Processing Systems_. NeurIPS, Barcelona, Spain, 3315–3323.
  * Jo et al. (2021) Nathanael Jo, Sina Aghaei, Andrés Gómez, and Phebe Vayanos. 2021.  Learning optimal prescriptive trees from observational data.
  * Justin et al. (2021) Nathan Justin, Sina Aghaei, Andres Gomez, and Phebe Vayanos. 2021.  Optimal Robust Classification Trees. In _The AAAI-22 Workshop on Adversarial Machine Learning and Beyond_. AAAI Press, Vancouver, Canada.
  * Kamiran et al. (2010) Faisal Kamiran, Toon Calders, and Mykola Pechenizkiy. 2010\.  Discrimination aware decision tree learning. In _Proceedings - IEEE International Conference on Data Mining, ICDM_. IEEE, Sydney, Australia, 869–874.  https://doi.org/10.1109/ICDM.2010.50
  * Kamiran et al. (2013) Faisal Kamiran, Indrė Žliobaitė, and Toon Calders. 2013\.  Quantifying explainable discrimination and removing illegal discrimination in automated decision making.  _Knowledge and information systems_ 35, 3 (2013), 613–644.
  * Kamishima et al. (2012) Toshihiro Kamishima, Shotaro Akaho, Hideki Asoh, and Jun Sakuma. 2012.  Fairness-aware classifier with prejudice remover regularizer. In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_. Springer, ECML PKDD, Prague, Czech Republic, 35–50.
  * Kanamori and Arimura (2019) Kentaro Kanamori and Hiroki Arimura. 2019.  Fairness-aware Edit of Thresholds in a Learned Decision Tree Using a Mixed Integer Programming Formulation. In _33rd Annual Conference of the Japanese Society for Artificial Intelligence, (2019)_. Japanese Society for Artificial Intelligence, Japanese Society for Artificial Intelligence, Niigata City, Japan, 3Rin211–3Rin211.
  * Mehrabi et al. (2021) Ninareh Mehrabi, Fred Morstatter, Nripsuta Saxena, Kristina Lerman, and Aram Galstyan. 2021\.  A survey on bias and fairness in machine learning.  _ACM Computing Surveys (CSUR)_ 54, 6 (2021), 1–35.
  * Mišić (2020) Velibor V Mišić. 2020\.  Optimization of tree ensembles.  _Operations Research_ 68, 5 (2020), 1605–1624.
  * Pedreshi et al. (2008) Dino Pedreshi, Salvatore Ruggieri, and Franco Turini. 2008\.  Discrimination-aware data mining. In _Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. Association for Computing Machinery, Las vegas, USA, 560–568.
  * Pleiss et al. (2017) Geoff Pleiss, Manish Raghavan, Felix Wu, Jon Kleinberg, and Kilian Q. Weinberger. 2017.  On fairness and calibration. In _31st International Conference on Neural Information Processing Systems_ , Vol. 2017-December. NeurIPS, Long Beach, USA, 5684––5693.
  * Ranzato et al. (2021) Francesco Ranzato, Caterina Urban, and Marco Zanella. 2021\.  Fair Training of Decision Tree Classifiers.
  * Rudin (2019) Cynthia Rudin. 2019\.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature Machine Intelligence_ 1, 5 (2019), 206–215.
  * Verwer and Zhang (2019) Sicco Verwer and Yingqian Zhang. 2019.  Learning optimal classification trees using a binary linear program formulation. In _33rd AAAI Conference on Artificial Intelligence_ , Vol. 33. Association for the Advancement of Artificial Intelligence, Honolulu, USA, 1625–1632.
  * Zemel et al. (2013) Rich Zemel, Yu Wu, Kevin Swersky, Toni Pitassi, and Cynthia Dwork. 2013.  Learning fair representations. In _30th International Conference on Machine Learning_. PMLR, ICML, Atlanta, USA, 325–333.
  * Zhang et al. (2020) Jiang Zhang, Ivan Beschastnikh, Sergey Mechtaev, and Abhik Roychoudhury. 2020.  Fairness-guided SMT-based Rectification of Decision Trees and Random Forests.
  * Zhang and Ntoutsi (2019) Wenbin Zhang and Eirini Ntoutsi. 2019.  FaHT: An adaptive fairness-aware decision tree classifier. In _IJCAI International Joint Conference on Artificial Intelligence_ , Vol. 2019-August. IJCAI, Macao, China, 1480–1486.  https://doi.org/10.24963/ijcai.2019/205
  * Žliobaitė (2017) Indrė Žliobaitė. 2017\.  Measuring discrimination in algorithmic decision making.  _Data Mining and Knowledge Discovery_ 31, 4 (2017), 1060–1089.