# Fair Streaming Principal Component Analysis: Statistical and Algorithmic
Viewpoint

Junghyun Lee Hanseul Cho∗ Se-Young Yun Chulhee Yun
Kim Jaechul Graduate School of AI, KAIST
{jh_lee00, jhs4015, yunseyoung<EMAIL_ADDRESS>Equal contributions

###### Abstract

Fair Principal Component Analysis (PCA) is a problem setting where we aim to
perform PCA while making the resulting representation fair in that the
projected distributions, conditional on the sensitive attributes, match one
another. However, existing approaches to fair PCA have two main problems:
theoretically, there has been no statistical foundation of fair PCA in terms
of learnability; practically, limited memory prevents us from using existing
approaches, as they explicitly rely on full access to the entire data. On the
theoretical side, we rigorously formulate fair PCA using a new notion called
_probably approximately fair and optimal_ (PAFO) learnability. On the
practical side, motivated by recent advances in streaming algorithms for
addressing memory limitation, we propose a new setting called _fair streaming
PCA_ along with a memory-efficient algorithm, fair noisy power method (FNPM).
We then provide its statistical guarantee in terms of PAFO-learnability, which
is the first of its kind in fair PCA literature. Lastly, we verify the
efficacy and memory efficiency of our algorithm on real-world datasets.

##  1 Introduction

Algorithmic fairness ensures that machine learning algorithms do not propagate
nor exacerbate bias, which may lead to discriminatory decision-making (Barocas
and Selbst, 2016) and thus has been a very active area of research. This has
direct implications in our everyday life, including but not limited to
criminal justice (Kirchner et al., 2016), education (Kizilcec and Lee, 2021),
and more. See Mehrabi et al. (2021) for a comprehensive survey of bias and
fairness in machine learning.

Often, one needs to consider fairness for a large number of high-dimensional
data points. One of the standard tools for dealing with such high-dimensional
data is PCA (Pearson, 1901; Hotelling, 1933), a classical yet still one of the
most popular algorithms for performing interpretable dimensionality reduction.
It has been adapted as a baseline and/or standard tool in exploratory data
analysis, whose application ranges from natural sciences, engineering (Abdi
and Williams, 2010; Jolliffe and Cadima, 2016), and even explainable AI (Tjoa
and Guan, 2021; Li et al., 2023). Due to its ubiquity and wide applicability,
several works study defining fairness in PCA and developing a fair variant of
it. A recent line of research (Lee et al., 2022; Olfat and Aswani, 2019;
Kleindessner et al., 2023) defines PCA fairness in the context of fair
representation (Zemel et al., 2013) in that the projected group conditional
distributions should match.

However, existing fair PCA approaches suffer from two problems. Theoretically,
they provide no statistical foundation of fair PCA or guarantees for their
algorithms. By statistical foundation, we mean the usual PAC-learnability
(Shalev-Schwartz and Ben-David, 2014) guarantees, e.g., sample complexity for
ensuring optimality in explained variance and fairness constraint with high
probability.

On top of that, the second problem arises from a practical viewpoint: memory
limitation. All the aforementioned fair PCA algorithms assume that the learner
can store the entire data points and incurs memory complexity of order at
least ${\mathcal{O}}(d\max(N,d))$, where $d$ is the dimensionality of the data
and $N$ is the number of data points. As memory limitation is often a critical
bottleneck in deploying machine learning algorithms (Mitliagkas et al., 2013),
as much as fairness is important, it is also paramount that imposing fairness
to PCA does not incur too much memory overhead. A popular approach to mitigate
such memory limitation for PCA is to consider the one-pass, streaming setting.
In this setting, each data point is revealed to the learner sequentially, each
point is irretrievably gone unless she explicitly stores it, and she can use
only ${\mathcal{O}}(dk)$ memory, with $k$ being the target dimension of
projection. Indeed, without the fairness constraint, streaming PCA has been
studied extensively; see Balzano et al. (2018) and references therein.

In this work, we address both problems in a principled manner. Our
contributions are as follows:

  * •

We provide an alternative formulation of fair PCA based on the “Null It Out”
approach (Section 3). Based on the new formulation, we introduce the concept
of probably approximately fair and optimal (PAFO)-learnability to formalize
the problem of fair PCA (Section 4).

  * •

To address the memory limitation, we propose a new problem setting called fair
streaming PCA, as well as fair noisy power method (FNPM), a simple yet memory-
efficient algorithm based on the noisy power method. (Section 5). We note that
our algorithm incurs a much lower memory complexity even compared to the most
efficient variant of fair PCA proposed by Kleindessner et al. (2023).

  * •

We then prove that our algorithm achieves the PAFO-learnability for fair
streaming PCA (Section 6). Such statistical guarantee is the first of its kind
in fair PCA literature.

  * •

Lastly, we empirically validate our algorithm on CelebA and UCI datasets.
Notably, we run FNPM on the original full-resolution CelebA dataset on which
existing fair PCA algorithms fail due to high memory requirements. It shows
turning such a non-streaming setting into a streaming setting and applying our
algorithm allows one to bypass the memory limitation (Section 7).

##  2 Preliminaries

#### Notation.

For $\ell\geq 1$, let ${\bm{I}}_{\ell}$ be the identity matrix of size
$\ell\times\ell$. For $k<d$, we bring the Stiefel manifold
$\operatorname{St}(d,k)=\\{{\bm{A}}\in{\mathbb{R}}^{d\times
k}:{\bm{A}}^{\intercal}{\bm{A}}={\bm{I}}_{k}\\}$, which is the collection of
all rank-$k$ semi-orthogonal matrices. We denote an orthonormal column basis
of a (full column rank) matrix ${\bm{M}}\in\mathbb{R}^{d\times k}$ obtained by
QR decomposition as $\operatorname{{\tt
QR}}({\bm{M}})\in\operatorname{St}(d,k)$ and denote its column space by
$\mathrm{col}({\bm{M}})$. Also, for ${\bm{A}}\in\operatorname{St}(d,k)$, we
denote the orthogonal projection matrix to
$\mathrm{col}({\bm{A}})^{\perp}=\mathrm{null}({\bm{A}}^{\intercal})$ as
$\bm{\Pi}_{\bm{A}}^{\perp}={\bm{I}}_{d}-{\bm{A}}{\bm{A}}^{\intercal}={\bm{I}}_{d}-\bm{\Pi}_{\bm{A}}$.
Moreover, we denote the collection of all possible $d$-dimensional probability
distributions as ${\mathcal{P}}_{d}$. For a zero-mean random matrix
${\bm{Z}}$, its (scalar-valued) variance is defined as
$\mathrm{Var}({\bm{Z}})=\max\left(\left\|\mathbb{E}\left[{\bm{Z}}{\bm{Z}}^{\intercal}\right]\right\|_{2},\left\|\mathbb{E}\left[{\bm{Z}}^{\intercal}{\bm{Z}}\right]\right\|_{2}\right).$
In general,
$\mathrm{Var}({\bm{Z}})=\mathrm{Var}({\bm{Z}}-\mathbb{E}[{\bm{Z}}])$. Lastly,
we use the usual ${\mathcal{O}}$, $\Omega$, and $\Theta$ notations for
asymptotic analyses, where tildes ($\tilde{{\mathcal{O}}}$, $\tilde{\Omega}$,
and $\tilde{\Theta}$, resp.) are used for hiding logarithmic factors.

#### Setup.

Assume that the sensitive attribute variable, which we will be imposing
fairness, is binary, denoted by $a\in\\{0,1\\}$. For each group $a$, let
${\mathcal{D}}_{a}$ be a $d$-dimensional distribution of mean $\bm{\mu}_{a}$
and covariance $\bm{\Sigma}_{a}$, both of which are assumed to be well-
defined. We often call them group-conditional mean and covariance,
respectively. With a fixed, unknown mixture parameter $p\in(0,1)$, let us
denote the total data distribution as
${\mathcal{D}}:=p{\mathcal{D}}_{0}+(1-p){\mathcal{D}}_{1}$. Equivalently, the
sensitive attribute follows $a\sim\mathrm{Bernoulli}(p)$, and the conditional
random variable ${\bm{x}}|a$ is sampled from ${\mathcal{D}}_{a}$. In that
case, $\bm{\mu}_{a}=\mathbb{E}[{\bm{x}}|a]$ and
$\bm{\Sigma}_{a}=\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}|a]-\bm{\mu}_{a}\bm{\mu}_{a}^{\intercal}$.
We often write $p_{0}=1-p$ and $p_{1}=p$ for brevity. We also define the mean
difference ${\bm{f}}:=\bm{\mu}_{1}-\bm{\mu}_{0}$ and the second moment
difference
${\bm{S}}:=\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}|a=1]-\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}|a=0]=\bm{\Sigma}_{1}-\bm{\Sigma}_{0}+\bm{\mu}_{1}\bm{\mu}_{1}^{\intercal}-\bm{\mu}_{0}\bm{\mu}_{0}^{\intercal}$.
Accordingly, denote the true mean and covariance of ${\mathcal{D}}$ as
$\bm{\mu}$ and $\bm{\Sigma}$, respectively. For simplicity, let us assume that
${\mathcal{D}}$ is centered, i.e., $\bm{\mu}={\bm{0}}$; note that this does
not mean that the group conditional distributions ${\mathcal{D}}_{a}$’s are
centered.

#### PCA.

In the _offline_ setting, the full covariance matrix $\bm{\Sigma}$ is given
which is often a sample covariance matrix
$\frac{1}{n}\sum_{i=1}^{n}{\bm{x}}_{i}{\bm{x}}_{i}^{\top}$ for $n$ data points
${\bm{x}}_{1},\dots,{\bm{x}}_{n}$. The goal of vanilla (offline) PCA
(Hotelling, 1933; Pearson, 1901) is to compute the loading matrix
${\bm{V}}\in{\mathbb{R}}^{d\times k}$ that preserves as much variance as
possible after projecting $\bm{\Sigma}$ via ${\bm{V}}$, i.e., maximize
$\mathrm{tr}(\bm{\Pi}_{\bm{V}}\bm{\Sigma})=\mathrm{tr}({\bm{V}}^{\intercal}\bm{\Sigma}{\bm{V}})$.
Here, $k<d$ is the target dimension to which the data’s dimensionality $d$ is
to be reduced and is chosen by the learner. We additionally consider the
constraint ${\bm{V}}\in\operatorname{St}(d,k)$ to ensure that the resulting
coordinate after the transformation is orthogonal and thus amenable to various
statistical interpretations (Johnson and Wichern, 2008). Without any fairness
constraint, Eckart-Young theorem (Eckart and Young, 1936) implies that the
solution is characterized as a matrix whose columns are the top-$k$
eigenvectors of $\bm{\Sigma}$.

#### Fair PCA.

Recently, it has been suggested that performing vanilla PCA on real-world
datasets may exhibit bias, making the final outputted projection “unfair”. As
is often the case, there can be multiple definitions of fairness in PCA, but
the following two are the most popular: equalizing reconstruction losses
(Samadi et al., 2018; Tantipongpipat et al., 2019; Vu et al., 2022; Kamani et
al., 2022), or equalizing the projected distributions (Olfat and Aswani, 2019;
Lee et al., 2022; Kleindessner et al., 2023) from the perspective of fair
representation (Zemel et al., 2013); we focus on the latter one.

##  3 An Alternative Approach to Fair PCA

###  3.1 “Null It Out” Formulation of Fair PCA

In this work, we consider fair PCA as learning fair representation (Zemel et
al., 2013). The goal is to preserve as much variance as possible while
obfuscating any information regarding the sensitive attribute. To this end, we
take the “Null It Out” approach as proposed in Ravfogel et al. (2020).
Intuitively, we want to nullify the directions in which the sensitive
attribute $a$ can be inferred, and in this work, we consider two such unfair
directions: mean difference ${\bm{f}}$ and eigenvectors of second moment
difference ${\bm{S}}$. To give the learner flexibility in choosing the trade-
off between fairness and performance (measured in explained variance), let
$m\geq 1$ be the number of top eigenvectors of ${\bm{S}}$ to nullify. Thus,
the learner is nullifying at most $(m+1)$-dimensional subspace that is unfair
with respect to $a$, which we refer to as the unfair subspace. Precisely, we
formulate our fair PCA as follows:

$\max_{{\bm{V}}\in\operatorname{St}(d,k)}\mathrm{tr}({\bm{V}}^{\intercal}\bm{\Sigma}{\bm{V}}),\quad\text{subject
to
}\mathrm{col}({\bm{V}})\subset\mathrm{col}([{\bm{P}}_{m}|{\bm{f}}])^{\perp},$
(1)

where $d$ is the data dimensionality, $k$ is the target dimension, and the
columns of ${\bm{P}}_{m}\in\operatorname{St}(d,m)$ is top-$m$ orthonormal
eigenvectors of ${\bm{S}}$.

###  3.2 An Explicit Characterization for Solution of Fair PCA

To first construct the unfair subspace that is spanned by ${\bm{f}}$ as well
as ${\bm{P}}_{m}$, let us define ${\bm{U}}\in\operatorname{St}(d,m^{\prime})$
to be the orthogonal matrix whose columns form a basis of
$\mathrm{col}([{\bm{P}}_{m}|{\bm{f}}])$. Then, ${\bm{U}}$ has a closed form as
follows: $m^{\prime}=m$ if ${\bm{f}}\in\mathrm{col}({\bm{P}}_{m})$ and
$m^{\prime}=m+1$ otherwise, and

$\displaystyle{\bm{U}}$
$\displaystyle=\begin{cases}{\bm{P}}_{m},\quad&\text{if
${\bm{f}}\in\mathrm{col}({\bm{P}}_{m})$,}\\\ \operatorname{{\tt
QR}}([{\bm{P}}_{m}|{\bm{f}}])=\left[{\bm{P}}_{m}\,\middle|\,\frac{{\bm{g}}}{\left\|{\bm{g}}\right\|_{2}}\right],\quad&\text{otherwise,}\end{cases}$
(2)

where
${\bm{g}}=\bm{\Pi}_{{\bm{P}}_{m}}^{\perp}{\bm{f}}\in\mathrm{col}({\bm{P}}_{m})^{\perp}$.
Note that ${\bm{g}}$ is a vector in a direction that ${\bm{f}}$ is projected
onto
$\mathrm{col}({\bm{P}}_{m})^{\perp}=\mathrm{null}({\bm{P}}_{m}^{\intercal})$.
For this ${\bm{U}}$, our constraint in (1) can be interpreted as an equivalent
nullity constraint $\bm{\Pi}_{{\bm{U}}}{\bm{V}}=\bm{0}$:

$\max_{{\bm{V}}\in\operatorname{St}(d,k)}\mathrm{tr}({\bm{V}}^{\intercal}\bm{\Sigma}{\bm{V}}),\quad\text{subject
to }\bm{\Pi}_{{\bm{U}}}{\bm{V}}=\bm{0}.$ (3)

The above is equivalent to the following problem without any constraint other
than semi-orthgonality:

$\max_{{\bm{V}}\in\operatorname{St}(d,k)}\mathrm{tr}\left({\bm{V}}^{\intercal}\bm{\Pi}_{\bm{U}}^{\perp}\bm{\Sigma}\bm{\Pi}_{\bm{U}}^{\perp}{\bm{V}}\right),$
(4)

which is basically the vanilla $k$-PCA problem of a matrix
$\bm{\Pi}_{\bm{U}}^{\perp}\bm{\Sigma}\bm{\Pi}_{\bm{U}}^{\perp}$. Therefore, a
top-$k$ orthonormal column basis of this matrix is indeed a solution of our
problem (4).

###  3.3 Comparison to the Existing Covariance Matching Constraint

Previous works on fair PCA (Kleindessner et al., 2023; Olfat and Aswani, 2019)
consider an exact covariance-matching constraint (namely,
${\bm{V}}^{\intercal}(\bm{\Sigma}_{1}-\bm{\Sigma}_{0}){\bm{V}}=\bm{0}$). In
fact, this is equivalent to the condition
${\bm{V}}^{\intercal}{\bm{S}}{\bm{V}}=\bm{0}$ under the mean-matching
constraint ${\bm{f}}^{\intercal}{\bm{V}}={\bm{0}}$, which can be derived as

$\displaystyle{\bm{V}}^{\intercal}(\bm{\Sigma}_{1}-\bm{\Sigma}_{0}){\bm{V}}$
$\displaystyle={\bm{V}}^{\intercal}\left(\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}|a=1]-\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}|a=0]\right){\bm{V}}-{\bm{V}}^{\intercal}\left(\bm{\mu}_{1}\bm{\mu}_{1}^{\intercal}-\bm{\mu}_{0}\bm{\mu}_{0}^{\intercal}\right){\bm{V}}$
$\displaystyle={\bm{V}}^{\intercal}{\bm{S}}{\bm{V}}-{\bm{V}}^{\intercal}\left({\bm{f}}\bm{\mu}_{1}^{\intercal}+\bm{\mu}_{0}{\bm{f}}^{\intercal}\right){\bm{V}}={\bm{V}}^{\intercal}{\bm{S}}{\bm{V}}=\bm{0}.$

One immediate problem with this is that the constraint may be infeasible
depending on the choice of $\bm{\Sigma}_{0}$, $\bm{\Sigma}_{1}$, or
${\bm{S}}$; e.g., when $\bm{\Sigma}_{1}-\bm{\Sigma}_{0}$ is positive definite.
For this reason, Kleindessner et al. (2023); Olfat and Aswani (2019) propose
relaxations of the fairness constraints but provide no further discussions on
its impact on statistical guarantees. On the contrary, our formulation is
always feasible without any need for relaxation, allowing us to consider a
rigorous definition of fair PCA (Definition 2) for the first time in fair PCA
literature.

##  4 Statistical Viewpoint: PAFO-Learnability of PCA

As all the distribution statistics ($\bm{\Sigma},p,\cdots$) are unknown, the
learner, given some finite number of samples, must learn all of them and solve
fair PCA. In supervised learning, such a problem is often formalized in a PAC-
learnability framework (Shalev-Schwartz and Ben-David, 2014). In the context
of PAC-learnability for unsupervised learning settings, TV-learning, which is
the task of learning distribution, has been mainly considered so far
(Ananthakrishnan et al., 2021; Hopkins et al., 2023). However, unlike TV-
learning, it is unnecessary to learn the whole distribution in fair PCA;
moreover, fair PCA has the fairness constraint
$\bm{\Pi}_{{\bm{U}}}{\bm{V}}=\bm{0}$ to be satisfied. Inspired by the
unsupervised PAC-learnability as well as constrained PAC-learnability (Chamon
and Ribeiro, 2020), we propose a new notion of learnability for fair PCA,
called PAFO (Probably Approximately Fair and Optimal) learnability, as
follows:

######  Definition 1 (Projection Learner).

A projection learner is a function that takes $k\geq 1$ and $d$-dimensional
samples as input and outputs a loading matrix
${\bm{V}}\in\operatorname{St}(d,k)$.

######  Definition 2 (PAFO-Learnability of PCA).

Let $d,k,m$ be integers such that $1\leq k<d$ and $m<d$. We say that
${\mathcal{F}}_{d}\subset{\mathcal{P}}_{d}\times{\mathcal{P}}_{d}\times(0,1)$
is PAFO-learnable for PCA if there exists a function
$N_{{\mathcal{F}}_{d}}:(0,1)^{3}\rightarrow{\mathbb{N}}$ and a projection
learner ${\mathcal{A}}$ satisfying the following:

RGB]230,230,230 For every $\varepsilon_{\rm o},\varepsilon_{\rm
f},\delta\in(0,1)$ and
$({\mathcal{D}}_{0},{\mathcal{D}}_{1},p)\in{\mathcal{F}}_{d}$, when running
${\mathcal{A}}$ on $N\geq N_{{\mathcal{F}}_{d}}(\varepsilon_{\rm
o},\varepsilon_{\rm f},\delta)$ i.i.d. samples from
${\mathcal{D}}:=p{\mathcal{D}}_{1}+(1-p){\mathcal{D}}_{0}$ of the form
$(a,{\bm{x}})$, ${\mathcal{A}}$ returns ${\bm{V}}$ s.t., with probability at
least $1-\delta$ (over the draws of the $N$ samples),
$\mathrm{tr}\left({\bm{V}}^{\intercal}\bm{\Sigma}{\bm{V}}\right)\geq\mathrm{tr}\left({{\bm{V}}_{\star}}^{\intercal}\bm{\Sigma}{\bm{V}}_{\star}\right)-\varepsilon_{\rm
o},\quad\lVert\bm{\Pi}_{{\bm{U}}}{\bm{V}}\rVert_{2}\leq\varepsilon_{\rm f},$
(5) where ${\bm{U}}$ is as defined in Eqn. (2) and ${\bm{V}}_{\star}$ is any
solution to Eqn. (4) (with prescribed $k$ and $m$).

Like in the usual PAC-learnability, $N_{{\mathcal{F}}_{d}}$ is referred to as
the sample complexity of fair PCA. Observe how the optimality is measured
w.r.t. the optimal solution of fair PCA, not the vanilla PCA. Also, the two
conditions are not overlapping: vanilla PCA (overly) satisfies
$\varepsilon_{\rm o}$-optimality in explained variance but does not satisfy
$\varepsilon_{\rm f}$-optimality in fairness, and vice versa for a poorly
chosen ${\bm{V}}\in St(d,k)$ with
$\mathrm{col}({\bm{V}})\subseteq\mathrm{col}([{\bm{P}}_{m}|{\bm{f}}])^{\perp}$.

##  5 Algorithmic Viewpoint: Fair Streaming PCA

We now introduce a new problem setting, fair streaming PCA. In this setting,
the learner receives a stream of pairs
$(a_{t},{\bm{x}}_{t})\in\\{0,1\\}\times{\mathbb{R}}^{d}$ sequentially. Note
that the sensitive attribute information $a_{t}$ is also available at each
time-step; this is commonly assumed when considering fairness in streaming
setting (El Halabi et al., 2020; Bera et al., 2022). Precisely, we assume the
following model of the data generation process: at each time-step $t$, a
sensitive attribute is chosen as $a_{t}\sim\mathrm{Bernoulli}(p)$, then the
data is sampled from the corresponding sensitive group’s conditional
distribution ${\bm{x}}_{t}\mid a_{t}\sim{\mathcal{D}}_{s_{t}}$. Importantly,
as done in previous streaming PCA literature (Mitliagkas et al., 2013), we
assume that the learner has only ${\mathcal{O}}(dk)$ memory, where $d$ is the
data dimension and $k$ is the target dimension. We can formally define the
PAFO-learnability in this streaming setting:

######  Definition 3.

We say that
${\mathcal{F}}_{d}\subseteq{\mathcal{P}}_{d}\times{\mathcal{P}}_{d}\times(0,1)$
is PAFO-learnable for streaming PCA if the projection learner ${\mathcal{A}}$
for which Definition 2 holds uses only ${\mathcal{O}}(dk)$ memory for
streaming data.

###  5.1 Our Algorithm: Fair Noisy Power Method (FNPM)

One only needs to estimate ${\bm{U}}$ to use the off-the-shelf streaming PCA
algorithm. As ${\bm{U}}$ is of size $d\times m$, storing its estimate is no
problem for the memory constraint as long as $m={\mathcal{O}}(k)$. Naturally,
we proceed via a two-stage approach; first, estimate ${\bm{U}}$ sufficiently
well, then with the fixed estimate of ${\bm{U}}$, apply the noisy power method
(Mitliagkas et al., 2013; Hardt and Price, 2014) for ${\bm{V}}$.

For estimating ${\bm{U}}$, one needs to estimate ${\bm{f}}$ and
${\bm{P}}_{m}$. Estimating ${\bm{f}}$ can be done using the usual cumulative
averaging. As for ${\bm{P}}_{m}$, we can consider the two main approaches for
streaming PCA: Oja’s method (Oja, 1982; Oja and Karhunen, 1985; Huang et al.,
2021) and noisy power method (NPM) (Mitliagkas et al., 2013; Hardt and Price,
2014). We first show that Oja’s method is inapplicable for our purpose, as it
may ignore some eigenvectors corresponding to negative (but large in
magnitude) eigenvalues of ${\bm{S}}$. For instance, if
${\bm{S}}=-2{\bm{e}}_{1}{\bm{e}}_{1}^{\intercal}+{\bm{e}}_{2}{\bm{e}}_{2}^{\intercal}+4{\bm{e}}_{3}{\bm{e}}_{3}^{\intercal}$
with ${\bm{e}}_{i}$ being the standard basis vectors, then Oja’s method with
$m=2$ would yield $[{\bm{e}}_{2}|{\bm{e}}_{3}]$ when we actually want
$[{\bm{e}}_{1}|{\bm{e}}_{3}]$. For the same reason, simply shifting the
eigenvalue spectrum by considering ${\bm{S}}+\lVert{\bm{S}}\rVert_{2}{\bm{I}}$
does not work. Thus we apply NPM for estimating ${\bm{P}}_{m}$ in our case,
which is known to converge as long as the singular value gap of ${\bm{P}}_{m}$
is large enough and norms of the noise matrices at each iterate are properly
bounded (Hardt and Price, 2014; Balcan et al., 2016).

1 Input: $m$, Block size $b$, Number of iterations $T$;

2 Output: A matrix $\widehat{\bm{U}}$ with orthonormal columns;

3 ${\bm{W}}_{0}=\operatorname{{\tt QR}}({\mathcal{N}}(0,1)^{d\times m})$;

4
$(\overline{{\bm{m}}}^{(0)},\overline{{\bm{m}}}^{(1)},B^{(0)},B^{(1)})=(\bm{0}_{d},\bm{0}_{d},0,0)$;

5 for  _$t\in[T]$_ do

6 Receive $\\{(a_{i},{\bm{x}}_{i})\\}_{i=(t-1)b+1}^{tb}$;

7 foreach  _$a\in\\{0,1\\}$_ do

8 Compute $b_{t}^{(a)},{\bm{m}}_{t}^{(a)},{\bm{C}}_{t}^{(a)}$ as Eqn. (6);

9
$\overline{{\bm{m}}}^{(a)}\\!\leftarrow\\!\frac{B^{(a)}}{B^{(a)}+b_{t}^{(a)}}\overline{{\bm{m}}}^{(a)}+\frac{b_{t}^{(a)}}{B^{(a)}+b_{t}^{(a)}}{\bm{m}}_{t}^{(a)}$;

10 $B^{(a)}\\!\leftarrow B^{(a)}\\!+b_{t}^{(a)}$;

11

12 end foreach

13 ${\bm{W}}_{t}=\operatorname{{\tt
QR}}\left({\bm{C}}_{t}^{(1)}-{\bm{C}}_{t}^{(0)}\right)$;

14

15 end for

16$\widehat{{\bm{f}}}\leftarrow\overline{{\bm{m}}}^{(1)}-\overline{{\bm{m}}}^{(0)}$;

17
$\widehat{{\bm{g}}}\leftarrow\widehat{{\bm{f}}}-{\bm{W}}_{T}{\bm{W}}_{T}^{\intercal}\widehat{{\bm{f}}}$;

18 if  _$\|\widehat{{\bm{g}}}\|_{2}=0$_ then

19 $\widehat{\bm{U}}={\bm{W}}_{T}$

20 else

21
$\widehat{\bm{U}}=\left[{\bm{W}}_{T}\,\big{|}\,\frac{\widehat{{\bm{g}}}}{\|\widehat{{\bm{g}}}\|_{2}}\right]$

22 end if

return $\widehat{\bm{U}}$

Algorithm 1 UnfairSubspace

1 Input: $k$, Block sizes $\mathcal{B}$, $b$, Numbers of iterations
${\mathcal{T}}$, $T$;

2 Output: ${\bm{V}}_{\mathcal{T}}\in\operatorname{St}(d,k)$;

3

$\widehat{{\bm{U}}}\leftarrow{\tt UnfairSubspace}(b,T)$;

// Algorithm 1

4 ${\bm{V}}_{0}\leftarrow\operatorname{{\tt QR}}({\mathcal{N}}(0,1)^{d\times
k})$;

5 for  _$\tau\in[{\mathcal{T}}]$_ do

6
${\bm{V}}_{\tau}\leftarrow{\bm{V}}_{\tau-1}-\widehat{{\bm{U}}}\widehat{{\bm{U}}}^{\intercal}{\bm{V}}_{\tau-1}$;

7 Receive
$\\{(\ast,\tilde{\bm{x}}_{j})\\}_{j=(\tau-1){\mathcal{B}}+1}^{\tau{\mathcal{B}}}$;

8
${\bm{V}}_{\tau}\leftarrow\frac{1}{{\mathcal{B}}}\sum_{j=1}^{\mathcal{B}}\tilde{\bm{x}}_{j}\tilde{\bm{x}}_{j}^{\intercal}{\bm{V}}_{\tau}$;

9 ${\bm{V}}_{\tau}\leftarrow\operatorname{{\tt
QR}}\left({\bm{V}}_{\tau}-\widehat{{\bm{U}}}\widehat{{\bm{U}}}^{\intercal}{\bm{V}}_{\tau}\right)$;

10

11 end for

return ${\bm{V}}_{\mathcal{T}}$

Algorithm 2 Fair NPM

#### Description of the algorithms.

The pseudocode of our algorithm is shown in Algorithms 1 and 2. The goal of
Algorithm 1 is to estimate ${\bm{U}}=\left[{\bm{P}}_{m}\ |\
\frac{{\bm{g}}}{\left\|{\bm{g}}\right\|_{2}}\right]$ in Eqn. (2) as accurately
as possible. Lines 1–1 do the estimation of ${\bm{P}}_{m}$ and
${\bm{g}}=\bm{\Pi}_{{\bm{P}}_{m}}^{\perp}{\bm{f}}$; line 1 is the NPM to find
${\bm{P}}_{m}$, lines 1 and 1 are the estimation of ${\bm{f}}$ and ${\bm{g}}$
respectively, and line 1 is the concatenation of the estimates of
${\bm{P}}_{m}$ and ${{\bm{g}}}/\\!\left\|{\bm{g}}\right\|_{2}$. Especially at
line 1, the algorithm determines whether to incorporate mean difference by
checking $\hat{{\bm{g}}}$, which can be proved to be correct, _i.e._ ,
${\bm{g}}={\bm{0}}$ if and only if $\hat{{\bm{g}}}={\bm{0}}$ with high
probability. With the estimated ${\bm{U}}$ from Algorithm 1, Algorithm 2
performs the usual NPM on
$\bm{\Pi}_{\bm{U}}^{\perp}\bm{\Sigma}\bm{\Pi}_{\bm{U}}^{\perp}$, as in Eqn.
(4). The memory complexity of Algorithm 2 is ${\mathcal{O}}(d\max(m,k))$,
since we do not have to store all $b$ or ${\mathcal{B}}$ data points at each
time, and all the operations used can be implemented in a manner that conforms
to the memory limitation; the full pseudocodes are provided in Appendix B.

At time step $t$ of Algorithm 1, for each $a\in\\{0,1\\}$, $b_{t}^{(a)}$ is
the number of data points ${\bm{x}}_{i}$’s such that $a_{i}=a$;
${\bm{m}}_{t}^{(a)}$ is the term used for estimation of the group-wise sample
mean of ${\bm{x}}_{i}$’s; ${\bm{C}}_{t}^{(a)}$ is used for the group-wise
sample second moment. Their forms are as follows and can be computed
incrementally in the streaming setup:

$\displaystyle
b_{t}^{(a)}\\!=\\!\\!\sum_{i=b(t-1)+1}^{bt}\\!\\!\mathbbm{1}_{[a_{i}=a]},\quad{\bm{m}}_{t}^{(a)}\\!=\\!\\!\sum_{i=b(t-1)+1}^{bt}\\!\\!\frac{\mathbbm{1}_{[a_{i}=a]}}{b_{t}^{(a)}}{\bm{x}}_{i},\quad{\bm{C}}_{t}^{(a)}\\!=\\!\\!\sum_{i=b(t-1)+1}^{bt}\\!\\!\frac{\mathbbm{1}_{[a_{i}=a]}}{b_{t}^{(a)}}{\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}{\bm{W}}_{t-1},$
(6)

where we set the last two quantities to ${\bm{0}}$ when $b_{t}^{(a)}=0$.

Note that as $b_{t}^{(a)}$ itself is random, this presents some technical
challenges in the proofs of the theoretical guarantees. For instance, the
above estimators for the mean and covariance are biased. Still, by properly
using peeling argument and matrix concentration inequalities as well as
perturbation theories (Golub and Loan, 2013; Tropp, 2015), we could
sufficiently bound their errors. Informally speaking, we show that line 11
corresponds to the noisy power method for the matrix ${\bm{S}}$, which incurs
the memory complexity ${\mathcal{O}}(dm)$.

###  5.2 Previous Approaches are not Suitable for Streaming Setup

All the existing approaches to fair PCA require the full knowledge of
$\bm{\mu}_{1}-\bm{\mu}_{0}$ and $\bm{\Sigma}_{1}-\bm{\Sigma}_{0}$, or even the
full data matrix ${\bm{X}}$. Olfat and Aswani (2019) need ${\bm{f}}$ and
${\bm{S}}$ to formulate the convex matrix constraints for their semi-positive
definite programming (SDP), which is then solved with commercial SDP solver;
Lee et al. (2022) need ${\bm{X}}$ to compute the derivative of their maximum
mean discrepancy (MMD) penalty term, which requires ${\mathcal{O}}(d^{2})$
memory to compute kernel Gram matrices. One may hope that the PCA-type
approach taken by Kleindessner et al. (2023) may be easily extendable to our
streaming setting by just using the standard techniques (e.g., using matrix-
vector products) from streaming PCA (Mitliagkas et al., 2013). Indeed, they
also propose a similar relaxation of the covariance constraint, leading to a
closed-form solution. However, their formulation is not memory-efficient and,
more importantly, is not as applicable to streaming settings as our
formulation; see Appendix D for more discussions.

###  5.3 Extension to Multiple Sensitive Groups/Attributes

So far, we have considered a single and binary-sensitive attribute for
simplicity. Here, we briefly discuss extending our framework to multiple non-
binary sensitive attributes. Suppose that there are $\ell$ different
attributes (e.g., gender and race), and for each $r\in[\ell]$, the $r$-th
sensitive attribute has $g_{r}$ groups (e.g., male, female, and non-binary:
$g_{\rm gender}=3$ in this case). We first describe the new data generation
process: at each time $t$, we obtain
$((a_{1,t},\cdots,a_{\ell,t}),{\bm{x}}_{t})$, where for each $r$,
$a_{r,t}\sim\operatorname{Categorical}(p_{r,1},\cdots,p_{r,g_{r}})$ with
$\sum_{a=1}^{g_{r}}p_{r,a}=1$. Then, our formulation of fair PCA easily
extends via a one-versus-all comparison approach. That is to say, for each
$r\in[\ell]$, define

$\displaystyle{\bm{f}}_{r,a}=\mathbb{E}[{\bm{x}}\mid
a_{r}=a]-\mathbb{E}[{\bm{x}}\mid a_{r}\neq
a],\quad{\bm{S}}_{r,a}=\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}\mid
a_{r}=a]-\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}\mid a_{r}\neq a].$

Denote the top-$m_{r}$ column basis of eigenspace of ${\bm{S}}_{r,a}$ by
${\bm{P}}_{m_{r}}({\bm{S}}_{r,a})$. Then our notion of fair PCA can be
extended as follows: where $d,k,m_{1},\cdots,m_{n}$ are given as
$d>k+\ell+\sum_{r}m_{r}$,

$\displaystyle\max_{{\bm{V}}\in\operatorname{St}(d,k)}\mathrm{tr}({\bm{V}}^{\intercal}\bm{\Sigma}{\bm{V}}),\quad\text{subject
to
}\mathrm{col}({\bm{V}})\subset\bigcap_{r\in[\ell]}\bigcap_{a\in[g_{r}]}\mathrm{col}\left(\left[{\bm{P}}_{m_{r}}({\bm{S}}_{r,a})\mid{\bm{f}}_{r,a}\right]\right)^{\perp}.$
(7)

Our algorithm can also be naturally adapted to this scenario via one-versus-
all manner. In Appendix C, we provide the full pseudocode for the general case
of multiple and non-binary sensitive attributes.

##  6 FNPM is a PAFO-Learning Algorithm

We now show that our proposed memory-efficient algorithm, FNPM, is actually a
PAFO-learning algorithm in that with a certain sample complexity, it satisfies
the definition of PAFO-learnability. The proofs of all the theoretical results
stated here are deferred to Appendix E.

To use proper matrix concentration inequalities for our error term analysis,
various streaming PCA literature adopt some assumption on the underlying data
distribution, e.g., sub-Gaussianity (Jain et al., 2016; Yang et al., 2018;
Bienstock et al., 2022). Throughout the paper and for our theoretical
discussions, we consider a collection of distributions
${\mathcal{F}}_{d}\subset{\mathcal{P}}_{d}\times{\mathcal{P}}_{d}\times(0,1)$
satisfying the following assumptions to our data generation process in terms
of the data point conditioned on the sensitive attribute and its
eigenspectrum:

######  Assumption 1.

Consider our data generation process $a\sim\mathrm{Bernoulli}(p)$ and
${\bm{x}}\mid a\sim{\mathcal{D}}_{a}$. There exists
$\sigma,V,{\mathcal{M}},{\mathcal{V}},g_{\min}>0$,
$f_{\max}\in(g_{\min},\infty)$, and $p_{\min}\in(0,0.5]$ such that the
followings hold for any
$({\mathcal{D}}_{0},{\mathcal{D}}_{1},p)\in{\mathcal{F}}_{d}$: for
$a\in\\{0,1\\}$, ${\mathcal{D}}_{a}\in\mathrm{nSG}(\sigma)$,111${\bm{y}}$ is
norm-sub-Gaussian, denoted as $\mathrm{nSG}(\sigma)$, if
${\mathbb{P}}[\lVert{\bm{y}}-\mathbb{E}{\bm{y}}\rVert\geq t]\leq
2e^{-\frac{t^{2}}{2\sigma^{2}}}$. We refer interested readers to Jin et al.
(2019, Section 2) for more discussions on norm-sub-Gaussianity.

$\displaystyle{\mathbb{P}}\left[\left\|{\bm{x}}{\bm{x}}^{\intercal}-(\bm{\Sigma}_{a}+\bm{\mu}_{a}\bm{\mu}_{a}^{\intercal})\right\|_{2}\leq{\mathcal{M}}\
\middle|\ a\right]=1,\ \
\left\|\bm{\Sigma}_{a}+\bm{\mu}_{a}\bm{\mu}_{a}^{\intercal}\right\|_{2}\leq
V,\ \ \mathrm{Var}\left({\bm{x}}{\bm{x}}^{\intercal}\mid
a\right)\leq{\mathcal{V}},$

and

$\displaystyle\left\|{\bm{g}}\right\|_{2}=\left\|\bm{\Pi}_{{\bm{P}}_{m}}^{\perp}{\bm{f}}\right\|_{2}\in\\{0\\}\cup[g_{\min},f_{\max}],\
\ \left\|{\bm{f}}\right\|_{2}\leq f_{\max},\ \ \text{and}\ \
p\in[p_{\min},1-p_{\min}].$

######  Assumption 2.

Fix $m,k\in{\mathbb{N}}$. There exist
$\Delta_{m,\nu},\Delta_{k,\kappa},K_{m,\nu},K_{k,\kappa}\in(0,\infty)$ such
that for any
$({\mathcal{D}}_{0},{\mathcal{D}}_{1},\cdot)\in{\mathcal{F}}_{d}$, the
followings hold:

$\displaystyle\nu_{m}-\nu_{m+1}\geq\Delta_{m,\nu},\
\kappa_{k}-\kappa_{k+1}\geq\Delta_{k,\kappa},\ \nu_{m}\leq K_{m,\nu},\text{
and }\kappa_{k}\leq K_{k,\kappa},$

where $\nu_{1}\geq\cdots\geq\nu_{d}\geq 0$ and
$\kappa_{1}\geq\cdots\geq\kappa_{d}\geq 0$ are the singular values of
${\bm{S}}$ and
$\bm{\Pi}_{\bm{U}}^{\perp}\bm{\Sigma}\bm{\Pi}_{\bm{U}}^{\perp}$, respectively.

We start by establishing the sample complexity bounds of Algorithms 1 and 2
based on the convergence bound for NPM by Hardt and Price (2014). Recall that
NPM is an algorithm for finding top-$r$ eigenspace (in magnitude) of a
symmetric (but not necessarily PSD) matrix ${\bm{A}}$ under a random noise
${\bm{Z}}\in\mathbb{R}^{d\times k}$, by update ${\bm{V}}_{t}\leftarrow{\tt
QR}({\bm{A}}{\bm{V}}_{t-1}+{\bm{Z}}_{t})$. We start by recalling their meta-
sample complexity result for NPM, which we have slightly reformulated for our
convenience:

######  Lemma 1 (Corollary 1.1 of Hardt and Price (2014)).

Let $1\leq r<d$, $\epsilon\in(0,1/2)$ and $\delta\in(0,2e^{-cd})$, where $c$
is an absolute constant.222It depends polynomially only in the sub-Gaussian
moment of the data distribution ${\mathcal{D}}$; see Theorem 1.1 of Rudelson
and Vershynin (2009). Let ${\bm{L}}_{r}$ be the $d\times r$ matrix, whose
columns correspond to the top-$r$ eigenvectors (in magnitude) of the symmetric
(not necessarily PSD) matrix ${\bm{A}}\in\mathbb{R}^{d\times d}$, and let
$\xi_{1}\geq\cdots\geq\xi_{d}\geq 0$ be the singular values of ${\bm{A}}$.
Assume that the noise matrices ${\bm{Z}}_{t}\in\mathbb{R}^{d\times r}$ satisfy

$5\left\|{\bm{Z}}_{t}\right\|_{2}\leq\epsilon(\xi_{r}-\xi_{r+1})\quad\text{and}\quad
5\left\|{\bm{L}}_{r}^{\intercal}{\bm{Z}}_{t}\right\|_{2}\leq\frac{\delta(\xi_{r}-\xi_{r+1})}{2\sqrt{dr}},\quad\forall
t\geq 1.$ (8)

Then, after
$T={\Theta}\left(\frac{\xi_{r}}{\xi_{r}-\xi_{r+1}}\log\left(\frac{d}{\epsilon\delta}\right)\right)$
steps of NPM,
$\left\|\bm{\Pi}_{{\bm{V}}_{T}}^{\perp}{\bm{L}}_{r}\right\|_{2}\leq\epsilon$
with probability at least $1-\delta$.

First, we prove that the ${\bm{W}}_{T}$ resulting from Algorithm 1 (NPM for
the second moment gap) converges to the true value. We can view Algorithm 1’s
updates as ${\bm{W}}_{t}\leftarrow{\tt
QR}({\bm{S}}{\bm{W}}_{t-1}+{\bm{Z}}_{t,1})$, where the noise matrix in this
case is
${\bm{Z}}_{t,1}:=({\bm{C}}_{t}^{(1)}-{\bm{C}}_{t}^{(0)})-{\bm{S}}{\bm{W}}_{t-1}$
and ${\bm{C}}_{t}^{(a)}$ is as defined in Eqn. (6). The following lemma
asserts that with large enough block size $b$, the error matrices are
sufficiently bounded such that the NPM iterates converge:

######  Lemma 2.

Let $\epsilon,\delta\in(0,1)$. It is sufficient to choose the block size $b$
in Algorithm 1 as

$b=\Omega\left(\frac{{\mathcal{V}}}{\Delta_{m,\nu}^{2}p_{\min}}\left(\frac{dm}{\delta^{2}}\log\frac{m}{\delta}+\frac{1}{\epsilon^{2}}\log\frac{d}{\delta}\right)+\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}p_{\min}}\log\frac{d}{\delta}\right)$
(9)

to make the following hold with probability at least $1-\delta$:

$5\left\|{\bm{Z}}_{t,1}\right\|_{2}\leq\epsilon\Delta_{m,\nu}\quad\text{and}\quad
5\left\|{\bm{P}}_{m}^{\intercal}{\bm{Z}}_{t,1}\right\|_{2}\leq\frac{\delta\Delta_{m,\nu}}{2\sqrt{dm}},\quad\forall
t\geq 1,$ (10)

where we recall that the columns of ${\bm{P}}_{m}$ are the top-$m$ (in
magnitude) orthonormal eigenvectors of ${\bm{S}}$.

Let $\widehat{{\bm{U}}}$ be the final estimate of the true ${\bm{U}}$
outputted by Algorithm 1. For Algorithm 2, the noise matrix is
${\bm{Z}}_{\tau,2}:=\left(\bm{\Pi}_{\widehat{{\bm{U}}}}^{\perp}\widehat{\bm{\Sigma}}_{\tau}\bm{\Pi}_{\widehat{{\bm{U}}}}^{\perp}-\bm{\Pi}_{\bm{U}}^{\perp}\bm{\Sigma}\bm{\Pi}_{\bm{U}}^{\perp}\right){\bm{V}}_{\tau-1}$
, where
$\widehat{\bm{\Sigma}}_{\tau}:=\frac{1}{{\mathcal{B}}}\sum_{j}{\bm{x}}_{j}{\bm{x}}_{j}^{\intercal}$
is the sample covariance at time step $\tau$ of Algorithm 2. Similarly, with
large enough block size ${\mathcal{B}}$, we have the following lemma:

######  Lemma 3.

Let $\epsilon,\delta\in(0,1)$. Suppose that
$\left\|\bm{\Pi}_{\widehat{{\bm{U}}}}-\bm{\Pi}_{{\bm{U}}}\right\|_{2}\leq\frac{\Delta_{k,\kappa}}{20V}\min\left(\epsilon,\frac{\delta}{2\sqrt{dk}}\right)$.
Then, it is sufficient to choose the block size ${\mathcal{B}}$ in Algorithm 2
as

${\mathcal{B}}=\Omega\left(\frac{{\mathcal{V}}+V^{2}}{\Delta_{k,\kappa}^{2}}\left(\frac{dk}{\delta^{2}}\log\frac{k}{\delta}+\frac{1}{\epsilon^{2}}\log\frac{d}{\delta}\right)+\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}+V^{2}}\log\frac{d}{\delta}\right),$
(11)

to make the following hold with probability at least $1-\delta$:

$5\left\|{\bm{Z}}_{\tau,2}\right\|_{2}\leq\epsilon\Delta_{k,\kappa}\quad\text{and}\quad
5\left\|{\bm{Q}}_{k}^{\intercal}{\bm{Z}}_{\tau,2}\right\|_{2}\leq\frac{\delta\Delta_{k,\kappa}}{2\sqrt{dk}},\quad\forall\tau\geq
1,$ (12)

where the columns of ${\bm{Q}}_{k}$ is the top-$k$ (in magnitude) orthonormal
eigenvectors of
$\bm{\Pi}_{\bm{U}}^{\perp}\bm{\Sigma}\bm{\Pi}_{\bm{U}}^{\perp}$. (We note that
${\bm{Q}}_{k}$ is a solution for Eqn. (4).)

######  Remark 1.

We emphasize that, regardless of the block size requirement in the above
lemmas, all the operations can be implemented within $o(d^{2})$ space
requirement.

Combining the convergence results, we can prove the following PAFO-
learnability guarantee in the memory-limited, streaming setting:

######  Theorem 1.

Let $d,m,k\in{\mathbb{N}}$ be fixed. Consider a collection
${\mathcal{F}}_{d}\subset{\mathcal{P}}_{d}\times{\mathcal{P}}_{d}\times(0,1)$
satisfying Assumptions 1 and 2. Then, ${\mathcal{F}}_{d}$ is PAFO-learnable
for streaming PCA with FNPM, where the sufficient number of samples is given
as $N_{{\mathcal{F}}_{d}}(\varepsilon_{\rm o},\varepsilon_{\rm
f},\delta)=N_{1}+N_{2}$, with

$\displaystyle N_{1}$
$\displaystyle=\tilde{\Omega}\left(\frac{K_{m,\nu}}{p_{\min}}\left\\{\frac{{\mathcal{V}}}{\Delta_{m,\nu}^{3}}\left(\\!\frac{dm}{\delta^{2}}\\!+\\!\frac{\alpha}{\eta_{k}^{2}}\right)+\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}\Delta_{m,\nu}}\right\\}+\frac{\mathbbm{1}[{\bm{g}}\neq\bm{0}]\sigma^{2}}{p_{\min}g_{\min}^{2}\eta_{k}^{2}}\right),$
(Algorithm 1) $\displaystyle N_{2}$
$\displaystyle=\tilde{\Omega}\left(\frac{K_{k,\kappa}({\mathcal{V}}+V^{2})}{\Delta_{k,\kappa}^{3}}\left(\frac{dk}{\delta^{2}}+\frac{k^{2}V^{2}}{\varepsilon_{\rm
o}^{2}}\right)+\frac{K_{k,\kappa}{\mathcal{M}}^{2}}{\Delta_{k,\kappa}({\mathcal{V}}+V^{2})}\right),$
(Algorithm 2)

where $\eta_{k}=\Theta\left(\min\left\\{\varepsilon_{\rm
f},\frac{\Delta_{k,\kappa}}{kV^{2}}\varepsilon_{\rm
o},\frac{\Delta_{k,\kappa}}{V\sqrt{dk}}\delta\right\\}\right)$ and
$\alpha=1+\frac{\mathbbm{1}[{\bm{g}}\neq\bm{0}]f_{\max}^{2}}{g_{\min}^{2}}$.

Let us take a moment to digest the sample complexity in Theorem 1. The second
term $N_{2}$ is determined by the $\varepsilon_{\rm o}$-optimality requirement
in our PAFO-learnability. Note that $N_{2}$ has no dependencies on fairness-
related quantities such as $\varepsilon_{\rm f}$, $p_{\min}$, and
$\Delta_{m,\nu}$. On the other hand, the first term $N_{1}$ is not only
determined by the $\varepsilon_{\rm f}$-fairness requirement, but also the
$\varepsilon_{\rm o}$-optimality as well. This is the “price” of pursuing
fairness; ${\bm{U}}$, which encodes the unfair subspace needed to be
nullified, is required to be accurately estimated as it impacts not only the
level of fairness but also the resulting solution’s optimality. This is clear
in our formulation of fair PCA; in Eqn. (4), note how the optimal solution
depends heavily on $\bm{\Pi}_{\bm{U}}$.

We further elaborate on the dependencies of $N_{1}$ on fairness-related
quantities, namely $p_{\min}$ and ${\bm{g}}$. First, if $p_{\min}\rightarrow
0$, i.e., if one of the two groups is never sampled, then the sample
complexity tends to infinity, and the learnability does not hold; this aligns
with our intuition, as we need samples from both of the sensitive groups. Its
dependency is also quite natural, as the minimum expected number of samples
from either group depends linearly on $\frac{1}{p_{\min}}$. Also, when
${\bm{g}}\neq\bm{0}$, $N_{1}$ has an additional additive term scaling linearly
with $\frac{1}{g_{\min}^{2}}$. This is because when ${\bm{g}}\neq\bm{0}$, we
must explicitly account for the approximation error of
$\frac{{\bm{g}}}{\left\|{\bm{g}}\right\|_{2}}$ due to the QR decomposition at
the end of Algorithm 1.

##  7 Experiments on real-world datasets

The code for all experiments is available at github.com/HanseulJo/fair-
streaming-pca.

###  7.1 CelebA dataset

Figure 1: Experimental results on full-resolution CelebA dataset. Original
image v.s. FNPM output. (Sensitive attribute: “Eyeglasses")

We evaluate the efficacy of our proposed FNPM on the CelebA dataset (Liu et
al., 2015b). It has been considered by Kleindessner et al. (2023) to show the
superior efficiency of their fair PCA algorithm compared to previous
approaches (Ravfogel et al., 2022a; Lee et al., 2022; Olfat and Aswani, 2019).
However, even in Kleindessner et al. (2023), the images were resized and grey-
scaled, reducing the dimension from the original 218$\times$178$=$38,804 to
80$\times$80$=$6,400. Indeed, for a modest-sized computer, it is impossible to
load all 162,770 original images in training set to the memory at once, while
they require a full dataset to run each step of their algorithm. Here, we use
the original resolution, full-color CelebA dataset. We implement our FNPM
using Python JAX NumPy Module (Harris et al., 2020; Bradbury et al., 2023) and
Pytorch (Paszke et al., 2017). All experiments were performed on Apple 2020
Mac mini M1 with 16GB RAM.

Although CelebA is not streaming in nature, we intend to show that
transforming it to one and using our memory-efficient approach allow us to
scale up fair PCA. Since there are three channels of color, we run FNPM
channel-wise but in _parallel_ as usual in vision tasks (Priorov et al.,
2013). For each channel of colors, we project the data onto a
$k=1000$-dimensional subspace while nullifying $m=2$ leading eigenvectors of
covariance difference.

The resulting images are displayed in Figure 1. Here, we consider ‘Eyeglasses’
as a sensitive attribute to divide groups. We adopt the predefined train-
validation split and run our algorithm only on the training set for 5
iterations with block sizes of $b={\mathcal{B}}=32,000$. Then, using the
output ${\bm{V}}$ of FNPM, we project images selected from the validation set.
We observe that we have images of faces wearing colorful glasses by nullifying
some of the leading eigenvectors of covariance difference. For the images
originally with sunglasses, their glasses get blurred, and “virtual” eyes are
added to them. We provide more results on the other attributes and ablation
studies on varying $m$ and $b$ in Appendix F.2.

###  7.2 UCI datasets

Table 1: Dataset = Adult Income [feature dim=102, #(train data)=31,655], $k=2$ Method | %Var($\uparrow$) | MMD2($\downarrow$) |  %Acc($\uparrow$) |  $\Delta_{\rm DP}$($\downarrow$)
---|---|---|---|---
kernel SVM
PCA |  6.88(0.14) |  0.374(0.006) |  82.4(0.2) |  0.19(0.01)
Olfat and Aswani (2019) (0.1) | Memory Out
Olfat and Aswani (2019) (0.0) | Memory Out
Lee et al. (2022) (1e-3) |  5.68(0.11) |  0.0(0.0) |  80.34(0.24) |  0.05(0.01)
Lee et al. (2022) (1e-6) |  5.42(0.11) |  0.0(0.0) |  79.41(0.23) |  0.02(0.01)
Kleindessner et al. (2023) (mean) |  5.74(0.11) |  0.002(0.0) |  80.6(0.2) |  0.07(0.01)
Kleindessner et al. (2023) (0.85) |  4.09(0.17) |  0.001(0.001) |  75.52(0.21) |  0.0(0.0)
Kleindessner et al. (2023) (0.5) |  2.63(0.07) |  0.0(0.0) |  75.38(0.18) |  0.0(0.0)
Kleindessner et al. (2023) (kernel) | Takes too long time
Ravfogel et al. (2020) |  1.91(0.08) |  0.001(0.001) |  75.67(0.31) |  0.0(0.0)
Ravfogel et al. (2022a) |  1.91(0.09) |  0.006(0.011) |  75.59(0.34) |  0.0(0.0)
Samadi et al. (2018) | N/A | N/A |  82.63(0.18) |  0.15(0.01)
Ours (offline, mean) |  5.74(0.11) |  0.002(0.0) |  80.6(0.2) |  0.07(0.01)
Ours (FNPM, mean) |  5.74(0.11) |  0.002(0.0) |  80.6(0.2) |  0.07(0.01)
Ours (offline, $m$=15) |  4.04(0.14) |  0.001(0.001) |  75.51(0.23) |  0.0(0.0)
Ours (FNPM, $m$=15) |  4.07(0.13) |  0.001(0.0) |  75.54(0.25) |  0.0(0.0)
Ours (offline, $m$=50) |  2.63(0.07) |  0.0(0.0) |  75.38(0.18) |  0.0(0.0)
Ours (FNPM, $m$=50) |  2.64(0.06) |  0.0(0.0) |  75.44(0.21) |  0.0(0.0)

For the sake of completeness, we conduct a quantitative evaluation of our
algorithm on UCI datasets (Adult Income, COMPAS, German Credit) and compare it
with six previous works (Olfat and Aswani, 2019; Lee et al., 2022;
Kleindessner et al., 2023; Ravfogel et al., 2020, 2022a; Samadi et al., 2018).
The results for the Adult Income dataset are shown in Table 1. We assess
several variants of our methods: in Table 1, “mean” is when we match the means
and not second moments, “FNPM” is when we run our Algorithm 2 with a block
size of full-batch, and “offline” is when we directly solve our “Null It Out”
formulation of fair PCA (Eqn. (4)) via offline eigen-decomposition. In the
table, we report the explained variance (%Var), representation fairness
measured by maximum mean discrepancy (MMD2), downstream task accuracy (%Acc),
and downstream task fairness in demographic parity (%$\Delta_{\textrm{DP}}$).
The result showcases that our method yields competitive quantitative
performance even for the common tabular datasets while being much more memory
efficient. We defer the full results for other UCI datasets to Appendix F.3.

##  8 Other Related Works

#### Fairness in ML.

There are largely two directions in the literature of algorithmic fairness.
One direction is to propose a suitable and meaningful fairness definition
(Feldman et al., 2015; Hardt et al., 2016; Dwork et al., 2012). The other
direction is to develop efficient fair algorithms, although often the fairness
constraint forces the algorithm to be much more inefficient than its unfair
counterpart, or it calls for a need for a completely new algorithmic approach.
There are also different ways of imposing fairness in an ML pipeline, such as
learning fair pre-processing (Biswas and Rajan, 2021), fair in-processing
(Zafar et al., 2019; Roh et al., 2021; Wan et al., 2023), and more. The reader
is encouraged to check Barocas et al. (2019) for a more comprehensive
treatment of this subject.

#### Fair online/streaming Learning.

Gillen et al. (2018); Bechavod et al. (2020) study individual fairness in
online learning in a learning theoretic framework, even when the underlying
metric is unavailable. Stemming from the concept of fair clustering as
proposed in Chierichetti et al. (2017), Bera et al. (2022); Schmidt et al.
(2020) study imposing demographic parity on clustering in the streaming
setting. Such fairness has been considered in various other streaming problems
such as online selection (Correa et al., 2021), streaming submodular
optimization (El Halabi et al., 2020), and diversity maximization (Wang et
al., 2022). Quite surprisingly, demographic parity (or any other concept of
fairness) has never been considered in the setting of streaming PCA.

#### Streaming PCA.

Without the fairness constraint, streaming PCA has been studied much from
statistics and the machine learning community. Two prominent algorithms have
been studied; the noisy power method (Mitliagkas et al., 2013) and Oja’s
method (Oja, 1982). Much work has been done in improving the theoretical
guarantees of streaming PCA (Hardt and Price, 2014; Balcan et al., 2016; Jain
et al., 2016; Liang, 2023), improving the algorithm itself (Xu, 2023; Yun,
2018), or extending the guarantees to various different settings (Balzano et
al., 2018; Bienstock et al., 2022; Kumar and Sarkar, 2023). Memory-limited,
streaming versions of somewhat related problems, such as community detection
(Yun et al., 2014) and low-rank matrix completion (Yun et al., 2015), have
been tackled as well using similar spectral techniques as PCA (e.g., power
method). However, to the best of our knowledge, fairness (regardless of the
definition) has never been considered in this context of streaming PCA, which
we tackle in this work and which we believe is of great importance.

##  9 Conclusion

In this work, we tackled the two outstanding problems of the existing fair PCA
literature. From the theoretical side, we illustrated a new formulation of
fair PCA based on the “Null It Out” approach and then provided a novel
statistical framework called PAFO-learnability of fair PCA. From the practical
side, we addressed the memory-limited scenarios by proposing a new problem
setting called fair streaming PCA and a memory-efficient two-stage algorithm
called FNPM. Based on these, we established a statistical guarantee that our
algorithm achieves the PAFO-learnability for fair streaming PCA. Lastly, we
ran experiments on the CelebA and UCI datasets to certify the scalability of
our method.

## Acknowledgments and Disclosure of Funding

We thank the anonymous reviewers for their helpful comments and suggestions.
We also thank Gwangsu Kim (Jeonbuk National University) for the helpful
discussions in the initial phase of the research. This work was supported by
the Institute of Information & Communications Technology Planning & Evaluation
(IITP) grants funded by the Korean government(MSIT) (No.2019-0-00075,
Artificial Intelligence Graduate School Program(KAIST); No. 2022-0-00184,
Development and Study of AI Technologies to Inexpensively Conform to Evolving
Policy on Ethics).

## References

  * Abdi and Williams (2010) Hervé Abdi and Lynne J. Williams.  Principal component analysis.  _WIREs Computational Statistics_ , 2(4):433–459, 2010.
  * Ananthakrishnan et al. (2021) Nivasini Ananthakrishnan, Shai Ben-David, Tosca Lechner, and Ruth Urner.  Identifying regions of trusted predictions.  In _Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence_ , volume 161 of _Proceedings of Machine Learning Research_ , pages 2125–2134. PMLR, 27–30 Jul 2021.
  * Balcan et al. (2016) Maria-Florina Balcan, Simon Shaolei Du, Yining Wang, and Adams Wei Yu.  An Improved Gap-Dependency Analysis of the Noisy Power Method.  In _29th Annual Conference on Learning Theory_ , volume 49 of _Proceedings of Machine Learning Research_ , pages 284–309, Columbia University, New York, New York, USA, 23–26 Jun 2016. PMLR.
  * Balzano et al. (2018) Laura Balzano, Yuejie Chi, and Yue M. Lu.  Streaming PCA and Subspace Tracking: The Missing Data Case.  _Proceedings of the IEEE_ , 106(8):1293–1310, 2018.
  * Barocas and Selbst (2016) Solon Barocas and Andrew D. Selbst.  Big Data’s Disparate Impact.  _104 California Law Review 671_ , 2016.
  * Barocas et al. (2019) Solon Barocas, Moritz Hardt, and Arvind Narayanan.  _Fairness and Machine Learning: Limitations and Opportunities_.  fairmlbook.org, 2019.
  * Bechavod et al. (2020) Yahav Bechavod, Christopher Jung, and Steven Z. Wu.  Metric-Free Individual Fairness in Online Learning.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 11214–11225. Curran Associates, Inc., 2020.
  * Bera et al. (2022) Suman K. Bera, Syamantak Das, Sainyam Galhotra, and Sagar Sudhir Kale.  Fair K-Center Clustering in MapReduce and Streaming Settings.  In _Proceedings of the ACM Web Conference 2022_ , WWW ’22, page 1414–1422, New York, NY, USA, 2022. Association for Computing Machinery.
  * Bienstock et al. (2022) Daniel Bienstock, Minchan Jeong, Apurv Shukla, and Se-Young Yun.  Robust Streaming PCA.  In _Advances in Neural Information Processing Systems_ , volume 35, pages 4231–4243. Curran Associates, Inc., 2022.
  * Biswas and Rajan (2021) Sumon Biswas and Hridesh Rajan.  Fair Preprocessing: Towards Understanding Compositional Fairness of Data Transformers in Machine Learning Pipeline.  In _Proceedings of the 29th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering_ , ESEC/FSE 2021, page 981–993, New York, NY, USA, 2021. Association for Computing Machinery.
  * Bradbury et al. (2023) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2023\.  URL http://github.com/google/jax.
  * Chamon and Ribeiro (2020) Luiz Chamon and Alejandro Ribeiro.  Probably Approximately Correct Constrained Learning.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 16722–16735. Curran Associates, Inc., 2020.
  * Chierichetti et al. (2017) Flavio Chierichetti, Ravi Kumar, Silvio Lattanzi, and Sergei Vassilvitskii.  Fair Clustering Through Fairlets.  In _Advances in Neural Information Processing Systems_ , volume 30, pages 5036–5044. Curran Associates, Inc., 2017.
  * Correa et al. (2021) Jose Correa, Andres Cristi, Paul Duetting, and Ashkan Norouzi-Fard.  Fairness and Bias in Online Selection.  In _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pages 2112–2121. PMLR, 18–24 Jul 2021.
  * Davis and Kahan (1969) Chandler Davis and W. M. Kahan.  Some new bounds on perturbation of subspaces.  _Bulletin of the American Mathematical Society_ , 75(4):863 – 868, 1969.
  * Dwork et al. (2012) Cynthia Dwork, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Richard Zemel.  Fairness through Awareness.  In _Proceedings of the 3rd Innovations in Theoretical Computer Science Conference_ , ITCS ’12, page 214–226, New York, NY, USA, 2012. Association for Computing Machinery.
  * Eckart and Young (1936) Carl Eckart and Gale Young.  The approximation of one matrix by another of lower rank.  _Psychometrika_ , 1:211—218, 1936.
  * El Halabi et al. (2020) Marwa El Halabi, Slobodan Mitrović, Ashkan Norouzi-Fard, Jakab Tardos, and Jakub M Tarnawski.  Fairness in Streaming Submodular Maximization: Algorithms and Hardness.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 13609–13622. Curran Associates, Inc., 2020.
  * Feldman et al. (2015) Michael Feldman, Sorelle A. Friedler, John Moeller, Carlos Scheidegger, and Suresh Venkatasubramanian.  Certifying and Removing Disparate Impact.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’15, page 259–268, New York, NY, USA, 2015. Association for Computing Machinery.
  * Ghashami et al. (2016) Mina Ghashami, Daniel J. Perry, and Jeff Phillips.  Streaming Kernel Principal Component Analysis.  In _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ , volume 51 of _Proceedings of Machine Learning Research_ , pages 1365–1374, Cadiz, Spain, 09–11 May 2016. PMLR.
  * Gillen et al. (2018) Stephen Gillen, Christopher Jung, Michael Kearns, and Aaron Roth.  Online Learning with an Unknown Fairness Metric.  In _Advances in Neural Information Processing Systems_ , volume 31, pages 2605–2614. Curran Associates, Inc., 2018.
  * Golub and Loan (2013) Gene H. Golub and Charles F. Van Loan.  _Matrix Computations_.  Johns Hopkins Studies in the Mathematical Sciences. The Johns Hopkins University Press, 4 edition, 2013.
  * Hardt and Price (2014) Moritz Hardt and Eric Price.  The Noisy Power Method: A Meta Algorithm with Applications.  In _Advances in Neural Information Processing Systems_ , volume 27, pages 2861–2869. Curran Associates, Inc., 2014.
  * Hardt et al. (2016) Moritz Hardt, Eric Price, and Nati Srebro.  Equality of Opportunity in Supervised Learning.  In _Advances in Neural Information Processing Systems_ , volume 29, pages 3323–3331. Curran Associates, Inc., 2016.
  * Harris et al. (2020) Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.
  * Hopkins et al. (2023) Max Hopkins, Daniel M. Kane, Shachar Lovett, and Gaurav Mahajan.  Do PAC-Learners Learn the Marginal Distribution?  _arXiv preprint arXiv:2302.06285_ , 2023.
  * Hotelling (1933) Harold Hotelling.  Analysis of a complex of statistical variables into principal components.  _Journal of Educational Psychology_ , 24(6):417–441, 1933.
  * Huang et al. (2021) De Huang, Jonathan Niles-Weed, and Rachel Ward.  Streaming k-PCA: Efficient guarantees for Oja’s algorithm, beyond rank-one updates.  In _Proceedings of Thirty Fourth Conference on Learning Theory_ , volume 134 of _Proceedings of Machine Learning Research_ , pages 2463–2498. PMLR, 15–19 Aug 2021.
  * Jain et al. (2016) Prateek Jain, Chi Jin, Sham M. Kakade, Praneeth Netrapalli, and Aaron Sidford.  Streaming PCA: Matching Matrix Bernstein and Near-Optimal Finite Sample Guarantees for Oja’s Algorithm.  In _29th Annual Conference on Learning Theory_ , volume 49 of _Proceedings of Machine Learning Research_ , pages 1147–1164, Columbia University, New York, New York, USA, 23–26 Jun 2016. PMLR.
  * Jin et al. (2019) Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M. Kakade, and Michael I. Jordan.  A Short Note on Concentration Inequalities for Random Vectors with SubGaussian Norm.  _arXiv preprint arXiv:1902.03736_ , 2019.
  * Johnson and Wichern (2008) Richard A. Johnson and Dean W. Wichern.  _Applied Multivariate Statistical Analysis_.  Pearson, 6 edition, 2008.
  * Jolliffe and Cadima (2016) Ian T. Jolliffe and Jorge Cadima.  Principal component analysis: a review and recent developments.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 374(2065), 2016.
  * Kamani et al. (2022) Mohammad Mahdi Kamani, Farzin Haddadpour, Rana Forsati, and Mehrdad Mahdavi.  Efficient fair principal component analysis.  _Machine Learning_ , 111(10):3671–3702, 2022\.
  * Kirchner et al. (2016) L Kirchner, J. Larson, S. Mattu, and J. Angwin.  Machine bias.  ProPublica, 2016.
  * Kizilcec and Lee (2021) René F. Kizilcec and Hansol Lee.  Algorithmic Fairness in Education.  _arXiv preprint arXiv:2007.05443_ , 2021.
  * Kleindessner et al. (2019) Matthäus Kleindessner, Samira Samadi, Pranjal Awasthi, and Jamie Morgenstern.  Guarantees for Spectral Clustering with Fairness Constraints.  In _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 3458–3467. PMLR, 09–15 Jun 2019.
  * Kleindessner et al. (2023) Matthäus Kleindessner, Michele Donini, Chris Russell, and Muhammad Bilal Zafar.  Efficient fair PCA for fair representation learning.  In _Proceedings of The 26th International Conference on Artificial Intelligence and Statistics_ , volume 206 of _Proceedings of Machine Learning Research_ , pages 5250–5270. PMLR, 25–27 Apr 2023.
  * Kohler and Lucchi (2017) Jonas Moritz Kohler and Aurelien Lucchi.  Sub-sampled Cubic Regularization for Non-convex Optimization.  In _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 1895–1904. PMLR, 06–11 Aug 2017.
  * Kumar and Sarkar (2023) Syamantak Kumar and Purnamrita Sarkar.  Streaming PCA for Markovian Data.  _arXiv preprint arXiv:2305.02456_ , 2023.
  * Lahoti et al. (2020) Preethi Lahoti, Alex Beutel, Jilin Chen, Kang Lee, Flavien Prost, Nithum Thain, Xuezhi Wang, and Ed Chi.  Fairness without Demographics through Adversarially Reweighted Learning.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 728–740. Curran Associates, Inc., 2020.
  * Lee et al. (2022) Junghyun Lee, Gwangsu Kim, Mahbod Olfat, Mark Hasegawa-Johnson, and Chang D. Yoo.  Fast and Efficient MMD-Based Fair PCA via Optimization over Stiefel Manifold.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, pages 7363–7371, Jun. 2022.
  * Li et al. (2023) Xun Li, Dongsheng Chen, Weipan Xu, Haohui Chen, Junjun Li, and Fan Mo.  Explainable dimensionality reduction (XDR) to unbox AI ‘black box’ models: A study of AI perspectives on the ethnic styles of village dwellings.  _Humanities and Social Sciences Communications_ , 10(1):35, Jan 2023.
  * Liang (2023) Xin Liang.  On the optimality of the Oja’s algorithm for online PCA.  _Statistics and Computing_ , 33(3):62, Mar 2023\.
  * Liberty (2013) Edo Liberty.  Simple and Deterministic Matrix Sketching.  In _Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’13, page 581–588, New York, NY, USA, 2013. Association for Computing Machinery.
  * Liu et al. (2015a) Ziqi Liu, Yu-Xiang Wang, and Alexander Smola.  Fast Differentially Private Matrix Factorization.  In _Proceedings of the 9th ACM Conference on Recommender Systems_ , RecSys ’15, page 171–178, New York, NY, USA, 2015a. Association for Computing Machinery.
  * Liu et al. (2015b) Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.  Deep Learning Face Attributes in the Wild.  In _2015 IEEE International Conference on Computer Vision (ICCV)_ , pages 3730–3738, 2015b.
  * Marshall et al. (2011) Albert W. Marshall, Ingram Olkin, and Barry C. Arnold.  _Inequalities: Theory of Majorization and Its Applications_.  Springer Series in Statistics. Springer New York, 2 edition, 2011.
  * Mehrabi et al. (2021) Ninareh Mehrabi, Fred Morstatter, Nripsuta Saxena, Kristina Lerman, and Aram Galstyan.  A Survey on Bias and Fairness in Machine Learning.  _ACM Computing Surveys_ , 54(6), jul 2021.
  * Mirsky (1975) L. Mirsky.  A trace inequality of John von Neumann.  _Monatshefte für Mathematik_ , 79(4):303–306, Dec 1975.
  * Mitliagkas et al. (2013) Ioannis Mitliagkas, Constantine Caramanis, and Prateek Jain.  Memory Limited, Streaming PCA.  In _Advances in Neural Information Processing Systems_ , volume 26, pages 2886–2894. Curran Associates, Inc., 2013.
  * Oja (1982) Erkki Oja.  Simplified neuron model as a principal component analyzer.  _Journal of Mathematical Biology_ , 15(3):267–273, Nov 1982.
  * Oja and Karhunen (1985) Erkki Oja and Juha Karhunen.  On Stochastic Approximation of the Eigenvectors and Eigenvalues of the Expectation of a Random Matrix.  _Journal of Mathematical Analysis and Applications_ , 106:69–84, 1985.
  * Olfat and Aswani (2019) Matt Olfat and Anil Aswani.  Convex Formulations for Fair Principal Component Analysis.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 663–670, Jul. 2019.
  * Paszke et al. (2017) Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in PyTorch.  In _NIPS 2017 Workshop Autodiff_ , 2017.
  * Pearson (1901) Karl Pearson.  Liii. on lines and planes of closest fit to systems of points in space.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 2(11):559–572, 1901.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher Manning.  GloVe: Global Vectors for Word Representation.  In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1532–1543, Doha, Qatar, October 2014. Association for Computational Linguistics.
  * Priorov et al. (2013) Andrey Priorov, Kirill Tumanov, Vladimir Volokhov, Evgeny Sergeev, and Ivan Mochalov.  Applications of image filtration based on principal component analysis and nonlocal image processing.  _IAENG International Journal of Computer Science_ , 40(2):62–80, 2013.
  * Ravfogel et al. (2020) Shauli Ravfogel, Yanai Elazar, Hila Gonen, Michael Twiton, and Yoav Goldberg.  Null It Out: Guarding Protected Attributes by Iterative Nullspace Projection.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 7237–7256, Online, July 2020. Association for Computational Linguistics.
  * Ravfogel et al. (2022a) Shauli Ravfogel, Michael Twiton, Yoav Goldberg, and Ryan D Cotterell.  Linear Adversarial Concept Erasure.  In _Proceedings of the 39th International Conference on Machine Learning_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 18400–18421. PMLR, 17–23 Jul 2022a.
  * Ravfogel et al. (2022b) Shauli Ravfogel, Francisco Vargas, Yoav Goldberg, and Ryan Cotterell.  Adversarial Concept Erasure in Kernel Space.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , pages 6034–6055, Abu Dhabi, United Arab Emirates, December 2022b. Association for Computational Linguistics.
  * Roh et al. (2021) Yuji Roh, Kangwook Lee, Steven Euijong Whang, and Changho Suh.  FairBatch: Batch Selection for Model Fairness.  In _International Conference on Learning Representations_ , 2021.
  * Rudelson and Vershynin (2009) Mark Rudelson and Roman Vershynin.  Smallest singular value of a random rectangular matrix.  _Communications on Pure and Applied Mathematics_ , 62(12):1707–1739, 2009.
  * Samadi et al. (2018) Samira Samadi, Uthaipon Tantipongpipat, Jamie H Morgenstern, Mohit Singh, and Santosh Vempala.  The Price of Fair PCA: One Extra dimension.  In _Advances in Neural Information Processing Systems_ , volume 31, pages 10999–11010. Curran Associates, Inc., 2018.
  * Schmidt et al. (2020) Melanie Schmidt, Chris Schwiegelshohn, and Christian Sohler.  Fair coresets and streaming algorithms for fair k-means.  In _Approximation and Online Algorithms_ , pages 232–251, Cham, 2020\. Springer International Publishing.
  * Schölkopf et al. (1998) Bernhard Schölkopf, Alexander Smola, and Klaus-Robert Müller.  Nonlinear Component Analysis as a Kernel Eigenvalue Problem.  _Neural Computation_ , 10(5):1299–1319, 07 1998\.
  * Shalev-Schwartz and Ben-David (2014) S. Shalev-Schwartz and S. Ben-David.  _Understanding Machine Learning: From Theory to Algorithms_.  Cambridge University Press, 2014.
  * Tantipongpipat et al. (2019) Uthaipon Tantipongpipat, Samira Samadi, Mohit Singh, Jamie H Morgenstern, and Santosh Vempala.  Multi-Criteria Dimensionality Reduction with Applications to Fairness.  In _Advances in Neural Information Processing Systems_ , volume 32, pages 15161–15171. Curran Associates, Inc., 2019.
  * Tjoa and Guan (2021) Erico Tjoa and Cuntai Guan.  A Survey on Explainable Artificial Intelligence (XAI): Toward Medical XAI.  _IEEE Transactions on Neural Networks and Learning Systems_ , 32(11):4793–4813, 2021.
  * Tropp (2015) Joel A. Tropp.  An Introduction to Matrix Concentration Inequalities.  _Foundations and Trends® in Machine Learning_ , 8(1-2):1–230, 2015.
  * Ullah et al. (2018) Enayat Ullah, Poorya Mianjy, Teodor Vanislavov Marinov, and Raman Arora.  Streaming Kernel PCA with $\tilde{{\mathcal{O}}}(\sqrt{n})$ Random Features.  In _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.
  * von Neumann (1937) John von Neumann.  Some Matrix-Inequalities and Metrization of Matrix-Space.  _Tomsk University Review_ , 1:286–300, 1937.
  * Vu et al. (2022) Hieu Vu, Toan Tran, Man-Chung Yue, and Viet Anh Nguyen.  Distributionally Robust Fair Principal Components via Geodesic Descents.  In _International Conference on Learning Representations_ , 2022.
  * Wan et al. (2023) Mingyang Wan, Daochen Zha, Ninghao Liu, and Na Zou.  In-Processing Modeling Techniques for Machine Learning Fairness: A Survey.  _ACM Transactions on Knowledge Discovery from Data_ , 17(3), mar 2023.
  * Wang and Lu (2016) Chuang Wang and Yue M. Lu.  Online learning for sparse PCA in high dimensions: Exact dynamics and phase transitions.  In _2016 IEEE Information Theory Workshop (ITW)_ , pages 186–190, 2016.
  * Wang et al. (2023) Ji Wang, Ding Lu, Ian Davidson, and Zhaojun Bai.  Scalable Spectral Clustering with Group Fairness Constraints.  In _Proceedings of The 26th International Conference on Artificial Intelligence and Statistics_ , volume 206 of _Proceedings of Machine Learning Research_ , pages 6613–6629. PMLR, 25–27 Apr 2023.
  * Wang et al. (2022) Yanhao Wang, Francesco Fabbri, and Michael Mathioudakis.  Streaming Algorithms for Diversity Maximization with Fairness Constraints.  In _2022 IEEE 38th International Conference on Data Engineering (ICDE)_ , pages 41–53, Los Alamitos, CA, USA, 2022. IEEE Computer Society.
  * Xu (2023) Zhiqiang Xu.  On the Accelerated Noise-Tolerant Power Method.  In _Proceedings of The 26th International Conference on Artificial Intelligence and Statistics_ , volume 206 of _Proceedings of Machine Learning Research_ , pages 7147–7175. PMLR, 25–27 Apr 2023.
  * Xu and Li (2022) Zhiqiang Xu and Ping Li.  Faster Noisy Power Method.  In _Proceedings of The 33rd International Conference on Algorithmic Learning Theory_ , volume 167 of _Proceedings of Machine Learning Research_ , pages 1138–1164. PMLR, 29 Mar–01 Apr 2022.
  * Yang et al. (2018) Puyudi Yang, Cho-Jui Hsieh, and Jane-Ling Wang.  History PCA: A new algorithm for streaming PCA.  _arXiv preprint arXiv:1802.05447_ , 2018.
  * Yang and Xu (2015) Wenzhuo Yang and Huan Xu.  Streaming Sparse Principal Component Analysis.  In _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pages 494–503, Lille, France, 07–09 Jul 2015. PMLR.
  * Yun (2018) Se-Young Yun.  Noisy Power Method with Grassmann Average.  In _2018 IEEE International Conference on Big Data and Smart Computing (BigComp)_ , pages 709–712, 2018.
  * Yun et al. (2014) Se-Young Yun, Marc Lelarge, and Alexandre Proutière.  Streaming, Memory Limited Algorithms for Community Detection.  In _Advances in Neural Information Processing Systems_ , volume 27, pages 3167–3175. Curran Associates, Inc., 2014.
  * Yun et al. (2015) Se-Young Yun, Marc Lelarge, and Alexandre Proutière.  Fast and Memory Optimal Low-Rank Matrix Approximation.  In _Advances in Neural Information Processing Systems_ , volume 28, pages 3166–3185. Curran Associates, Inc., 2015.
  * Zafar et al. (2019) Muhammad Bilal Zafar, Isabel Valera, Manuel Gomez-Rodriguez, and Krishna P. Gummadi.  Fairness Constraints: A Flexible Approach for Fair Classification .  _Journal of Machine Learning Research_ , 20(75):1–42, 2019.
  * Zemel et al. (2013) Rich Zemel, Yu Wu, Kevin Swersky, Toni Pitassi, and Cynthia Dwork.  Learning Fair Representations.  In _Proceedings of the 30th International Conference on Machine Learning_ , volume 28 of _Proceedings of Machine Learning Research_ , pages 325–333, Atlanta, Georgia, USA, 17–19 Jun 2013. PMLR.
  * Zhao et al. (2022) Tianxiang Zhao, Enyan Dai, Kai Shu, and Suhang Wang.  Towards Fair Classifiers Without Sensitive Attributes: Exploring Biases in Related Features.  In _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_ , WSDM ’22, page 1433–1442, New York, NY, USA, 2022\. Association for Computing Machinery.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 An Alternative Approach to Fair PCA
    1. 3.1 “Null It Out” Formulation of Fair PCA
    2. 3.2 An Explicit Characterization for Solution of Fair PCA
    3. 3.3 Comparison to the Existing Covariance Matching Constraint
  4. 4 Statistical Viewpoint: PAFO-Learnability of PCA
  5. 5 Algorithmic Viewpoint: Fair Streaming PCA
    1. 5.1 Our Algorithm: Fair Noisy Power Method (FNPM)
    2. 5.2 Previous Approaches are not Suitable for Streaming Setup
    3. 5.3 Extension to Multiple Sensitive Groups/Attributes
  6. 6 FNPM is a PAFO-Learning Algorithm
  7. 7 Experiments on real-world datasets
    1. 7.1 CelebA dataset
    2. 7.2 UCI datasets
  8. 8 Other Related Works
  9. 9 Conclusion
  10. A Broader Impacts, Limitations, and Future Directions
    1. A.1 Broader Impacts
    2. A.2 Limitations and Future Directions
  11. B Full Pseudocodes of Algorithms 1 and 2 for Fair “Streaming” PCA
  12. C Extension of Algorithms 3 to Multiple and Non-binary Attributes
  13. D More Detailed Comparison to Kleindessner et al. (2023)
    1. D.1 Their Approach
    2. D.2 Unsuitability for the Streaming Setting
  14. E Proofs of Lemma 2, 3, and Theorem 1
    1. E.1 Notations and Assumptions
    2. E.2 Matrix/Vector Concentration Inequalities
    3. E.3 Proof of Lemma 2 \- Bounding Error in Second Moment Gap
    4. E.4 Proof of Lemma 3 \- Bounding the Final Error
    5. E.5 Bounding the Estimation Error of Mean Difference
    6. E.6 Proof of Theorem 1 \- Sample Complexity for PAFO-learnability
  15. F More Experiments
    1. F.1 Synthetic Example
    2. F.2 Additional Results on CelebA Dataset
    3. F.3 Full Results on UCI Dataset

## Appendix

##  Appendix A Broader Impacts, Limitations, and Future Directions

###  A.1 Broader Impacts

This work proposes a dimensionality reduction method that addresses memory
efficiency _and_ fairness while providing a statistical guarantee. By
identifying and nullifying the “unfair direction” inherent in the data
distribution, we offer an alternative approach to fair PCA. We anticipate that
our approach will motivate researchers to explore other dimensionality
reduction techniques (e.g., auto-encoder) with fairness constraints.
Significantly, our contribution includes a rigorous theoretical framework,
PAFO learnability, which enables sample complexity analysis of fair PCA. We
envision the potential for our theory to be further generalized to broader
contexts, including alternative definitions of fairness and optimality of
different algorithms for fair machine learning.

On the application side, the memory efficiency of our method can facilitate
the scalability of fair PCA, making it viable for processing high-dimensional
datasets, even in scenarios where data points arrive in a streaming fashion.
One possible application of our approach is data pre-/post-processing to
alleviate unfairness, such as generating fair word embeddings by eliminating
sensitive attribute information through orthogonal projection. For more
detailed discussions on potential future directions, please refer to below.

###  A.2 Limitations and Future Directions

Here, we list some of our work’s limitations and possible extensions/future
directions.

#### Individual Fairness for PCA.

Our definition of fairness in PCA only covers group fairness via fair
presentation learning, which was also the case for the prior works (Olfat and
Aswani, 2019; Lee et al., 2022; Kleindessner et al., 2023). None of the works,
including ours, have yet to consider the notion of individual fairness (Dwork
et al., 2012) in the context of PCA, for which we do not have a definitive
answer.

#### Knowledge of sensitive attribute.

Our framework requires the full knowledge of the sensitive attribute $s$ for
all data points ${\bm{x}}$’s, which was also the case for all the previous
works (Kleindessner et al., 2023; Lee et al., 2022; Olfat and Aswani, 2019).
But, the assumption of such knowledge may not be feasible in the real world
due to privacy or legal reasons (Zhao et al., 2022; Lahoti et al., 2020), or
even just due to some extrinsic noises. Considering the case where the dataset
may lack sensitive attributes for some or all of the data points is an
important future direction.

#### Making the algorithm anytime.

In our formulation of fair PCA, our algorithm is two-phase, with the first
phase as a “burn-in” period for estimating the unfair subspace. Thus, it is
not an anytime algorithm in that if the algorithm stops whilst in the first
phase, then the resulting ${\bm{V}}$ is random and completely uninformative.
Designing an anytime variant of our algorithm is an important future
direction. One possible way to achieve that is to consider other streaming PCA
algorithms such as Oja’s method (Oja and Karhunen, 1985) or accelerated NPM
(Xu and Li, 2022; Xu, 2023).

#### Improving gap dependence.

Our algorithm’s current sample complexity analysis relies on the analysis of
NPM by Hardt and Price (2014), which relies on the immediate singular value
gap, $\sigma_{k}-\sigma_{k+1}$. Balcan et al. (2016) showed that considering
greater iteration rank $q\geq k$, i.e., by considering optimization variable
of greater size, leads to a better gap dependency: from
$\sigma_{k}-\sigma_{k+1}$ to $\sigma_{k}-\sigma_{q+1}$. However, this is
incompatible with our current definition of PAFO-learnability (Definition 2)
because with greater iteration rank, the solution ${\bm{V}}^{\star}$ to which
the explained variance should be compared against is not clear. The problem is
that the sample complexity guarantee of Hardt and Price (2014); Balcan et al.
(2016) is derived in terms of the sine angle between the top $k$-eigenspace of
the true covariance and the estimated $q$-dimensional subspace. Clearing this
up would allow for a better theoretical guarantee in sample complexity and
thus an important future direction.

#### Kernelizing our framework

Kernel PCA (Schölkopf et al., 1998) is a cornerstone in modern machine
learning that has lent itself to be an inspiration to many applications, both
theory and practice-wise. Unlike previous fair PCA works (Olfat and Aswani,
2019; Kleindessner et al., 2023), in which the authors provided a kernelized
version of their fair PCA algorithm, both our statistical framework and
memory-efficient algorithm for streaming setting do not readily extend to the
kernelized version. Algorithmically, to tackle the streaming setting, one may
take inspiration from streaming kernel PCA (Ghashami et al., 2016; Ullah et
al., 2018), which was in turn inspired by matrix sketching (Liberty, 2013).

#### Extending to other settings.

In a similar spirit, extending our formulation of fair PCA to other streaming
settings such as sparse (Wang and Lu, 2016; Yang and Xu, 2015), nonstationary
(Bienstock et al., 2022), or even distributionally robust settings (Vu et al.,
2022) would also be interesting.

#### More experiments

Last but not least, we’ve only performed experiments on synthetic, CelebA, and
UCI datasets. It would be interesting to try our framework (either the
statistical formulation, the streaming setting, or both) on datasets from
other domains, such as NLP and graphs. GloVe vectors (Pennington et al., 2014)
has been used as a benchmark in a similar task called “concept erasure”
(Ravfogel et al., 2020, 2022b); group fairness in spectral clustering of
graphs has been studied as well (Kleindessner et al., 2019; Wang et al.,
2023).

##  Appendix B Full Pseudocodes of Algorithms 1 and 2 for Fair “Streaming” PCA

Throughout the appendices, from hereon, we will refer to our algorithms by
Algorithms 3 and 4 instead of Algorithms 1 and 2.

1 Input: $d$, $m$, Block size $b$, Number of iterations $T$;

2 Output: A matrix $\widehat{\bm{U}}$ with orthonormal columns;

3 ${\bm{W}}_{0}\leftarrow\operatorname{{\tt QR}}({\mathcal{N}}(0,1)^{d\times
m})$;

4
$(\overline{{\bm{m}}}^{(0)},\overline{{\bm{m}}}^{(1)},B^{(0)},B^{(1)})\leftarrow(\bm{0}_{d},\bm{0}_{d},0,0)$;

5 for  _$t\in[T]$_ do

6 $(b_{t}^{(0)},b_{t}^{(1)})\leftarrow(0,0)$;

7
$({\bm{m}}_{t}^{(0)},{\bm{m}}_{t}^{(1)},{\bm{C}}_{t}^{(0)}\\!,{\bm{C}}_{t}^{(1)})\leftarrow(\bm{0}_{d},\bm{0}_{d},\bm{0}_{d\times
m},\bm{0}_{d\times m})$;

8 for  _$i\in[b]$_ do

9 Receive $(a,{\bm{x}})=(a_{(t-1)b+i},{\bm{x}}_{(t-1)b+i})$;

10 $b_{t}^{(a)}\leftarrow b_{t}^{(a)}+1$;

11 ${\bm{m}}_{t}^{(a)}\leftarrow{\bm{m}}_{t}^{(a)}+{\bm{x}}$;

12
${\bm{C}}_{t}^{(a)}\leftarrow{\bm{C}}_{t}^{(a)}+{\bm{x}}{\bm{x}}^{\intercal}{\bm{W}}_{t-1}$;

13

14 end for

15 foreach  _$a\in\\{0,1\\}$_ do

16 if  _$b_{t}^{(a)} >0$_ then

17 ${\bm{m}}_{t}^{(a)}\leftarrow\frac{1}{b_{t}^{(a)}}{\bm{m}}_{t}^{(a)}$;

18 ${\bm{C}}_{t}^{(a)}\leftarrow\frac{1}{b_{t}^{(a)}}{\bm{C}}_{t}^{(a)}$;

19

20 end if

21
$\overline{{\bm{m}}}^{(a)}\\!\leftarrow\\!\frac{B^{(a)}}{B^{(a)}+b_{t}^{(a)}}\overline{{\bm{m}}}^{(a)}+\frac{b_{t}^{(a)}}{B^{(a)}+b_{t}^{(a)}}{\bm{m}}_{t}^{(a)}$;

22 $B^{(a)}\\!\leftarrow B^{(a)}\\!+b_{t}^{(a)}$;

23

24 end foreach

25 ${\bm{W}}_{t}=\operatorname{{\tt
QR}}\left({\bm{C}}^{(1)}-{\bm{C}}^{(0)}\right)$;

26

27 end for

28

29$\widehat{{\bm{f}}}\leftarrow\overline{{\bm{m}}}^{(1)}-\overline{{\bm{m}}}^{(0)}$;

30
$\widetilde{{\bm{g}}}\leftarrow\widehat{{\bm{f}}}-{\bm{W}}_{T}{\bm{W}}_{T}^{\intercal}\widehat{{\bm{f}}}$;

31 if  _$\|\widehat{{\bm{g}}}\|_{2}=0$_ then

32 $\widehat{\bm{U}}={\bm{W}}_{T}$

33 else

34
$\widehat{\bm{U}}=\left[{\bm{W}}_{T}\,\big{|}\,\frac{\widehat{{\bm{g}}}}{\|\widehat{{\bm{g}}}\|_{2}}\right]$

35 end if

return $\widehat{\bm{U}}$

Algorithm 3 UnfairSubspace (= Algorithm 1)

1 Input: $d$, $k$, Block sizes $\mathcal{B}$, $b$, Numbers of iterations
${\mathcal{T}}$, $T$;

2 Output: ${\bm{V}}_{\mathcal{T}}\in\operatorname{St}(d,k)$;

3

4$\widehat{{\bm{U}}}\leftarrow{\tt UnfairSubspace}(b,T)$;

5 ${\bm{V}}_{0}\leftarrow\operatorname{{\tt QR}}({\mathcal{N}}(0,1)^{d\times
k})$;

6 for  _$\tau\in[{\mathcal{T}}]$_ do

7
${\bm{V}}_{\tau}\leftarrow{\bm{V}}_{\tau-1}-\widehat{{\bm{U}}}\widehat{{\bm{U}}}^{\intercal}{\bm{V}}_{\tau-1}$;

8 ${\bm{C}}\leftarrow\bm{0}_{d\times k}$;

9 for  _$j\in[{\mathcal{B}}]$_ do

10 Receive
$(\ast,\tilde{\bm{x}})=(\ast,\tilde{\bm{x}}_{(\tau-1){\mathcal{B}}+j})$;

11
${\bm{C}}\leftarrow{\bm{C}}+\frac{1}{{\mathcal{B}}}\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}{\bm{V}}_{\tau}$;

12

13 end for

14 ${\bm{V}}_{\tau}\leftarrow\operatorname{{\tt
QR}}\left({\bm{C}}-\widehat{{\bm{U}}}\widehat{{\bm{U}}}^{\intercal}{\bm{C}}\right)$;

15

16 end for

return ${\bm{V}}_{\mathcal{T}}$

Algorithm 4 Fair NPM (= Algorithm 2)

##  Appendix C Extension of Algorithms 3 to Multiple and Non-binary Attributes

We remark that Algorithm 3 is a special case of the pseudocode below. By
applying Algorithm 4 after applying this algorithm, we can handle the multiple
non-binary attribute case and approximately and memory-efficiently solve the
extended notion of fair PCA elaborated in Section 5.3.

1 Input: $d$, $m_{1},\cdots,m_{\ell}$, Block size $b$, Number of iterations
$T$;

2 Output: A matrix $\widehat{\bm{U}}$ with orthonormal columns;

3 for  _$r\in[\ell]$_ do

4 for  _$a\in[g_{r}]$_ do

5 ${\bm{W}}_{0}^{(r,a)}=\operatorname{{\tt QR}}({\mathcal{N}}(0,1)^{d\times
m_{r}})$;

6 $(\overline{{\bm{m}}}^{(r,a,{\rm in})},\overline{{\bm{m}}}^{(r,a,{\rm
out})},B^{(r,a)})\leftarrow(\bm{0}_{d},\bm{0}_{d},0)$

7 end for

8

9 end for

10for  _$t\in[T]$_ do

11 Receive
$\\{((a_{1,i},\cdots,a_{\ell,i}),{\bm{x}}_{i})\\}_{i=b(t-1)+1}^{bt}$;

12

13 for  _$r\in[\ell]$_ do

14 for  _$a\in[g_{r}]$_ do

15 $b_{t}^{(r,a)}=\sum_{i=b(t-1)+1}^{bt}\mathbbm{1}_{[a_{r,i}=a]}$;

16 $B^{(r,a)}\leftarrow B^{(r,a)}+b_{t}^{(r,a)}$;

17 if  _$b_{t}^{(a)} >0$_ then

18 ${\bm{m}}_{t}^{(r,a,{\rm
in})}=\frac{1}{b_{t}^{(a)}}\sum_{i=b(t-1)+1}^{bt}\mathbbm{1}_{[a_{r,i}=a]}{\bm{x}}_{i}$;

19 ${\bm{C}}_{t}^{(r,a,{\rm
in})}=\frac{1}{b_{t}^{(a)}}\sum_{i=b(t-1)+1}^{bt}\mathbbm{1}_{[a_{r,i}=a]}{\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}{\bm{W}}_{t-1}^{(r,a)}$;

20

21 else

22 $({\bm{m}}_{t}^{(r,a,{\rm in})},{\bm{C}}_{t}^{(r,a,{\rm
in})})\leftarrow(\bm{0}_{d},\bm{0}_{d\times m_{r}})$

23 end if

24 if  _$b_{t}^{(a)} <b$_ then

25 ${\bm{m}}_{t}^{(r,a,{\rm
out})}=\frac{1}{b-b_{t}^{(a)}}\sum_{i=b(t-1)+1}^{bt}\mathbbm{1}_{[a_{r,i}\neq
a]}{\bm{x}}_{i}$;

26 ${\bm{C}}_{t}^{(r,a,{\rm
out})}=\frac{1}{b-b_{t}^{(a)}}\sum_{i=b(t-1)+1}^{bt}\mathbbm{1}_{[a_{r,i}\neq
a]}{\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}{\bm{W}}_{t-1}^{(r,a)}$;

27

28 else

29 $({\bm{m}}_{t}^{(r,a,{\rm out})},{\bm{C}}_{t}^{(r,a,{\rm
out})})\leftarrow(\bm{0}_{d},\bm{0}_{d\times m_{r}})$

30 end if

31 $\overline{{\bm{m}}}^{(r,a,{\rm
in})}\leftarrow\frac{B^{(r,a)}-b_{t}^{(r,a)}}{B^{(a)}}\overline{{\bm{m}}}^{(r,a,{\rm
in})}+\frac{b_{t}^{(a)}}{B^{(a)}}{\bm{m}}_{t}^{(r,a,{\rm in})}$;

32 $\overline{{\bm{m}}}^{(r,a,{\rm
out})}\leftarrow\frac{(t-1)b-B^{(r,a)}+b_{t}^{(r,a)}}{tb-B^{(a)}}\overline{{\bm{m}}}^{(r,a,{\rm
out})}+\frac{b-b_{t}^{(a)}}{tb-B^{(a)}}{\bm{m}}_{t}^{(r,a,{\rm out})}$;

33 ${\bm{W}}_{t}^{(r,a)}=\operatorname{{\tt QR}}\left({\bm{C}}_{t}^{(r,a,{\rm
in})}-{\bm{C}}_{t}^{(r,a,{\rm out})}\right)$;

34

35 end for

36

37 end for

38

39 end for

40for  _$r\in[\ell]$_ do

41 for  _$a\in[g_{r}]$_ do

42 $\widehat{\bm{f}}^{(r,a)}=\overline{{\bm{m}}}_{t}^{(r,a,{\rm
in})}-\overline{{\bm{m}}}_{t}^{(r,a,{\rm out})}$;

43

44 end for

45

46 end for

47$\widehat{\bm{U}}=\operatorname{{\tt
QR}}\left(\left[{\bm{W}}_{t}^{(1,1)}\Bigm{|}\cdots\Bigm{|}{\bm{W}}_{t}^{(\ell,g_{\ell})}\Bigm{|}\widehat{\bm{f}}^{(1,1)}\Bigm{|}\cdots\Bigm{|}\widehat{\bm{f}}^{(\ell,g_{\ell})}\right]\right)$;

return $\widehat{\bm{U}}$

Algorithm 5 UnfairSubspace_Extended

##  Appendix D More Detailed Comparison to Kleindessner et al. (2023)

###  D.1 Their Approach

Kleindessner et al. (2023) considered the following formulation of fair PCA:

$\max_{{\bm{V}}\in\operatorname{St}(d,k)}\mathrm{tr}({\bm{V}}^{\intercal}{\bm{X}}^{\intercal}{\bm{X}}{\bm{V}}),\quad\text{subject
to }{\bm{f}}^{\intercal}{\bm{V}}={\bm{0}}\ \land\
{\bm{V}}^{\intercal}(\bm{\Sigma}_{1}-\bm{\Sigma}_{0}){\bm{V}}={\bm{0}},$ (13)

and proposed a reasonable approximation of the covariance constraint, which we
briefly describe here. Let us write
${\bm{S}}^{\prime}=\bm{\Sigma}_{1}-\bm{\Sigma}_{0}$ for brevity. As implicitly
assumed in Kleindessner et al. (2023), let us assume that
${\bm{f}}\not={\bm{0}}$. The mean constraint is dealt with first. Denoting
${\bm{U}}_{\bm{f}}\in\operatorname{St}(d,d-1)$ to be the matrix whose columns
form a basis $(d-1)$-dimensional nullspace of ${\bm{f}}$, the mean constraint
is satisfied if and only if ${\bm{V}}$ is of the form
${\bm{U}}_{\bm{f}}{\bm{U}}$ with ${\bm{U}}\in\operatorname{St}(d-1,k)$ as the
intermediate optimization variable. With this first reparametrization, the
optimization now becomes

$\max_{{\bm{U}}\in\operatorname{St}(d-1,k)}\mathrm{tr}({\bm{U}}^{\intercal}{\bm{U}}_{\bm{f}}^{\intercal}{\bm{X}}^{\intercal}{\bm{X}}{\bm{U}}_{\bm{f}}{\bm{U}}),\quad\text{subject
to
}{\bm{U}}^{\intercal}{\bm{U}}_{\bm{f}}^{\intercal}{\bm{S}}^{\prime}{\bm{U}}_{\bm{f}}{\bm{U}}={\bm{0}},$
(14)

which now has only the covariance constraint.

To deal with the possible infeasibility of the covariance constraint,
Kleindessner et al. (2023) proposed the following approach: letting
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}\in\operatorname{St}(d-1,l)$ be the
matrix whose columns are the $l$ smallest eigenvectors of
${\bm{U}}_{\bm{f}}^{\intercal}{\bm{S}}^{\prime}{\bm{U}}_{\bm{f}}$ (in
magnitude), ${\bm{U}}$ only needs to nullify the eigenspace spanned by the
remaining $d-1-l$ eigenvectors. To see which terms are ignored with the
relaxation, let
$\sum_{i=1}^{d}q_{i}^{\prime}{\bm{v}}^{\prime}{{\bm{v}}^{\prime}_{i}}^{\intercal}$
be the eigenvalue decomposition of
${\bm{U}}_{\bm{f}}^{\intercal}{\bm{S}}^{\prime}{\bm{U}}_{\bm{f}}$ with
$|q_{1}^{\prime}|\geq|q_{2}^{\prime}|\geq\cdots|q_{d}^{\prime}|\geq 0$. This
approach essentially ignores the constraint
${\bm{U}}^{\intercal}{\bm{E}}_{l}{\bm{U}}={\bm{0}}$, where

${\bm{E}}_{l}={\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}^{\intercal}{\bm{U}}_{\bm{f}}^{\intercal}{\bm{S}}^{\prime}{\bm{U}}_{\bm{f}}{\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}=\sum_{i=d-l}^{d-1}q_{i}^{\prime}{\bm{v}}^{\prime}{{\bm{v}}^{\prime}_{i}}^{\intercal}.$
(15)

The number $l\in\\{k,\cdots,d-1\\}$ controls how much the group-conditional
covariances will be equalized; a smaller $l$ means that the covariance
constraint is enforced more stringently and vice versa. Ultimately, the
learner can control $l$ as a hyperparameter to create a trade-off between
fairness and the explained variance. Moreover, if $|q_{i}^{\prime}|$’s are
negligible for all $i\geq d-1-l$, then
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}^{\intercal}{\bm{U}}_{\bm{f}}^{\intercal}{\bm{S}}^{\prime}{\bm{U}}_{\bm{f}}{\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}\approx{\bm{0}}$,
and the relaxation becomes tighter.

As any ${\bm{U}}$ of the form
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}\bm{\Lambda}$ satisfies the relaxed
covariance constraint, a second reparameterization in terms of the new
optimization variable $\bm{\Lambda}\in\operatorname{St}(l,k)$ gives us

$\max_{\bm{\Lambda}\in\operatorname{St}(l,k)}\mathrm{tr}(\bm{\Lambda}^{\intercal}{\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}^{\intercal}{\bm{U}}_{\bm{f}}^{\intercal}{\bm{X}}^{\intercal}{\bm{X}}{\bm{U}}_{\bm{f}}{\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}\bm{\Lambda}),$
(16)

which can be solved via the standard SVD-based approach for vanilla PCA. Then,
the final solution is obtained as
${\bm{V}}^{*}={\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}{\bm{U}}_{\bm{f}}\bm{\Lambda}^{*}$.

###  D.2 Unsuitability for the Streaming Setting

In Section 5.2, we provided a rough overview of why existing approaches to
fair PCA (Olfat and Aswani, 2019; Lee et al., 2022; Kleindessner et al., 2023)
are not amenable to our streaming setting. Especially as the approach of
Kleindessner et al. (2023) (described above) is almost like a PCA, one may
wonder if standard techniques used in streaming PCA (Mitliagkas et al., 2013)
can be used. Here, we argue in detail why that is not the case, which is in
sharp contrast to our approach and our reformulation of fair PCA that led to a
memory-efficient fair streaming PCA algorithm.

#### Memory constraint.

For streaming PCA without fairness constraints, the main objective could be
written as
$\mathrm{tr}({\bm{V}}^{\intercal}{\bm{X}}^{\intercal}{\bm{X}}{\bm{V}})=\sum_{i=1}^{N}\mathrm{tr}\left({\bm{V}}^{\intercal}{\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}{\bm{V}}\right)$,
which is easily amenable to memory-limited algorithm such as noisy power
method (Mitliagkas et al., 2013) or stochastic optimization (Oja, 1982). Both
approaches utilize the fact that instead of storing $d\times d$ matrices, one
only needs to store matrix-vector product of size $d\times 1$, e.g.,
${\bm{V}}^{\intercal}{\bm{x}}_{i}$. This is not the case for the approach of
Kleindessner et al. (2023). To see this, consider Eqn. (14) without the
covariance constraint, i.e., fair PCA with only the mean constraint. Even
here, although the objective can be written as
$\sum_{i=1}^{N}\mathrm{tr}\left({\bm{U}}^{\intercal}{\bm{U}}_{\bm{f}}^{\intercal}{\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}{\bm{U}}_{\bm{f}}{\bm{U}}\right)$,
we must know ${\bm{U}}_{\bm{f}}$ in order to proceed further. Moreover, as
${\bm{U}}_{\bm{f}}$ is of size ${\mathcal{O}}(d^{2})$, it cannot be stored nor
estimated explicitly. When the covariance constraint is also taken into
account, although the matrix in question
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}$ is of dimension $(d-1)\times l$, which
is within the memory constraint if $l={\mathcal{O}}(1)$, the computation of
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}$ requires the knowledge of
${\bm{U}}_{\bm{f}}$. This is because Kleindessner et al. (2023) dealt with the
two constraints sequentially (mean first, then covariance), which forced the
reparametrization to be done twice and, more importantly, coupling the memory
requirement for the computation of ${\bm{U}}_{\bm{f}}$ and
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}$.

#### Statistical consideration.

The statistical guarantee (global convergence, sample complexity) of streaming
PCA algorithms is often obtained by appropriately bounding the error term and
using proper matrix concentration inequalities (Tropp, 2015). For example, for
the sample complexity guarantee of noisy power method (Hardt and Price, 2014;
Balcan et al., 2016), as the learner performs power method on the empirical
covariance $\sum_{i}{\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}$ instead of the true
covariance $\bm{\Sigma}$, the proof proceeds by first showing that the
empirical covariance is close enough to the true covariance (i.e., the norm of
their error term is sufficiently bounded), which then implies that the
variance of the estimation isn’t too high. Thus to analyze the error bound of
the final iterates ${\bm{V}}$ when the approach of Kleindessner et al. (2023)
is extended to a streaming setting, one would have to bound the estimation
error of
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}^{\intercal}{\bm{U}}_{\bm{f}}^{\intercal}\bm{\Sigma}{\bm{U}}_{\bm{f}}{\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}$.
There are three sources of estimation error:
${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}},{\bm{U}}_{\bm{f}}^{\intercal},\bm{\Sigma}$.
Recalling that ${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}$ consists of the
eigenvectors of
${\bm{U}}_{\bm{f}}^{\intercal}{\bm{S}}^{\prime}{\bm{U}}_{\bm{f}}$, one can see
that the estimation error of ${\bm{M}}_{{\bm{f}},{\bm{S}}^{\prime}}$ is
actually nontrivially dependent on the estimation quality of both ${\bm{U}}$
and ${\bm{S}}^{\prime}$. Here, we say nontrivially because the error isn’t
simply bounded in a linear sense via the usual triangle inequality; it
requires rather intricate techniques involving eigenvector perturbation theory
(Golub and Loan, 2013; Davis and Kahan, 1969), which may require additional
assumptions on eigenvalue gaps. As one can see later in the proof, our
approach considers both constraints simultaneously and thus allows for a quite
simple theoretical analysis.

##  Appendix E Proofs of Lemma 2, 3, and Theorem 1

###  E.1 Notations and Assumptions

We first recall the data generation process of our interest. At every time
step, a pair of a sensitive attribute and a data point is sampled as
$a\sim\mathrm{Bernoulli}(p)$ and ${\bm{x}}|a\sim{\mathcal{D}}_{a}$,
respectively, where ${\mathcal{D}}_{a}$ has mean $\bm{\mu}_{a}$ and covariance
$\bm{\Sigma}_{a}$. This can be written more compactly as
${\bm{x}}\sim{\mathcal{D}}:=p_{0}{\mathcal{D}}_{0}+p_{1}{\mathcal{D}}_{1}$,
where we denote $p_{0}:=1-p$ and $p_{1}:=p$.

We recall a set of notation needed for the proof. The mean difference is
denoted by ${\bm{f}}=\bm{\mu}_{1}-\bm{\mu}_{0}$, while the second moment
difference is
${\bm{S}}=\mathbb{E}[{\bm{x}}{\bm{x}}^{\intercal}|a=1]-E[{\bm{x}}{\bm{x}}^{\intercal}|a=0]=\bm{\Sigma}_{1}-\bm{\Sigma}_{0}+\bm{\mu}_{1}\bm{\mu}_{1}^{\intercal}-\bm{\mu}_{0}\bm{\mu}_{0}^{\intercal}$.
Let ${\bm{P}}_{m}\in\mathbb{R}^{d\times m}$ be the matrix whose columns are
top $m$ eigenvectors of ${\bm{S}}$ in magnitude of eigenvalues. Note that
$t\in[T]$ is the time variable for Algorithm 3, which we often omit if the
context is clear. Moreover, $b$ (without any sub-/superscript) is the block
size of Algorithm 3. The remaining notation is listed below; most of them came
from the pseudocode of our algorithms.

  * •

Counts of data points for each sensitive attribute $a$:

$\displaystyle
b^{(a)}_{t}=\sum_{i=(t-1)b+1}^{tb}\mathbbm{1}_{[a_{i}=a]},\qquad
B^{(a)}=\sum_{t=1}^{T}b_{t}^{(a)}=\sum_{i=1}^{Tb}\mathbbm{1}_{[a_{i}=a]}$

  * •

Estimates of group-conditional means:

$\displaystyle\overline{{\bm{m}}}^{(a)}=\begin{cases}\displaystyle\frac{1}{TB^{(a)}}\sum_{i=1}^{Tb}\mathbbm{1}_{[a_{i}=a]}{\bm{x}}_{i}&\text{
if }B^{(a)}>0,\\\ \bm{0}&\text{ if }B^{(a)}=0\end{cases}$

  * •

Estimate of the group-conditional mean gap:
$\widehat{\bm{f}}=\overline{{\bm{m}}}^{(1)}-\overline{{\bm{m}}}^{(0)}$

  * •

Estimate of the group-conditional second moments:

$\displaystyle{\bm{A}}^{(a)}_{t}=\begin{cases}\displaystyle\frac{1}{b^{(a)}_{t}}\sum_{i=(t-1)b+1}^{tb}\mathbbm{1}_{[a_{i}=a]}{\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}&\text{
if }b^{(a)}_{t}>0,\\\ \bm{0}&\text{ if }b^{(a)}_{t}=0\end{cases}$

  * •

Estimate of the group-conditional second moment gap:
$\widehat{\bm{S}}_{t}={\bm{A}}^{(1)}_{t}-{\bm{A}}^{(0)}_{t}$

  * •

${\bm{W}}_{t}\in\operatorname{St}(d,n)$ (decision variable of Algorithm 1 and
3)

  * •

${\bm{C}}^{(a)}_{t}={\bm{A}}^{(a)}_{t}{\bm{W}}_{t-1}$

  * •

${\bm{F}}_{t}={\bm{C}}^{(1)}_{t}-{\bm{C}}^{(0)}_{t}=\widehat{\bm{S}}_{t}{\bm{W}}_{t-1}$

  * •

${\bm{W}}_{t}=\operatorname{{\tt QR}}({\bm{F}}_{t})$ (i.e.,
${\bm{F}}_{t}\stackrel{{\scriptstyle\operatorname{{\tt
QR}}}}{{=}}{\bm{W}}_{t}{\bm{R}}_{t}$, ${\bm{R}}_{t}\in\mathbb{R}^{m\times m}$
is a lower triangular matrix)

Recall the assumptions on the data distribution, which we will assume
throughout the proof: See 1 See 2

We provide some intuition on the assumptions that we impose here. Assumption 1
consists of three parts. The first part, which involves
$\sigma,V,{\mathcal{M}}$ and ${\mathcal{V}}$, ensures that the maximum
deviation in mean and covariance of each ${\mathcal{D}}_{a}$ are well bounded;
this is critical in allowing for us to use proper matrix concentration
inequalities (to be described in the next subsection) and has been used in
various streaming PCA literature (Huang et al., 2021; Bienstock et al., 2022;
Jain et al., 2016). The second part, which involves $g_{\min}$ and $f_{\max}$,
imposes a bound on the maximum mean separation,
${\bm{f}}=\bm{\mu}_{1}-\bm{\mu}_{0}$, $\ell_{2}$-wise and angle-wise,
respectively. If the mean difference can be arbitrarily large, then the
$\ell_{2}$-estimation error that one has to achieve becomes arbitrarily small;
precisely speaking, $\lVert{\bm{f}}\rVert_{2}$ acts as a Lipschitz constant.
On the other hand, if the mean difference can be arbitrarily small, then the
angle-wise estimation error becomes arbitrarily large. The last part, which
involves $p_{\min}$, ensures that both groups are selected with some positive,
nonvanishing probability. Assumption 2 is standard in streaming PCA literature
to ensure convergence; indeed, if the singular value gap is zero, a definitive
convergence result can never be obtained, as the ground-truth solution becomes
vague.

###  E.2 Matrix/Vector Concentration Inequalities

Before moving on to our proof, we review some useful concentration
inequalities for our theoretical analysis.

######  Definition 4 (Variance of random matrix).

For a zero-mean random matrix ${\bm{Z}}$, its variance is defined as

$\mathrm{Var}({\bm{Z}})=\max\left(\left\|\mathbb{E}\left[{\bm{Z}}{\bm{Z}}^{\intercal}\right]\right\|_{2},\left\|\mathbb{E}\left[{\bm{Z}}^{\intercal}{\bm{Z}}\right]\right\|_{2}\right).$

In general,
$\mathrm{Var}({\bm{Z}})=\mathrm{Var}({\bm{Z}}-\mathbb{E}[{\bm{Z}}])$.

######  Proposition 1 (Matrix Bernstein inequality (Theorem 6.6.1 of Tropp
(2015))).

Consider a finite collection $\\{{\bm{Y}}_{j}\\}_{j=1}^{b}$ of independent
matrices with the same size ($d_{1}\times d_{2}$). Suppose they are zero mean
and they have uniformly bounded singular values, i.e.,

$\mathbb{E}[{\bm{Y}}_{j}]=\bm{0}\quad\text{ and
}\quad{\mathbb{P}}\left(\left\|{\bm{Y}}_{j}\right\|\leq{\mathcal{M}}\right)=1\quad\text{for
each }j\in[b].$

Let ${\mathcal{V}}$ be an upper bound of matrix variance,
${\mathcal{V}}\geq\mathrm{Var}\left({\bm{Y}}_{j}\right)$ for all $j\in[b]$.
Then, for all $x\geq 0$,

${\mathbb{P}}\left(\left\|\frac{1}{b}\sum_{j=1}^{b}{\bm{Y}}_{j}\right\|\geq
x\right)\leq(d_{1}+d_{2})\exp\left(\frac{-bx^{2}}{2({\mathcal{V}}+{\mathcal{M}}x/3)}\right).$

In particular, if $0\leq x\leq\frac{3{\mathcal{V}}}{{\mathcal{M}}}$,

${\mathbb{P}}\left(\left\|\frac{1}{b}\sum_{j=1}^{b}{\bm{Y}}_{j}\right\|\geq
x\right)\leq(d_{1}+d_{2})\exp\left(\frac{-bx^{2}}{4{\mathcal{V}}}\right).$

######  Proposition 2 (Vector Hoeffding inequality (Corollary 7 of Jin et al.
(2019))).

There exists an absolute constant $\mathfrak{c}$ such that if
${\bm{y}}_{1},\cdots,{\bm{y}}_{b}$ be independent random vectors with common
dimension $d$, and assume the following:

$\mathbb{E}[{\bm{y}}_{i}]=0,\quad{\bm{y}}_{i}\in\mathrm{nSG}(\sigma).$

Then, for any $x\geq 0$,

${\mathbb{P}}\left(\left\lVert\frac{1}{b}\sum_{j=1}^{b}{\bm{y}}_{j}\right\rVert_{2}\geq
x\right)\leq 2d\exp\left(-\frac{bx^{2}}{\mathfrak{c}^{2}\sigma^{2}}\right).$

######  Remark 2.

With an additional assumption that the random variables are bounded, we can
also consider using vector Bernstein inequality (Lemma 18 of Kohler and Lucchi
(2017)), which removes the factor of dimension $d$ from the RHS of the
concentration inequality.

###  E.3 Proof of Lemma 2 \- Bounding Error in Second Moment Gap

Recall that Algorithm 3 performs NPM on the second moment difference matrix
${\bm{S}}$ so that an iteration can be written as

${\bm{W}}_{t}=\operatorname{{\tt
QR}}(\widehat{\bm{S}}_{t}{\bm{W}}_{t-1})=\operatorname{{\tt
QR}}({\bm{S}}{\bm{W}}_{t-1}+{\bm{Z}}_{t,1}),$

where ${\bm{Z}}_{t,1}=(\widehat{\bm{S}}_{t}-{\bm{S}}){\bm{W}}_{t-1}$ is the
noise matrix.

The lemma below provides a general statement on sufficient size $b$ of blocks
to bound the second moment difference error $\widehat{\bm{S}}_{t}-{\bm{S}}$
multiplied by some fixed matrices.

######  Lemma 4.

Let any $\delta>0$ and matrices ${\bm{M}}\in\mathbb{R}^{m\times d}$ and
${\bm{N}}\in\mathbb{R}^{d\times n}$ be given. If we choose
$b\geq\frac{4}{p_{\min}}\max\left(1,\frac{8{\mathcal{M}}^{2}}{9{\mathcal{V}}}\right)\log\frac{2(m+n)}{\delta}$,
then with probability at least $1-\delta$,

$\displaystyle\left\|{\bm{M}}(\widehat{\bm{S}}_{t}-{\bm{S}}){\bm{N}}\right\|_{2}\leq{\mathcal{E}}_{m+n}^{(\rm
S)}\triangleq\sqrt{\frac{32{\mathcal{V}}}{bp_{\min}}\log\frac{2(m+n)}{\delta}}.$

We now proceed the proof of Lemma 2. To this end, recall that
${\bm{Z}}_{t,1}=(\widehat{\bm{S}}_{t}-{\bm{S}}){\bm{W}}_{t-1}$. Firstly,
applying Lemma 4 with ${\bm{M}}={\bm{I}}_{d}$ and
${\bm{N}}={\bm{W}}_{t-1}\in\mathbb{R}^{d\times m}$, a sufficient condition for
having $5\left\|{\bm{Z}}_{t,1}\right\|_{2}\leq\epsilon\Delta_{m,\nu}$ with
probability at least $1-\delta$ is that

$\displaystyle
b\geq\frac{400{\mathcal{V}}}{\epsilon^{2}\Delta_{m,\nu}^{2}p_{\min}}\log\frac{2(d+m)}{\delta}+\frac{4}{p_{\min}}\max\left(1,\frac{8{\mathcal{M}}^{2}}{9{\mathcal{V}}}\right)\log\frac{2(d+m)}{\delta}.$

Secondly, applying Lemma 4 with
${\bm{M}}={\bm{P}}_{m}^{\intercal}\in\mathbb{R}^{m\times d}$ and
${\bm{N}}={\bm{W}}_{t-1}\in\mathbb{R}^{d\times m}$, a sufficient condition for
having
$5\left\|{\bm{P}}_{m}^{\intercal}{\bm{Z}}_{t,1}\right\|_{2}\leq\frac{\delta\Delta_{m,\nu}}{2\sqrt{dm}}$
with probability at least $1-\delta$ is that

$\displaystyle
b\geq\frac{3200dm{\mathcal{V}}}{\delta^{2}\Delta_{m,\nu}^{2}p_{\min}}\log\frac{4m}{\delta}+\frac{4}{p_{\min}}\max\left(1,\frac{8{\mathcal{M}}^{2}}{9{\mathcal{V}}}\right)\log\frac{4m}{\delta}.$

Combining all bounds for $b$, we conclude that a sufficient block size $b$ is

$b=\Omega\left(\frac{{\mathcal{V}}}{\Delta_{m,\nu}^{2}p_{\min}}\left(\frac{dm}{\delta^{2}}\log\frac{m}{\delta}+\frac{1}{\epsilon^{2}}\log\frac{d}{\delta}\right)+\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}p_{\min}}\log\frac{d}{\delta}\right).$

###### Proof of Lemma 4.

Note that $(a_{i},{\bm{x}}_{i})$’s are i.i.d. samples
($i\in\\{(t-1)b+1,\dots,tb\\}$). Consider independent random matrices

${\bm{Y}}_{i}^{(a)}={\bm{M}}\left({\bm{x}}_{i}{\bm{x}}_{i}^{\intercal}-(\bm{\Sigma}_{a}+\bm{\mu}_{a}\bm{\mu}_{a}^{\intercal})\right){\bm{N}}.$

For each set $A\subseteq[b]$, define an event
$E_{A}:=\left\\{a_{i}=\mathbbm{1}_{[i\in A]}\ \forall i\in[b]\right\\}$. Note
that ${\mathbb{P}}[E_{A}]=p_{0}^{b-|A|}p_{1}^{|A|}$. To exploit this, we first
apply a peeling argument as follows: for any $c_{0},c_{1}\in(0,1)$ with
$c_{0}+c_{1}=1$,

$\displaystyle{\mathbb{P}}\left(\left\|{\bm{M}}(\widehat{\bm{S}}_{t}-{\bm{S}}){\bm{N}}\right\|_{2}\geq
x\right)$
$\displaystyle=\sum_{r=0}^{b}\sum_{A\in\binom{[b]}{r}}{\mathbb{P}}\left(\left.\left\|{\bm{M}}(\widehat{\bm{S}}_{t}-{\bm{S}}){\bm{N}}\right\|_{2}\geq
x\right|E_{A}\right){\mathbb{P}}\left(E_{A}\right)$
$\displaystyle=\sum_{r=1}^{b-1}\sum_{A\in\binom{[b]}{r}}{\mathbb{P}}\left(\left.\left\|\frac{1}{r}\sum_{i\in
A}{\bm{Y}}_{i}^{(1)}-\frac{1}{b-r}\sum_{i\in
A^{\mathsf{c}}}{\bm{Y}}_{i}^{(0)}\right\|_{2}\geq
x\right|E_{A}\right){\mathbb{P}}\left(E_{A}\right)$
$\displaystyle\quad+{\mathbb{P}}\left(\left.\left\|\frac{1}{b}\sum_{i=1}^{b}{\bm{Y}}_{i}^{(0)}\right\|_{2}\geq
x\right|E_{\emptyset}\right){\mathbb{P}}\left(E_{\emptyset}\right)+{\mathbb{P}}\left(\left.\left\|\frac{1}{b}\sum_{i=1}^{b}{\bm{Y}}_{i}^{(1)}\right\|_{2}\geq
x\right|E_{[b]}\right){\mathbb{P}}\left(E_{[b]}\right)$
$\displaystyle\leq\sum_{r=1}^{b-1}\sum_{A\in\binom{[b]}{r}}\left\\{{\mathbb{P}}\left(\left.\left\|\frac{1}{b-r}\sum_{i\in
A^{\mathsf{c}}}{\bm{Y}}_{i}^{(0)}\right\|_{2}\geq
c_{0}x\right|E_{A}\right)+{\mathbb{P}}\left(\left.\left\|\frac{1}{r}\sum_{i\in
A}{\bm{Y}}_{i}^{(1)}\right\|_{2}\geq
c_{1}x\right|E_{A}\right)\right\\}p_{0}^{b-r}p_{1}^{r}$
$\displaystyle\quad+{\mathbb{P}}\left(\left.\left\|\frac{1}{b}\sum_{i=1}^{b}{\bm{Y}}_{i}^{(0)}\right\|_{2}\geq
c_{0}x\right|E_{\emptyset}\right)p_{0}^{b}+{\mathbb{P}}\left(\left.\left\|\frac{1}{b}\sum_{i=1}^{b}{\bm{Y}}_{i}^{(1)}\right\|_{2}\geq
c_{1}x\right|E_{[b]}\right)p_{1}^{b},$

where the last inequality holds for the following fact: if
$\left\|{\bm{a}}-{\bm{b}}\right\|\geq x$ then
$\left\|{\bm{a}}\right\|\geq\lambda x$ or
$\left\|{\bm{b}}\right\|\geq(1-\lambda)x$ for any $\lambda\in(0,1)$.

Now, we make the crucial observation that conditioned on the event $E_{A}$,
both $\left\\{{\bm{Y}}_{i}^{(1)}\right\\}_{i\in A}$ and
$\left\\{{\bm{Y}}_{i}^{(0)}\right\\}_{i\in A^{\mathsf{c}}}$ are sets of i.i.d.
random zero-mean $m\times n$ matrices. In addition, with Assumption 1, we have
that for each $i\in A$,

${\mathbb{P}}\left[\left.\left\|{\bm{Y}}_{i}^{(1)}\right\|_{2}\leq{\mathcal{M}}\right|i\in
A\right]=1,\quad\mathrm{Var}\left(\left.{\bm{Y}}_{i}^{(1)}\right|i\in
A\right)\leq{\mathcal{V}}.$

(and analogously for each $i\in A^{\mathsf{c}}$). Thus, applying the matrix
Bernstein inequality333Precisely, we use its conditional version where the
means and the variances are replaced with their conditional counterparts.
(Proposition 1), we obtain the following: when $0\leq
x\leq\frac{3{\mathcal{V}}}{{\mathcal{M}}}$,

$\displaystyle{\mathbb{P}}\left(\left\|(\widehat{\bm{S}}_{t}-{\bm{S}}){\bm{W}}_{t-1}\right\|_{2}\geq
x\right)$
$\displaystyle\leq(m+n)\sum_{r=1}^{b-1}p_{0}^{b-r}p_{1}^{r}\sum_{A\in\binom{[b]}{r}}\left\\{\exp\left(-\frac{(b-r)c_{0}^{2}x^{2}}{4{\mathcal{V}}}\right)+\exp\left(-\frac{rc_{1}^{2}x^{2}}{4{\mathcal{V}}}\right)\right\\}$
$\displaystyle\quad+(m+n)p_{0}^{b}\exp\left(-\frac{bc_{0}^{2}x^{2}}{4{\mathcal{V}}}\right)+(m+n)p_{1}^{b}\exp\left(-\frac{bc_{1}^{2}x^{2}}{4{\mathcal{V}}}\right)$
$\displaystyle\leq(m+n)\sum_{r=0}^{b}\binom{b}{r}p_{0}^{b-r}p_{1}^{r}\left\\{\exp\left(-\frac{(b-r)c_{0}^{2}x^{2}}{4{\mathcal{V}}}\right)+\exp\left(-\frac{rc_{1}^{2}x^{2}}{4{\mathcal{V}}}\right)\right\\}$
$\displaystyle=(m+n)\left\\{\left(p_{0}+p_{1}\exp\left({-\frac{c_{1}^{2}x^{2}}{4{\mathcal{V}}}}\right)\right)^{b}+\left(p_{1}+p_{0}\exp\left({-\frac{c_{0}^{2}x^{2}}{4{\mathcal{V}}}}\right)\right)^{b}\right\\}<\delta.$

For this to occur, it suffices for both terms to be bounded by
$\frac{\delta}{2}$. Let us consider only the first term, as the second term
follows from symmetry. Taking the log, we have

$\log\left((1-p_{1})+p_{1}\exp\left({-\frac{c_{1}^{2}x^{2}}{4{\mathcal{V}}}}\right)\right)\leq\frac{1}{b}\log\frac{\delta}{2(m+n)}.$

Using $\log(1+y)\leq y$, it suffices to have

$\frac{1}{b}\log\frac{2(m+n)}{\delta}\leq
p_{1}\left(1-\exp\left({-\frac{c_{1}^{2}x^{2}}{4{\mathcal{V}}}}\right)\right).$

Using $y\leq 1-e^{-2y}$ for $y\in[0,1/2]$, it now suffices to have

$\frac{8{\mathcal{V}}}{bp_{1}c_{1}^{2}}\log\frac{2(m+n)}{\delta}\leq
x^{2}\leq\min\left(\frac{2{\mathcal{V}}}{c_{1}^{2}},\frac{9{\mathcal{V}}^{2}}{{\mathcal{M}}^{2}}\right).$

Combining this with the other term and, for simplicity444One could try to
optimize for $c_{0},c_{1}$ to obtain an “optimal” sample complexities.
However, from some preliminary computations, this seems to be not worth
pursuing, and we conjecture that it would yield the same asymptotic dependency
as when $c_{0}=c_{1}=1/2$. choosing $c_{0}=c_{1}=1/2$, we have our desired
statement. ∎

###  E.4 Proof of Lemma 3 \- Bounding the Final Error

Recall that the iterates are ${\bm{V}}_{\tau+1}=\operatorname{{\tt
QR}}(\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\widehat{\bm{\Sigma}}_{\tau}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}_{\tau})$,
where $\widehat{{\bm{U}}}$ is the output of Algorithm 3,
$\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}:={\bm{I}}_{d}-\widehat{{\bm{U}}}\widehat{{\bm{U}}}^{\intercal}$,
and
$\widehat{\bm{\Sigma}}_{\tau}:=\frac{1}{{\mathcal{B}}}\sum_{j=(\tau-1){\mathcal{B}}+1}^{\tau{\mathcal{B}}}\tilde{\bm{x}}_{j}\tilde{\bm{x}}_{j}^{\intercal}$
is the sample covariance at time step $\tau$ of Algorithm 4. This can be
rewritten as

${\bm{V}}_{\tau+1}=\operatorname{{\tt
QR}}(\bm{\Pi}^{\perp}_{{\bm{U}}}\bm{\Sigma}\bm{\Pi}^{\perp}_{{\bm{U}}}{\bm{V}}_{\tau}+{\bm{Z}}_{\tau,2}),$

where
${\bm{Z}}_{\tau,2}=(\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\widehat{\bm{\Sigma}}_{\tau}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}-\bm{\Pi}^{\perp}_{{\bm{U}}}\bm{\Sigma}\bm{\Pi}^{\perp}_{{\bm{U}}}){\bm{V}}_{\tau}$
is the noise matrix and
$\bm{\Sigma}=\sum_{a\in\\{0,1\\}}p_{a}(\bm{\Sigma}_{a}+\bm{\mu}_{a}\bm{\mu}_{a}^{\intercal})$.
By Assumption 1, we have that $\lVert\bm{\Sigma}\rVert_{2}\leq V$.

The lemma below provides a general statement on a sufficient block size
${\mathcal{B}}$ to bound the covariance error
$\widehat{\bm{\Sigma}}-\bm{\Sigma}$ multiplied by some fixed matrices.

######  Lemma 5.

Let any $\delta>0$ and matrices ${\bm{M}}\in\mathbb{R}^{m\times d}$ and
${\bm{N}}\in\mathbb{R}^{d\times n}$ be given. If we choose
${\mathcal{B}}\geq\frac{4}{9}\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}+V^{2}}\log\frac{4(m+n)}{\delta}$,
then the following holds with probability at least $1-\delta$:

$\displaystyle\left\|{\bm{M}}(\widehat{\bm{\Sigma}}-\bm{\Sigma}){\bm{N}}\right\|_{2}$
$\displaystyle\leq{\mathcal{E}}_{m+n}^{(\rm\Sigma)}\triangleq\sqrt{\frac{4({\mathcal{V}}+V^{2})}{{\mathcal{B}}}\log\frac{4(m+n)}{\delta}}.$

We now proceed the proof of Lemma 3. Fix some
${\bm{V}}\in\operatorname{St}(d,k)$, and from the design of our Algorithm 4,
$\widehat{{\bm{U}}}$ is also fixed. By our given assumption on upper bound of
$\left\|\bm{\Pi}_{\widehat{{\bm{U}}}}-\bm{\Pi}_{{\bm{U}}}\right\|_{2}$, and
the fact that $\lVert{\bm{V}}\rVert_{2}=\left\|{\bm{Q}}_{k}\right\|_{2}=1$, we
have

$\displaystyle\left\|{\bm{Z}}_{\tau,2}\right\|_{2}$
$\displaystyle\leq\left\|\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}(\widehat{\bm{\Sigma}}-\bm{\Sigma})\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\right\|_{2}+2\lVert\bm{\Sigma}\rVert_{2}\lVert\bm{\Pi}_{\widehat{{\bm{U}}}}-\bm{\Pi}_{{\bm{U}}}\rVert_{2}$
$\displaystyle\leq\left\|\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}(\widehat{\bm{\Sigma}}-\bm{\Sigma})\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\right\|_{2}+\frac{\epsilon\Delta_{k,\kappa}}{10}$

and, with a similar logic,

$\displaystyle\left\|{\bm{Q}}_{k}^{\intercal}{\bm{Z}}_{\tau,2}\right\|_{2}\leq\left\|{\bm{Q}}_{k}^{\intercal}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}(\widehat{\bm{\Sigma}}-\bm{\Sigma})\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\right\|_{2}+\frac{\delta\Delta_{k,\kappa}}{20\sqrt{dk}}.$

Applying Lemma 5 with
${\bm{M}}=\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\in\mathbb{R}^{d\times d}$ and
${\bm{N}}=\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\in\mathbb{R}^{d\times
k}$, a sufficient condition for having
$\left\|\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}(\widehat{\bm{\Sigma}}-\bm{\Sigma})\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\right\|_{2}\leq\frac{\epsilon\Delta_{k,\kappa}}{10}$
with probability at least $1-\delta$ is that

$\displaystyle{\mathcal{B}}\geq\frac{400({\mathcal{V}}+V^{2})}{\epsilon^{2}\Delta_{k,\kappa}^{2}}\log\frac{4(d+k)}{\delta}+\frac{4}{9}\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}+V^{2}}\log\frac{4(d+k)}{\delta}.$

Furthermore, applying Lemma 5 with
${\bm{M}}={\bm{Q}}_{k}^{\intercal}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\in\mathbb{R}^{k\times
d}$ and
${\bm{N}}=\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\in\mathbb{R}^{d\times
k}$, a sufficient condition for having
$\left\|{\bm{Q}}_{k}^{\intercal}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}(\widehat{\bm{\Sigma}}-\bm{\Sigma})\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\right\|_{2}\leq\frac{\delta\Delta_{k,\kappa}}{20\sqrt{dk}}$
with probability at least $1-\delta$ is that

$\displaystyle{\mathcal{B}}\geq\frac{1600dk({\mathcal{V}}+V^{2})}{\delta^{2}\Delta_{k,\kappa}^{2}}\log\frac{8k}{\delta}+\frac{4}{9}\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}+V^{2}}\log\frac{8k}{\delta}.$

Combining all bounds for ${\mathcal{B}}$, we conclude that a sufficient block
size ${\mathcal{B}}$ to have our desired statement is

${\mathcal{B}}=\Omega\left(\frac{{\mathcal{V}}+V^{2}}{\Delta_{k,\kappa}^{2}}\left(\frac{dk}{\delta^{2}}\log\frac{k}{\delta}+\frac{1}{\epsilon^{2}}\log\frac{d}{\delta}\right)+\frac{{\mathcal{M}}^{2}}{{\mathcal{V}}+V^{2}}\log\frac{d}{\delta}\right).$

###### Proof of Lemma 5.

Note that $\tilde{\bm{x}}_{j}$’s are i.i.d. samples from ${\mathcal{D}}$
$(j\in\\{(\tau-1){\mathcal{B}}+1,\dots,\tau{\mathcal{B}}\\})$. Let

${\bm{Y}}_{j}={\bm{M}}\left(\tilde{\bm{x}}_{j}\tilde{\bm{x}}_{j}^{\intercal}-\bm{\Sigma}\right){\bm{N}}.$

Then, ${\bm{Y}}_{i}$ are i.i.d. random zero-mean $d\times k$ matrices across
$j$. In addition, with Assumption 1,

${\mathbb{P}}\left[\left\|\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}-\bm{\Sigma}\right\|_{2}\leq{\mathcal{M}}\right]=\sum_{a^{\prime}\in\\{0,1\\}}{\mathbb{P}}\left[\left.\left\|\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}-\bm{\Sigma}\right\|_{2}\leq{\mathcal{M}}\right|a=a^{\prime}\right]{\mathbb{P}}[a=a^{\prime}]=1,$

and

$\displaystyle\mathrm{Var}(\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal})$
$\displaystyle=\left\|\mathbb{E}\left[\mathbb{E}[(\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal})^{2}|a]\right]-\mathbb{E}\left[\mathbb{E}[\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}|a\right]^{2}\right\|_{2}$
$\displaystyle\leq\left\|\mathbb{E}\left[\mathbb{E}\left[(\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal})^{2}|a\right]-\mathbb{E}[\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}|a]^{2}\right]\right\|_{2}+\left\|\mathbb{E}\left[\mathbb{E}[\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}|a]^{2}\right]-\mathbb{E}[\mathbb{E}[\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}|a]]^{2}\right\|_{2}$
$\displaystyle\leq\sum_{a\in[0,1]}p_{a}\underbrace{\left\|\mathbb{E}[(\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal})^{2}|a]-\mathbb{E}[\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}|a]^{2}\right\|_{2}}_{\mathrm{Var}(\tilde{\bm{x}}\tilde{\bm{x}}^{\intercal}|a)}+\left\|\mathbb{E}\left[(\bm{\Sigma}_{a}+\bm{\mu}_{a}\bm{\mu}_{a}^{\intercal})^{2}\right]-\mathbb{E}\left[\bm{\Sigma}_{a}+\bm{\mu}_{a}\bm{\mu}_{a}^{\intercal}\right]^{2}\right\|_{2}$
$\displaystyle\leq{\mathcal{V}}+p_{0}p_{1}\left\|\bm{\Sigma}_{0}+\bm{\mu}_{0}\bm{\mu}_{0}^{\intercal}-\bm{\Sigma}_{1}+\bm{\mu}_{1}\bm{\mu}_{1}^{\intercal}\right\|_{2}^{2}$
$\displaystyle\leq{\mathcal{V}}+4p_{0}(1-p_{0})V^{2}$
$\displaystyle\leq{\mathcal{V}}+V^{2}.$

Applying the matrix Bernstein inequality (Proposition 1), we obtain the
following: when $0\leq x\leq\frac{3({\mathcal{V}}+V^{2})}{{\mathcal{M}}}$,

${\mathbb{P}}\left(\left\|\frac{1}{{\mathcal{B}}}\sum_{i=1}^{{\mathcal{B}}}{\bm{Y}}_{i}\right\|_{2}\geq
x\right)\leq(m+n)\exp\left(-\frac{{\mathcal{B}}x^{2}}{4({\mathcal{V}}+V^{2})}\right)\leq\frac{\delta}{4}.$

Solving for $x$, we have

$\frac{4({\mathcal{V}}+V^{2})}{{\mathcal{B}}}\log\frac{4(m+n)}{\delta}\leq
x^{2}\leq 9\left(\frac{({\mathcal{V}}+V^{2})}{{\mathcal{M}}}\right)^{2},$

from which the statement naturally follows. ∎

###  E.5 Bounding the Estimation Error of Mean Difference

######  Lemma 6.

For any $\delta>0$, choose $b\geq\frac{2}{Tp_{\min}}\log\frac{4d}{\delta}$.
Then, the following holds with probability at least $1-\delta$:

$\displaystyle\left\|\widehat{\bm{f}}-{\bm{f}}\right\|_{2}$
$\displaystyle\leq{\mathcal{E}}^{\rm(f)}\triangleq\sqrt{\frac{8\mathfrak{c}^{2}\sigma^{2}}{Tbp_{\min}}\log\frac{4d}{\delta}},$

where $\mathfrak{c}$

###### Proof.

Again, note that $(a_{i},{\bm{x}}_{i})$’s are i.i.d. samples ($i\in[Tb]$).
Consider independent random vectors

${\bm{y}}_{i}^{(a)}={\bm{x}}_{i}-\bm{\mu}_{a}.$

For each $A\subseteq[Tb]$, $E_{A}=\left\\{a_{i}=\mathbbm{1}[i\in A]\ \forall
i\in[b]\right\\}$ is an event satisfying
${\mathbb{P}}(E_{A})=p_{0}^{Tb-|A|}p_{1}^{|A|}$. To exploit this, we now apply
the peeling argument as follows: for any $c_{0},c_{1}\in(0,1)$ with
$c_{0}+c_{1}=1$,

$\displaystyle{\mathbb{P}}\left(\left\|\widehat{\bm{f}}-{\bm{f}}\right\|_{2}\geq
x\right)$
$\displaystyle=\sum_{r=0}^{Tb}\sum_{A\in\binom{[Tb]}{r}}{\mathbb{P}}\left(\left.\left\|\widehat{\bm{f}}-{\bm{f}}\right\|_{2}\geq
x\right|E_{A}\right){\mathbb{P}}\left(E_{A}\right)$
$\displaystyle=\sum_{r=1}^{Tb-1}\sum_{A\in\binom{[Tb]}{r}}{\mathbb{P}}\left(\left.\left\|\frac{1}{r}\sum_{i\in
A}{\bm{y}}_{i}^{(1)}-\frac{1}{Tb-r}\sum_{i\in
A^{\mathsf{c}}}{\bm{y}}_{i}^{(0)}\right\|_{2}\geq
x\right|E_{A}\right){\mathbb{P}}\left(E_{A}\right)$
$\displaystyle\quad+{\mathbb{P}}\left(\left.\left\|\frac{1}{Tb}\sum_{i=1}^{Tb}{\bm{y}}_{i}^{(0)}\right\|_{2}\geq
x\right|E_{\emptyset}\right){\mathbb{P}}\left(E_{\emptyset}\right)+{\mathbb{P}}\left(\left.\left\|\frac{1}{Tb}\sum_{i=1}^{Tb}{\bm{y}}_{i}^{(1)}\right\|_{2}\geq
x\right|E_{[Tb]}\right){\mathbb{P}}\left(E_{[Tb]}\right)$
$\displaystyle\leq\sum_{r=1}^{Tb-1}\sum_{A\in\binom{[Tb]}{r}}\left\\{{\mathbb{P}}\left(\left.\left\|\frac{1}{Tb-r}\sum_{i\in
A^{\mathsf{c}}}{\bm{y}}_{i}^{(0)}\right\|_{2}\geq
c_{0}x\right|E_{A}\right)+{\mathbb{P}}\left(\left.\left\|\frac{1}{r}\sum_{i\in
A}{\bm{y}}_{i}^{(1)}\right\|_{2}\geq
c_{1}x\right|E_{A}\right)\right\\}p_{0}^{Tb-r}p_{1}^{r}$
$\displaystyle\quad+{\mathbb{P}}\left(\left.\left\|\frac{1}{Tb}\sum_{i=1}^{Tb}{\bm{y}}_{i}^{(0)}\right\|_{2}\geq
x\right|E_{\emptyset}\right)p_{0}^{Tb}+{\mathbb{P}}\left(\left.\left\|\frac{1}{Tb}\sum_{i=1}^{Tb}{\bm{y}}_{i}^{(1)}\right\|_{2}\geq
x\right|E_{[Tb]}\right)p_{1}^{Tb}$

Observe that conditioned on the event $E_{A}$, both
$\left\\{{\bm{y}}_{i}^{(1)}\right\\}_{i\in A}$ and
$\left\\{{\bm{y}}_{i}^{(0)}\right\\}_{i\in A^{\mathsf{c}}}$ are sets of i.i.d.
random zero-mean $d$-dimensional vectors. In addition, with Assumption 1, we
have that for each $i\in A$, ${\bm{y}}_{i}^{(1)}\in\mathrm{nSG}(\sigma_{1})$
(and analogously for each $i\in A^{\mathsf{c}}$). Thus, applying the vector
Bernstein inequality (Proposition 2), we obtain the following: for $x>0$,

$\displaystyle{\mathbb{P}}\left(\left\|\widehat{\bm{f}}-{\bm{f}}\right\|_{2}\geq
x\right)$ $\displaystyle\leq
2d\sum_{r=1}^{Tb-1}p_{0}^{Tb-r}p_{1}^{r}\sum_{A\in\binom{[Tb]}{r}}\left\\{\exp\left(-\frac{(Tb-r)c_{0}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{0}^{2}}\right)+\exp\left(-\frac{rc_{1}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{1}^{2}}\right)\right\\}$
$\displaystyle\quad+2dp_{0}^{Tb}\exp\left(-\frac{Tbx^{2}}{\mathfrak{c}^{2}\sigma_{0}^{2}}\right)+2dp_{1}^{Tb}\exp\left(-\frac{Tbx^{2}}{\mathfrak{c}^{2}\sigma_{1}^{2}}\right)$
$\displaystyle\leq
2d\sum_{r=0}^{Tb}\binom{Tb}{r}p_{0}^{Tb-r}p_{1}^{r}\left\\{\exp\left(-\frac{(Tb-r)c_{0}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{0}^{2}}\right)+\exp\left(-\frac{rc_{1}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{1}^{2}}\right)\right\\}$
$\displaystyle=2d\left(p_{0}+p_{1}\exp\left({-\frac{c_{1}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{1}^{2}}}\right)\right)^{Tb}+2d\left(p_{1}+p_{0}\exp\left({-\frac{c_{0}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{1}^{2}}}\right)\right)^{Tb}<\delta.$

It suffices for each of both terms to be bounded by $\frac{\delta}{2}$. Let us
consider only the first term, as the second term follows from symmetry. Taking
the log, we have

$\log\left((1-p_{1})+p_{1}\exp\left({-\frac{c_{1}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{1}^{2}}}\right)\right)\leq\frac{1}{Tb}\log\frac{\delta}{4d}.$

Using $\log(1+y)\leq y$, it suffices to have

$\frac{1}{Tb}\log\frac{4d}{\delta}\leq
p_{1}\left(1-\exp\left({-\frac{c_{1}^{2}x^{2}}{\mathfrak{c}^{2}\sigma_{1}^{2}}}\right)\right).$

Using $y\leq 1-e^{-2y}$ for $y\leq 1/2$, it now suffices to have

$\frac{2\mathfrak{c}^{2}\sigma_{1}^{2}}{Tbp_{1}c_{1}^{2}}\log\frac{4d}{\delta}\leq
x^{2}\leq\frac{\mathfrak{c}^{2}\sigma_{1}^{2}}{c_{1}^{2}}.$

Combining this with the other term and, for simplicity (see the footnote on
pg. 25), choosing $c_{0}=c_{1}=1/2$, we have our desired statement. ∎

###  E.6 Proof of Theorem 1 \- Sample Complexity for PAFO-learnability

Let us fix $\varepsilon_{\rm f},\varepsilon_{\rm o},\delta\in(0,1)$. From the
definition of PAFO-learnability (Definition 2), with a large enough sample
size, we must guarantee the following with probability at least $1-\delta$:

$\mathrm{tr}\left({\bm{V}}^{\intercal}\bm{\Sigma}{\bm{V}}\right)\geq\mathrm{tr}\left({{\bm{V}}^{\star}}^{\intercal}\bm{\Sigma}{\bm{V}}^{\star}\right)-\varepsilon_{\rm
f},\quad\lVert\bm{\Pi}_{{\bm{U}}}{\bm{V}}\rVert_{2}\leq\varepsilon_{\rm o}.$

Let $\widehat{{\bm{U}}}$ be the final estimate of the unfair subspace
${\bm{U}}$ from Algorithm 1, and let ${\bm{V}}$ be the final estimate of
${\bm{Q}}_{k}={\bm{V}}^{\star}$, whose columns are top $k$ eigenvectors of
$\bm{\Pi}^{\perp}_{\bm{U}}\bm{\Sigma}\bm{\Pi}^{\perp}_{\bm{U}}$.

We first recall the von Neumann trace inequality (Mirsky, 1975; von Neumann,
1937):

######  Proposition 3 (Theorem H.1.g of Marshall et al. (2011)).

If ${\bm{A}},{\bm{B}}$ are two $n\times n$ Hermitian matrices, then

$\mathrm{tr}({\bm{A}}{\bm{B}})\leq\sum_{i=1}^{n}\lambda_{i}({\bm{A}})\lambda_{i}({\bm{B}}),$

where $\lambda_{i}(\cdot)$ is the $i$-th smallest eigenvalue.

Also, the following lemma implies that the sine angle between subspaces (with
the same dimension) and the corresponding Grassmannian (projection) distance
are the same.

######  Lemma 7.

Let ${\bm{A}}_{1}\in\mathbb{R}^{d\times k}$ and
${\bm{B}}_{1}\in\mathbb{R}^{d\times\ell}$ be semi-orthogonal matrices:
${\bm{A}}_{1}^{\intercal}{\bm{A}}_{1}={\bm{I}}_{k}$ and
${\bm{B}}_{1}^{\intercal}{\bm{B}}_{1}={\bm{I}}_{\ell}$. Define

$\displaystyle\operatorname{dist}(\mathrm{col}({\bm{A}}_{1}),\mathrm{col}({\bm{B}}_{1}))$
$\displaystyle\triangleq\left\|\bm{\Pi}_{{\bm{A}}_{1}}-\bm{\Pi}_{{\bm{B}}_{1}}\right\|_{2},$
$\displaystyle\sin{(\mathrm{col}({\bm{A}}_{1}),\mathrm{col}({\bm{B}}_{1}))}$
$\displaystyle\triangleq\left\|\bm{\Pi}_{{\bm{A}}_{1}}^{\perp}{\bm{B}}_{1}\right\|_{2}.$

Then, (i) if $k=\ell$,

$\displaystyle\operatorname{dist}(\mathrm{col}({\bm{A}}_{1}),\mathrm{col}({\bm{B}}_{1}))$
$\displaystyle=\sin{(\mathrm{col}({\bm{A}}_{1}),\mathrm{col}({\bm{B}}_{1}))}$
$\displaystyle=\sin{(\mathrm{col}({\bm{B}}_{1}),\mathrm{col}({\bm{A}}_{1}))}=\sqrt{1-\sigma_{\min}({\bm{A}}_{1}^{\intercal}{\bm{B}}_{1})^{2}},$

where $\sigma_{\min}({\bm{M}})$ is the minimum singular value of ${\bm{M}}$.
On the other hand, (ii) if $k<\ell$,

$\displaystyle\operatorname{dist}(\mathrm{col}({\bm{A}}_{1}),\mathrm{col}({\bm{B}}_{1}))$
$\displaystyle=\sin{(\mathrm{col}({\bm{A}}_{1}),\mathrm{col}({\bm{B}}_{1}))}=1$
$\displaystyle\geq\sin{(\mathrm{col}({\bm{B}}_{1}),\mathrm{col}({\bm{A}}_{1}))}=\sqrt{1-\sigma_{\min}({\bm{A}}_{1}^{\intercal}{\bm{B}}_{1})^{2}}.$

Now we can write the optimality (the first inequality) as follows:

$\displaystyle\mathrm{tr}\left({\bm{Q}}_{k}^{\intercal}\bm{\Sigma}{\bm{Q}}_{k}\right)-\mathrm{tr}\left({\bm{V}}^{\intercal}\bm{\Sigma}{\bm{V}}\right)$
$\displaystyle=\mathrm{tr}\left(\bm{\Sigma}({\bm{Q}}_{k}{\bm{Q}}_{k}^{\intercal}-{\bm{V}}{\bm{V}}^{\intercal})\right)$
$\displaystyle\overset{(a)}{\leq}\sum_{i=1}^{d}\lambda_{i}(\bm{\Sigma})\lambda_{i}({\bm{Q}}_{k}{\bm{Q}}_{k}^{\intercal}-{\bm{V}}{\bm{V}}^{\intercal})$
$\displaystyle\overset{(b)}{\leq}2k\lVert\bm{\Sigma}\rVert_{2}\lVert{\bm{Q}}_{k}{\bm{Q}}_{k}^{\intercal}-{\bm{V}}{\bm{V}}^{\intercal}\rVert_{2}$
$\displaystyle\overset{(c)}{\leq}2kV\lVert\bm{\Pi}_{{\bm{V}}}^{\perp}{\bm{Q}}_{k}\rVert_{2}\leq\varepsilon_{\rm
f}.$

Here, $(a)$ follows from the von Neumann trace inequality, $(b)$ follows from
the fact that
${\bm{Q}}_{k}{\bm{Q}}_{k}^{\intercal}-{\bm{V}}{\bm{V}}^{\intercal}$ is a
symmetric matrix of rank at most $2k$, and $(c)$ follows from Lemma 7. Thus,
for the optimality, it suffices to ensure that
$\lVert\bm{\Pi}_{{\bm{V}}}^{\perp}{\bm{Q}}_{k}\rVert_{2}\leq\frac{\varepsilon_{\rm
f}}{2kV}$.

Similarly, for the fairness (the second inequality), it suffices to ensure
that
$\lVert\bm{\Pi}_{{\bm{U}}}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\rVert_{2}\leq\varepsilon_{\rm
o}$ or
$\lVert\bm{\Pi}_{\widehat{{\bm{U}}}}-\bm{\Pi}_{\bm{U}}\rVert_{2}\leq\varepsilon_{\rm
o}$, as

$\displaystyle\left\|\bm{\Pi}_{{\bm{U}}}{\bm{V}}\right\|_{2}\stackrel{{\scriptstyle(\ast)}}{{=}}\left\|\bm{\Pi}_{{\bm{U}}}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}\right\|_{2}\leq\left\|\bm{\Pi}_{{\bm{U}}}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\right\|_{2}\stackrel{{\scriptstyle(\star)}}{{\leq}}\left\|\bm{\Pi}_{\widehat{{\bm{U}}}}-\bm{\Pi}_{\bm{U}}\right\|_{2}\leq\varepsilon_{\rm
o},$

where $(\ast)$ follows from
${\bm{V}}=\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}{\bm{V}}$ (due to the design of
Algorithm 4) and $(\star)$ follows from Lemma 7.

From Lemma 1 (letting
${\mathcal{T}}=\Theta\left(\frac{K_{k,\kappa}}{\Delta_{k,\kappa}}\log\frac{d}{\epsilon\delta}\right)$)
and 3 (substituting $\epsilon\mapsto\frac{\varepsilon_{\rm f}}{2kV}$), given
that

$\left\|\bm{\Pi}_{\widehat{{\bm{U}}}}-\bm{\Pi}_{{\bm{U}}}\right\|_{2}\leq\eta_{k}=\eta_{k}(\varepsilon_{\rm
f},\varepsilon_{\rm o},\delta)\triangleq\min\left(\varepsilon_{\rm
o},\frac{\Delta_{k,\kappa}\varepsilon_{\rm
f}}{40V^{2}k},\frac{\Delta_{k,\kappa}\delta}{160V\sqrt{dk}}\right),$ (17)

a total of $N_{2}={\mathcal{T}}{\mathcal{B}}$ samples are sufficient in
Algorithm 4 for ensuring $\varepsilon_{\rm f}$-optimality with probability at
least $1-\frac{\delta}{2}$, where

$N_{2}\gtrsim\left(\frac{K_{k,\kappa}({\mathcal{V}}+V^{2})}{\Delta_{k,\kappa}^{3}}\left(\frac{dk}{\delta^{2}}\log\frac{k}{\delta}+\frac{k^{2}V^{2}}{\varepsilon_{\rm
f}^{2}}\log\frac{d}{\delta}\right)+\frac{K_{k,\kappa}{\mathcal{M}}^{2}}{\Delta_{k,\kappa}({\mathcal{V}}+V^{2})}\log\frac{d}{\delta}\right)\log\frac{dkV}{\varepsilon_{\rm
f}\delta}.$ (18)

We now focus on obtaining the sample complexity for Algorithm 3 to satisfy
Eqn. (17) with probability at least $1-\frac{\delta}{2}$. Then combining those
with a union bound gives us the desired statement.

From hereon and forth, we write

$\displaystyle{\bm{g}}$
$\displaystyle\triangleq\bm{\Pi}_{{\bm{P}}_{m}}^{\perp}{\bm{f}},\quad{\bm{h}}\triangleq\frac{1}{\left\|{\bm{g}}\right\|_{2}}{\bm{g}},$
$\displaystyle\widehat{\bm{g}}$
$\displaystyle\triangleq\bm{\Pi}_{{\bm{W}}_{T}}^{\perp}\widehat{\bm{f}},\quad\widehat{\bm{h}}\triangleq\frac{1}{\left\|\widehat{\bm{g}}\right\|_{2}}\widehat{\bm{g}}.$

Also, we introduce the following lemma that provides a general bound of sine
angle between a pair of vectors with a small $\ell_{2}$-distance.

######  Lemma 8.

Consider two vectors ${\bm{a}},{\bm{b}}\in\mathbb{R}^{d}$. If we denote the
(acute) angle between ${\bm{a}}$ and ${\bm{b}}$ as
$\theta({\bm{a}},{\bm{b}})$,

$\sin\theta({\bm{a}},{\bm{b}})=\sin(\operatorname{span}({\bm{a}}),\operatorname{span}({\bm{b}}))=\left\|\frac{1}{\left\|{\bm{a}}\right\|_{2}^{2}}{\bm{a}}{\bm{a}}^{\intercal}-\frac{1}{\left\|{\bm{b}}\right\|_{2}^{2}}{\bm{b}}{\bm{b}}^{\intercal}\right\|_{2}.$

Consider any
${\epsilon}\in\left(0,\frac{\left\|{\bm{a}}\right\|_{2}}{2}\right]$ and
suppose $\left\|{\bm{a}}-{\bm{b}}\right\|_{2}\leq{\epsilon}$. Then,

$\sin\theta({\bm{a}},{\bm{b}})\leq\frac{\sqrt{2}{\epsilon}}{\left\|{\bm{a}}\right\|_{2}}.$

#### Case 1.

Assume ${\bm{g}}\neq\bm{0}$. In this case, ${\bm{U}}=[{\bm{P}}_{m}|{\bm{h}}]$,
$\left\|{\bm{g}}\right\|_{2}\geq g_{\min}$, and
$\left\|{\bm{f}}\right\|_{2}\leq f_{\max}$ by Assumption 1. We want to upper-
bound the following probability to be $<\delta/2$:

$\displaystyle{\mathbb{P}}\left[\left\|\bm{\Pi}_{{\bm{U}}}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\right\|_{2}>\eta_{k}\right]$
$\displaystyle={\mathbb{P}}\left[\left(\widehat{\bm{g}}\neq\bm{0}\text{ \bf
and
}\left\|\bm{\Pi}_{{\bm{U}}}-\bm{\Pi}_{\widehat{{\bm{U}}}}\right\|_{2}>\eta_{k}\right)\text{
\bf or }\left(\widehat{\bm{g}}=\bm{0}\text{ \bf and
}\left\|\bm{\Pi}_{{\bm{U}}}\bm{\Pi}^{\perp}_{\widehat{{\bm{U}}}}\right\|_{2}>\eta_{k}\right)\right]$