[a]Debasish Banerjee

# Introducing Fermionic Link Models

Emilie Huffman  Lukas Rammelmüller

###### Abstract

Quantum link models (QLMs) are extensions of Wilson-type lattice gauge
theories, and show rich physics beyond the phenomena of conventional Wilson
gauge theories. Here we explore the physics of $U(1)$ symmetric QLMs, both
using a more conventional quantum spin-1/2 representation, as well as a
fermionic representation. In 2D, we show that both bosonic and fermionic QLMs
display the same physics. We then explore the models in 3D and find different
behavior for the two QLMs. For the bosons, we see evidence for a quantum phase
transition from a symmetry broken phase to a potential quantum spin liquid
phase. For the fermions, we identify not one but two distinct phases in
addition to a symmetry broken phase. We explore the symmetries of the ground
state in the strong coupling limit, which breaks lattice symmetries and
examine the spectrum for both models.

##  1 Introduction

The Hamiltonian formulation of lattice gauge theories [1] (LGTs) has seen
renewed interest in the context of recent theoretical and experimental
developments. On the theoretical side, the introduction and the exploration of
quantum field theories in finite-dimensional Hilbert spaces [2, 3, 4, 5] have
gathered a lot of momentum given the potential of realizing such models on
quantum computers in the near future (see [6] for an overview). On the other
hand, the last decades have seen astonishing developments in precise control
of quantum systems, and engineering a variety of interactions between them
[7]. This has brought to life many toy models in tabletop quantum experiments
[8, 9, 10, 11]. Moreover, such experiments have opened up the possibility of
actually addressing physical problems, such as phases of frustrated magnets,
quantum chromodynamics (QCD) at finite density, or real-time evolution of
quantum many-body systems, where the sign problem presents an obstacle to
simulations but not to experiments.

One challenging aspect in this direction is the formulation of gauge theories
on finite-dimensional Hilbert spaces. Compact gauge groups used to formulate
LGTs non-perturbatively, as was first done by Wilson [12] to describe QCD,
give rise to infinite-dimensional Hilbert spaces even locally. Naive
truncation risks breaking gauge invariance. A theoretically appealing approach
is to regulate this local infinity by gauge-invariant finite-dimensional
operators, called quantum links [13, 14, 15]. For the simpler case of an
$U(1)$ Abelian lattice gauge theory, these operators are quantum spins [16],
while for non-Abelian gauge theories, such as QCD, they can be chosen to be
appropriate finite-dimensional operators [17, 18]. One particular
representation using rishons (fermion bilinears) [17] is particularly suited
for realizing these models in analog quantum simulators [19]. Furthermore, it
is possible to realize the physics of Wilson formulation in the quantum link
approach by using large representations of the link operators [20, 21]. One
should note that different truncation schemes of the Hamiltonian formulation
for gauge theories are also under active investigation [22, 23, 24].

Here we introduce a new class of models, generalizing the quantum link
construction. After reviewing the formulation of the $U(1)$ quantum link model
(QLM) with spin-$\frac{1}{2}$ quantum links, we show it is possible to realize
a different gauge theory by choosing fermionic operators. For the square
lattice, the theories realized with spin and with fermionic operators are
identical to each other. However, for the cubic lattice the two models lead to
differing physics, and give rise to rich phase diagrams. Various implications
of this construction are discussed.

##  2 Plaquette operator in Abelian $U(1)$ LGTs

We first review the Hamiltonian formulation of Wilson LGTs. For an Abelian
$U(1)$ gauge theory, the degree of freedom is a quantum rotor on the bond
joining two lattice sites $x$ and $x+\hat{i}$, where the unit vector $\hat{i}$
is a spatial direction. The rotor is characterized by an angle
$\phi\in[0,2\pi[$, and the gauge field operators are
$U_{x,\hat{i}}=\exp(i\phi_{x,\hat{i}})$ and
$U_{x,\hat{i}}^{\dagger}=\exp(-i\phi_{x,\hat{i}})$. The corresponding electric
field is the angular momentum, $E_{x,\hat{i}}=-i\partial_{\phi_{x,\hat{i}}}$
and satisfies the following commutation relations:

$[E_{x,\hat{i}},U_{y,\hat{j}}]=U_{x,\hat{i}}\delta_{xy}\delta_{ij},~{}~{}[E_{x,\hat{i}},U_{y,\hat{j}}^{\dagger}]=-U_{x,\hat{i}}^{\dagger}\delta_{xy}\delta_{ij};~{}~{}[U_{x,\hat{i}},U_{y,\hat{j}}^{\dagger}]=0.$
(1)

It is convenient to work in the electric flux basis where $E_{x,\hat{i}}$ is
integer-valued $m_{x,\hat{i}}=0,\pm 1,\pm 2,\cdots$. Thus the Hilbert space of
even a single degree of freedom is infinite.

For an extended lattice system, a configuration can be specified by the value
of the electric flux at each link. The gauge field operators
$U_{x,\hat{i}}=L_{x,\hat{i}}^{+}$ and
$~{}U_{x,\hat{i}}^{\dagger}=L_{x,\hat{i}}^{-}$ raise and lower the electric
flux, respectively. The Hamiltonian for this system consists of the electric
and the magnetic energy terms, the latter constructed from plaquette operators
$U_{\Box}$ (with $e$ as the gauge coupling):

$H=\frac{e^{2}}{2}\sum_{x,i}E_{x,\hat{i}}^{2}-\frac{1}{2e^{2}}\sum_{\Box}(U_{\Box}+U^{\dagger}_{\Box});~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}U_{\Box}=U_{x,i}U_{x+\hat{i},\hat{j}}U^{\dagger}_{x+\hat{j},\hat{i}}U^{\dagger}_{x,\hat{j}}.$
(2)

The gauge invariance of the theory is due to the presence of local conserved
charges

$G_{x}=\sum_{i}(E_{x,\hat{i}}-E_{x-\hat{i},\hat{i}}),$ (3)

the lattice Gauss law operators, which commute with the Hamiltonian
$[G_{x},H]=0$. $G_{x}$ are the generators of gauge transformation
$V=\prod_{x}\exp(i\theta_{x}G_{x})$, with $\theta_{x}\in[0,2\pi[$. The
geometry of the plaquette and the Gauss Law operator are shown in Fig. 1
(left) and (middle) respectively.

It is now interesting to consider how to realize the same local $U(1)$
symmetry using a finite-dimensional Hilbert space. One needs operators with
the same commutations as Eq.(1), but nevertheless operating on finite-
dimensional Hilbert spaces. The simplest example is a quantum spin-$S$, having
a a $(2S+1)-$d Hilbert space, represented via $\vec{S}=(S^{+},S^{-},S^{z})$.
From elementary quantum mechanics we know that
$[S^{z},S^{+}]=S^{+},~{}[S^{z},S^{-}]=-S^{-}$. Identifying the electric flux
with the $S^{z}$ component, the gauge fields are the raising and lowering
operators of the electric flux, as before. The Hamiltonian of this QLM is
given by Eq. (2), and the Gauss law by Eq. (3), with the respective spin
operators for the gauge and electric fields. Therefore, using quantum spin-$S$
as a degree of freedom on each link and identifying
$U_{x,\hat{i}}=S^{+}_{x,\hat{i}},U^{\dagger}_{x,\hat{i}}=S^{-}_{x,\hat{i}}$,
and $E_{x,\hat{i}}=S^{z}_{x,\hat{i}}$ realizes a $U(1)$ Abelian LGT on
$(2S+1)-$d local Hilbert space.

Figure 1: (Left) The plaquette operator on a cubic lattice. The dashed line
represents the third direction, perpendicular to the plane of the paper.
(Middle) Six links touch a single site of the 3d spatial lattice. Gauss’ law
thus involves 6 number operators, 3 in the forward and 3 in the backward
links, each of which can be either empty or occupied. There are 20 states
which satisfy the Gauss law. (Right) The plaquette using fermionic operators
can be understood a simultaneous hop of the fermionic particles along the
lines indicated, see text for details.

It turns out that even the simplest non-trivial QLM has very rich physics.
Realized at its most "quantum" limit, the spin-$\frac{1}{2}$ representation of
the $U(1)$ group has a local $2$-d Hilbert space on the links. The gauge
symmetry is identical to that in the Wilson theory, since we only require the
first two expressions of Eq. (1) to prove the commutation $[H,G_{x}]=0$. The
most interesting consequence of the QLM construction is the non-unitarity of
the gauge link operator:
$[U_{x,\hat{i}},U_{y,\hat{j}}^{\dagger}]=2E_{x,\hat{i}}\delta_{xy}\delta_{ij}$.
In fact, this causes the QLMs to give rise to a wider class of physical
phenomena than possible in Wilson-type gauge theories. This theory in
$(2+1)-$d gives rise to crystalline confined phases with fractionalized
electric fluxes [25, 26]. It also has anomalous excited states, called quantum
scars, which exhibit slow thermalization [27].

The rich physics of the QLM as compared to the Wilson theory motivates us to
further our research into the class of models which gives rise to newer phases
of quantum gauge matter [28]. A conceptually new model can be constructed
using the $2$-d local Hilbert space of fermionic operators instead of spins.
The presence or the absence of the fermion exhausts the two possibilities.
Similar to the raising and lowering operators, the gauge fields are
represented by creation or annihilation operators:
$U_{x,\hat{i}}=c^{\dagger}_{x,\hat{i}},~{}U^{\dagger}_{x,\hat{i}}=c_{x,\hat{i}}$.
Moreover, the generator of the gauge transformations is the number operator,
which is also identified with the electric flux:
$E_{x,\hat{i}}=(n_{x,\hat{i}}-\frac{1}{2})$. Of course, this entire
construction only works since these operators also satisfy the commutation
relations amongst each other as the quantum spins:

$\displaystyle[n_{x,\hat{i}},c^{\dagger}_{y,\hat{j}}]$
$\displaystyle=c^{\dagger}_{x,\hat{i}}\delta_{xy}\delta_{ij};~{}[n_{x,\hat{i}},c_{y,\hat{j}}]=c_{x,\hat{i}}\delta_{xy}\delta_{ij};~{}[c^{\dagger}_{x,\hat{i}},c_{y,\hat{j}}]=2E_{x,\hat{i}}=2(n_{x,\hat{i}}-\frac{1}{2});$
(4) $\displaystyle\\{c_{x,\hat{i}},c_{y,\hat{j}}\\}$
$\displaystyle=\\{c^{\dagger}_{x,\hat{i}},c^{\dagger}_{y,\hat{j}}\\}=0;~{}~{}\\{c^{\dagger}_{x,\hat{i}},c_{y,\hat{j}}\\}=\delta_{xy}\delta_{ij}.$
(5)

Because of the anti-commutation relations, naively one can expect that the
physics of this gauge theory will be very different from the one formed using
quantum spins (i.e. bosonic links).

##  3 The plaquette model with fermionic links

In terms of the fermionic operators, the plaquette and Gauss’ law are:

$U_{\Box}=c^{\dagger}_{x,\hat{i}}c^{\dagger}_{x+\hat{i},\hat{j}}c_{x+\hat{j},\hat{i}}c_{x,\hat{j}};~{}~{}G_{x}=\sum_{i}\left(n_{x,\hat{i}}-n_{x-\hat{i},\hat{i}}\right).$
(6)

Even though it might look onerous at first sight, the plaquette operator can
be interpreted as a simultaneous hop of two fermions as shown in the right
panel of Fig. 1. However, the hopping can only occur if the shaded sites are
occupied and their directed neighbors are empty. Thus, among the six hoppings
one can imagine, only two are actually allowed.

As long as we construct the Hamiltonian from gauge-invariant operators, the
resulting model will retain the gauge invariance. Moreover, the electric flux
energy term contributes a constant, and hence we omit it from the Hamiltonian.
However, it is always possible to add a background electric field term which
couples to the total electric flux. This can be considered as a self-adjoint
extension, which we disregard for the moment. Instead, we consider a Rokhsar-
Kivelson-like term, also considered for the usual QLM realized with spins:

$H_{\rm
F}=-J\sum_{\Box}(U_{\Box}+U^{\dagger}_{\Box})+\lambda\sum_{\Box}(U_{\Box}+U^{\dagger}_{\Box})^{2}.$
(7)

The Hamiltonian satisfies the usual lattice point group symmetries, such as
lattice translation, rotation, parity, as well as charge conjugation, a
discrete $\mathbb{Z}(2)$ symmetry. A detailed discussion of the lattice
symmetries has implications on the understanding of its phase diagram and is
discussed in depth in [28]. The only other global internal symmetry of the
model is the $\mathbb{U}(1)^{d}$ generated by large gauge transformation
(winding numbers). In 2 spatial dimensions, these operators are:

${W}_{x}=\frac{1}{{L}_{x}}\sum_{x}({n}_{r=(x,y_{0}),\hat{j}}-\frac{1}{2});~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{W}_{y}=\frac{1}{{L}_{y}}\sum_{y}({n}_{r=(x_{0},y),\hat{i}}-\frac{1}{2}).$
(8)

Every state in the electric flux basis can be characterized by a specific
winding number, and two such examples are shown in Fig. 2. The Hamiltonian
commutes with ${\rm W}_{\rm x},{\rm W}_{\rm y}$, and they can be used to
block-diagonalize the Hamiltonian.

Figure 2:  Examples of configurations with zero (left) and finite (right)
winding number states. The winding number equals half of the difference
between the occupied and empty links touched by the vertical (horizontal) line
in the respective directions. These are topological quantities, and the exact
location of the line is irrelevant as long as the line winds across the
lattice.

##  4 From two to three spatial dimensions

In this section, we argue that the proposed model has identical physics with
spin-$\frac{1}{2}$ model in two spatial dimensions, but not in three. To
understand the difference, we ask if it is possible for the fermions to
interchange their spatial positions. If this is allowed, then the fermion
anti-commutation properties will lead to additional negative sign factors and
potentially change the physics. Otherwise, the physics will be the same, since
the Hamiltonian, the Gauss law and the Hilbert space are identical for both
formulations.

Figure 3:  From left to right, top to bottom: A set of plaquette hops to make
the fermions (shaded in red) exchange their positions with their neighbours
(shaded in blue) on the same diagonal lines (blue and cyan). The nature of the
plaquette interactions moves both the red fermions, and at the end of the
sequence _both_ pairs of fermions exchange positions, even in finite volume.
There are no unpaired negative signs due to Pauli exclusion principle, and the
physics is the same as in the bosonic version.

We will present a geometrical proof of the statement. Note from Fig. 1 (right)
that the plaquette interaction causes the fermions to move diagonally along
the dotted lines, which pass through the mid-point of the neighbouring links.
As indicated in Fig. 3, the motion of the fermions is restricted to 1-d non-
overlapping lines, which wrap across the boundary, but do not intersect. As a
result, it is impossible for two fermions to interchange places, except by
going across the boundary. However, another key feature of the interaction is
that it involves a _simultaneous_ hop of a fermion on an adjacent line.
Therefore, for every instance when two fermions on a line interchange their
positions, another exchange happens between two other fermions in one of the
two adjacent lines. Consequently, there is never a net negative sign
encountered. Thus the physics of the model with spins and fermions in two
dimensions are identical.

In three dimensions, this argument does not work. Drawing the paths on which
the fermions can move, one traces out a set of non-intersecting planes, as
shown in Fig. 4. On the plane, it is easily possible for a pair of fermions to
interchange their positions (shown in the figure as red and yellow) keeping
the position of all other fermions unchanged (shown in grey). We have verified
the existence of paths which have an extra negative sign due to this
interchange. Since such negative signs are absent in the corresponding model
formulated with quantum spin-$\frac{1}{2}$, it is reasonable to expect that
the two models show different physics.

Figure 4:  From left to right: Two fermions colored in red and yellow have
been chosen among the inactive (grey) fermions, whose positions are exchanged
through a series of flips. Finally, the red and yellow have switched places,
keeping all the grey ones at their original positions. This is an example of a
configuration which has an overall negative sign compared to the spin version.

##  5 Physics of the 3d model and Outlook

Motivated by the prospect of investigating novel phases in this class of
models, we have carried out large scale exact diagonalization (ED) studies for
both the bosonic and the fermionic model in $(3+1)-$d [28]. The complete
Hamiltonian we consider is the one in Eq. (7), with the Gauss law as in Eq.
(6). ED calculations are exponentially expensive, but nevertheless by
restricting to the Gauss law sector and by exploiting the global symmetries
extensively, we have been able to diagonalize systems of up to 72 particles
(spins) and gain some insight into the phase diagram using finite size scaling
(FSS).

On very general grounds, we expect that the phase of both the models will be
similar for large and negative $\lambda$. The $\lambda$-term counts whether
the plaquettes are flippable (if the position of the fermions in the plaquette
are in the preferred two of the six total possibilities). For large negative
$\lambda$ therefore, this term dominates over the fermion-hopping. In the
absence of fermionic hopping, the anti-commutations cannot play a role, and
the two models would have similar physics.

Figure 5: Gap of the lowest energy eigenstates realized with spin (left) and
fermion (right) operators.

In Fig. 5 we have plotted the lowest 15 eigenvalues of the Hamiltonian as a
function of $\lambda$ both the bosonic (left) and the fermionic (right)
models. In region of large negative $\lambda$, the spectra are similar, even
quantitatively. A closer analysis of the energy levels indicate a spontaneous
symmetry breaking (SSB) of the rotational symmetry, discussed in detail in
[28]. In between $-1<\lambda<0$, both models show different behaviour of the
first excited state. For the fermionic model there is a level crossing, while
for the bosonic model the sign of the SSB disappears. The FSS is consistent
with the existence of a Coulomb phase in this region [28]. On increasing
$\lambda\rightarrow 1$, both models flow to the Rokhsar-Kivelson point, which
brings in the fluxes from the excited states. The fermionic model seems to
have a narrow region for $\lambda<1$ where the flux is condensed [28].

In addition to the physics issues explored in the [28], there are several very
intriguing conceptual ideas about this particular model which branch out from
this work. An almost obvious question is how to explore the model on larger
lattices in three spatial dimensions. Quantum Monte Carlo is not efficient for
the bosonic model, and for the fermionic model there is even a sign problem.
On the experimental side, our reinterpretation of the plaquette interaction
could be an important way of realizing this on cold atom quantum simulators,
an issue which has remained challenging. We may also ask if dynamical
signatures about the thermalization of this model in two dimensions manifest
in three dimensions.

Another completely different set of questions could result in further
development of the ideas presented here. Since we are working with fermions,
how does one reach the Kogut-Susskind-Wilson limit? If there is a phase with
large correlation length, can we take the continuum limit using dimensional
reduction? More importantly, what would the low-energy effective theory with
fermionic photons look like? Would this idea generalize to non-Abelian gauge
theories? We believe that answering even some of these questions would
increase our understanding of quantum field theories at the frontiers.

## References

  * [1] J.B. Kogut and L. Susskind, _Hamiltonian Formulation of Wilson’s Lattice Gauge Theories_ , _Phys. Rev. D_ 11 (1975) 395.
  * [2] R. Brower, S. Chandrasekharan, S. Riederer and U.J. Wiese, _D theory: Field quantization by dimensional reduction of discrete variables_ , _Nucl. Phys. B_ 693 (2004) 149 [hep-lat/0309182].
  * [3] H. Singh and S. Chandrasekharan, _Qubit regularization of the $O(3)$ sigma model_, _Phys. Rev. D_ 100 (2019) 054505 [1905.13204].
  * [4] NuQS collaboration, _$\sigma$ Models on Quantum Computers_, _Phys. Rev. Lett._ 123 (2019) 090501 [1903.06577].
  * [5] N. Klco and M.J. Savage, _Digitization of scalar fields for quantum computing_ , _Phys. Rev. A_ 99 (2019) 052335 [1808.10378].
  * [6] J. Preskill, _Quantum computing in the nisq era and beyond_ , _Quantum_ 2 (2018) 79.
  * [7] I.M. Georgescu, S. Ashhab and F. Nori, _Quantum simulation_ , _Rev. Mod. Phys._ 86 (2014) 153.
  * [8] E.A. Martinez, C.A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl et al., _Real-time dynamics of lattice gauge theories with a few-qubit quantum computer_ , _Nature_ 534 (2016) 516.
  * [9] N. Klco, E.F. Dumitrescu, A.J. McCaskey, T.D. Morris, R.C. Pooser, M. Sanz et al., _Quantum-classical computation of schwinger model dynamics using quantum computers_ , _Physical Review A_ 98 (2018) 032331.
  * [10] A. Mil, T.V. Zache, A. Hegde, A. Xia, R.P. Bhatt, M.K. Oberthaler et al., _A scalable realization of local u (1) gauge invariance in cold atomic mixtures_ , _Science_ 367 (2020) 1128.
  * [11] B. Yang, H. Sun, R. Ott, H.-Y. Wang, T.V. Zache, J.C. Halimeh et al., _Observation of gauge invariance in a 71-site quantum simulator_ , _arXiv preprint arXiv:2003.08945_ (2020) .
  * [12] K.G. Wilson, _Confinement of Quarks_ , _Phys. Rev. D_ 10 (1974) 2445.
  * [13] D. Horn, _Finite Matrix Models With Continuous Local Gauge Invariance_ , _Phys. Lett. B_ 100 (1981) 149.
  * [14] P. Orland and D. Rohrlich, _Lattice Gauge Magnets: Local Isospin From Spin_ , _Nucl. Phys. B_ 338 (1990) 647.
  * [15] S. Chandrasekharan and U.-J. Wiese, _Quantum link models: A discrete approach to gauge theories_ , _Nuclear Physics B_ 492 (1997) 455\.
  * [16] U.-J. Wiese, _From Quantum Link Models to D-Theory: A Resource Efficient Framework for the Quantum Simulation and Computation of Gauge Theories_ , 2107.09335.
  * [17] R. Brower, S. Chandrasekharan and U.J. Wiese, _QCD as a quantum link model_ , _Phys. Rev. D_ 60 (1999) 094502 [hep-th/9704106].
  * [18] D. Banerjee, F.J. Jiang, T.Z. Olesen, P. Orland and U.J. Wiese, _From the $SU(2)$ quantum link model on the honeycomb lattice to the quantum dimer model on the kagomé lattice: Phase transition and fractionalized flux strings_, _Phys. Rev. B_ 97 (2018) 205108 [1712.08300].
  * [19] D. Banerjee, M. Bögli, M. Dalmonte, E. Rico, P. Stebler, U.J. Wiese et al., _Atomic Quantum Simulation of U(N) and SU(N) Non-Abelian Lattice Gauge Theories_ , _Phys. Rev. Lett._ 110 (2013) 125303 [1211.2242].
  * [20] B. Schlittgen and U.J. Wiese, _Low-energy effective theories of quantum spin and quantum link models_ , _Phys. Rev. D_ 63 (2001) 085007 [hep-lat/0012014].
  * [21] T.V. Zache, M. Van Damme, J.C. Halimeh, P. Hauke and D. Banerjee, _Achieving the continuum limit of quantum link lattice gauge theories on quantum devices_ , 2104.00025.
  * [22] E. Zohar and J.I. Cirac, _Eliminating fermionic matter fields in lattice gauge theories_ , _Phys. Rev. B_ 98 (2018) 075119 [1805.05347].
  * [23] I. Raychowdhury and J.R. Stryker, _Loop, string, and hadron dynamics in SU(2) Hamiltonian lattice gauge theories_ , _Phys. Rev. D_ 101 (2020) 114502 [1912.06133].
  * [24] Z. Davoudi, I. Raychowdhury and A. Shaw, _Search for efficient formulations for Hamiltonian simulation of non-Abelian lattice gauge theories_ , _Phys. Rev. D_ 104 (2021) 074505 [2009.11802].
  * [25] D. Banerjee, P. Widmer, F.J. Jiang and U.J. Wiese, _Crystalline Confinement_ , _PoS_ LATTICE2013 (2014) 333 [1311.2459].
  * [26] D. Banerjee, S. Caspar, F.J. Jiang, J.H. Peng and U.J. Wiese, _Nematic confined phases in the $u(1)$ quantum link model on a triangular lattice: An opportunity for near-term quantum computations of string dynamics on a chip_, 2021\.
  * [27] D. Banerjee and A. Sen, _Quantum Scars from Zero Modes in an Abelian Lattice Gauge Theory on Ladders_ , _Phys. Rev. Lett._ 126 (2021) 220601 [2012.08540].
  * [28] D. Banerjee, E. Huffman and L. Rammelmüller, _Exploring bosonic and fermionic link models on $(3+1)-$d tubes_, 2021.


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

RecSys in HR’22: The 2nd Workshop on Recommender Systems for Human Resources,
in conjunction with the 16th ACM Conference on Recommender Systems, September
18–23, 2022, Seattle, USA.

<EMAIL_ADDRESS>]

[orcid=0000-0002-0465-535X<EMAIL_ADDRESS>]

[ orcid= 0000-0002-0458-9233<EMAIL_ADDRESS>]

[orcid=0000-0003-3946-928X<EMAIL_ADDRESS>]

# Closing the Gender Wage Gap: Adversarial Fairness in Job Recommendation

Clara Rus  Radboud University, Houtlaan 4, 6525XZ, Nijmegen, the Netherlands
DPG Media Online Services, Jacob Bontiusplaats 9, 1018LL, Amsterdam, the
Netherlands  Jeffrey Luppes  Harrie Oosterhuis  Gido H. Schoenmacker

(2022)

###### Abstract

The goal of this work is to help mitigate the already existing gender wage gap
by supplying unbiased job recommendations based on resumes from job seekers.
We employ a generative adversarial network to remove gender bias from word2vec
representations of 12M job vacancy texts and 900k resumes. Our results show
that representations created from recruitment texts contain algorithmic bias
and that this bias results in real-world consequences for recommendation
systems. Without controlling for bias, women are recommended jobs with
significantly lower salary in our data. With adversarially fair
representations, this wage gap disappears, meaning that our debiased job
recommendations reduce wage discrimination. We conclude that adversarial
debiasing of word representations can increase real-world fairness of systems
and thus may be part of the solution for creating fairness-aware
recommendation systems.

###### keywords:

Generative adversarial networks Fairness-aware machine learning Recruitment
Gender bias

##  1 Introduction

The recruitment industry relies more and more on automation for processing,
searching, and matching job vacancies to job seekers. However, automation of
the recruitment process can lead to discriminatory results with respect to
certain groups, based on gender, ethnicity or age [1]. Inequality in
employment and remuneration still exists between for example ethnic groups [2,
3, 4] and gender groups [5, 6], thus naive implementations of AI recruitment
systems are at risk of copying and perpetuating these inequalities.

One reason for an algorithm to show discriminatory behaviour is the input data
[7]. If the data is under–representative or if historical bias is present,
then the system can propagate this in its predictions [1]. Ignoring the
presence of bias in the data, can perpetuate existing (gender) stereotypes and
inequalities in employment.

Examples of systems that have shown biased behaviour with respect to gender
include the Amazon recruitment
system111https://www.reuters.com/article/idUSKCN1MK08G and the Facebook Add
algorithm [8]. Also widely used models, such as BERT [9] and word2vec [10],
have been shown to create biased representations [11, 12]. Obtaining fair
representations could eliminate the bias present in the data and help a system
achieve fairer predictions [13].

One way to learn debiased representations is through adversarial learning.
State-of-the-art adversarial debiasing methods [14, 15, 16, 17, 18] rely on
the same general approach as generative adversarial networks [19]. A generator
model is trained to produce new data representations, that are critiqued by an
adversary neural network. The adversary tries to predict the sensitive
variable (in our case, gender) from the produced representation. By training
the representations together with an adversary and classifier, they are aimed
to be both fair and useful for the task.

This work is motivated by the desire to supply unbiased job recommendations to
job seekers. We focus specifically on mitigating gender bias in word
embeddings obtained from recruitment texts using adversarial learning. Our
work adds to existing research by applying state-of-the-art debiasing [14, 20]
to industrial sized free-format recruitment textual data. Firstly, we
investigate gender bias in the existing representations and the unfairness it
results in. Secondly, we apply two debiasing methods to create new
representations. These methods balance multi-label classification to ensure
that task-relevant information has been preserved, with an adversarial setup
that attempts to remove the effects of gender bias. The resulting new
representations are tested in a job recommendation setting where the
difference in wage between jobs recommended based on female/male resumes is
evaluated.

To summarize, our contributions are three-fold: (i) we measure whether
adversarial learning can mitigate gender bias in representations of industrial
sized free-format recruitment textual data; (ii) we show whether debiased
representations help achieve fairness and performance on a multi-label
classification task; and (iii) to the authors’ best knowledge, we are the
first to successfully apply debiased representations to help solve the gender
wage-gap in a job recommendation setting. Moreover, our implementation of the
adversarial debiasing method is publicly available.

In the next section, our data and methods are described in detail. After that,
the results are presented. Lastly, these results are discussed together with
our final conclusions and suggestions for future directions.

##  2 Data and Method

###  2.1 Data

The recruitment data set used throughout this research consists of job
vacancies and job seeker information provided by DPG Recruitment. Job vacancy
information included (i) salary ranges, (ii) working hours, and (iii)
anonymised free-format job vacancy texts. In total there are 12 millions
vacancies.

Job seeker information consisted of (i) one or more industry group(s) that the
job seeker expressed interest in (out of a total of 21 pre-defined groups),
(ii) inferred dichotomous gender, and (iii) anonymised free-format resume
texts. Gender of the job seeker was inferred based on first name. From the
total of available resumes, entries with missing data (65%) or ambiguous first
name (3%) were excluded, leaving 904,576 (32%) complete resumes with a female
to male ratio of 0.93. Anonymisation included removal of all names (including
company names), dates, addresses, telephone numbers, email addresses,
websites, and other contact information. A more complete overview of this data
is given in Appendix A.

Both vacancy and resume texts were embedding into 300-dimensional word vector
using a word2vec [10] model trained on all vacancy texts. Finally, each text
was represented as the mean over the embeddings of the words composing the
text.

###  2.2 Bias and debiasing

Previous research has shown that popular models such as BERT [9] and word2vec
[10] can create biased representations [11, 12, 21]. In this work, two
debiasing methods were employed to combat this bias.

Firstly, to create a simple baseline, we attempt to debias the representations
by replacing gendered words with neutral words. For example, gendered pronouns
“she”/“he”, “her”/“his” are replaced with neutral pronouns “they” and
“theirs”. Gendered words such as: “woman”/“man”, “girl”/“boy” are replaced
with the word “person”. The full list of substitutions can be found in
Appendix B. A new word2vec model was trained on this augmented corpus,
resulting in new representations for both the resumes and the vacancies. In
the remaining text, “original representations” will refer to the
representations trained on the original texts, whereas “word-substitution
representations” will refer to the representations trained on the altered
texts.

Secondly, we applied the adversarial approach as proposed by Edwards and
Storkey [14]. This method consists of three neural network components: a
generator, a classifier, and an adversary. Inspired by Özdenizci et al. [20],
we chose the following the architecture: The generator is a multilayer
perceptron with three hidden layers of $128$ neurons that outputs a $300$
dimensional vector representing the new representations. The classifier and
the adversary have one hidden layer of $128$ neurons. The output dimension of
the classifier is $21$ (industry group classes), and the output dimension of
the adversary is one (gender). An architecture schematic is included as Figure
1.

The generator creates new representations for the classification task, while
the adversary attempts to predict a sensitive variable gender from these new
representations. The goal of the generator is to create representations that
can fool the adversary in such a way that the sensitive variable can no longer
be predicted, while also obtaining a good performance on the classification
task. The classification task is considered to be a multi-label task of 21
classes, predicting the industry group(s) for each job seeker. This means that
the classification loss should be minimized while the adversarial loss should
be maximized. The final loss (Equation 1) of the model is a weighted sum of
the adversarial loss and the classification loss, where Z are the newly
generated representations, Y’ are the predictions of the classifier and S’ are
the predictions of the adversary:

$L=\alpha L_{cla}(Z,Y^{\prime})+\beta L_{adv}(Z,S^{\prime}).$ (1)

We will call representations created by this method “adversarial
representations”. Because the adversarial process could be unstable, all
results pertaining to these are the mean of $5$ independent complete training
runs.

Figure 1: Architecture of the adversarial setup. The left section (green)
represents the generator, consisting of an input layer (d=300) for the
word2vec representations, three hidden layers (d=128), and an output layer
(d=300) for the debiased representations. The top section (blue) represents
the classifier consisting of a hidden layer (d=128) and an output layer
$\hat{Y}$ (d=21) encoding the industry groups. The bottom section (red)
represents the adversary and consists of a hidden layer (d=128) and an output
neuron $\hat{S}$ (d=1) encoding the sensitive variable gender.

###  2.3 Evaluation

Classifiers for both industry groups and sensitive variable are evaluated in
terms of accuracy and area under the receiver operating characteristic curve
(AUC). Fairness was evaluated using statistical parity [22]:

$P(cla(Z)=1|S=1)-P(cla(Z)=1|S=0)<\epsilon.$ (2)

In the recruitment industry, if a system designed to match resumes and
vacancies perpetuates biased associations, it could lead to a wage gap between
salaries of women and men [23]. To specifically test differences in salary, a
salary association test was performed between the representations of the
resumes and of the vacancies. Using the embeddings of the resumes and the
vacancies the L2 distance matrix was computed and each resume was matched to
the closest vacancy. The salary distribution of the matched vacancies of the
female-inferred group were compared with the male-inferred group.

###  2.4 Experimental Setup

The train split was created by taking 30% of random samples for the validation
split, and the rest of the full data is used for training. The full data set
was not used for the salary association due to computational limitations.
Instead, 10,000 resumes were associated with all jobs from the time period
June 2020–June 2021 that provided salary information. This resulted in 23,501
total vacancies. All experiments were conducted using a fixed 70-30% split and
the Adam optimizer with a learning rate of $1\mathrm{e}{-5}$. For all
components the binary cross-entropy loss was used. Parameters of the final
loss (Equation 1) are set in the following way: $\alpha$ = 1, $\beta$ = 1. The
implementation of the adversarial debiasing method can be found at:
https://github.com/ClaraRus/Debias-Embeddings-Recruitment.

##  3 Results

###  3.1 Prediction of sensitive variable

Firstly, the discriminatory power to predict the sensitive variable gender was
tested using the original, word-substitution, and adversarial representations.
Using the original representations, gender is predicted with 94% AUC and an
accuracy of 86%; the word-substitution representations result in an 93% AUC
and an accuracy of 85%; lastly, the adversarial representations lowered both
the accuracy and the AUC to 82%.

###  3.2 Prediction of industry group

Secondly, the information contents and statistical parity of the three
representation types were tested by attempting to predict the function group
based on resume representation. Table 1 shows the result obtained in terms of
performance and statistical parity.

Training a classifier with the original representations of the resumes
obtained a statistical parity of $0.076$. The word-substitution
representations obtained similar results. Using an adversarial approach
improved the statistical parity by $21$%, at the cost of lowering the accuracy
by $2$ percentage point and the true positive rate by $16$ percentage point.

Table 1: Statistical parity and performance of multi-label classification of
21 industry groups from three different types of representations. Original
representations were obtained using word2vec. Word-substitution
representations were obtained using a word-substitution debiasing method.
Adversarial representations were obtained using the adversarial debiasing
method. The “Overall” row represents the weighted mean. Parity: statistical
parity (Eq. 2), closer to zero is better. TPR: True positive rate.

| Original | Word-substitution | Adversarial
---|---|---|---
| Parity | Accuracy | TPR | Parity | Accuracy | TPR | Parity | Accuracy | TPR
Overall | 0.076 | 0.90 | 0.37 | 0.760 | 0.90 | 0.38 | 0.060 | 0.89 | 0.21
Administration/Secretarial | 0.267 | 0.74 | 0.52 | 0.271 | 0.74 | 0.53 | 0.260 | 0.72 | 0.45
Automation/Internet | 0.066 | 0.83 | 0.46 | 0.069 | 0.83 | 0.47 | 0.045 | 0.82 | 0.32
Policy/Executive | 0.000 | 0.79 | 0.19 | 0.001 | 0.79 | 0.21 | 0.005 | 0.79 | 0.08
Security/Defence/Police | 0.010 | 0.83 | 0.16 | 0.009 | 0.83 | 0.16 | 0.000 | 0.83 | 0.03
Commercial/Sales | 0.074 | 0.74 | 0.36 | 0.070 | 0.74 | 0.35 | 0.059 | 0.72 | 0.19
Consultancy/Advice | 0.026 | 0.75 | 0.20 | 0.033 | 0.75 | 0.22 | 0.012 | 0.74 | 0.07
Design/Creative/Journalism | 0.004 | 0.82 | 0.26 | 0.005 | 0.82 | 0.30 | 0.001 | 0.82 | 0.09
Management | 0.070 | 0.78 | 0.32 | 0.063 | 0.78 | 0.29 | 0.052 | 0.77 | 0.22
Financial/Accounting | 0.021 | 0.81 | 0.47 | 0.021 | 0.81 | 0.48 | 0.020 | 0.80 | 0.29
Financial services | 0.012 | 0.79 | 0.22 | 0.012 | 0.79 | 0.28 | 0.007 | 0.79 | 0.10
HR/Training | 0.041 | 0.80 | 0.32 | 0.043 | 0.80 | 0.34 | 0.014 | 0.78 | 0.09
Catering/Retail | 0.037 | 0.76 | 0.33 | 0.023 | 0.76 | 0.27 | 0.018 | 0.75 | 0.14
Procurement/Logistics/Transport | 0.115 | 0.77 | 0.38 | 0.102 | 0.77 | 0.35 | 0.087 | 0.76 | 0.24
Legal | 0.015 | 0.85 | 0.45 | 0.015 | 0.85 | 0.44 | 0.002 | 0.84 | 0.09
Customer service/Call centre/Front office | 0.039 | 0.76 | 0.12 | 0.031 | 0.76 | 0.10 | 0.001 | 0.76 | 0.01
Marketing/PR/Communications | 0.031 | 0.77 | 0.41 | 0.031 | 0.77 | 0.45 | 0.028 | 0.76 | 0.27
Medical/Healthcare | 0.115 | 0.76 | 0.40 | 0.116 | 0.77 | 0.40 | 0.100 | 0.75 | 0.27
Education/Research/Science | 0.045 | 0.77 | 0.32 | 0.057 | 0.77 | 0.39 | 0.031 | 0.75 | 0.16
Other | 0.005 | 0.68 | 0.04 | 0.009 | 0.68 | 0.05 | 0.000 | 0.68 | 0.00
Production/Operational | 0.063 | 0.78 | 0.27 | 0.062 | 0.78 | 0.27 | 0.043 | 0.77 | 0.15
Technology | 0.165 | 0.79 | 0.51 | 0.169 | 0.79 | 0.52 | 0.153 | 0.78 | 0.43

Table 2: Salary Association Test between resumes and vacancies. For each
resume the most similar vacancy was assigned based on Euclidean distance in
the representation space. The values represent the salary per hour in Euros
(€). Original representations were obtained using word2vec. Word-substitution
representations were obtained using the word-substitution debiasing method.
Adversarial representations were obtained using the adversarial debiasing
method. The top three rows represent the weighted summary statistics. The
industries names with an asterisk (*) are the ones for which the adversarial
method reduced the wage gap.

| Original | Word-substitution | Adversarial
---|---|---|---
| Female | Male | Wage gap | Female | Male | Wage gap | Female | Male | Wage gap
Mean | 25.28 | 26.09 | 0.81 | 25.19 | 26.14 | 0.95 | 27.06 | 27.15 | 0.09
Standard deviation | 9.43 | 9.90 | 0.47 | 9.54 | 10.07 | 0.53 | 10.14 | 9.94 | -0.20
Median | 23.40 | 23.62 | 0.22 | 22.95 | 23.62 | 0.67 | 23.97 | 24.30 | 0.33
Administration/Secretarial* | 23.50 | 24.94 | 1.44 | 23.44 | 24.95 | 1.51 | 26.45 | 26.41 | -0.04
Automation/Internet | 28.34 | 28.58 | 0.24 | 28.05 | 29.02 | 0.97 | 29.94 | 28.44 | -1.50
Policy/Executive* | 29.90 | 31.23 | 1.33 | 30.16 | 31.53 | 1.37 | 30.35 | 31.18 | 0.83
Security/Defence/Police* | 24.81 | 22.78 | -2.03 | 24.81 | 23.09 | -1.72 | 25.51 | 26.44 | 0.93
Commercial/Sales* | 23.76 | 25.77 | 2.01 | 23.88 | 25.53 | 1.65 | 26.66 | 27.39 | 0.73
Consultancy/Advice* | 29.27 | 30.49 | 1.22 | 29.25 | 30.42 | 1.17 | 29.92 | 30.63 | 0.71
Design/Creative/Journalism | 26.39 | 26.33 | -0.06 | 26.13 | 26.12 | -0.01 | 28.22 | 28.02 | -0.20
Management* | 29.49 | 31.22 | 1.73 | 29.49 | 31.47 | 1.98 | 30.66 | 30.31 | -0.35
Financial/Accounting* | 24.30 | 27.94 | 3.64 | 24.43 | 28.07 | 3.64 | 27.20 | 28.62 | 1.42
Financial services* | 24.33 | 27.85 | 3.52 | 24.19 | 27.76 | 3.57 | 26.80 | 28.67 | 1.87
HR/Training | 28.59 | 28.87 | 0.28 | 28.80 | 29.10 | 0.30 | 29.52 | 29.15 | -0.37
Catering/Retail* | 22.80 | 23.76 | 0.96 | 22.76 | 23.49 | 0.73 | 25.15 | 24.50 | 0.65
Procurement/Logistics/Transport | 23.70 | 23.28 | -0.42 | 23.46 | 23.30 | -0.16 | 25.96 | 25.10 | -0.86
Legal* | 24.89 | 28.79 | 3.90 | 25.52 | 28.91 | 3.39 | 28.82 | 29.01 | 0.19
Customer service/Call centre/Front office* | 22.89 | 23.78 | 0.89 | 23.01 | 23.85 | 0.84 | 25.35 | 25.96 | 0.61
Marketing/PR/Communications* | 26.64 | 27.65 | 1.01 | 26.71 | 27.55 | 0.84 | 28.86 | 29.22 | 0.36
Medical/Healthcare* | 26.30 | 27.51 | 1.21 | 26.11 | 27.19 | 1.08 | 27.17 | 28.07 | 0.90
Education/Research/Science | 28.82 | 27.43 | -1.39 | 28.30 | 27.65 | -0.65 | 27.66 | 29.07 | 1.41
Other* | 24.91 | 24.58 | -0.33 | 24.79 | 24.84 | 0.05 | 26.07 | 26.32 | 0.25
Production/Operational* | 21.69 | 22.61 | 0.92 | 21.15 | 22.49 | 1.34 | 24.39 | 23.94 | -0.45
Technology* | 25.51 | 24.09 | -1.42 | 24.57 | 24.07 | -0.50 | 25.79 | 25.79 | 0.00

###  3.3 Salary Association Test

Thirdly, a salary association test was performed using the three
representation types. Table 2 describes the salary distribution of the female
and male groups for each debiasing method. In the female group there are 4827
samples and in the male group 5173.

Using the original representations, female-inferred resumes were associated
with a mean salary of €25.28 per hour, whereas male-inferred resumes were
associated with a mean salary of €26.09 per hour, which is significantly
(p<$1\mathrm{e}{-5}$) higher. This results in an estimated average annual wage
gap of €1680.

Using the word-substitution representations, female-inferred resumes were
associated with a mean salary of €25.19 per hour, whereas male-inferred
resumes were associated with a mean salary of €26.14 per hour. The difference
between the means of the female group and the male group increased, broadening
the annual wage gap to €1900 (with a significant difference between groups,
p<$1\mathrm{e}{-7}$).

With the adversarial representations, female-inferred resumes were associated
with a mean salary of €27.06 an hour, whereas male-inferred resumes were
associated with a mean salary of €27.15 an hour. Using the adversarial method
to generate fair representations for both the resumes and vacancies decreased
the mean gap, lowering the annual wage gap to €180. This resulted in the
female/male difference now being non-significant (p=$0.47$).

Table 2 shows the mean salary per hour for each industry group. Ideally
females and males belonging to the same industry group should have similar
salaries. The word-substitution representations lowered the wage gap in 13 of
the industry groups by €620 per year, while increasing the gap for 7 with an
average of €460 per year. For “Financial/Accounting” there is no change in the
salary association. The adversarial method lowered the wage gap in 16 out of
the 21 industries by an average of €2160 per year but it increased the gap in
the rest of the industries by an average of €780 per year.

##  4 Discussion and Conclusion

This work focused on removing gender bias from word embeddings of vacancy
texts and resumes with the goal of creating debiased job recommendations. It
showed that gender can be predicted extremely well from anonymised resume
embeddings and that naive resume-to-job recommendations based on these
embeddings can perpetuate the “wage gap” that exists between women and men.
Adversarial debiasing improved statistical parity for industry classification
based on resume and eliminated the female/male salary difference in job
recommendations. This suggests that adversarial debiasing can help make fairer
recommendations in realistic scenarios.

Our results indicate that anonymisation alone is not enough to remove indirect
information about the gender of the job seeker. Namely, from our $900k$
anonymised resumes, gender could be predicted with an AUC of $0.94$. This
exceeds similar results that have been shown in a smaller data set
(AUC=$0.81$) [24]. This is a common problem in fairness-aware machine
learning, where removal of directly sensitive information is undermined by
correlated features that allow the sensitive information to be inferred [22].

The difficultly of removing gender bias from language was further illustrated
by our data augmentation attempt to substitute a selection of gendered words
by neutral words before word2vec training. The resulting embeddings did not
effect much change in any of our tests. Previous work on word substitution
data augmentation has been shown effective [25, 26], so it may be that our
results are limited by the quantity and/or selection of our word substitution
pairs (Table 4), which were taken from [21]. While it is possible to improve
upon our substitution pairs, creating a complete list of gendered words as
used in vacancies and resumes is challenging if not unfeasible, especially in
multiple languages.

In contrast, the adversarial approach improved both statistical parity and the
wage gap in our data. Using the adversarial representations, prediction of
gender dropped from an AUC of $0.94$ to $0.82$ while performance of industry
group prediction, in terms of accuracy, dropped only minimally (Table 1).
However, the true positive rate was decreased, indicating that performance was
affected. These results are linked and can be adapted by changing the $\alpha$
and $\beta$ parameters in Equation 1: more gender-neutral embeddings will
likely lead to improved statistical parity but decreased industry prediction
performance. Since statistical parity balances for equal true positive rate,
the false positive and negative rates are likely to be affected.

Our analysis reveals that ignoring the presence of bias in recruitment texts,
that are used to match resumes and vacancies, could lead to severe unwanted
discriminatory behaviour.

The original representations produced a wage gap of €1680 per year between the
female group and the male group.

The adversarial representations eliminated this wage gap to a statistically
insignificant difference.

This result is especially important, because it shows that the adversarial
representations did not just perform better on selected in-vitro metrics, but
also improved fairness in a real application. This suggests that the
adversarial representations do not remove bias only “cosmetically” [27], but
instead are effective for improving fairness in job recommendation. The
adversarial method increased the mean salary for both the female group and the
male group, with a higher increase for the female group to balance the gap.
This is a positive outcome as the method did not sacrifice the salaries of one
of the groups in order to reduce the wage gap.

This work was limited by several factors. Firstly, while the fairness of job
recommendations was assessed, due to unavailability of data, the quality of
recommendations could not be evaluated. This was mitigated by performing a
related classification task: predicting which industry groups a job seeker is
interested in based on resume. The accuracy of $0.89$ on this task suggests
that salient information relevant to job placement has been preserved.
However, since the true positive rate was impacted, it seems likely that the
recall of the recommendations would be impacted too. Secondly, the recommender
system to suggest jobs based on representation distances was relatively
simple; if job-to-resume association data were available, a more complex
solution might be preferable. Thirdly, because gender was inferred for this
research, it was not possible to include non-binary gender identities [28].
Since this group is vulnerable to employment discrimination [29, 30, 31], it
should not be overlooked and more research here is needed. Fourthly, results
reported in this research use only word2vec document embeddings; other types
of embeddings are not considered. Lastly, training of the models was performed
using a fixed split instead of cross-validation, which was infeasible due to
time and costs. However, the results are likely to be representative given the
large size of the data set.

The strengths of our work include the application of adversarial debiasing for
fairness-aware machine learning on real and large industry data. While
adversarial debiasing for fairness is not novel [14, 17, 18, 32, 33],
applications generally extend to publicly available benchmark data sets that
make it difficult to assess its applicability to real-world recommendation
systems. Our work is one of the first to show the results of adversarial
fairness in a real, industrial-scale system. In addition, this research
obtained an acceptable trade-off between fairness and performance for a
complex multi-label classification task. Finally, this work showed that the
adversarial approach eliminated the female/male wage gap in our job
recommendations, even though it was not trained for this task.

In conclusion, this work identified gender bias in word representations and
salary associations based on recruitment industry texts and successfully
applied adversarial debiasing to combat gender bias in job recommendation.
With adversarial representations, the mean female/male wage gap was no longer
statistically significant due to being reduced by 89% from €1680 to €180
annually. Our results show that adversarial debiasing of word representations
can increase real-world fairness of recommendation systems and thus may
contribute to creating fairness-aware machine learning systems.

## References

  * Köchling and Wehner [2020] A. Köchling, M. C. Wehner,  Discriminated by an algorithm: A systematic review of discrimination and fairness by algorithmic decision-making in the context of HR recruitment and HR development,  Business Research 13 (2020) 795–848. doi:10.1007/s40685-020-00134-w.
  * Thijssen et al. [2021] L. Thijssen, B. Lancee, S. Veit, R. Yemane,  Discrimination against Turkish minorities in Germany and the Netherlands: field experimental evidence on the effect of diagnostic information on labour market outcomes,  Journal of Ethnic and Migration Studies 47 (2021) 1222–1239. doi:10.1080/1369183X.2019.1622793.
  * Bisschop et al. [2020] P. Bisschop, B. ter Weel, J. Zwetsloot,  Ethnic employment gaps of graduates in the Netherlands,  De Economist 168 (2020) 577–598. doi:10.1007/s10645-020-09375-w.
  * Ramos et al. [2021] M. Ramos, L. Thijssen, M. Coenders,  Labour market discrimination against Moroccan minorities in the Netherlands and Spain: A cross-national and cross-regional comparison,  Journal of Ethnic and Migration Studies 47 (2021) 1261–1284. doi:10.1080/1369183X.2019.1622824.
  * Matteazzi et al. [2018] E. Matteazzi, A. Pailhé, A. Solaz,  Part-time employment, the gender wage gap and the role of wage-setting institutions: Evidence from 11 European countries,  European Journal of Industrial Relations 24 (2018) 221–241. doi:10.1177/0959680117738857.
  * Ciminelli et al. [2021] G. Ciminelli, C. Schwellnus, B. Stadler,  Sticky floors or glass ceilings? The role of human capital, working time flexibility and discrimination in the gender wage gap,  OECD Economics Department Working Papers (2021). doi:10.1787/02ef3235-en.
  * Chouldechova and Roth [2020] A. Chouldechova, A. Roth,  A snapshot of the frontiers of fairness in machine learning,  Communications of the ACM 63 (2020) 82–89. doi:10.1145/3376898.
  * Ali et al. [2019] M. Ali, P. Sapiezynski, M. Bogen, A. Korolova, A. Mislove, A. Rieke,  Discrimination through optimization: How Facebook’s ad delivery can lead to biased outcomes,  Proceedings of the ACM on human-computer interaction 3 (2019) 1–30. doi:10.1145/3359301.
  * Devlin et al. [2018] J. Devlin, M.-W. Chang, K. Lee, K. Toutanova,  BERT: Pre-training of deep bidirectional transformers for language understanding,  arXiv preprint (2018). doi:10.48550/arXiv.1810.04805.
  * Mikolov et al. [2013] T. Mikolov, K. Chen, G. Corrado, J. Dean,  Efficient estimation of word representations in vector space,  Proceedings of Workshop at ICLR 2013 (2013). doi:10.48550/arXiv.1301.3781.
  * Kurita et al. [2019] K. Kurita, N. Vyas, A. Pareek, A. W. Black, Y. Tsvetkov,  Measuring bias in contextualized word representations,  arXiv preprint (2019). doi:10.48550/arXiv.1906.07337.
  * Bolukbasi et al. [2016] T. Bolukbasi, K.-W. Chang, J. Y. Zou, V. Saligrama, A. T. Kalai,  Man is to computer programmer as woman is to homemaker? Debiasing word embeddings,  Advances in neural information processing systems 29 (2016). doi:10.48550/arXiv.1607.06520.
  * Beutel et al. [2017] A. Beutel, J. Chen, Z. Zhao, E. H. Chi,  Data decisions and theoretical implications when adversarially learning fair representations,  arXiv preprint (2017). doi:10.48550/arXiv.1707.00075.
  * Edwards and Storkey [2016] H. Edwards, A. J. Storkey,  Censoring representations with an adversary,  in: Y. Bengio, Y. LeCun (Eds.), 4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016, pp. 1–14. doi:10.48550/arXiv.1511.05897.
  * Xu et al. [2019] D. Xu, S. Yuan, L. Zhang, X. Wu,  FairGAN+: Achieving fair data generation and classification through generative adversarial nets,  in: 2019 IEEE International Conference on Big Data (Big Data), 2019, pp. 1401–1406. doi:10.1109/BigData47090.2019.9006322.
  * Madras et al. [2018] D. Madras, E. Creager, T. Pitassi, R. Zemel,  Learning adversarially fair and transferable representations,  in: J. Dy, A. Krause (Eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, PMLR, 2018, pp. 3384–3393. URL: https://proceedings.mlr.press/v80/madras18a.html. doi:10.48550/arXiv.1802.06309.
  * Wu et al. [2021] C. Wu, F. Wu, X. Wang, Y. Huang, X. Xie,  Fairness–aware news recommendation with decomposed adversarial learning,  in: Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, 2021, pp. 4462–4469. doi:10.48550/arXiv.2006.16742.
  * Liu et al. [2022] H. Liu, N. Zhao, X. Zhang, H. Lin, L. Yang, B. Xu, Y. Lin, W. Fan,  Dual constraints and adversarial learning for fair recommenders,  Knowledge-Based Systems 239 (2022) 108058. doi:10.1016/j.knosys.2021.108058.
  * Goodfellow et al. [2014] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio,  Generative adversarial networks,  arXiv preprint (2014). doi:10.48550/arXiv.1406.2661.
  * Özdenizci et al. [2020] O. Özdenizci, Y. Wang, T. Koike-Akino, D. Erdoğmuş,  Learning invariant representations from EEG via adversarial inference,  IEEE access 8 (2020) 27074–27085. doi:10.1109/ACCESS.2020.2971600.
  * Caliskan et al. [2017] A. Caliskan, J. Bryson, A. Narayanan,  Semantics derived automatically from language corpora contain human-like biases,  Science 356 (2017) 183–186. doi:10.1126/science.aal4230.
  * Mehrabi et al. [2021] N. Mehrabi, F. Morstatter, N. Saxena, K. Lerman, A. Galstyan,  A survey on bias and fairness in machine learning,  ACM Comput. Surv. 54 (2021). doi:10.1145/3457607.
  * Calanca et al. [2019] F. Calanca, L. Sayfullina, L. Minkus, C. Wagner, E. Malmi,  Responsible team players wanted: an analysis of soft skill requirements in job advertisements,  EPJ Data Science 8 (2019) 1–20. doi:10.1140/epjds/s13688-019-0190-z.
  * Parasurama and Sedoc [2022] P. Parasurama, J. Sedoc,  Gendered information in resumes and its role in algorithmic and human hiring bias,  in: Academy of Management Proceedings, volume 2022, Academy of Management Briarcliff Manor, NY 10510, 2022, p. 17133. doi:https://doi.org/10.5465/AMBPP.2022.285.
  * Hall Maudslay et al. [2019] R. Hall Maudslay, H. Gonen, R. Cotterell, S. Teufel,  It’s all in the name: Mitigating gender bias with name–based counterfactual data substitution,  in: Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP), Association for Computational Linguistics, Hong Kong, China, 2019, pp. 5267–5275. doi:10.18653/v1/D19-1530.
  * Pruksachatkun et al. [2021] Y. Pruksachatkun, S. Krishna, J. Dhamala, R. Gupta, K.-W. Chang,  Does robustness improve fairness? Approaching fairness with word substitution robustness methods for text classification,  in: Findings of the Association for Computational Linguistics: ACL-IJCNLP 2021, 2021, pp. 3320–3331. doi:10.48550/arXiv.2106.10826.
  * Gonen and Goldberg [2019] H. Gonen, Y. Goldberg,  Lipstick on a pig: Debiasing methods cover up systematic gender biases in word embeddings but do not remove them,  in: Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), Association for Computational Linguistics, Minneapolis, Minnesota, 2019, pp. 609–614. doi:10.18653/v1/N19-1061.
  * Richards et al. [2016] C. Richards, W. P. Bouman, L. Seal, M. J. Barker, T. O. Nieder, G. T’Sjoen,  Non-binary or genderqueer genders,  International Review of Psychiatry 28 (2016) 95–102. doi:10.3109/09540261.2015.1106446.
  * Harrison et al. [2012] J. Harrison, J. Grant, J. L. Herman,  A gender not listed here: Genderqueers, gender rebels, and otherwise in the national transgender discrimination survey,  LGBTQ Public Policy Journal at the Harvard Kennedy School 2 (2012) 13. URL: https://escholarship.org/uc/item/2zj46213.
  * Davidson [2016] S. Davidson,  Gender inequality: Nonbinary transgender people in the workplace,  Cogent Social Sciences 2 (2016) 1236511. doi:10.1080/23311886.2016.1236511.
  * Fogarty and Zheng [2018] A. A. Fogarty, L. Zheng, Gender ambiguity in the workplace: Transgender and gender-diverse discrimination, ABC-CLIO, 2018\. URL: https://publisher.abc-clio.com/9781440863233/.
  * Zhang et al. [2018] B. H. Zhang, B. Lemoine, M. Mitchell,  Mitigating unwanted biases with adversarial learning,  in: Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society, 2018, pp. 335–340. doi:10.1145/3278721.3278779.
  * Sattigeri et al. [2019] P. Sattigeri, S. C. Hoffman, V. Chenthamarakshan, K. R. Varshney,  Fairness GAN: Generating datasets with fairness properties using a generative adversarial network,  IBM Journal of Research and Development 63 (2019) 3–1. doi:10.1147/JRD.2019.2945519.

##  Appendix A Industry and resumes distributions

Table 3 shows the distribution of the samples in each industry group over the
whole data set. In more than half of the industries the resumes from the
female group are under-represented. The “Technology” industry has the fewest
samples from the female group. The are industries where the resumes from the
female group are over-represented, such as the “Administration/Secretarial”
and the “Customer service/Call centre/Front office” industries, where the
resumes from the male group are more than three times less present.

Table 3: Distribution of samples over each industry group. Counts and percentages per industry group do not sum to the expected totals, because job seekers were free to select multiple groups. “F-M Ratio” represent the ratio between the number of females within an industry group and the number of males. | Male | Female | Total | F-M Ratio
---|---|---|---|---
Overall | 467173 | 437403 | 904576 | 0.93
Administration/Secretarial | 45293 (9%) | 167585 (38%) | 212878 (23%) | 3.7
Automation/Internet | 49527 (10%) | 8547 (1%) | 58074 (6%) | 0.17
Policy/Executive | 40086 (8%) | 33541 (7%) | 73627 (8%) | 0.83
Security/Defence/Police | 23134 (4%) | 8821 (2%) | 31955 (3%) | 0.38
Commercial/Sales | 92801 (19%) | 66461 (15%) | 159262 (17%) | 0.71
Consultancy/Advice | 69914 (14%) | 42245 (9%) | 112159 (12%) | 0.6
Design/Creative/Journalism | 19279 (4%) | 24839 (5%) | 44118 (4%) | 1.28
Management | 67412 (14%) | 32153 (7%) | 99565 (11%) | 0.48
Financial/Accounting | 34233 (7%) | 25523 (5%) | 59756 (6%) | 0.74
Financial services | 34342 (7%) | 29882 (6%) | 64224 (7%) | 0.87
Catering/Retail | 44647 (9%) | 60588 (13%) | 105235 (11%) | 1.35
HR/Training | 26852 (5%) | 53679 (12%) | 80531 (8%) | 1.99
Procurement/Logistics/Transport | 99429 (21%) | 29677 (6%) | 129106 (14%) | 0.29
Legal | 8638 (1%) | 18488 (4%) | 27126 (2%) | 2.14
Customer service/Call centre/Front office | 20000 (4%) | 71090 (16%) | 91090 (10%) | 3.55
Marketing/PR/Communications | 46832 (10%) | 58598 (13%) | 105430 (11%) | 1.25
Medical/Healthcare | 24018 (5%) | 85414 (19%) | 109432 (12%) | 1.25
Education/Research/Science | 38430 (8%) | 66318 (15%) | 104748 (11%) | 1.72
Other | 86749 (18%) | 82728 (18%) | 169477 (18%) | 0.95
Production/Operational | 77790 (5%) | 25452 (25%) | 103242 (11%) | 0.32
Technology | 102798 (22%) | 9097 (2%) | 111895 (12%) | 0.08

##  Appendix B Word-substitution debiasing method

Table 4 shows the substitutions of the gendered words with the neutral words
for both English and Dutch.

Table 4: Substitutions of gendered words with neutral words used in the word-substitution debiasing method in both English (top) and Dutch (bottom). Male Word | Female Word | Neutral Word
---|---|---
he | she | they
his | hers | theirs
himself | herself | themselves
male | female | person
boy | girl | person
man | woman | person
hij | zij/ze | u
zijn | haar | uw
hijzelf | zijzelf | uzelf
jongen | meisje | persoon
man | vrouw | persoon


* [2] Peter Aczel.  An introduction to inductive definitions.  In J. Barwise, editor, Handbook of Mathematical Logic, pages 739–782. North-Holland, 1977.
  * [3] Ross Anderson and Roger Needham.  Programming Satan’s computer.  In Jan van Leeuwen, editor, Computer Science Today: Recent Trends and Developments, LNCS 1000, pages 426–440. Springer, 1995.
  * [4] Giampaolo Bella and Lawrence C. Paulson.  Using Isabelle to prove properties of the Kerberos authentication system.  In Orman and Meadows [30].
  * [5] Dominique Bolignano.  An approach to the formal verification of cryptographic protocols.  In Third ACM Conference on Computer and Communications Security, pages 106–118. ACM Press, 1996.
  * [6] Dominique Bolignano and Valérie Ménissier-Morain.  Formal verification of cryptographic protocols using Coq.  Technical report, INRIA-Rocquencourt, 1996.
  * [7] Stephen H. Brackin.  A HOL extension of GNY for automatically analyzing cryptographic protocols.  In Computer Security Foundations Workshop [13], pages 62–75.
  * [8] Jeremy Bryans and Steve Schneider.  CSP, PVS and a recursive authentication protocol.  In Orman and Meadows [30].
  * [9] John A. Bull and David J. Otway.  The authentication protocol.  Technical Report DRA/CIS3/PROJ/CORBA/SC/1/CSM/436-04/0.5b, Defence Research Agency, Malvern, UK, 1997.  In press??
  * [10] M. Burrows, M. Abadi, and R. M. Needham.  A logic of authentication.  Proceedings of the Royal Society of London, 426:233–271, 1989.
  * [11] John Clark and Jeremy Jacob.  On the security of recent protocols.  Information Processing Letters, 56(3):151–155, 1995.
  * [12] 8th Computer Security Foundations Workshop. IEEE Computer Society Press, 1995\.
  * [13] 9th Computer Security Foundations Workshop. IEEE Computer Society Press, 1996\.
  * [14] 10th Computer Security Foundations Workshop. IEEE Computer Society Press, 1997\.
  * [15] Tim Dierks and Christopher Allen.  The TLS protocol: Version 1.0, January 1999.  Request for Comments: 2246, on the Internet at ftp://ftp.isi.edu/in-notes/rfc2246.txt.
  * [16] Bruno Dutertre and Steve Schneider.  Using a PVS embedding of CSP to verify authentication protocols.  In Elsa L. Gunter and Amy Felty, editors, Theorem Proving in Higher Order Logics: TPHOLs ’97, LNCS 1275, pages 121–136. Springer, 1997\.
  * [17] Dieter Gollmann.  What do we mean by entity authentication?  In Symposium on Security and Privacy, pages 46–54. IEEE Computer Society, 1996.
  * [18] M. J. C. Gordon and T. F. Melham.  Introduction to HOL: A Theorem Proving Environment for Higher Order Logic.  Cambridge University Press, 1993.
  * [19] James W. Gray, III and John McLean.  Using temporal logic to specify and verify cryptographic protocols.  In Computer Security Foundations Workshop [12], pages 108–116.
  * [20] Matthew Hennessy.  The Semantics of Programming Languages: An Elementary Introduction Using Structural Operational Semantics.  Wiley, 1990.
  * [21] C. A. R. Hoare.  Communicating Sequential Processes.  Prentice-Hall, 1985.
  * [22] Richard Kemmerer, Catherine Meadows, and Jonathan Millen.  Three systems for cryptographic protocol analysis.  Journal of Cryptology, 7(2):79–130, 1994.
  * [23] Richard A. Kemmerer.  Analyzing encryption protocols using formal verification techniques.  IEEE Journal on Selected Areas in Communications, 7(4):448–457, May 1989.
  * [24] Gavin Lowe.  Breaking and fixing the Needham-Schroeder public-key protocol using CSP and FDR.  In T. Margaria and B. Steffen, editors, Tools and Algorithms for the Construction and Analysis of Systems: second international workshop, TACAS ’96, LNCS 1055, pages 147–166. Springer, 1996.
  * [25] Gavin Lowe.  Some new attacks upon security protocols.  In Computer Security Foundations Workshop [13], pages 162–169.
  * [26] Gavin Lowe.  SPLICE/AS: A case study in using CSP to detect errors in security protocols.  Technical report, Oxford University Computing Laboratory, 1996.
  * [27] Gavin Lowe.  Casper: A compiler for the analysis of security protocols.  Journal of Computer Security, 6:53–84, 1998.
  * [28] Wenbo Mao and Colin Boyd.  Towards formal analysis of security protocols.  In Computer Security Foundations Workshop VI, pages 147–158. IEEE Computer Society Press, 1993.
  * [29] Catherine A. Meadows.  Analyzing the Needham-Schroeder public-key protocol: A comparison of two approaches.  In E. Bertino, H. Kurth, G. Martella, and E. Montolivo, editors, Computer Security — ESORICS 96, LNCS 1146, pages 351–364. Springer, 1996\.
  * [30] Hilarie Orman and Catherine Meadows, editors.  Workshop on Design and Formal Verification of Security Protocols. DIMACS, September 1997.
  * [31] Sam Owre, John Rushby, Natarajan Shankar, and Friedrich von Henke.  Formal verification for fault-tolerant architectures: Prolegomena to the design of PVS.  IEEE Transactions on Software Engineering, 21(2):107–125, February 1995.
  * [32] Lawrence C. Paulson.  Isabelle: A Generic Theorem Prover.  Springer, 1994.  LNCS 828.
  * [33] Lawrence C. Paulson.  ML for the Working Programmer.  Cambridge University Press, 2nd edition, 1996.
  * [34] Lawrence C. Paulson.  Generic automatic proof tools.  In Robert Veroff, editor, Automated Reasoning and its Applications: Essays in Honor of Larry Wos, chapter 3. MIT Press, 1997.
  * [35] Lawrence C. Paulson.  Inductive analysis of the Internet protocol TLS.  ACM Transactions on Information and System Security, 2(3):332–351, August 1999.
  * [36] Lawrence C. Paulson.  Relations between secrets: Two formal analyses of the Yahalom protocol.  Journal of Computer Security, 9(3):197–216, 2001.
  * [37] Lawrence C. Paulson.  Tool support for logics of programs.  In Manfred Broy, editor, Mathematical Methods in Program Development: Summer School Marktoberdorf 1996, NATO ASI Series F, pages 461–498. Springer, Published 1997.
  * [38] R. L. Rivest, A. Shamir, and L. Adleman.  A method for obtaining digital signatures and public-key cryptosystems.  Communications of the ACM, 21(2):120–126, February 1978.
  * [39] A. W. Roscoe.  Modelling and verifying key-exchange protocols using CSP and FDR.  In Computer Security Foundations Workshop [12], pages 98–107.
  * [40] Peter Y. A. Ryan.  The design and verification of security protocols.  Technical Report DRA/CIS3/SISG/CR/96/1.0, Defence Research Agency, May 1996.
  * [41] Peter Y. A. Ryan and Steve A. Schneider.  An attack on a recursive authentication protocol: A cautionary tale.  Information Processing Letters, 65(1):7–10, January 1998.
  * [42] Steve Schneider.  Verifying authentication protocols with CSP.  In Computer Security Foundations Workshop [14], pages 3–17.


# An Improvement on the Hotelling $T^{2}$ Test Using the Ledoit-Wolf Nonlinear
Shrinkage Estimator

Benjamin Robinson  Sensors Directorate
US Air Force Research Laboratory
Dayton, Ohio
ORCID: 0000-0002-9391-4543  Robert Malinas  Van Latimer  EECS
University of Michigan
Ann Arbor, MI, USA
ORCID: 0000-0002-1318-5564  Mathematics
UCLA
Los Angeles, California, USA
ORCID: 0000-0001-9789-7739  Beth Bjorkman Morrison  Sensors Directorate
US Air Force Research Laboratory
Dayton, Ohio
ORCID: 0000-0001-7406-0311  Alfred O. Hero, III  EECS
University of Michigan
Ann Arbor, MI, USA
ORCID: 0000-0002-2531-9670

###### Abstract

Hotelling’s $T^{2}$ test is a classical approach for discriminating the means
of two multivariate normal samples that share a population covariance matrix.
Hotelling’s test is not ideal for high-dimensional samples because the
eigenvalues of the estimated sample covariance matrix are inconsistent
estimators for their population counterparts. We replace the sample covariance
matrix with the nonlinear shrinkage estimator of Ledoit and Wolf 2020. We
observe empirically for sub-Gaussian data that the resulting algorithm
dominates past methods (Bai and Saranadasa 1996, Chen and Qin 2010, and Li et
al. 2020) for a family of population covariance matrices that includes
matrices with high or low condition number and many or few nontrivial—i.e.,
spiked—eigenvalues.

###### Index Terms:

Two-sample testing, high-dimensional limit, shrinkage covariance estimation,
Ledoit-Wolf estimator, Hotelling $T^{2}$ test

##  I Introduction

A fundamental problem in statistics and signal processing is determining
whether two independent samples have the same mean. For multivariate samples
that are Gaussian and have a shared population covariance matrix, classical
methods like the Hotelling $T^{2}$ test apply [1]. However, there is no
standard technique if the sample dimension is substantial compared to the
sample sizes: the so-called large-dimensional regime.

Hotelling’s $T^{2}$ test relies on estimating the shared population covariance
matrix using the sample covariance matrix. In the large-dimensional regime,
the sample covariance matrix’s eigenvalues are inconsistent estimators for
their population counterparts [2], leading to poor performance of Hotelling’s
test. Additionally, the sample covariance matrix can be considerably more ill-
conditioned than the population covariance matrix, resulting in numerical
instabilities [3]. As a result, several authors have proposed alternatives to
the standard $T^{2}$ test. Bai and Saranadasa [4] and Chen and Qin [5] propose
tests (BS96 and CQ10) that are an improvement for a well-conditioned
population covariance matrix, and Li et al. [6] propose a test (LAPPW20) that
is an improvement if the population covariance follows the spiked model of
Johnstone [3]. To the best of our knowledge, until now no test improves upon
of all these methods under more general assumptions.

In this paper, we propose a replacement for Hotelling’s $T^{2}$ test and
present some simulations in which it dominates BS96, CQ10, and LAPPW20. We do
not restrict the condition number or the number of spiked eigenvalues of the
population covariance matrix in our testing. Our method, similar to LAPPW20,
is to replace the sample covariance matrix in Hotelling’s $T^{2}$ with a
matrix that has a smaller condition number. More precisely, we use the
covariance matrix estimator of Ledoit and Wolf [7], which belongs to Stein’s
shrinkage class [8, 9]. By contrast, LAPPW20 uses a diagonal-loading
estimator, which also belongs to Stein’s class, but is simply a sum of the
sample covariance matrix and a scalar multiple of the identity. We argue for
our test’s optimality within Stein’s class using ideas similar to those some
of us have applied to other detection problems [10]. We note that we are not
the first to have thought of applying Ledoit-Wolf-type estimators to two-
sample testing, but are developing the idea toward maturity [11].

In Section II, we explore past improvements on Hotelling’s $T^{2}$ test. In
Section III, we give the Ledoit-Wolf estimator utilized in our approach,
present the proposed test, and argue for the test’s asymptotic optimality and
predict a precise asymptotically constant false-alarm rate. In Section IV, we
provide simulations that show empirical improvement over past methods.
Finally, in Section V, we present our conclusions.

##  II Background

Suppose we have independent random $p$-dimensional column vectors
$\mathbf{x}_{ij}\sim\mathcal{N}(\boldsymbol{\mu}_{i},\mathbf{R})$ for
$j=1,2,\dots,n_{i}$ and $i\in\\{1,2\\}$, where $\boldsymbol{\mu}_{i}$ are
unknown means and $\mathbf{R}$ is an unknown $p\times p$ symmetric positive-
definite population covariance matrix. We want to test the following
hypotheses

$\begin{cases}\mathcal{H}_{0}:&\boldsymbol{\mu}_{1}=\boldsymbol{\mu}_{2},\\\
\mathcal{H}_{1}:&\boldsymbol{\mu}_{1}\neq\boldsymbol{\mu}_{2}.\end{cases}$

If the covariance matrix $\mathbf{R}$ is known, a reasonable detector is the
Mahalanobis-distance detector

$(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})^{\prime}\mathbf{R}^{-1}(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\tau,$
(1)

where $(\,\cdot\,)^{\prime}$ denotes the transpose, $\tau\in(0,\infty)$, and
$\overline{\mathbf{x}}_{i}$ is the sample mean of the $\mathbf{x}_{ij}$’s. For
unknown $\mathbf{R}$, the classical replacement of (1) is Hotelling’s $T^{2}$
test:

$(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})^{\prime}\mathbf{S}_{n}^{-1}(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\tau,$
(2)

where $n=n_{1}+n_{2}-2$ and $\mathbf{S}_{n}$ is the “pooled” sample covariance
matrix, given by

$\mathbf{S}_{n}=\frac{1}{n}\sum_{i=1}^{2}\sum_{j=1}^{n_{i}}(\mathbf{x}_{ij}-\overline{\mathbf{x}}_{i})(\mathbf{x}_{ij}-\overline{\mathbf{x}}_{i})^{\prime}.$
(3)

The performance of the above test is well-characterized when
$n_{1},n_{2}\to\infty$ and $p$ is fixed [1]; however, this test can become
inadmissible in the regime where $p,n_{1},n_{2}\to\infty$ and
$p/n_{1}\rightarrow\gamma_{1}$ and $p/n_{2}\rightarrow\gamma_{2}$ for some
$\gamma_{1},\gamma_{2}\in(0,\infty)$ [4]. This latter limit is known as the
_high-dimensional asymptotic regime_.

Bai and Saranadasa [4] suggest an improvement, BS96, under the assumption that
$\left\|\mathbf{R}\right\|=o\left(\sqrt{\mathrm{tr}\left(\mathbf{R}^{2}\right)}\right)$.
This assumption yields that the matrix $\mathbf{R}$ is well-conditioned, for
otherwise the two referenced quantities would be similar in size. Their
proposed test statistic is

$\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}}\left\|\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2}\right\|^{2}-\mathrm{tr}\mathbf{S}_{n}}{\sqrt{\frac{2(n+1)}{n}}B_{n}},$
(4)

where

$B_{n}=\frac{n^{2}}{(n+2)(n-1)}\left(\mathrm{tr}\left(\mathbf{S}_{n}^{2}\right)-\frac{1}{n}\left(\mathrm{tr}\mathbf{S}_{n}\right)^{2}\right).$

The corresponding detector is desirable in that, asymptotically, it has known
constant false-alarm rate and predictable detection power.

In response to the need to analyze ultra-high dimensional data, Chen and Qin
[5] propose a test, CQ10, that is similar to BS96 but is additionally well
suited even to values of $n$ much smaller than $p$. Their test possesses
desirable provable optimality properties, provided again an assumption that
implies a low condition number:
$\mathrm{tr}(\mathbf{R}^{4})=o\\{\mathrm{tr}^{2}(\mathbf{R}^{2})\\}$. Their
test statistic is

$\frac{\sum_{i\neq
j}^{n_{1}}\mathbf{x}_{1i}^{\prime}\mathbf{x}_{1j}}{n_{1}(n_{1}-1)}+\frac{\sum_{i\neq
j}^{n_{2}}\mathbf{x}_{2i}^{\prime}\mathbf{x}_{2j}}{n_{2}(n_{2}-1)}-\frac{2\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}\mathbf{x}_{1i}^{\prime}\mathbf{x}_{2j}}{n_{1}n_{2}}.$

Li et al. [6] assume instead that $\mathbf{R}$ follows the well-known _spiked
covariance model_ of Johnstone [3], i.e., that all the eigenvalues of
$\mathbf{R}$ are equal to one except for a fixed finite number that are larger
than one. This assumption potentially contrasts with the well-conditionedness
assumption of BS96 and CQ10. Li et al. propose a test, LAPPW20, that replaces
$\mathbf{S}_{n}$ in (2) by the _diagonal-loading_ estimator
$\mathbf{S}_{n}+\lambda\mathbf{I}$, for some optimal _loading factor_
$\lambda>0$. Like BS96 and CQ10, LAPPW20 has a known asymptotically constant
false-alarm rate and predictable detection power. Further, Li et al. remark
that LAPPW20 reduces to BS96 for population covariance matrices that are both
well-conditioned and spiked.

The covariance estimator of LAPPW20 belongs to the shrinkage class of Stein
[8, 9]. These are estimators that differ from $\mathbf{S}_{n}$ only in their
eigenvalues. That is, all of the eigenvectors of the estimator are
eigenvectors of $\mathbf{S}_{n}$, but the eigenvalues may be different from
$\mathbf{S}_{n}$’s. Shrinkage covariance estimators have been studied in the
high-dimensional regime in the spiked model by Donoho et al. [12], and in a
more general model by Ledoit and Wolf [13, 14, 15]. In particular, Ledoit and
Wolf have devised closed-form expressions for shrinkage eigenvalues that are
asymptotically optimal with respect to many criteria, including Stein’s loss,
inverse Stein’s loss, Frobenius loss, inverse Frobenius loss, and so-called
minimum-variance loss [7]. In what follows, we will propose again replacing
$\mathbf{S}_{n}$ in Hotelling’s $T^{2}$ test, but rather than a diagonal-
loading estimator, we utilize a shrinkage estimator of Ledoit and Wolf.

##  III Proposed Test

In this section, we define our proposed test, give an indication of its
asymptotic false-alarm rate, and argue for a type of asymptotic detection-
theoretic optimality.

###  III-A Definition

In order to define our proposed test, we must first describe Ledoit and Wolf’s
nonlinear shrinkage estimator from [7]. Throughout this paper, we will
consider the high-dimensional limit in which $n,p\to\infty$ and
$p/n=p/(n_{1}+n_{2}-2)\to\gamma\in(0,\infty)$. The notation $\to$ will always
refer to the limit in which $n,p\to\infty$ and $p/n\to\gamma$. For notational
convenience, we will follow the convention of [16] that quantities that are
not explicitly deemed constant are varying with $n$ and $p$. Identifying
$\mathbf{S}_{n}$ with $\mathbf{S}$, then, we let
$\mathbf{S}=\mathbf{U}\mathrm{diag}(\boldsymbol{\lambda})\mathbf{U}^{\prime}$
be an eigen-decomposition, where for emphasis $\mathbf{U}$ and
$\boldsymbol{\lambda}$ depend on $n$ and $p$. We list the eigenvalues
$\boldsymbol{\lambda}=(\lambda_{1},\dots,\lambda_{p})$ in non-increasing
order, with corresponding eigenvectors
$\mathbf{U}=[\mathbf{u}_{1},\mathbf{u}_{2},\dots,\mathbf{u}_{p}]$.

The Ledoit-Wolf nonlinear shrinkage estimator can be described as follows. Let
$h_{j}=n^{-1/3}\lambda_{j}$ and $[y]^{+}=\max\\{y,0\\}$. Define

$\displaystyle a(\lambda,\boldsymbol{\lambda})$
$\displaystyle:=\sum_{j=[p-n]^{+}+1}^{p}\left\\{-\frac{3(\lambda-\lambda_{j})}{10\pi
h_{j}^{2}}+\right.\left.\frac{3}{4\sqrt{5}\pi h_{j}}\times\right.$
$\displaystyle\qquad\left.\left[1-\frac{1}{5}\left(\frac{\lambda-\lambda_{j}}{h_{j}}\right)^{2}\right]\log\left|\frac{\sqrt{5}h_{j}-\lambda+\lambda_{j}}{\sqrt{5}h_{j}+\lambda-\lambda_{j}}\right|\right\\},$

and

$b(\lambda,\boldsymbol{\lambda}):=\sum_{j=[p-n]^{+}+1}^{p}\frac{3}{4\sqrt{5}h_{j}}\left[1-\frac{1}{5}\left(\frac{\lambda-\lambda_{j}}{h_{j}}\right)^{2}\right]^{+},$

and

$s(\lambda,\boldsymbol{\lambda})=\pi(a(\lambda,\boldsymbol{\lambda})+ib(\lambda,\boldsymbol{\lambda}))/\min\\{n,p\\}.$

Next, the shrunken eigenvalues $\hat{d}_{i}$ are defined as

$\hat{d}_{i}:=\begin{cases}\frac{\lambda_{i}}{\left|1-p/n-(p/n)\lambda_{i}s(\lambda_{i},\boldsymbol{\lambda})\right|^{2}},&\lambda_{i}>0\\\
\frac{1}{(p/n-1)a(0,\boldsymbol{\lambda})/n},&\lambda_{i}=0.\end{cases}$

(We assume $p\neq n$.) Finally, the Ledoit-Wolf estimator, denoted
$\hat{\mathbf{R}}_{\mathrm{LW}}$, is given by
$\hat{\mathbf{R}}_{\mathrm{LW}}=\mathbf{U}\mathrm{diag}(\hat{d}_{1},\hat{d}_{2},\dots,\hat{d}_{p})\mathbf{U}^{\prime}$.

We define our test to be the Hotelling $T^{2}$ test with
$\hat{\mathbf{R}}_{\mathrm{LW}}$ substituted for $\mathbf{S}$:

$T_{\mathrm{LW}}^{2}:=\frac{n_{1}n_{2}}{n_{1}+n_{2}}(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})^{\prime}\hat{\mathbf{R}}_{\mathrm{LW}}^{-1}(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\tau.$

For the false-alarm rate analysis that follows, we shift and re-scale
$T_{\mathrm{LW}}^{2}$ for the equivalent test

$Z=\frac{1}{\sqrt{2p}}(T_{\mathrm{LW}}^{2}-p)\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\tau.$
(5)

###  III-B False-Alarm Rate

A key consideration about any test statistic is whether its null distribution
can be characterized. In the following, we argue intuitively that the null
distribution of $Z$ should be asymptotically standard normal for Gaussian data
and provide empirical support for this assertion. We expect a similar result
holds for sub-Gaussian data.

Let $\sigma_{i}^{2}=\mathbf{u}_{i}^{\prime}\mathbf{R}\mathbf{u}_{i}$. The test
in subsection III-A is related to the test statistic $\tilde{Z}$, defined as
follows. Let $\hat{\mathbf{R}}$ be a shrinkage estimator
$\mathbf{U}\mathbf{D}\mathbf{U}^{\prime}$, where without loss of generality,

$\mathbf{D}=\mathrm{diag}\left(\frac{\sigma_{1}^{2}}{c_{1}},\frac{\sigma_{2}^{2}}{c_{2}},\dots,\frac{\sigma_{p}^{2}}{c_{p}}\right),$

for some positive coefficients $\\{c_{i}$}. Define

$\tilde{Z}:=\frac{1}{\sqrt{2\sum_{i}c_{i}^{2}}}\left(\tilde{T}^{2}-\sum_{i=1}^{p}c_{i}\right),$
(6)

where $\tilde{T}$ is given by

$\tilde{T}^{2}:=\frac{n_{1}n_{2}}{n_{1}+n_{2}}(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2})^{\prime}\hat{\mathbf{R}}^{-1}(\overline{\mathbf{x}}_{1}-\overline{\mathbf{x}}_{2}).$
(7)

By the Berry-Esseen theorem, the conditional distribution of $\tilde{Z}$ given
$\mathbf{S}$ is almost surely asymptotically standard normal.

The shrunken eigenvalues $\hat{d}_{i}$ of Ledoit and Wolf are designed to
approximate $\sigma_{i}^{2}$. A weak form of this approximation that holds
under a general matrix model laid out [13, 7] is that for all fixed intervals
$I\subset[0,\infty)$

$\epsilon:=\frac{1}{p}\sum_{\lambda_{j}\in
I}(\hat{d}_{j}-\sigma_{j}^{2})\overset{a.s.}{\to}0.$ (8)

To borrow notation from [17], we conjecture further that, under the conditions
of [13, 7], $\epsilon=O_{\prec}(p^{-d})$ uniformly in intervals $I$ for some
$d>0$, where, $a_{p}=O_{\prec}(b_{p})$ means that for any $\delta,D>0$, we
have for sufficiently large $p$ that

$\Pr\left[|a_{p}|>p^{\delta}|b_{p}|\right]<p^{-D}.$ (9)

In other words, with polynomially high confidence, $a_{p}$ does not grow much
more quickly than $b_{p}$, if at all. We further conjecture that

$\max_{i,j}\left|\mathbf{u}_{i}^{\prime}\mathbf{R}\mathbf{u}_{j}-\delta_{ij}\hat{d}_{j}\right|=O_{\prec}(p^{-d})$
(10)

for some $d>0$. This conjecture is motivated by an analogous result involving
quadratic forms of Wigner-matrix eigenvectors [18, Equation 1].

Taking $c_{i}=\sigma_{i}^{2}/\hat{d}_{i}$, the test statistic in (6) becomes

$\frac{1}{\sqrt{2\sum_{i}\sigma_{i}^{4}/\hat{d}_{i}^{2}}}\left(T_{\mathrm{LW}}^{2}-\sum_{i=1}^{p}\frac{\sigma_{i}^{2}}{\hat{d}_{i}}\right).$

Using the fact that the above is asymptotically almost surely standard normal
given $\mathbf{S}$, and given the approximations of the last paragraph, we
expect $Z$ to be asymptotically almost surely standard normal as well. In
particular, we expect $Z\gtrless\tau$ has an asymptotically constant false-
alarm rate.

We present a simulation comparing a finite-sample distribution of $Z$ to a
standard normal in Figure 1. Using 1000 Monte-Carlo values of $Z$ with
$\mathbf{R}$ being the 200$\times$200 matrix $\mathbf{R}_{4}$ defined in the
Section IV and $n_{1}=n_{2}=200$, one can see that the finite-sample
distribution of $Z$ approximates a standard normal, as predicted.

Figure 1: A plot showing $Z$ from (5) to be roughly standard normal. Here, the
population covariance matrix follows the model $\mathbf{R}_{4}$ in Section IV.

###  III-C Maximizing Detection Power

We now provide an intuitive argument that $\hat{\mathbf{R}}_{\mathrm{LW}}$ in
(7) optimizes conditional detection power given $\mathbf{S}$ among shrinkage
estimators $\hat{\mathbf{R}}$, for Gaussian data. It can be shown using (8)
that $t=\mathrm{tr}(\mathbf{R}\hat{\mathbf{R}}^{-1})/p$ depends in the limit
only on the limiting sample spectrum, so $\hat{\mathbf{R}}\leftarrow
t\hat{\mathbf{R}}$ can be considered to be a shrinkage estimator, and we
assume without loss of generality that $t$ converges to 1.

Consider the conditional false-alarm rate of the proposed detector given
$\mathbf{S}$:
$\Pr\left[\mathbf{d}^{\prime}\hat{\mathbf{R}}^{-1}\mathbf{d}>\tau\mid\mathbf{S}\right]$,
where $\mathbf{d}=\mathbf{R}^{1/2}\mathbf{w}$ and
$\mathbf{w}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. Using the result [19,
Theorem 5.1.4] regarding concentration of Lipschitz functions on the sphere
and the fact that $p\to\infty$, this probability can be approximated by

$\Pr\left[\left\|\mathbf{w}\right\|^{2}\mathrm{tr}(\mathbf{R}^{1/2}\hat{\mathbf{R}}^{-1}\mathbf{R}^{1/2})>\tau\mid\mathbf{S}\right].$
(11)

Since $\hat{\mathbf{R}}$ depends only on $\mathbf{S}$, the probability above
is a function of $t$, $\tau$, and $p$ alone. Since by assumption $t$ converges
to 1, (11) is asymptotically dependent on $\tau$ and $p$ alone, and the same
can be said of the corresponding unconditional probability. Thus, we wish to
maximize detection power for each $\tau$ and large $p$.

Consider the conditional probability of detection given $\mathbf{S}$:

$\Pr\left[(\boldsymbol{\mu}+\mathbf{d})^{\prime}\hat{\mathbf{R}}^{-1}(\boldsymbol{\mu}+\mathbf{d})>\tau\mid\mathbf{S}\right],$
(12)

where
$\boldsymbol{\mu}=\sqrt{n_{1}n_{2}/(n_{1}+n_{2})}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{2})$.
Expanding (12), we get

$\Pr\left[\alpha+\boldsymbol{\beta}^{\prime}\mathbf{v}+\left\|\mathbf{w}\right\|^{2}>\tau\mid\mathbf{S}\right],$
(13)

where $\alpha=\boldsymbol{\mu}^{\prime}\hat{\mathbf{R}}^{-1}\boldsymbol{\mu}$,
$\boldsymbol{\beta}=2\boldsymbol{\mu}^{\prime}\hat{\mathbf{R}}^{-1}\mathbf{R}^{1/2}\left\|\mathbf{w}\right\|,$
and $\mathbf{v}=\mathbf{w}/\left\|\mathbf{w}\right\|$. If we further condition
on $\left\|\mathbf{w}\right\|$, the argument of the probability in (13)
depends only on the randomness of $\mathbf{v}$, and the probability itself
depends only on $\alpha$ and $\left\|\boldsymbol{\beta}\right\|$. More
precisely, the conditional probability of detection given $\mathbf{S}$ and
$\left\|\mathbf{w}\right\|$ is maximized when
$\alpha/\left\|\boldsymbol{\beta}\right\|$ is maximized.

For any $\left\|\mathbf{w}\right\|$ and $\mathbf{S}$, maximizing
$\alpha/\left\|\boldsymbol{\beta}\right\|$ is equivalent to maximizing the
_signal-to-noise ratio_ :

$\frac{(\boldsymbol{\mu}^{\prime}\hat{\mathbf{R}}^{-1}\boldsymbol{\mu})^{2}}{\boldsymbol{\mu}^{\prime}\hat{\mathbf{R}}^{-1}\mathbf{R}\hat{\mathbf{R}}^{-1}\boldsymbol{\mu}}.$
(14)

Using [19, Theorem 5.1.4] again, (14) is well-approximated in probability for
almost all $\boldsymbol{\mu}$ by

$\frac{\mathrm{tr}(\hat{\mathbf{R}}^{-1})^{2}}{p\mathrm{tr}(\hat{\mathbf{R}}^{-1}\mathbf{R}\hat{\mathbf{R}}^{-1})},$
(15)

so that the ideal choice of $\hat{\mathbf{R}}$ subject to the constraint
$t\approx 1$ is the Ledoit-Wolf estimator discussed in Section III-A. This
justifies our choice of detector.

An additional feature of equation (15) is that it can be used to numerically
find the optimum in the class of diagonal-loading estimators, which is a class
of shrinkage estimators that includes LAPPW20. This can be done by simply
plugging $\hat{\mathbf{R}}=\mathbf{S}+\lambda\mathbf{I}$ in (15) and
iteratively solving for the optimal $\lambda>0$ using, for example, fminsearch
in Matlab. Although this estimator is not technically the same as LAPPW20, we
refer to it henceforth as LAPPW20 since we can calculate it and its
performance is an upper bound for LAPPW20’s.

##  IV Simulations

A key attribute of our proposed algorithm is its performance in simulation. In
this section, we show performance that is consistent with average-case
dominance over BS96, CQ10, and LAPPW20, under the model that the difference in
population means of the two samples is drawn uniformly from the sphere.

In this Section, the dimension $p$ is $200$, $n_{1}$ and $n_{2}$ are taken to
be 150, and we consider diagonal population covariance matrices
$\mathbf{R}=\mathbf{R}_{P}$ parametrized by an integer $P\in\\{0,2,4\\}$. The
eigenvalues of $\mathbf{R}_{P}$ are defined for $1\leq j\leq 40$ by

$\left(\mathbf{R}_{P}\right)_{jj}=10^{(41-j)P/40}+\epsilon_{j},\qquad(1\leq
j\leq 40),$

where $\epsilon_{j}$ is chosen i.i.d. uniformly at random from $[0,1]$. On the
other hand, for $j>40$, we set $(\mathbf{R}_{P})_{jj}=1$. Thus, largest
eigenvalue of $\mathbf{R}_{P}$ is $10^{P}$, the next 39 decrease
exponentially, and the remaining 160 are equal to unity. This choice of
spectrum mirrors the generalized spiked structure widely encountered in
sensing applications, such as radar, and $P$ corresponds to the order of the
$\mathbf{R}$’s condition number.

For each detector and each $P\in\\{0,2,4\\}$, we generate 100,000 detection
scores as follows. We generate a sub-Gaussian $p\times n_{1}$ and a $p\times
n_{2}$ data matrix $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$, both by coloring a
matrix of i.i.d. uniform random variables with mean zero and variance one.
(Simulations appear to be similar for Gaussian data matrices.) We perturb the
columns of one of these data matrices by a random vector uniformly drawn from
a sphere of radius $1$ and treat the mean of the other as $\mathbf{0}$. We
then plug the data into the Proposed, LAPPW20, BS96, CQ10, and Hotelling
detectors. The resulting 100,000 detection scores for each detector are used
to generate the ROC curves in Figures 2, 3, and 4, which correspond to
$P=0,2,4$, respectively. Figure 2 shows the case $P=0$, where BS96, LAPPW20,
CQ10, and the proposed method all perform similarly, and all dominate
Hotelling. In Figures 3 and 4, the proposed method dominates BS96, LAPPW20,
CQ10, and Hotelling, with the advantage being more pronounced for increasing
values of $P$. This is likely due to the fact that the number of large
eigenvalues in these cases challenges both the spiked assumption of LAPPW20
and the well-conditioned assumption of BS96 and CQ10, mentioned in Section II.

Figure 2: For $\mathbf{R}=\mathbf{R}_{0}$ with $O(10^{0})$ condition number,
all methods perform similarly, except Hotelling. Figure 3: For
$\mathbf{R}=\mathbf{R}_{2}$ with $O(10^{2})$ condition number, Proposed method
performs similarly to LAPPW20, which outperforms Hotelling, which outperforms
BS96 and CQ10, which lie on the chance line. Figure 4: For
$\mathbf{R}=\mathbf{R}_{4}$ with $O(10^{4})$ condition number, Proposed method
outperforms LAPPW20 and Hotelling, which outperform BS96 and CQ10, which lie
on the chance line.

##  V Conclusion

In this paper, we have proposed an alternative to Hotelling’s $T^{2}$ test
with desirable theoretical and empirical properties in the high-dimensional
setting. In simulation, the proposed test appears to dominate the state-of-
the-art alternatives LAPPW20, CQ10, and BS96, likely due in part to the
comparatively less forgiving assumptions they require. We provided an
intuitive argument that our test has an asymptotically constant false-alarm
rate and maximal detection power among shrinkage-based Hotelling-type tests
under mild assumptions. We conjecture that a more rigorous foundation for this
theory can be established by finding conditions under which (10) and
strengthenings of (8) hold, which we intend to investigate in future work.

## Acknowledgments

This work was supported by the United States Air Force Sensors Directorate,
AFOSR grant 19RYCOR036, ARO grant W911NF-15-1-0479, and the DoD SMART
Scholarship SEED grant. However, the views and opinions expressed in this
article are those of the authors and do not necessarily reflect the official
policy or position of any agency of the U.S. government. Examples of analysis
performed within this article are only examples. Assumptions made within the
analysis are also not reflective of the position of any U.S. Government
entity. The Public Affairs approval number of this document is AFRL-2022-2605.

## References

  * [1] T. W. Anderson, “Asymptotic theory for principal component analysis,” _The Annals of Mathematical Statistics_ , vol. 34, no. 1, pp. 122–148, 1963\.
  * [2] D. Paul, “Asymptotics of sample eigenstructure for a large dimensional spiked covariance model,” _Statistica Sinica_ , pp. 1617–1642, 2007.
  * [3] I. M. Johnstone, “On the distribution of the largest eigenvalue in principal components analysis,” _Annals of Statistics_ , pp. 295–327, 2001.
  * [4] Z. Bai and H. Saranadasa, “Effect of high dimension: By an example of a two sample problem,” _Statistica Sinica_ , pp. 311–329, 1996.
  * [5] S. X. Chen and Y.-L. Qin, “A two-sample test for high-dimensional data with applications to gene-set testing,” _The Annals of Statistics_ , vol. 38, no. 2, pp. 808–835, 2010.
  * [6] H. Li, A. Aue, D. Paul, J. Peng, and P. Wang, “An adaptable generalization of Hotelling’s T2 test in high dimension,” _The Annals of Statistics_ , vol. 48, no. 3, pp. 1815–1847, 2020.
  * [7] O. Ledoit and M. Wolf, “Analytical nonlinear shrinkage of large-dimensional covariance matrices,” _The Annals of Statistics_ , vol. 48, no. 5, pp. 3043–3065, 2020.
  * [8] C. Stein, “Estimation of a covariance matrix, Rietz lecture,” in _39th Annual Meeting IMS, Atlanta, GA, 1975_ , 1975.
  * [9] ——, “Lectures on the theory of estimation of many parameters,” _Journal of Soviet Mathematics_ , vol. 34, no. 1, pp. 1373–1403, 1986.
  * [10] B. D. Robinson, R. Malinas, and A. O. Hero, “Space-time adaptive detection at low sample support,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 2939–2954, 2021.
  * [11] J. Namdari, D. Paul, and L. Wang, “High-dimensional linear models: A random matrix perspective,” _Sankhya A_ , vol. 83, no. 2, pp. 645–695, 2021.
  * [12] D. L. Donoho, M. Gavish, and I. M. Johnstone, “Optimal shrinkage of eigenvalues in the spiked covariance model,” _Annals of Statistics_ , vol. 46, no. 4, p. 1742, 2018.
  * [13] O. Ledoit and S. Péché, “Eigenvectors of some large sample covariance matrix ensembles,” _Probability Theory and Related Fields_ , vol. 151, no. 1-2, pp. 233–264, 2011.
  * [14] O. Ledoit and M. Wolf, “Nonlinear shrinkage estimation of large-dimensional covariance matrices,” _The Annals of Statistics_ , vol. 40, no. 2, pp. 1024–1060, 2012.
  * [15] ——, “Optimal estimation of a large-dimensional covariance matrix under Stein’s loss,” _Bernoulli_ , vol. 24, no. 4B, pp. 3791–3832, 2018.
  * [16] F. Benaych-Georges and A. Knowles, “Lectures on the local semicircle law for Wigner matrices,” _arXiv preprint arXiv:1601.04055_ , 2016.
  * [17] L. Erdős, A. Knowles, and H.-T. Yau, “Averaging fluctuations in resolvents of random band matrices,” in _Annales Henri Poincaré_ , vol. 14, no. 8. Springer, 2013, pp. 1837–1926.
  * [18] G. Cipolloni, L. Erdős, and D. Schröder, “Eigenstate thermalization hypothesis for Wigner matrices,” _Communications in Mathematical Physics_ , vol. 388, no. 2, pp. 1005–1048, 2021.
  * [19] R. Vershynin, _High-Dimensional Probability: An Introduction with Applications in Data Science_. Cambridge University Press, 2018, vol. 47.