# Cold damping of levitated optically coupled nanoparticles

Vojtěch Liška  Tereza Zemánková  Vojtěch Svak  Petr Jákl  Jan Ježek  Martin
Bránecký  Stephen H. Simpson  Pavel Zemánek  Oto Brzobohatý
<EMAIL_ADDRESS>The Czech Academy of Sciences, Institute of Scientific
Instruments, Královopolská 147, 612 64 Brno, Czech Republic

###### Abstract

Methods for controlling the motion of single particles, optically levitated in
vacuum, have developed rapidly in recent years. The technique of cold damping
makes use of feedback-controlled, electrostatic forces to increase dissipation
without introducing additional thermal fluctuations. This process has been
instrumental in the ground-state cooling of individual electrically charged
nanoparticles. Here we show that the same method can be applied to a pair of
nanoparticles, coupled by optical binding forces. These optical binding forces
are about three orders of magnitude stronger than typical Coulombic inter-
particle force and result in a coupled motion of both nanoparticles
characterized by a pair of normal modes. We demonstrate cold damping of these
normal modes, either independently or simultaneously, to sub-Kelvin
temperatures at pressures of $5\times 10^{-3}$ mbar. Experimental observations
are captured by a theoretical model which we use to survey the parameter space
more widely and to quantify the limits imposed by measurement noise and time
delays. Our work paves the way for the study of quantum interactions between
meso-scale particles and the exploration of multiparticle entanglement in
levitated optomechanical systems.

##  I Introduction

Due to its isolation from the environment, a single nanoparticle, optically
levitated in an ultrahigh vacuum, provides a promising experimental platform
for weak force sensing Ranjit  _et al._ (2016); Hempston  _et al._ (2017);
Hebestreit  _et al._ (2018a) and for testing fundamental physics at the
boundary between the classical and quantum regimes Svak  _et al._ (2018);
Rondin  _et al._ (2017). With the recent achievement of ground state cooling
for individual nanoparticles Delić  _et al._ (2020); Magrini  _et al._ (2021);
Tebbenjohanns  _et al._ (2021), levitational optomechanics opens up
opportunities for experimentally exploring quantum effects at previously
unattainable length scales and, potentially, for designing sensors with
quantum-enhanced sensitivity. Extending levitational optomechanics to arrays
of multiple, interacting particles is an exciting new development, still in
its infancy Rieser  _et al._ (2022), which promises to open new research
directions in quantum gravityMarletto and Vedral (2017), quantum friction
measurements Zhao  _et al._ (2012), dark matter detection Moore and Geraci
(2021) or probing quantum correlations and entanglement Rudolph  _et al._
(2022).

To fully exploit the potential of these optomechanical arrays requires both
understanding and control of the interaction forces acting between the
particles. It also requires new protocols capable of cooling the multiple
degrees of freedom defining these complex systems. The number of experimental
demonstrations of optical trapping of interacting multiple particles is
currently very limited. Bykov et al. reported long-range optical binding of
multiple levitated microparticles mediated by intermodal scattering and
interference in the evacuated core of a hollow-core photonic crystal fibre
Bykov  _et al._ (2018), while Arita et al. demonstrated optical binding
between two rotating chiral microparticles confined in vacuum in independent
circularly polarised optical traps Arita  _et al._ (2018).

The tunable longitudinal Svak  _et al._ (2021) and lateral Rieser  _et al._
(2022) optical interactions between levitated nanoparticles have been
demonstrated only very recently. An all-optical cold damping cooling approach
Simmons (1996); Kamba  _et al._ (2022) has been very recently demonstrated for
non-interacting nanoparticles Vijayan  _et al._ (2022) while a sympathetic
cooling scheme, inspired by experiments of cooling motion of atoms and ions
Myatt  _et al._ (1997); Larson  _et al._ (1986), has been applied for
optically interacting rotating birefringent microparticles Arita  _et al._
(2022). Finally, the sympathetic cooling of two particles, coupled by a
Coulombic interaction, and levitated in a Paul trap has been shown recently
Bykov  _et al._ (2023); Penny  _et al._ (2023). In this latter case, the
mechanical oscillatory frequencies were typically several orders of magnitude
smaller than in the case of optical levitation.

Here, we demonstrate cooling of the motion of two optically bound Wei  _et
al._ (2016) nanoparticles levitated in two parallel optical traps in a vacuum
to the sub-Kelvin level. This all-optical interactionBurns  _et al._ (1990);
Tatarkova  _et al._ (2002); Singer  _et al._ (2003); Dholakia and Zemánek
(2010); Čižmár  _et al._ (2010); Demergis and Florin (2012); Yan  _et al._
(2014); Li  _et al._ (2018); Čižmár  _et al._ (2011); Brzobohatý  _et al._
(2011) is mediated by light scattering and, in our case, offers the ability to
precisely control the interactions between particles Rieser  _et al._ (2022);
Wei  _et al._ (2016). For cooling the motion of the particles, we employ a
cold-damping scheme Tebbenjohanns  _et al._ (2019); Magrini  _et al._ (2021)

Figure 1:  Experimental geometry and observations. a) Two silica nanoparticles
levitated in two parallel optical traps created by two pairs of interfering
counter-propagating Gaussian beams. The position of each particle is monitored
by the corresponding QPD. Independently, for QPD’s calibration, the $x-z$
plane was imaged on an ultra-fast CMOS camera. The box illustrates the
feedback loop using FPGA processing. b) The coupling rate can be controlled by
varying the trap separation, $d_{0}$. The black arrow marks the separation
between the traps used in our experiments. c) PSDs of the particle
coordinates, $x_{1}$ and $x_{2}$, plotted together with PSDs of the center of
mass $a_{\Sigma}$ and the breathing modes $a_{\Delta}$ for the coupled
harmonic oscillators.

which makes use of an external electric field to induce controllable,
electrostatic forces on the pair of charged nanoparticles. This approach has
previously been shown to be the most efficient feedback-based cooling method
currently available Penny  _et al._ (2021); Iwasaki  _et al._ (2019) and is,
in theory, capable of preparing entangled states of interacting particles or
acting as an ultra-sensitive sensor of the gradient of a probed external force
field Rudolph  _et al._ (2022).

##  II Experimental geometry

In our experiments, two silica nanoparticles are levitated in two parallel
optical traps created by two pairs of interfering counter-propagating Gaussian
beams, see Fig. 1a). We employ a nebulizer to spray nanoparticles directly
into the trapping region in the vacuum chamber. After successfully trapping
two nanoparticles, we evacuated the vacuum chamber to sub-millibar pressure
and monitored the motion of the particles using two custom-made quadrant
photodiodes (QPDs, Hamamatsu G6849) with homodyne detection of the light
transmitted through the optical traps, which has the particle motion encoded
in its interference pattern, see Fig. 1a). In addition, we illuminated the
particles along $y$ axis with a weak and wide laser beam ($\lambda=532$ nm,
beam waist $w_{0}\approx 50\,\mu$m) and imaged the scattering patterns of both
nanoparticles

with a microscope, onto a fast carefully calibrated CMOS camera (Vision
Research Phantom V611), see the colour inverted image of nanoparticles in Fig.
1a). The frame rate was set to 200 kHz and, typically, we recorded at least
100 000 frames, tracking the motion of the nanoparticles, to obtain
sufficiently long trajectories for the analysis of their stochastic dynamics.
This provided us with information about the in-plane $x$ and $z$ coordinates.
The record gave us quantitative, properly calibrated information about the
amplitude of particles’ motion without any prior knowledge of the particles’
properties Svak  _et al._ (2018). And by comparing parallel records from the
camera and QPDs we calibrated the QPDs signals too Svak  _et al._ (2018).

The radius of levitated nanoparticles ($\rho=305$ nm) was significantly
smaller than the wavelength of trapping beams and their beam waists
($\lambda=1550$ nm, $w_{0}=1.5\,\mu$m, $\rho/\lambda\approx 0.2$) and
therefore the nanoparticles behave as induced dipoles. Due to the particular
symmetry of our system and our use of relatively wide counter-propagating
beams the optical binding force acts mainly along the vector connecting the
particles ($x$ axis) which is perpendicular to the beam axes [i.e. $z$ axis in
Fig. 1a)]. We note that our transversal direction of the binding force differs
from that which acts prominently along the optical axis in parallel optical
tweezers Rieser  _et al._ (2022).

The phase-coherent optical traps were generated by first-order diffraction
from a digital micromirror device (DMD, Vialux). The total trapping power
$2P=140$ mW was split into two independent optical traps in a vacuum chamber,
which allows us to independently set the trapping stiffness of both traps. The
direction of polarization of the beams was controlled with a pair of half-wave
plates and set along $y$ axis to maximize the scattering along $x$ axis and
thus maximize the light-induced particle-particle interaction (optical
binding). In the dipole approximation, the resulting optical binding forces
oscillate periodically Wei  _et al._ (2016); Burns  _et al._ (1989) and decay
as $F_{1,2}\propto\pm\sin(kd\pm\Delta\phi)/kd,$ where $k=2\pi/\lambda$, $d$ is
the inter-particle distance and $\Delta\phi$ denotes optical phase difference
between the trapping laser beams incident on particle 1 and 2. Whenever
$\Delta\phi\neq 0$, both binding forces have a component pointing in the same
direction and forming a non-reciprocal and non-conservative interaction Rieser
_et al._ (2022). Using the DMD we set $\Delta\phi=0$ to ensure that the non-
conservative part of the interaction Rieser  _et al._ (2022) is negligible and
set the distance between the optical traps to be approximately $d_{0}=5\,\mu$m
to maximize the conservative part of the optical binding force, see Fig. 1b).

Since the non-conservative part of the interaction is negligible, and the
particles remain within the linear range of the traps, our system can be
described as a pair of coupled oscillators Svak  _et al._ (2021) and thus
their motion can be described in terms of a linear combination of two normal
modes (Supplemental II A) i.e. the centre-of-mass (CoM) mode, with amplitude
$a_{\Sigma}=(x_{1}+x_{2})/2$ and the breathing (BR) mode, with amplitude
$a_{\Delta}=(x_{1}-x_{2})/2$. This can be clearly seen in the power spectral
density of the particles motion, Fig. 1c, which features two distinct resonant
frequencies. The lower frequency resonance corresponds to the CoM mode,
$\Omega_{\Sigma}/2\pi=(\Omega_{0}-b)/2\pi\approx 21$ kHz, and the higher
corresponds to BR, $\Omega_{\Delta}/2\pi=(\Omega_{0}+b)/2\pi\approx 26$ kHz,
where $\Omega_{0}$ is the resonant frequency of an individual trap. The
parameter, $b$, is the coupling rate of the optical interaction, given by the
derivative of the binding force evaluated at the stable positions of the
particles (i.e. the derivative of $F_{1,2}\propto\cos(kx)/kx$ with respect to
$x$, evaluated at $d$), Fig. 1b. We note that the separation of optical traps
$d_{0}$ differs slightly from inter-particle separation $d$ due to the
displacement of the particles in the traps, caused by the optical binding
force.

Each particle is randomly charged. To estimate the magnitude of the Coulomb
interaction between the particles we performed charge calibrationMagrini  _et
al._ (2021). A sinusoidal voltage
$V^{\mathrm{d}}(t)=V^{\mathrm{d}}_{0}\sin\Omega_{\mathrm{d}}t$ with the
driving frequency $\Omega_{\mathrm{d}}/2\pi=50$ kHz relatively close to the
resonant frequencies of the normal modes ($\Omega_{\Sigma}$ and
$\Omega_{\Delta}$). From the recorded trajectories of the particles, we
calculated the power spectral densities for both normal modes and, by
integrating them around driving frequency, we get variances $\langle
a_{\Sigma}^{2}\rangle$ and $\langle a_{\Delta}^{2}\rangle$. Assuming a model
for a point charge in a parallel-plate capacitor, we get driving force acting
on both modes in the form
$F^{\mathrm{d}}_{\Sigma/\Delta}(t)=q_{\Sigma/\Delta}V^{\mathrm{d}}(t)/D,$
where $D=(198\pm 5)\,\mu$m is the distance between electrodes and
$q_{\Sigma}=(q_{1}+q_{2})/2$ and $q_{\Delta}=(q_{1}-q_{2})/2$ is the effective
charge for the CoM and BR normal modes, respectively. These charges can be
then expressed as

$q_{\Sigma/\Delta}=\sqrt{2\langle
a_{\Sigma/\Delta}^{2}\rangle}mD/V^{\mathrm{d}}_{0}\sqrt{(\Omega_{\Sigma/\Delta}^{2}-\Omega_{\mathrm{d}}^{2})^{2}+\xi_{\Sigma/\Delta}^{2}\Omega_{\mathrm{d}}^{2}},$
(1)

where $\xi_{\Sigma/\Delta}$ are the effective damping coefficients for the
normal modes Svak  _et al._ (2021), $m$ is the mass of the particle and
$V^{\mathrm{d}}_{0}=750$ mV is amplitude of the applied voltage. The charges
on the particles can then be written as $q_{1,2}=(q_{\Sigma}\pm q_{\Delta})$,
and the typical number of elementary charges on the particles was determined
to be $N<100$ and thus the magnitude of coupling rate
$b/\Omega_{0}=-q_{1}q_{2}/(8\pi\varepsilon_{0}m\Omega_{0}^{2}d^{3})<10^{-4}$,
for Coulomb interaction Rudolph  _et al._ (2022); Rieser  _et al._ (2022) is 3
order of magnitude smaller than that one obtained for optical binding
interaction, see Supplemental II C. This conclusion we confirmed
experimentally when we suppressed the optical binding interaction by setting
the polarization of the laser beams parallel to the $x$ axis and leaving only
the Coulomb interaction. In this geometry, no mode splitting was detectable,
showing that the Coulomb interaction is indeed negligible.

Figure 2:  Experimental cold damping of normal modes. With increasing feedback
gain, the PSD broadens in width and simultaneously shrinks in area. a) and b)
shows PSDs of both modes when the center of mass ($\Sigma$) and the breathing
($\Delta$) mode, respectively, are cooled independently. c) Both modes are
cooled simultaneously. The dashed line denotes the fit of PSD using Eqs. (4a)
and (4b). The external electric field necessarily couples the normal modes
causing the height of the second peak to increase as well. The shift of the
resonant oscillatory frequencies is caused by the finite time delay in the
feedback loop $\tau=5.6\,\mu$s. Experimental parameters: ambient pressure
$p=0.2$ mbar, particle charges $q_{1}=88$e and $q_{2}=6$e.

##  III Results and discussion

We applied a controllable, time-dependent electrostatic force to the net
charges carried by the optically trapped nano-particles by applying a voltage
to a pair of electrodes enclosing the trap [see Fig 1a)], to implement the
cold damping feedback cooling scheme Tebbenjohanns  _et al._ (2019). This
allowed us to control the effective temperatures of the normal modes
independently. A complete and general analysis for feedback cooling of coupled
conservative oscillators is provided in the Supplementary Information. To
start with, we focus on the case where the system is perfectly symmetrical
i.e. the traps and particles are identical. The equation of motion of our
system of the optically bound particles can be expressed in the normal mode
basis as

$m\ddot{\mathbf{a}}(t)=-\mathbf{K}\mathbf{a}(t)-\mathbf{Z}\dot{\mathbf{a}}(t)-\mathbf{G}\dot{\mathbf{a}}(t-\tau)+\mathbf{f}^{L}(t),$
(2)

where $\mathbf{a}=(a_{\Sigma},a_{\Delta})$ are the amplitudes of the centre of
mass and breathing modes, with $\dot{\mathbf{a}}$ and $\ddot{\mathbf{a}}$ the
first and second derivatives with respect to time,
$\mathbf{K}=\mathrm{Diag}(K_{\Sigma},K_{\Delta})$ are the modal stiffness,
$\mathbf{Z}=\mathrm{Diag}(\xi_{\Sigma},\xi_{\Delta})$ are the effective
damping coefficients ($\xi_{\Sigma}\neq\xi_{\Delta}$ due to hydrodynamic
coupling), $\mathbf{f}^{L}$ is the Langevin force and finally,
$-\mathbf{G}\dot{\mathbf{a}}(t-\tau)$ is cold damping force, deriving from the
external electric force, and $\tau$ is a time delay incurred due to finite,
experimental response times. Since this electric force acts on both particles,
in proportion to their differing charges, it necessarily couples the normal
modes:

$\mathbf{f}^{\mathrm{fb}}(t)=\begin{bmatrix}f^{\mathrm{fb}}_{\Sigma}\\\
f^{\mathrm{fb}}_{\Delta}\end{bmatrix}=g\begin{bmatrix}c_{\Sigma}q_{\Sigma}&c_{\Delta}q_{\Sigma}\\\
c_{\Sigma}q_{\Delta}&c_{\Delta}q_{\Delta},\end{bmatrix}\dot{\mathbf{a}}(t-\tau)\equiv
g^{\prime}\begin{bmatrix}c_{\Sigma}&c_{\Delta}\\\
c_{\Sigma}r&c_{\Delta}r\end{bmatrix}\dot{\mathbf{a}}(t-\tau)\equiv-\mathbf{G}\dot{\mathbf{a}}(t-\tau),$
(3)

where $g$ is global gain, and $c_{\Sigma}$, $c_{\Delta}$ are modal gains, $r$
is the charge ratio, $r=q_{\Delta}/q_{\Sigma}$ and $g^{\prime}=gq_{\Sigma}$ is
the effective scalar gain.

We measured the delay of our feedback loop to be $\tau=5.6\,\mu$s which is
approximately $1/10$ of the time period of the oscillations in the system. Our
linear feedback signal was implemented on an FPGA card, see Fig. 1a). In the
majority of the experiments, we use a low-pass filter whose cut-off frequency
was set to above 100 kHz which is far from the normal mode frequencies, such
that we achieve a flat phase response. In the case of the lowest pressure
($p=5\times 10^{-3}$ mbar) we applied bandpass filters centred around modal
frequencies to get rid of low-frequency noise and detection crosstalks with
$y$ axis. We calculate instantaneous mode velocities numerically and, for the
products, we multiplied modal gains $c_{\Sigma}$ and $c_{\Delta}$. The final
feedback voltage from the FPGA card was amplified and applied to one of the
electrodes (the second one was grounded).

We investigated the cooling performance as a function of the gas pressure and
feedback gain to explore the limitations of the method. In Fig. 2, we show the
single-sided PSDs for both modes (blue colour denotes the CoM i.e. $\Sigma$
mode, red colour denotes the BR i.e. $\Delta$ mode). Three different
combinations were tested. In particular, we applied the cold damping feedback
cooling scheme on each of the modes separately (e.g., $c_{\Sigma}=1$ or
$c_{\Delta}=0$), see Fig 2a and 2b, or we applied feedback for both modes
simultaneously ($c_{\Sigma}=c_{\Delta}=1/\sqrt{2}$), see Fig 2c. These
scenarios are described in Supplemental IV B.

The measured signal corresponds to the power spectral densities (PSD),
$S^{xx}_{\Sigma/\Delta}$, where the superscript $xx$ indicates the positional
PSD, as opposed to the PSD for velocities which we denote,
$S^{vv}_{\Sigma/\Delta}$. The PSDs are added to a spectrally flat noise floor,
$S_{\mathrm{nn}}$, associated with noise in our quadrant photodiode detectors
and to a further contribution, $S^{\mathrm{fb}}_{\Sigma/\Delta}$, which
quantifies the effective heating of the particles due to measurement noise
entering the feedback loop (see Supplemental III). We extracted the damping
coefficients $\mathbf{Z}$ and effective feedback-induced damping rates
$\mathbf{G}$ by fitting the PSDs of modes :

$\displaystyle S^{xx}_{\Sigma}$
$\displaystyle=\frac{2k_{\mathrm{B}}T\xi_{0}}{|P_{\Sigma}P_{\Delta}+\Omega^{2}G_{21}G_{12}e^{-2i\Omega\tau}|^{2}}(|P_{\Delta}|^{2}+\Omega^{2}G_{12}^{2})+S^{\mathrm{fb}}_{\Sigma}(S_{\mathrm{nn}})+S_{\mathrm{nn}},$
(4a) $\displaystyle S^{xx}_{\Delta}$
$\displaystyle=\frac{2k_{\mathrm{B}}T\xi_{0}}{|P_{\Sigma}P_{\Delta}+\Omega^{2}G_{21}G_{12}e^{-2i\Omega\tau}|^{2}}(|P_{\Sigma}|^{2}+\Omega^{2}G_{21}^{2})+S^{\mathrm{fb}}_{\Delta}(S_{\mathrm{nn}})+S_{\mathrm{nn}},$
(4b)

where $\Omega$ is the frequency and $\xi_{0}=6\pi\mu\rho$ is the usual Stokes
drag. The contributions from noise-induced heating are,

$\displaystyle S^{\mathrm{fb}}_{\Sigma}$
$\displaystyle=\Bigg{(}\frac{|P_{\Delta}|^{2}+2r\Omega
G_{12}\Im(P_{\Delta}e^{i\Omega\tau})+r^{2}\Omega^{2}G_{12}^{2}}{|P_{\Sigma}P_{\Delta}+\Omega^{2}G_{12}G_{21}e^{-2i\Omega\tau}|^{2}}\Bigg{)}\Omega^{2}g^{\prime
2}S_{\mathrm{nn}},$ (5a) $\displaystyle S^{\mathrm{fb}}_{\Delta}$
$\displaystyle=\Bigg{(}\frac{r^{2}|P_{\Sigma}|^{2}+2r\Omega
G_{21}\Im(P_{\Sigma}e^{i\Omega\tau})+\Omega^{2}G_{21}^{2}}{|P_{\Sigma}P_{\Delta}+\Omega^{2}G_{12}G_{21}e^{-2i\Omega\tau}|^{2}}\Bigg{)}\Omega^{2}g^{\prime
2}S_{\mathrm{nn}},$ (5b)

where $P_{\Sigma/\Omega}$ are quadratic in $\Omega$,

$P_{\Sigma}(\Omega)=m\left[\left(\Omega^{2}-\frac{G_{11}\sin(\Omega\tau)}{m}\Omega-\frac{K_{\Sigma}}{m}\right)-\frac{i\Omega}{m}\big{(}\xi_{\Sigma}+G_{11}\cos(\Omega\tau)\big{)}\right],$
(6)
$P_{\Delta}(\Omega)=m\left[\left(\Omega^{2}-\frac{G_{22}\sin(\Omega\tau)}{m}\Omega-\frac{K_{\Delta}}{m}\right)-\frac{i\Omega}{m}\big{(}\xi_{\Delta}+G_{22}\cos(\Omega\tau)\big{)}\right].$
(7)

The resonant peaks appearing in the PSD correspond approximately to the minima
of $|P_{\Sigma}(\Omega)|^{2}$ and $|P_{\Delta}(\Omega)|^{2}$. When $\tau=0$
this gives peaks at $\Omega_{\Sigma}\approx\sqrt{K_{\Sigma}/m}$ and
$\Omega_{\Delta}\approx\sqrt{K_{\Delta}/m}$. However, when the time delta,
$\tau$, is appreciable, the resonant frequencies shift significantly with
increasing gain, see Fig. 2 and the Supplemental IV C.

Figure 3:  Cooling performance visualized using the effective temperatures of
the modes. The values of the feedback damping $\mathbf{G}$ coefficients were
determined by fitting the experimental PSDs, $S^{xx}_{\Sigma/\Delta}$. a) and
b) shows the effective temperatures when the centre of mass and the breathing
mode, respectively, are cooled independently. c) Both modes cooled
simultaneously. Sub-Kelvin effective temperature was achieved for both modes
at $5\times 10^{-3}$ mbar. Full curves show temperatures determined from
theoretical PSDs [Eqs. (4a) and (4b)], dotted curves show the temperature when
the values of coupling terms were set to be zero $G_{12}=G_{21}=0$, dashed and
dash-dotted curves illustrate heating of the motion of the particles due to
the coupling induced by the external field and due to the measurement noise
fed back to control feedback loop, respectively. The filled areas in part a)
illustrate the increase of the effective temperatures caused by feedback loop
delay $\tau=5.6\,\mu$s, i.e. the lower values of temperatures are calculated
for zero delay. (for $p=1$ and 0.2 mbar, $q_{1}=88$e and $q_{2}=6$e, for
$p=5\times 10^{-3}$ mbar, $q_{1}=27$e and $q_{2}=1$e).

As expected, as we increase the value of the global feedback gain $g$, the PSD
for cooled mode broadens in width and simultaneously shrinks in area, see Figs
2. When the global gain, $g$, is zero, each modal PSDs contain a single peak,
corresponding to the resonant frequency of the associated mode. However, as
$g$ is increased, a small, subsidiary peak begins to appear at the resonant
frequency of the second mode, indicating a degree of feed-back induced inter-
modal coupling, see Figs. 2a) and 2b). This limits the effect of the feedback
cooling scheme. For independently cooled modes, this effect would be absent in
a perfectly symmetric system (with identical particles and traps), see
Supplemental information. In contrast, when cooling both modes simultaneously,
the inter-modal coupling is introduced via off-diagonal terms of the gain
matrix, $\mathbf{G}$, and is completely unavoidable here, see Fig. 2c).

To quantify the cooling performance we calculated the effective temperatures
$T_{\Sigma/\Delta}$ of the oscillation modes from the measured displacement
time traces $a_{\Sigma}(t)$ and $a_{\Delta}(t)$. To do this we used the
kinetic energy of the oscillators

$T_{\Sigma/\Delta}=2\langle
E^{\mathrm{kin}}_{\Sigma/\Delta}\rangle/k_{\mathrm{B}}=m\langle\dot{a}_{\Sigma/\Delta}^{2}\rangle/k_{\mathrm{B}},$
(8)

for which the equipartition theorem is valid even for anharmonic trapping
potentials Hebestreit  _et al._ (2018b). For obtained
$\langle\dot{a}_{\Sigma/\Delta}^{2}\rangle$ we numerically integrate the power
spectral density for mode velocity which can be expressed using power spectral
density for particle displacement,

$\langle\dot{a}_{\Sigma/\Delta}^{2}\rangle=\int_{0}^{\infty}S^{vv}_{\Sigma/\Delta}(\Omega)\mathrm{d}\Omega=\int_{0}^{\infty}\Omega^{2}[S^{xx}_{\Sigma/\Delta}(\Omega)-S_{\mathrm{nn}}]\mathrm{d}\Omega.$
(9)

It should be noted that we have subtracted the detection noise floor
$S_{\mathrm{nn}}$ from our experimentally determined power spectral density.

Figures 3a) – 3c) show the magnitude of the effective temperature for
different values of gas pressure ($p=1,\,0.2$ and $5\,\times 10^{-3}$ mbar)
determined for increasing values of feedback induced damping
$G_{11}=gc_{\Sigma}q_{\Sigma}$ and $G_{22}=gc_{\Delta}q_{\Delta}$ which was
controlled, during the experiment, by the magnitude of the global feedback
gain $g$. Complementary to the power spectral densities presented in Fig. 2,
we compare here three cases, i.e. independent cooling of the modes
[$c_{\Delta}=0$ in Fig. 3a) and $c_{\Sigma}=0$ in Fig. 3b)] and simultaneous
cooling of both modes Fig. 3c). At the lowest studied pressure ($5\,\times
10^{-3}$ mbar) we achieved sub-Kelvin temperature for both modes. For all the
studied cases, there is an optimal value of feedback damping rate $G_{ii}/m$
for which we reach minimal values of effective temperature for a particular
surrounding gas pressure.

Figure 4:  Extrapolated cooling performance visualized using the effective
temperatures of CoM mode a) and BR mode b). The dashed line illustrates
heating due to the noise in detection, which differs for both modes see Eqs.
(5a) and (5b). The coloured area shows the effect of coupling between modes
due to the electric feedback force. c) The minimal values of effective
temperatures as a function of the gas pressure determined from a) and b).

For independently cooled modes [Figs. 3a) and 3b)] the presence of an optimal
condition arises due to competition between two mechanisms: the cooling effect
to increasing feedback damping, and the inevitable heating effect caused by
feeding back measurement noise into the control loop,
$T^{\mathrm{det}}_{\Sigma/\Delta}\propto\mathbf{G}S_{\mathrm{nn}}$ (see dash-
dotted lines). In this respect, the limits in the efficiency of independently
cooling the normal modes are similar to those operating in the single particle
case Tebbenjohanns  _et al._ (2019).

When independently cooling the CoM mode ($c_{\Sigma}\neq 0$ and
$c_{\Delta}=0$), the electric force, Eq. (3), acting on the CoM mode is not
influenced by the breathing mode,
$f_{\Sigma}=g[c_{\Sigma}q_{\Sigma}v_{\Sigma}+c_{\Delta}q_{\Sigma}v_{\Delta}]=gc_{\Sigma}q_{\Sigma}v_{\Sigma}$,
on the other hand, there is a coupling between un-cooled breathing mode and
CoM mode, mediated by feedback electric force
$f_{\Delta}=g[c_{\Sigma}q_{\Delta}v_{\Sigma}+c_{\Delta}q_{\Delta}v_{\Delta}]=gc_{\Sigma}q_{\Delta}v_{\Sigma}$.
This resulting variance in the BR mode is,

$\langle
a^{2}_{\Delta}\rangle=\frac{k_{\mathrm{B}}T}{K_{\Delta}}+2k_{\mathrm{B}}T\xi_{\Delta}g^{2}q_{\Delta}^{2}\int^{\infty}_{\infty}d\Omega\frac{\Omega^{2}}{|P_{\Sigma}|^{2}|P_{\Delta}|^{2}}.$
(10)

The first term in Eq. (10) is constant, while the second term is necessarily
positive and dependent on the gain. This second term is connected with the
heating of the uncooled mode, due to feedback-induced coupling. Indeed, the
heating of the un-cooled breathing mode was observed experimentally for 1
mbar, see the red line for the total value of total effective temperature and
the red dashed line which illustrates the heating due to the coupling of modes
by the electric field, in Fig. 3a). However, for lower pressure, we observed
an opposite trend, which we observed also in our numerical model by assuming
very subtle asymmetry in our system, where our idealized modes are not fully
valid, see Supplemental information for a full description. The asymmetry in
our numerical model was induced, by considering only a very small difference
in particle sizes $<$ 0.5%, which is much smaller than the standard deviation
of particle size specified by a manufacturer (3%). The non-zero value of time
delay $\tau$ in our feedback loop caused the increase of values of the
effective temperatures $T$ (about $2\times$), as illustrated using colour-
filled areas in Fig. 3a).

Figure 3c) compares the values of the effective temperatures when both modes
were cooled simultaneously. For the sake of clarity, the heating mechanisms
are here illustrated only for the CoM mode. Since the heating due to coupling
between the modes, scales with damping of the surrounding gas [e.g., for CoM
it is $2k_{B}T\xi_{0}\Omega^{2}G_{12}^{2}$, see Eq. (4a)], it is the main
limiting mechanism at pressures higher than $10^{-3}$ mbar, on the other hand
the heating due to measurement noise which is proportional to
$\Omega^{2}g^{\prime 2}S_{\mathrm{nn}}$ [see Eqs. (5a) and (5b)] starts to be
important mainly at the lower pressure, see also Figs. 4a) and 4b) for values
of effective temperature for both simultaneously cooled modes extrapolated for
low pressure [other parameters are same as in Fig. 3c)]. Note, that the
heating due to the detection noise differs for both modes, see Eqs. (5a) and
(5b), and thus the minimal values of effective temperatures determined from
curves presented in Figs. 4a) and 4b) differs too, see Fig 4c) for pressure
dependence of minimal values of effective temperatures.

Finally, a numerical parametric survey of optimal cooling conditions is
provided in Supplemental IV B. The main conclusions are as follows. First, for
a fixed total charge, $|q_{1}|+|q_{2}|$, cooling is most efficient when all
the charge is located on one particle. This is true whether we are cooling the
modes separately or together. In the former case, the cooling mechanism is
purely sympathetic (c.f. Arita  _et al._ (2022)), with the uncharged particle
cooled indirectly, via optical binding forces. Second, the CoM mode is most
efficiently cooled with $c_{\Sigma}=1$ and $c_{\Delta}=0$ while BR mode is
most effectively cooled when $c_{\Sigma}=0$ and $c_{\Delta}=1$. Third, for
symmetric systems with $q_{1}=q_{2}$ the breathing mode cannot be cooled in
the manner discussed here. Similarly, when $q_{1}=-q_{2}$ the CoM mode cannot
be cooled.

##  IV Conclusion

In summary, we have demonstrated cooling of the motion of optically
interacting levitated nanoparticles to sub-Kelvin effective temperatures using
a cold damping method employing an external electric field. We have discussed
the limitations of the cooling method using a rigorous theoretical model. For
individually cooled motional modes of coupled oscillators, the most stringent
limit is due to the detection efficiency and its noise, which is analogous to
the case of single levitated particles. In the case of simultaneous cooling of
the modes, the coupling between the modes by the feedback electric force is an
additional source of feedback-induced heating which limits the performance of
the method at pressures higher than $10^{-3}$ mbar. At lower pressures, the
detection noise is the dominant heating mechanism. Thus, in principle, the
cold damping method could get closer to the quantum ground state with an
optimized detection efficiency, operating at lower pressures. Motional cooling
of several interacting levitated nanoparticles is a significant step towards
the generation of cooled arrays of optically levitated particles, which will
allow multi-particle studies at the boundary of classical and quantum physics.

## References

  * Ranjit  _et al._ (2016) G. Ranjit, M. Cunningham, K. Casey, and A. A. Geraci, Phys. Rev. A 93, 053801 (2016).
  * Hempston  _et al._ (2017) D. Hempston, J. Vovrosh, M. Toroš, G. Winstone, M. Rashid, and H. Ulbricht, Applied Physics Letters 111, 133111 (2017).
  * Hebestreit  _et al._ (2018a) E. Hebestreit, M. Frimmer, R. Reimann, and L. Novotny, Phys. Rev. Lett. 121, 063602 (2018a).
  * Svak  _et al._ (2018) V. Svak, O. Brzobohatý, M. Šiler, P. Jákl, J. Kaňka, P. Zemánek, and S. H. Simpson, Nature Commun. 9, 5453 (2018).
  * Rondin  _et al._ (2017) L. Rondin, J. Gieseler, F. Ricci, R. Quidant, C. Dellago, and L. Novotny, Nature Nanotech 12, 1130 (2017).
  * Delić  _et al._ (2020) U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Science 367, 892 (2020).
  * Magrini  _et al._ (2021) L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, Nature 595, 373 (2021).
  * Tebbenjohanns  _et al._ (2021) F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Nature 595, 378 (2021).
  * Rieser  _et al._ (2022) J. Rieser, M. A. Ciampini, H. Rudolph, N. Kiesel, K. Hornberger, B. A. Stickler, M. Aspelmeyer, and U. Delić, Science 377, 987 (2022).
  * Marletto and Vedral (2017) C. Marletto and V. Vedral, Phys. Rev. Lett. 119, 240402 (2017).
  * Zhao  _et al._ (2012) R. Zhao, A. Manjavacas, F. J. García de Abajo, and J. B. Pendry, Phys. Rev. Lett. 109, 123604 (2012).
  * Moore and Geraci (2021) D. C. Moore and A. A. Geraci, Quantum Sci. Technol. 6, 014008 (2021).
  * Rudolph  _et al._ (2022) H. Rudolph, U. c. v. Delić, M. Aspelmeyer, K. Hornberger, and B. A. Stickler, Phys. Rev. Lett. 129, 193602 (2022).
  * Bykov  _et al._ (2018) D. S. Bykov, S. Xie, R. Zeltner, A. Machnev, G. K. L. Wong, T. G. Euser, and P. S. J. Russell, Light: Science & Applications 7, 22 (2018).
  * Arita  _et al._ (2018) Y. Arita, E. M. Wright, and K. Dholakia, Optica 5, 910 (2018).
  * Svak  _et al._ (2021) V. Svak, J. Flajšmanová, L. Chvátal, M. Šiler, A. Jonáš, J. Ježek, S. H. Simpson, P. Zemánek, and O. Brzobohatý, Optica 8, 220–229 (2021).
  * Simmons (1996) R. M. Simmons, Biophys. J. 70, 1813 (1996).
  * Kamba  _et al._ (2022) M. Kamba, R. Shimizu, and K. Aikawa, Opt. Express 30, 26716 (2022).
  * Vijayan  _et al._ (2022) J. Vijayan, Z. Zhang, J. Piotrowski, D. Windey, F. van der Laan, M. Frimmer, and L. Novotny, Nat. Nanotechnol. , 1 (2022).
  * Myatt  _et al._ (1997) C. J. Myatt, E. A. Burt, R. W. Ghrist, E. A. Cornell, and C. E. Wieman, Phys. Rev. Lett. 78, 586 (1997).
  * Larson  _et al._ (1986) D. J. Larson, J. C. Bergquist, J. J. Bollinger, W. M. Itano, and D. J. Wineland, Phys. Rev. Lett. 57, 70 (1986).
  * Arita  _et al._ (2022) Y. Arita, G. D. Bruce, E. M. Wright, S. H. Simpson, P. Zemánek, and K. Dholakia, Optica 9, 1000 (2022).
  * Bykov  _et al._ (2023) D. S. Bykov, L. Dania, F. Goschin, and T. E. Northup, Optica 10, 438 (2023).
  * Penny  _et al._ (2023) T. W. Penny, A. Pontin, and P. F. Barker, Phys. Rev. Res. 5, 013070 (2023).
  * Wei  _et al._ (2016) M.-T. Wei, J. Ng, C. T. Chan, and H. D. Ou-Yang, Scientific Reports 6, 38883 (2016).
  * Burns  _et al._ (1990) M. M. Burns, J.-M. Fournier, and J. A. Golovchenko, Science 249, 749 (1990).
  * Tatarkova  _et al._ (2002) S. A. Tatarkova, A. E. Carruthers, and K. Dholakia, Phys. Rev. Lett. 89, 283901 (2002).
  * Singer  _et al._ (2003) W. Singer, M. Frick, S. Bernet, and M. Ritsch-Marte, J. Opt. Soc. Am. B 20, 1568 (2003).
  * Dholakia and Zemánek (2010) K. Dholakia and P. Zemánek, Rev. Mod. Phys. 82, 1767 (2010).
  * Čižmár  _et al._ (2010) T. Čižmár, L. C. D. Romero, K. Dholakia, and D. L. Andrews, J. Phys. B 43, 102001 (2010).
  * Demergis and Florin (2012) V. Demergis and E.-L. Florin, Nano Lett. 12, 5756 (2012).
  * Yan  _et al._ (2014) Z. Yan, S. K. Gray, and N. F. Scherer, Nat Commun 5, 3751 (2014).
  * Li  _et al._ (2018) Y. Li, H. Xin, Y. Zhang, H. Lei, T. Zhang, H. Ye, J. J. Saenz, C.-W. Qiu, and B. Li, ACS Nano 12, 10703 (2018).
  * Čižmár  _et al._ (2011) T. Čižmár, O. Brzobohatý, K. Dholakia, and P. Zemánek, Laser Phys. Lett. 8, 50 (2011).
  * Brzobohatý  _et al._ (2011) O. Brzobohatý, V. Karásek, T. Čižmár, and P. Zemánek, Appl. Phys. Lett. 99, 101105 (2011).
  * Tebbenjohanns  _et al._ (2019) F. Tebbenjohanns, M. Frimmer, A. Militaru, V. Jain, and L. Novotny, Phys. Rev. Lett. 122, 223601 (2019).
  * Penny  _et al._ (2021) T. W. Penny, A. Pontin, and P. F. Barker, Phys. Rev. A 104, 023502 (2021).
  * Iwasaki  _et al._ (2019) M. Iwasaki, T. Yotsuya, T. Naruki, Y. Matsuda, M. Yoneda, and K. Aikawa, Phys. Rev. A 99, 051401 (2019).
  * Burns  _et al._ (1989) M. M. Burns, J.-M. Fournier, and J. A. Golovchenko, Phys. Rev. Lett. 63, 1233 (1989).
  * Hebestreit  _et al._ (2018b) E. Hebestreit, M. Frimmer, R. Reimann, C. Dellago, F. Ricci, and L. Novotny, Review of Scientific Instruments 89, 033111 (2018b).


#  Coping with Mistreatment in Fair Algorithms

Ankit Kulshrestha
Computer and Information Sciences
University of Delaware
Newark, DE
<EMAIL_ADDRESS>&Ilya Safro
Computer and Information Sciences
University of Delaware
Newark, DE
<EMAIL_ADDRESS>
(February 2021)

###### Abstract

Machine learning actively impacts our everyday life in almost all endeavors
and domains such as healthcare, finance, and energy. As our dependence on the
machine learning increases, it is inevitable that these algorithms will be
used to make decisions that will have a direct impact on the society spanning
all resolutions from personal choices to world-wide policies. Hence, it is
crucial to ensure that (un)intentional bias does not affect the machine
learning algorithms especially when they are required to take decisions that
may have unintended consequences. Algorithmic fairness techniques have found
traction in the machine learning community and many methods and metrics have
been proposed to ensure and evaluate fairness in algorithms and data
collection.

In this paper, we study the algorithmic fairness in a supervised learning
setting and examine the effect of optimizing a classifier for the Equal
Opportunity metric. We demonstrate that such a classifier has an increased
false positive rate across sensitive groups and propose a conceptually simple
method to mitigate this bias. We rigorously analyze the proposed method and
evaluate it on several real world datasets demonstrating its efficacy.

Reproducibility: All source code, and experimental results are available at
https://anonymous.4open.science/r/b6c6653c-5f50-477c-bb13-7dfdeb39d4f4/

##  1 Introduction

Machine learning has permeated almost every sphere of human endeavours. The
algorithms are being applied on both macro and micro scale in various degrees
in diverse fields like space [Kothari et al., 2020], finance [Dixon et al.,
2020]. Moreover, our dependence on machine learning algorithms to take
decisions has created an urgent need to ensure that the algorithms are
_unbiased_ when taking decisions, more so if that decision affects society.
There has been a steady growth of work in algorithmic fairness in the recent
times, e.g., [Ahmadian et al., 2019, 2020, Chierichetti et al., 2017, Dwork et
al., 2017, Pessach and Shmueli, 2020]. The underlying assumption in most of
the algorithmic fairness works is that there exists a unique sensitive
attribute and it influences the decision made by the algorithm. However, in
real world especially healthcare the problem is not so simple. For instance,
in the healthcare data the number of sick people is usually significantly
lower than that of the healthy people Dua et al. [2014] in many data samples.
However, this statistical difference is not just influenced by that sensitive
label (a diagnosis code in the electronic health records) or any single other
such as the gender of the patient but also on whether a patient has access to
medical insurance, or whether a patient would want another surgery. In this
situation, it becomes hard to not just quantify a fairness metric but also
design algorithms that discover and optimize such domain specific metrics.

Figure 1: The True Positive Rates and False Positive Rates of Vanilla SVM vs
Fair-Empirical Risk Minimization SVM.

For machine learning algorithms to be completely fair, they need to
simultaneously satisfy the dual objectives of not preferring a statistical
majority over a statistical minority to make a final decision and at the same
time ensuring that the statistical majority is not unfairly treated to satisfy
some constraint. In other words, a fair machine learning algorithm should be
able to discriminate between different outcomes given a data sample without
overly relying on a single sensitive attribute or a combination thereof.

Many metrics such as Disparate Impact [Barocas and Selbst, 2016, Chierichetti
et al., 2017, Ahmadian et al., 2020, Bera et al., 2019], Demographic Parity
[Zafar et al., , Zliobaite, 2015, Calders et al., 2009], and Equal Odds/Equal
Opportunity [Hardt et al., 2016, Donini et al., 2020, Zafar et al., 2017] have
been proposed to ensure the first fairness requirement but they haven’t been
studied to see if they also satisfy the second requirement as well.

We focus on the supervised learning paradigm and assume that we have access to
an underlying “protected attribute” (we also refer to it as a “sensitive
attribute”). Two important fairness measures in this paradigm are the Equal
Odds and Equal Opportunity measures [Hardt et al., 2016]. The Equal Odds
requires that a classifier has the same true positive rate (TPR) and false
positive rate (FPR) across sensitive groups, i.e.,

$|P[\hat{y}=1|s=0,y=1]-P[\hat{y}=1|s=1,y=1]|\leq\epsilon$ (1)

and

$|P[\hat{y}=1|s=0,y=0]-P[\hat{y}=0|s=1,y=0]|\leq\epsilon$ (2)

where $\hat{y}$ is the output of the classifier, $y$ is the ground truth,
$\epsilon$ is a user defined fairness criterion threshold and $s$ is an
indicator variable defining the sensitive group. For our purposes,
$s\in\\{0,1\\}$. Similarly, Equal Opportunity requires that

$|P[\hat{y}=1|s=0,y=1]-P[\hat{y}=1|s=1,y=1]|\leq\epsilon.$ (3)

It has been shown in [Pleiss et al., 2017] that the only way to satisfy the
Equal Odds metric is to set $\epsilon=0$. This requirement severely limits the
areas in which the Equal Odds metric can be a useful measure for ensuring
fairness. The Equal Opportunity can be seen as a relaxed version of Equalized
Odds. When $\epsilon=0$, Equal Opportunity requires that the true positive
rates across sensitive groups are same. Figure 1 shows the true positive and
false positive rates of an SVM model trained without fairness constraint and
another SVM trained with Equal Opportunity (EO) constrain across the sensitive
groups on two datasets - COMPAS recidivism dataset [Larson et al., 2016]
released by ProPublica and the Drug dataset [Fehrman et al., 2017]. We can see
that in both datasets, the false positive rates increase across the sensitive
groups. Moreover, the increase in TPR of a statistical minority group comes at
the expense of a statistical majority group. We name this phenomenon as
_disparate mistreatment_ (different from Zafar et al. ). A classifier trying
to ensure equal opportunity can make myopic decisions by inadvertently
associating the sensitive group as the most important feature for final
prediction. For instance, in the healthcare domain the data is often extremely
imbalanced [Razzaghi et al., 2019, 2015] due to such reasons as much smaller
number of sick people than healthy, unequal access to medical insurance, and
incorrect sampling among the people seeking for medical help. A trained
classifier with Equal Opportunity metric may falsely output more positive
responses on the group forming a minority in the overall disease rate. This
bias would then violate our requirement on fairness.

###  1.1 Our Contributions

In this paper, we note the existence of bias exhibited by classifiers
optimizing EO metric _against_ majority groups belonging to a single protected
attribute. We make the following contributions with our work:

  * •

We concretely define the disparate mistreatment phenomenon presented above and
propose a conceptually simple criteria to prevent the disparate mistreatment.

  * •

We rigorously show that minimizing the optimization objective subject to our
constraint can generalize to unseen samples.

  * •

We quantify the performance of classifier in terms of the DFPR metric
(difference of False Positive Rates) and show that a classifier trained with
our approach exhibits a significantly less DFPR when compared against other
baselines.

##  2 Related Work

Work on algorithmic fairness has progressed on different areas of machine
learning and different algorithms have been proposed in both supervised and
unsupervised settings. In this paper we limit ourselves to work concerning
statistical or group fairness. For an excellent introduction to algorithmic
fairness in general, please consider [Pessach and Shmueli, 2020] and the
references therein.

In the supervised learning paradigm, the proposed algorithms achieve fairness
by either data-preprocessing, optimizing a model satisfying a particular
criteria or fine tuning an unoptimal classifier to meet fairness criteria. For
instance, [Hardt et al., 2016] propose the Equal Odds and Equal Opportunity
fairness metric and first train an unoptimal (in the fairness sense)
classifier. Then the authors fine tune the classifier keeping in mind the
protected attribute of the dataset.

In a different approach [Donini et al., 2020] propose the FERM framework to
optimize the Equal Opportunity metric by integrating the fairness constraint
in an Empirical Risk Management (ERM) framework. [Zafar et al., 2017] propose
a framework to mitigate disparate mistreatment by measuring the distance of
users feature vectors that are misclassified with the decision boundary. The
authors then propose to mitigate disparate impact(DI) [Zafar et al., ] by
measuring the distance between a users feature vector and the decision
boundary. We note that our definition of disparate mistreatment is different
from [Zafar et al., 2017] and also note that their method assumes linear
separability of data while we do not make any such assumptions about our data.
Our work is more closely related to Donini et al. [2020] since we also use the
ERM framework to integrate fairness in the algorithm.

Another direction to ensure fairness in supervised setting is proposed by
[Zemel et al., 2013] in which a _new_ data representation is derived from the
original existing data in a way that removes the protected attribute from
participating directly into the decision process. An ensemble learning
approach is proposed by [Dwork et al., 2017] in which decoupled classifiers
are trained and each classifier is fine tuned jointly by minimizing a joint
loss.

Fairness and calibration are closely related. [Pleiss et al., 2017] define
calibration and rigorously show that a classifier cannot simultaneously
satisfy Equalized Odds and output calibrated responses. Moreover, it has been
shown that the only way to ensure Equalized Odds is to ensure that both true
positive rates and false positive rates across two sensitive groups are
exactly the same.

In the unsupervised learning domain, fairness has been considered within the
framework of clustering. Following the work of [Chierichetti et al., 2017] on
fairlets, different methods have been proposed [Ahmadian et al., 2019, 2020,
Chhabra and Mohapatra, 2020, Backurs et al., 2019] that either solve the
clustering problem in a heirarchical agglomerative setting or ensuring
fairness in a k-means clustering setting. The central idea in all works in
this vein is to prevent an over-representation of a particular sensitive group
in a cluster.

##  3 Our Approach

We begin with the notation. Let $\mathcal{D}=\\{(x_{i},y_{i})\\}^{n}_{i=1}$ be
a labelled dataset of $n$ samples, where $x_{i}$ and $y_{i}$ represent the
data point and label, respectively. Let $\cal{F}$ be the space of all
functions that minimize a prescribed loss function
$l(f(x_{i}),y_{i}),~{}\forall f\in\cal{F}$. We define $f_{a}$ to be a
classifier that optimizes accuracy and $f_{eo}$ to be a classifier that
minimizes $l(\cdot,\cdot)$ taking into account fairness using the Equal
Oportunity (EO) constraint. Further, we constrain the target classes
$\mathcal{Y}\in\\{-1,+1\\}$. For any classifier, let

$\mathbb{P}_{g}(\hat{y},y)=P[\hat{y}|s=g,y]$ (4)

be the probability of predicting $\hat{y}$ for sensitive group $g$ and target
$y\in\mathcal{Y}$. We further assume that there are only two classes of
protected attribute hence $g\in\\{a,b\\}$. We are now ready to define the
disparate mistreatment phenomenon.

######  Definition 3.1 (Disparate Mistreatment).

A classifier $f$ is said to exhibit disparate mistreatment if and only if for
groups $a,b$ an increase in $\mathbb{P}_{b}(1,1)$ directly corresponds to a
decrease in $\mathbb{P}_{a}(1,1)$ and vice versa.

It follows from Definition 3.1 that any classifier $f$ exhibiting disparate
mistreatment will have an increased rate of false negatives for one
statistical group and at the same time also result in an increased false
positive rate in the other. For any $f\in\mathcal{F}$, optimizing for
maximizing accuracy will return $f_{a}$ which tries to find the decision
boundary that differentiates between positive and negative samples conditioned
on the outcome alone.

However, optimizing EO results in $f_{eo}$ which constrains the classifer to
find a supoptimal decision boundary so that true positives increase across the
different sensitive groups. Geometrically, this can be seen as pushing the
clusters of positively labeled samples in either group closer to each other by
minimizing the distance between their respective means. This optimization
method however completely ignores the effect on the cluster of negatively
labeled samples in either sensitive group. We hypothesize that this is the
root cause of $f_{eo}$ exhibiting disparate mistreatment.

Our approach is to introduce a _minimum separation parameter_ into the
Empirical Risk Minimization (ERM) framework. We build on the work of [Donini
et al., 2020] due its attractive guarantees on upper bounds of the empirical
risk.

######  Definition 3.2.

_Let_

$R^{-,g}(f)=\mathbb{E}([l(f(x_{i}),y_{i})\mid y=-1,s=g])$

_be the risk on negatively labeled samples for a particular sensitive group in
a dataset. Then $f$ is $\rho$-minimum separated if_

$|R^{-,a}(f)-R^{-,b}(f)|\geq\rho.$ (5)

Based on Definition 5, we can formulate two constrained ERM problems. In the
case where we want to obtain $f$ such that the decision boundary respects a
minimal separation between negative label clusters we optimize:

$\displaystyle\underset{f\in\mathcal{F}}{\operatorname{\mathop{\mathrm{argmin}}\limits}}$
$\displaystyle R(f)$ (6) such that
$\displaystyle|R^{-,a}(f)-R^{-,b}(f)|\geq\rho$

We call Problem 6 as _Equal Treatment_(ET). In this case, we do not constrain
the optimal function to ensure equal true positive rates across sensitive
groups. Instead, fairness is achieved by ensuring that negative samples from
any sensitive group do not get classified as positive up to the tolerance
level $\rho$. We additionally require that $0<\rho\leq 1$. However, there may
be situations where we would want both EO and ET constraints to be satisfied.
In those cases we optimize:

$\displaystyle\underset{f\in\mathcal{F}}{\operatorname{\mathop{\mathrm{argmin}}\limits}}$
$\displaystyle R(f)$ (7) such that
$\displaystyle|R^{-,a}(f)-R^{-,b}(f)|\geq\rho$
$\displaystyle|R^{+,a}(f)-R^{+,b}(f)|\leq\epsilon$

We have defined the optimization problems 7 and 6 in a general manner. We now
analyze the bounds on the performance of our algorithm. First, we define
$R(f)$ to be the _true risk_ of the unknown underlying data distribution $P$.
Further,

$\hat{R}(f)=\frac{1}{n}\sum_{i=1}^{n}l_{c}(f(x_{i},y_{i}))$

is the _empirical risk_ associated with samples $(x_{i},y_{i})$ drawn i.i.d
from $\mathcal{D}^{n}$. We call $f^{*}$ to be the Bayes optimal classifier
minimizing true risk and $\hat{f}$ to be optimal classifier minimizing the
empirical risk. Furthermore, we constrain the loss function $l_{c}$ to be
convex.

######  Lemma 3.1.

If $f\in\mathcal{F}$ be a function minimizing empirical risk subject to
constraints in Problem 6 and the loss $l_{c}$ is convex. If
$C(\mathcal{F},\delta,n^{-})$ is the upper bound function on empirical risk,
then with probability of at least $1-4\delta$,

$\displaystyle|R^{-,a}(f^{*})-R^{-,b}(f^{*})|<$
$\displaystyle|R^{-,a}(\hat{f})-R^{-,b}(\hat{f})|+$ (8)
$\displaystyle\rho+2C(\mathcal{F},\delta,n^{-}).$

###### Proof.

For a given upper bound function $C(\mathcal{F},\delta,n^{-})$, the true risk
and empirical risk are related as:

$R(\hat{f})<\hat{R}(f)+C(\mathcal{F},\delta,n)$ (9)

with probability greater or equal $(1-\delta)$. Using Equation 9, we can that
with probability of atleast $1-2\delta$, where $\delta$ is the confidence
parameter that

$|R^{-,a}(\hat{f})-R^{-,b}(\hat{f})|<|\hat{R}^{-,a}(f)-\hat{R}^{-,b}(f)|+C(\mathcal{F},\delta,n).$
(10)

But, we have defined $|\hat{R}^{-,a}(f)-\hat{R}^{-,b}(f)|\geq\rho$. Then it
follows that,

$\displaystyle|R^{-,a}(\hat{f})-R^{-,b}(\hat{f})|<$
$\displaystyle|\hat{R}^{-,a}(f)-\hat{R}^{-,b}(f)|+$ (11) $\displaystyle
C(\mathcal{F},\delta,n^{-})+\rho.$

Essentially, Equation 11 means that the true risk of a classifier that
minimizes ERM is upper bounded by the cost function and the minimal separation
parameter $\rho$. However, we are interested in bounding the risk of the
optimal classifier on the training data (observed samples) in terms of the
Bayes optimal classifier $f^{*}$, so the following relation is applied

$\hat{R}(f^{*})<R(f^{*})+C(\mathcal{F},\delta,n)$ (12)

with probability greater or equal $(1-\delta)$. Using Equation 12 and noticing
that the right hand side of Equation 11 is valid $\forall f\in\mathcal{F}$, we
can state that with probability at least $1-4\delta$

$\begin{split}|R^{-,a}(\hat{f})-R^{-,b}(\hat{f})|<|R^{-,a}(f^{*})-R^{-,b}(f^{*})|+\\\
\rho+2*C(\mathcal{F},\delta,n^{-}).\end{split}$ (13)

∎

The result of Lemma 3.1 upper bounds the empirical risk of the optimal and can
guarantee that with a small penalty of $\rho$ the empirical risk will be close
to the true risk on the underlying distribution. However, the degree of
closeness will depend on how fast $C(\mathcal{F},\delta,n)$ converges to a
uniform value. If $\mathcal{F}$ is finite we can use Union bound to derive the
upper bound function as:

$C(\mathcal{F},\delta,n)=\sqrt{\frac{log(\mathcal{|F|})+log(1/\delta)}{2n}}$
(14)

On the other hand, if $\mathcal{F}$ is infinite then we can make use of
Rademacher bound to find an upper bound for the empirical risk. The key point
to note is that $C(\mathcal{F},\delta,n)$ decays to $0$ as
$n\rightarrow\infty$ which implies that the empirical risk can approximate the
true risk within a small penalty of $\rho$.

###  3.1 Non-Linear Optimization with Kernels

In this section we specify our problem when the underlying space of functions
is the Reproducing Kernel Hilbert Space (RKHS). For any function $f$ in RKHS
we parametrize it as $f(x)=\langle\textbf{w},\phi(\textbf{x})\rangle$ where
$w\in\mathbb{H}$ and $\phi(\textbf{x})$ is a feature map in high dimensional
space. For our purpose, we assume $\phi(\textbf{x})$ is a _non-linear_ mapping
in high dimensional space, i.e., $\phi:\mathcal{X}\rightarrow\mathbb{H}$.
Since our underlying function space is RKHS, we can define a positive
semidefinite kernel matrix
$K:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$, such that
$K(x,x^{\prime})=\langle\phi(x),\phi(x^{\prime})\rangle$.

Following [Donini et al., 2020], we define $\textbf{u}_{g}$ to be the
barycenter of a negatively labeled cluster of points belonging to a sensitive
group $g$:

$\textbf{u}_{g}=\frac{1}{n^{-,g}}\sum_{i\in\mathcal{Y}^{-},g=g}^{n}\phi(x_{i})$
(15)

In general, we optimize the following problem

$\displaystyle\min_{\textbf{w}}$
$\displaystyle\sum_{i=1}^{m}l(\langle\textbf{w},\phi(x_{i})\rangle,y_{i})+||\textbf{w}||_{2}^{2}$
(16) such that
$\displaystyle\langle\textbf{w},\bar{\textbf{u}}\rangle\geq\rho$

In Equation 16, $\bar{\textbf{u}}=\textbf{u}_{a}-\textbf{u}_{b}$ . We further
specify this problem to a soft-margin SVM by optimizing

$\displaystyle\min_{\textbf{w}}$
$\displaystyle\textbf{w}^{T}\textbf{w}+C\sum_{i=1}^{n}\zeta_{i}$ (17) such
that $\displaystyle\zeta_{i}\geq 0$ $\displaystyle
y_{i}(1-\phi(x_{i})^{T}\textbf{w})\geq 1-\zeta_{i}$
$\displaystyle\langle\textbf{w},\bar{\textbf{u}}\rangle\geq\rho$

To simplify solving Equation 16, we use the Representer theorem [Kimeldorf and
Wahba, 1971] to note that $\textbf{w}=\sum_{i=1}^{n}\alpha_{i}\phi(x_{i})$,
then

$\langle\textbf{w},\phi(\textbf{x})\rangle=\sum_{i=1}^{n}\alpha_{i}\langle\phi(x),\phi(x^{\prime})\rangle,$
(18)

which can be simply expressed as

$\langle\textbf{w},\phi(\textbf{x})\rangle=\sum_{i=1}^{n}\alpha_{i}k(x_{i},x^{\prime}_{i}),$
(19)

where $k(x_{i},x^{\prime}_{i})$ is an element of the gram matrix $K$. We also
write,
$||\textbf{w}||_{2}^{2}=\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}K(x_{i},x_{j})$.
The optimization problem can be restated as:

$\displaystyle\min_{\alpha}$
$\displaystyle\sum_{i=1}^{n}l(\alpha_{i}*K(x_{i},x^{\prime}_{i}),y_{i})+\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}K(x_{i},x_{j})$
(20) such that
$\displaystyle\sum_{i=1}^{n}\alpha_{i}\left[\frac{1}{n^{-,a}}\sum_{y\in\mathcal{Y}^{-}}^{n}K_{iy}-\frac{1}{n^{-,b}}\sum_{y\in\mathcal{Y}^{-}}^{n}K_{iy}\right]\geq\rho$

Geometric Intuitions: We stated earlier that $\textbf{w}\in\mathbb{H}$ where
$\mathbb{H}$ is the Hilbert space. Since Hilbert space is a complete vector
space, we can interpret any $\langle\textbf{x},\textbf{y}\rangle$ as a
distance between two infinite dimensional vectors embedded in that space.
Moreover, since the underlying function space is a RKHS, any functional will
be upper bounded by a finite quantity. If we solve Equation 16 without the
constraint, then we simply are choosing a weight vector which is closest to
points $\phi(x_{i})$ that minimize the loss. However, our added constraint
forces the algorithm to search for a weight vector that minimizes the loss
_while_ maintaining a distance of at-least $\rho$ between the center of masses
of points that have negative labels associated with them.

##  4 Experimental results

Dataset | # Data points | # Features | Sensitive Attribute
---|---|---|---
Adult |  $32561$/$16281$ | $12$ | Gender
COMPAS | $6171$ | $10$ | Race
Drug | $1889$ | $11$ | Race
Arrhythmia | $452$ | $279$ | Gender

Table 1: Statistics of Real World Datasets.

We now describe our experimental results on several real-world datasets and
show the effectiveness of our approach in reducing the FPR across the
sensitive groups. Table 1 summarizes the datasets.

Adult is a dataset from UCI repository [Dua and Graff, 2017] containing 14
features with 48842 data points about different aspects of an individual (such
as age, gender, and marital status) to determine if that person earns more
than $50,000 a year. For this dataset we consider _gender_ to be the sensitive
attribute.

COMPAS Correctional Offender Management Profiling for Alternative Sanctions
(COMPAS) algorithm [Larson et al., 2016] used by Broward County Prison to
assess the re-arrest rate of different individuals was shown to be biased
against black defendants. We used the dataset [Larson et al., 2016] provided
by Propublica and specifically focused on the violent-recidivism subset. The
dataset contains 6,171 data points with 12 features. We chose _race_ as the
sensitive feature.

Drug is a dataset provided by UCI [Dua and Graff, 2017] containing information
about 1889 individuals. Each individual has 12 features associated with them
including personality type, ethnicity and whether they have used a drug in the
past. There are 16 possible classification problems concerning use of
different drugs. In this paper we select the use of cocaine as our target
variable and discretize it into “Never used” and “Used” classes. We further
discretize the sensitive attribute _ethnicity_ into “White” and “Other”
classes.

Arrhythmia is a dataset containing 452 samples and 279 features to detect the
absence or presence of cardiac arrhythmia in a patient. We also turn this into
a binary classification problem by considering “Normal” vs “Others” as the
outcome variable. We choose _gender_ as the protected attribute. Since this
data contains a lot of missing values we pre-processed the data by dropping
feature with number of unknowns greater than a threshold and substituting the
rest of the unknowns with the median value of the feature. This dataset is
also provided by the UCI [Dua and Graff, 2017] repository.

We use the RBF kernel

$K(\textbf{x},\textbf{y})=e^{\gamma||\textbf{x}-\textbf{y}||^{2}}$

for all experiments. All real valued features are scaled to have a zero mean
and unit standard deviation. We further vectorize all categorical variables as
well. In datasets where the test/train split is not provided we search for the
optimal hyperparameters ($C,\gamma$) using 5-fold cross validation for each
given value of $\rho$ and report the best result. For baseline, we use the
code provided by [Donini et al., 2020] with our pre-processed data. In all
tables containing the results, our method is labeled by SVM-MT, the original
FERM and SVM by FERM and SVM, respectively.

Method | Adult | COMPAS | Drug | Arrhythmia
---|---|---|---|---
| Precision | DFPR | Precision | DFPR | Precision | DFPR | Precision | DFPR
SVM | 0.74 | 0.086 | $0.668\pm 0.009$ | $\textbf{0.217}\pm 0.058$ | $0.672\pm 0.027$ | $0.098\pm 0.066$ | $0.747\pm 0.061$ | $\textit{0.158}\pm 0.103$
SVM-MT | 0.73 | 0.020 | $0.664\pm 0.023$ | $\textit{0.225}\pm 0.058$ | $0.625\pm 0.017$ | $\textbf{0.033}\pm 0.026$ | $0.703\pm 0.072$ | $\textbf{0.056}\pm 0.04$
FERM | 0.71 | 0.067 | $0.667\pm 0.015$ | $0.25\pm 0.053$ | $0.627\pm 0.023$ | $\textit{0.069}\pm 0.078$ | $0.742\pm 0.048$ | $0.21\pm 0.124$
Table 2: A comparison of results on four datasets with baseline FERM method
and optimizing for minimal separation. Best results are shown in bold and
second best are in italics.

In order to measure the performance of the classifier after optimization we
define the _DFPR_ metric to be the absolute difference between the FPR across
the sensitive groups. Our search for optimal $\rho$ was performed by first
optimizing the SVM specific hyperparameters via cross fold validation. When we
do not explicitly balance the dataset, then we sweep the values of $\rho$ in
$[0.1,\dots,0.5]$ in steps of $0.1$ to determine the optimal value that
minimizes the DFPR.

In our experiments, we did not notice a significant decrease of DFPR value for
$\rho>0.5$. We explain this as follows: when $\rho$ is higher then a maximum
margin classifier would effectively optimize for accuracy. Table 2 shows the
results of our experiments on four real world datasets. _We observe that our
method consistently outperforms the baseline FERM[Donini et al., 2020] method
achieving significantly lower DFPR on almost all datasets without sacrificing
too much precision_.

Model | Adult | COMPAS | Drug | Arrhythmia
---|---|---|---|---
SVM | $1\pm 1.5$ | $3\pm 1$ | $1\pm 1$ | $0.2\pm 0.5$
SVM-MT | $7\pm 1$ | $121\pm 1$ | $10\pm 1$ | $1\pm 0.6$
FERM | $10\pm 1$ | $136\pm 1$ | $16\pm 1$ | $2\pm 1$
Table 3: Model optimization times on different datasets. The results are an
average of 5 different runs. All times except for Adult dataset are in
seconds. Adult dataset time is measured in minutes.

For the completeness of comparison we also trained SVM models without any
constraints. As expected, the unconstrained SVM has a higher accuracy, but
also a higher rate of FPR than any other models. In particular, the FERM model
has a comparable precision with the baseline SVM but has a much higher DFPR
rate on the Arrhytmia dataset, while on the other hand our method sacrifices
some precision to achieve the lowest DFPR rate of $0.05$ with $\rho=0.1$.
Table 3 shows the average optimization time over 5 different runs for all
three models considered in our study. As expected, an unconstrained baseline
SVM is significantly faster than the other two methods in finding the decision
boundary. However, compared against the FERM model, our approach is
consistently faster. In the case of Adult dataset optimizing both FERM and
SVM-MT takes considerably longer than the baseline SVM. This indicates that
for larger datasets a more scalable algorithm needs to be designed.

Figure 2: Statistics about sensitive groups in real world datasets.

Model | Drug | Arrhythmia
---|---|---
| Precision | DFPR | Precision | DFPR
SVM | $0.672\pm 0.028$ | $0.12\pm 0.087$ | $0.744\pm 0.05$ | $0.251\pm 0.135$
SVM-MT | $0.661\pm 0.032$ | $0.074\pm 0.04$ | $0.703\pm 0.03$ | $0.09\pm 0.06$
FERM | $0.648\pm 0.017$ | $0.102\pm 0.05$ | $0.735\pm 0.03$ | $0.121\pm 0.07$

Table 4: Results of experiments with undersampling the majority label class on
Drug and Arrhythmia datasets. The results are an average of 5 different runs.

###  4.1 Handling Imbalanced Classes

The real world datasets have a class imbalance since data collection is not
always uniform. For instance, in the Drug dataset there are 784 people
belonging to sensitive group (“Other”) that are labeled as a positive class.
By contrast, only 63 people belonging to the other sensitive group(“White”)
are positively labeled. The same trend is also observed in other datasets. In
order to evaluate the algorithms more fairly, we wanted to remove the
imbalance in the class labels. Hence, we undersampled the numerically superior
class and ran the optimization against all three methods i.e. baseline SVM,
SVM-MT(our method) and FERM. Table 4 shows the results of our experiments. All
the results are the average of 5 independent runs. We found $\rho=0.1$ to be a
stable choice for all the experiments.

In the case of Drug dataset, FERM model has a significantly reduced precision,
while our method has almost comparable precision to the baseline. Moreover,
our DFPR is the lowest amongst all three methods. Given that the drug dataset
exhibits a strong class imbalance and FERM optimizes to minimize the TPR,
undersampling the majority decreases the precision. In the Arrhythmia dataset
we also note that our method achieves the best DFPR. However, for the FERM
model the precision drops less due to the fact that difference between
positively labeled samples is not as pronounced (160 in one sensitive group
and 85 in the other). We also observe that our method does exhibit a lower
precision on the Arrhythmia dataset which we interpret to be because of the
low number of samples in the data and our objective to keep the negatively
labeled samples separated by atleast $\rho$.

##  5 Conclusion

In this paper we have studied the effects of optimizing a classifier under the
Equal Opportunity constraint. We showed that a classifier trying to minimize
the difference between true positive rates of sensitive groups causes an
increase in the false positive rates within each group. Furthermore, the
increase in true positive rate of the statistical minority is achieved by the
decrease in the true positive rate of the statistical majority. While this may
be a valid method to ensure fairness with a single protected attribute, we
reasoned that it would be counterproductive if the sensitive attribute could
not be precisely defined. To mitigate this mistreatment we proposed a minimum
separation parameter and showed that our approach is successful in decreasing
the false positive rates across the sensitive groups without sacrificing too
much accuracy.

Our work opens up a discussion about how far the assumption of a unique
sensitive attribute in data is valid. A fruitful direction of work can be to
examine the existing fairness metrics when the sensitive attribute is not
precise. Another direction that opens up is to define metrics that are robust
to imprecise sensitive attributes and create algorithms that make minimal
assumptions about the underlying data. Another important direction is related
to introducing the Equal Opportunity constraint into the recent scalable
accelerators of nonlinear support vector machines such as [Wen et al., 2018,
Sadrfaridpour et al., 2019] which are expected to provide a notable
scalability improvement to the entire optimization process.

Reproducibility: All source code, and experimental results are available at
https://anonymous.4open.science/r/b6c6653c-5f50-477c-bb13-7dfdeb39d4f4/

## References

  * Ahmadian et al. [2019] Sara Ahmadian, Alessandro Epasto, Ravi Kumar, and Mohammad Mahdian.  Clustering without Over-Representation.  _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 267–275, July 2019.  doi: 10.1145/3292500.3330987.  URL http://arxiv.org/abs/1905.12753.  arXiv: 1905.12753.
  * Ahmadian et al. [2020] Sara Ahmadian, Alessandro Epasto, Marina Knittel, Ravi Kumar, Mohammad Mahdian, Benjamin Moseley, Philip Pham, Sergei Vassilvitskii, and Yuyan Wang.  Fair Hierarchical Clustering.  _arXiv:2006.10221 [cs, stat]_ , June 2020.  URL http://arxiv.org/abs/2006.10221.  arXiv: 2006.10221.
  * Backurs et al. [2019] Arturs Backurs, Piotr Indyk, Krzysztof Onak, Baruch Schieber, Ali Vakilian, and Tal Wagner.  Scalable Fair Clustering.  _arXiv:1902.03519 [cs]_ , June 2019.  URL http://arxiv.org/abs/1902.03519.  arXiv: 1902.03519.
  * Barocas and Selbst [2016] Solon Barocas and A. Selbst.  Big data’s disparate impact.  _California Law Review_ , 104:671, 2016.
  * Bera et al. [2019] Suman Bera, Deeparnab Chakrabarty, Nicolas Flores, and Maryam Negahbani.  Fair algorithms for clustering.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper/2019/file/fc192b0c0d270dbf41870a63a8c76c2f-Paper.pdf.
  * Calders et al. [2009] T. Calders, F. Kamiran, and M. Pechenizkiy.  Building classifiers with independency constraints.  In _2009 IEEE International Conference on Data Mining Workshops_ , pages 13–18, 2009.  doi: 10.1109/ICDMW.2009.83.
  * Chhabra and Mohapatra [2020] Anshuman Chhabra and Prasant Mohapatra.  Fair Algorithms for Hierarchical Agglomerative Clustering.  _arXiv:2005.03197 [cs, stat]_ , July 2020.  URL http://arxiv.org/abs/2005.03197.  arXiv: 2005.03197.
  * Chierichetti et al. [2017] Flavio Chierichetti, Ravi Kumar, Silvio Lattanzi, and Sergei Vassilvitskii.  Fair Clustering Through Fairlets.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 5029–5037. Curran Associates, Inc., 2017.  URL http://papers.nips.cc/paper/7088-fair-clustering-through-fairlets.pdf.
  * Dixon et al. [2020] Matthew F Dixon, Igor Halperin, and Paul Bilokon.  _Machine Learning in Finance_.  Springer, 2020.
  * Donini et al. [2020] Michele Donini, Luca Oneto, Shai Ben-David, John Shawe-Taylor, and Massimiliano Pontil.  Empirical Risk Minimization under Fairness Constraints.  _arXiv:1802.08626 [cs, stat]_ , January 2020.  URL http://arxiv.org/abs/1802.08626.  arXiv: 1802.08626.
  * Dua and Graff [2017] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Dua et al. [2014] Sumeet Dua, U Rajendra Acharya, and Prerna Dua.  _Machine learning in healthcare informatics_ , volume 56.  Springer, 2014.
  * Dwork et al. [2017] Cynthia Dwork, Nicole Immorlica, Adam Tauman Kalai, and Max Leiserson.  Decoupled classifiers for fair and efficient machine learning.  _arXiv:1707.06613 [cs]_ , July 2017.  URL http://arxiv.org/abs/1707.06613.  arXiv: 1707.06613.
  * Fehrman et al. [2017] E. Fehrman, A. K. Muhammad, E. M. Mirkes, V. Egan, and A. N. Gorban.  The five factor model of personality and evaluation of drug consumption risk, 2017.
  * Hardt et al. [2016] Moritz Hardt, Eric Price, and Nathan Srebro.  Equality of Opportunity in Supervised Learning.  _arXiv:1610.02413 [cs]_ , October 2016.  URL http://arxiv.org/abs/1610.02413.  arXiv: 1610.02413.
  * Kimeldorf and Wahba [1971] George Kimeldorf and Grace Wahba.  Some results on tchebycheffian spline functions.  _Journal of mathematical analysis and applications_ , 33(1):82–95, 1971.
  * Kothari et al. [2020] Vivek Kothari, Edgar Liberis, and Nicholas D Lane.  The final frontier: Deep learning in space.  In _Proceedings of the 21st International Workshop on Mobile Computing Systems and Applications_ , pages 45–49, 2020.
  * Larson et al. [2016] Jeff Larson, Marjorie Roswell, and Vaggelis Atlidakis.  Correctional Offender Management Profilingfor Alternative Sanctions (COMPAS) dataset.  https://github.com/propublica/compas-analysis/, 2016.
  * Pessach and Shmueli [2020] Dana Pessach and Erez Shmueli.  Algorithmic Fairness.  _arXiv:2001.09784 [cs, stat]_ , January 2020.  URL http://arxiv.org/abs/2001.09784.  arXiv: 2001.09784.
  * Pleiss et al. [2017] Geoff Pleiss, Manish Raghavan, Felix Wu, Jon Kleinberg, and Kilian Q Weinberger.  On Fairness and Calibration.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 5680–5689. Curran Associates, Inc., 2017.  URL http://papers.nips.cc/paper/7151-on-fairness-and-calibration.pdf.
  * Razzaghi et al. [2015] Talayeh Razzaghi, Oleg Roderick, Ilya Safro, and Nick Marko.  Fast imbalanced classification of healthcare data with missing values.  In _2015 18th International Conference on Information Fusion (Fusion)_ , pages 774–781. IEEE, 2015.
  * Razzaghi et al. [2019] Talayeh Razzaghi, Ilya Safro, Joseph Ewing, Ehsan Sadrfaridpour, and John D Scott.  Predictive models for bariatric surgery risks with imbalanced medical datasets.  _Annals of Operations Research_ , 280(1):1–18, 2019.
  * Sadrfaridpour et al. [2019] Ehsan Sadrfaridpour, Talayeh Razzaghi, and Ilya Safro.  Engineering fast multilevel support vector machines.  _Machine Learning_ , 108(11):1879–1917, 2019\.
  * Wen et al. [2018] Zeyi Wen, Jiashuai Shi, Qinbin Li, Bingsheng He, and Jian Chen.  ThunderSVM: A fast SVM library on GPUs and CPUs.  _The Journal of Machine Learning Research_ , 19(1):797–801, 2018.
  * [25] Muhammad Bilal Zafar, Isabel Valera, Manuel Gomez Rodriguez, and Krishna P Gummadi.  Fairness Constraints: Mechanisms for Fair Classiﬁcation.  page 9.
  * Zafar et al. [2017] Muhammad Bilal Zafar, Isabel Valera, Manuel Gomez Rodriguez, and Krishna P. Gummadi.  Fairness Beyond Disparate Treatment & Disparate Impact: Learning Classification without Disparate Mistreatment.  _Proceedings of the 26th International Conference on World Wide Web_ , pages 1171–1180, April 2017.  doi: 10.1145/3038912.3052660.  URL http://arxiv.org/abs/1610.08452.  arXiv: 1610.08452.
  * Zemel et al. [2013] Rich Zemel, Yu Wu, Kevin Swersky, Toni Pitassi, and Cynthia Dwork.  Learning fair representations.  In _International conference on machine learning_ , pages 325–333. PMLR, 2013.
  * Zliobaite [2015] Indre Zliobaite.  On the relation between accuracy and fairness in binary classification.  _arXiv preprint arXiv:1505.05723_ , 2015.