Year: 2021

Abstract:
Principal component analysis is a simple yet useful dimensionality reduction technique in modern machine learning pipelines. In consequential domains such as college admission, healthcare and credit approval, it is imperative to take into account emerging criteria such as the fairness and the robustness of the learned projection. In this paper, we propose a distributionally robust optimization problem for principal component analysis which internalizes a fairness criterion in the objective function. The learned projection thus balances the trade-off between the total reconstruction error and the reconstruction error gap between subgroups, taken in the min-max sense over all distributions in a moment-based ambiguity set. The resulting optimization problem over the Stiefel manifold can be efficiently solved by a Riemannian subgradient descent algorithm with a sub-linear convergence rate. Our experimental results on real-world datasets show the merits of our proposed method over state-of-the-art baselines.

heading: 1 INTRODUCTION
text: Machine learning models are ubiquitous in our daily lives and supporting the decision-making process in diverse domains. With their flourishing applications, there also surface numerous concerns regarding the fairness of the models’ outputs (Mehrabi et al., 2021). Indeed, these models are prone to biases due to various reasons (Barocas et al., 2018). First, the collected training data is likely to include some demographic disparities due to the bias in the data acquisition process (e.g., conducting surveys on a specific region instead of uniformly distributed places), or the imbalance of observed events at a specific period of time. Second, because machine learning methods only care about data statistics and are objective driven, groups that are under-represented in the data can be neglected in exchange for a better objective value. Finally, even human feedback to the predictive models can also be biased, e.g., click counts are human feedback to recommendation systems but they are highly correlated with the menu list suggested previously by a potentially biased system. Real-world examples of machine learning models that amplify biases and hence potentially cause unfairness are commonplace, ranging from recidivism prediction giving higher false positive rates for African-American1 to facial recognition systems having large error rate for women2.
To tackle the issue, various fairness criteria for supervised learning have been proposed in the machine learning literature, which encourage the (conditional) independence of the model’s predictions on a particular sensitive attribute (Dwork et al., 2012; Hardt et al., 2016b; Kusner et al., 2017; Chouldechova, 2017; Verma & Rubin, 2018; Berk et al., 2021). Strategies to mitigate algorithmic bias are also investigated for all stages of the machine learning pipelines (Berk et al., 2021). For the pre-processing steps, (Kamiran & Calders, 2012) proposed reweighting or resampling techniques to achieve statistical parity between subgroups; or in the training steps, fairness can be encouraged by adding constraints (Donini et al., 2018) or regularizing the original objective function (Kamishima et al., 2012; Zemel et al., 2013); and in the post-processing steps, adjusting classification threshold based on examining black-box models over a holdout dataset can be used (Hardt et al., 2016b; Wei et al., 2019).
Since biases may already exist in the raw data, it is reasonable to demand the machine learning pipeline to combat biases as early as possible. We focus in this paper on the Principal Component 1 https://www.propublica.org/article/machine-bias-risk-assessments-in-criminal-sentencing 2 https://news.mit.edu/2018/study-finds-gender-skin-type-bias-artificial-intelligence-systems-0212
Analysis (PCA), which is a fundamental dimensionality reduction technique in the early stage of the pipelines (Pearson, 1901; Hotelling, 1933). PCA finds a linear transformation that embeds the original data into a lower-dimensional subspace that maximizes the variance of the projected data. Thus, PCA is prone to amplify biases if the data variability is different between the majority and the minority subgroups, see a toy example in Figure 1. A naive approach to promote fairness is to train one independent transformation for each subgroup. However, this requires knowing the sensitive attribute of each sample at test time, which would raise disparity concerns. On the contrary, using a single transformation for all subgroups is “group-blinded” and faces no discrimination problem (Lipton et al., 2018).
Learning a fair PCA has attracted attention from many fields from machine learning, statistics to signal processing. Samadi et al. (2018) and Zalcberg & Wiesel (2021) propose to find the principal components that minimize the maximum subgroup reconstruction error; the min-max formulations can be relaxed and solved as semidefinite programs. Olfat & Aswani (2019) propose to learn a transformation that minimizes the possibility of predicting the sensitive attribute from the projected data. Apart from being a dimensionality reduction technique, PCA can also be thought of as a representation learning toolkit. Viewed in this way, we can also consider a more general family of fair representation learning methods that can be applied before any further analysis steps. There are a number of works develop towards this idea (Kamiran & Calders, 2012; Zemel et al., 2013; Calmon et al., 2017; Feldman et al., 2015; Beutel et al., 2017; Madras et al., 2018; Zhang et al., 2018; Tantipongpipat et al., 2019), which apply a multitude of fairness criteria.
In addition, we also focus on the robustness criteria for the linear transformation. Recently, it has been observed that machine learning models are susceptible to small perturbations of the data (Goodfellow et al., 2014; Madry et al., 2017; Carlini & Wagner, 2017). These observations have fuelled many defenses using adversarial training (Akhtar & Mian, 2018; Chakraborty et al., 2018) and distributionally robust optimization (Rahimian & Mehrotra, 2019; Kuhn et al., 2019).
Contributions. This paper blends the ideas from the field of fairness in artifical intelligence and distributionally robust optimization. Our contributions can be described as follows.
• We propose the fair principal components which balance between the total reconstruction error and the absolute gap of reconstruction error between subgroups. Moreover, we also add a layer of robustness to the principal components by considering a min-max formulation that hedges against all perturbations of the empirical distribution in a moment-based ambiguity set.
• We provide the reformulation of the distributionally robust fair PCA problem as a finitedimensional optimization problem over the Stiefel manifold. We provide a Riemannian gradient descent algorithm and show that it has a sub-linear convergence rate.
Figure 1 illustrates the qualitative comparison between (fair) PCA methods and our proposed method on a 2-dimensional toy example. The majority group (blue dots) spreads on the horizontal axis, while the minority group (yellow triangles) spreads on the slanted vertical axis. The nominal PCA (red) captures the majority direction to minimize the total error, while the fair PCA of Samadi et al. (2018) returns the diagonal direction to minimize the maximum subgroup error. Our fair PCA can probe the full spectrum in between these two extremes by sweeping through our penalization parameters appropriately. If we do not penalize the error gap between subgroups, we recover the PCA method; if we penalize heavily, we recover the fair PCA of Samadi et al. (2018). Extensive numerical results on real datasets are provided in Section 5. Proofs are relegated to the appendix.

heading: 2 FAIR PRINCIPAL COMPONENT ANALYSIS

heading: 2.1 PRINCIPAL COMPONENT ANALYSIS
text: We first briefly revisit the classical PCA. Suppose that we are given a collection of N i.i.d. samples {x̂i}Ni=1 generated by some underlying distribution P. For simplicity, we assume that both the empirical and population mean are zero vectors. The goal of PCA is to find a k-dimensional linear subspace of Rd that explains as much variance contained in the data {x̂i}Ni=1 as possible, where k < d is a given integer. More precisely, we parametrize k-dimensional linear subspaces by orthonormal matrices, i.e., matrices whose columns are orthogonal and have unit Euclidean norm. Given any such matrix V , the associated k-dimensional subspace is the one spanned by the columns of V . The projection matrix onto the subspace is V V >, and hence the variance of the projected data is given by tr ( V V >ΞΞ> ) , where Ξ = [x̂1, · · · , x̂N ] ∈ Rd×N is the data matrix. By a slight abuse of terminology, sometimes we refer to V as the projection matrix. The problem of PCA then reads
max V ∈Rd×k,V >V=Ik
tr ( V V >ΞΞ> ) . (1)
For any vector X ∈ Rd and orthonormal matrix V , denote by `(V,X) the reconstruction error, i.e.,
`(V,X) = ‖X − V V >X‖22 = X>(Id − V V >)X. The problem of PCA can alternatively be formulated as a stochastic optimization problem
min V ∈Rd×k,V >V=Ik EP̂[`(V,X)], (2)
where P̂ is the empirical distribution associated with the samples {x̂i}Ni=1 and X ∼ P̂. It is wellknown that PCA admits an analytical solution. In particular, the optimal solution to problem (2) (and also problem (1)) is given by any orthonormal matrix whose columns are the eigenvectors associated with the k largest eigenvalues of the sample covariance matrix ΞΞ>.

heading: 2.2 FAIR PRINCIPAL COMPONENT ANALYSIS
text: In the fair PCA setting, we are also given a discrete sensitive attributeA ∈ A, whereAmay represent features such as race, gender or education. We consider binary attribute A and let A = {0, 1}. A straightforward idea to define fairness is to require the (strict) balance of a certain objective between the two groups. For example, this is the strategy in Hardt et al. (2016a) for developing fair supervised learning algorithms. A natural objective to balance in the PCA context is the reconstruction error. It is therefore tempted to adopt the following definition. Definition 2.1 (Fair projection). Let Q be an arbitrary distribution of (X,A). A projection matrix V ∈ Rd×k is fair relative to Q if the conditional expected reconstruction error is equal between subgroups, i.e.,
EQ[`(V,X)|A = a] = EQ[`(V,X)|A = a′] ∀(a, a′) ∈ A×A.
Unfortunately, Definition 2.1 is too stringent: for a general probability distribution Q, it is possible that there exists no fair projection matrix V . Proposition 2.2 (Impossibility result). For any distribution Q on X × A, let S = EQ[XX>|A = 0] − EQ[XX>|A = 1]. Then, there exists a fair projection matrix V ∈ Rd×k relative to Q if and only if rank(S) ≤ k.
One way to circumvent the impossibility result is to relax the requirement of strict balance to approximate balance. In other words, an inequality constraint of the following form is imposed:
|EQ[`(V,X)|A = a]− EQ[`(V,X)|A = a′]| ≤ ∀(a, a′) ∈ A×A, where > 0 is some prescribed fairness threshold. This approach has been adopted in other fair machine learning settings, see Donini et al. (2018) and Agarwal et al. (2019) for example.
In this paper, instead of imposing the fairness requirement as a constraint, we penalize the unfairness in the objective function. Specifically, for any projection matrix V , we define the unfairness as the absolute difference between the conditional loss between two subgroups:
U(V,Q) , |EQ[`(V,X)|A = 0]− EQ[`(V,X)|A = 1]|.
We thus consider the following fairness-aware PCA problem
min V ∈Rd×k, V >V=Ik EP̂[`(V,X)] + λU(V, P̂), (3)
where λ ≥ 0 is a penalty parameter to encourage fairness. Note that for fair PCA, the dataset is {(x̂i, âi)}Ni=1 and hence the empirical distribution P̂ is given by P̂ = 1N ∑N i=1 δ(x̂i,âi).

heading: 3 DISTRIBUTIONALLY ROBUST FAIR PCA
text: The weakness of empirical distribution-based stochastic optimization has been well-documented, see (Smith & Winkler, 2006; Homem-de Mello & Bayraksan, 2014). In particular, due to overfitting, the out-of-sample performance of the decision, prediction, or estimation obtained from such a stochastic optimization model is unsatisfactory, especially in the low sample size regime. Ideally, we could improve the performance by using the underlying distribution P instead of the empirical distribution P̂. But the underlying distribution P is unavailable in most practical situations, if not all. Distributional robustification is an emerging approach to handle this issue and has been shown to deliver promising out-of-sample performance in many applications (Delage & Ye, 2010; Namkoong & Duchi, 2017; Kuhn et al., 2019; Rahimian & Mehrotra, 2019). Motivated by the success of distributional robustification, we propose a robustified version of model (3), called the distributionally robust fairness-aware PCA:
min V ∈Rd×k,V >V=Ik sup Q∈B(P̂)
EQ[`(V,X)] + λU(V,Q), (4)
where B(P̂) is a set of probability distributions similar to the empirical distribution P̂ in a certain sense, called the ambiguity set. The empirical distribution P̂ is also called the nominal distribution. Many different ambiguity sets have been developed and studied in the optimization literature, see Rahimian & Mehrotra (2019) for an extensive overview.

heading: 3.1 THE WASSERSTEIN-TYPE AMBIGUITY SET
text: To present our ambiguity set and main results, we need to introduce some definitions and notations. Definition 3.1 (Wasserstein-type divergence). The divergenceW between two probability distributions Q1 ∼ (µ1,Σ1) ∈ Rd × Sd+ and Q2 ∼ (µ2,Σ2) ∈ Rd × Sd+ is defined as
W ( Q1 ‖ Q2 ) , ‖µ1 − µ2‖22 + tr ( Σ1 + Σ2 − 2 ( Σ 1 2 2 Σ1Σ 1 2 2 ) 1 2 ) .
The divergence W coincides with the squared type-2 Wasserstein distance between two Gaussian distributions N (µ1,Σ1) and N (µ2,Σ2) (Givens & Shortt, 1984). One can readily show that W is non-negative, and it vanishes if and only if (µ1,Σ1) = (µ2,Σ2), which implies that Q1 and Q2 have the same first- and second-moments.
Recall that the nominal distribution is P̂ = 1N ∑N i=1 δ(x̂i,âi). For any a ∈ A, its conditional distribution given A = a is given by
P̂a = 1 |Ia| ∑ i∈Ia δxi , where Ia , {i ∈ {1, . . . , N} : ai = a}.
We also use (µ̂a, Σ̂a) to denote the empirical mean vector and covariance matrix of X given A = a:
µ̂a = EP̂a [X] = EP̂[X|A = a] and Σ̂a + µ̂aµ̂ > a = EP̂a [XX >] = EP̂[XX >|A = a].
For any a ∈ A, the empirical marginal distribution of A is denoted by p̂a = |Ia|/N . Finally, for any set S, we use P(S) to denote the set of all probability distributions supported on S. For any integer k, the k-by-k identity matrix is denoted Ik. We then define our ambiguity set as
B(P̂) , Q ∈ P(X ×A) : ∃Qa ∈ P(X ) such that: Q(X× {a}) = p̂aQa(X) ∀X ⊆ Rd, a ∈ A W(Qa, P̂a) ≤ εa ∀a ∈ A  , (5)
where Qa is the conditional distribution of X|A = a. Intuitively, each Q ∈ B(P̂) is a joint distribution of the random vector (X,A), formed by taking a mixture of conditional distributions Qa with mixture weight p̂a. Each conditional distribution Qa is constrained in an εa-neighborhood from the nominal conditional distribution P̂a with respect to the W divergence. Notice that because the loss function ` is a quadratic function of X , the (conditional) expected losses only involve the first two moments of X , and thus prescribing the ambiguity set using W would suffice for the purpose of robustification.

heading: 3.2 REFORMULATION
text: We now present the reformulation of problem (4) under the ambiguity set B(P̂). Theorem 3.2 (Reformulation). Suppose that for any a ∈ A, either of the following two conditions holds:
(i) 0 ≤ λ ≤ p̂a, (ii) the empirical second moment matrix M̂a = 1Na ∑ i∈Ia x̂ix̂ > i satisfies ∑d−k j=1 σj(M̂a) ≥ εa,
where σj(M̂a) is the j-th smallest eigenvalues of M̂a.
Then problem (4) is equivalent to
min V ∈Rd×k,V >V=Ik
max{J0(V ), J1(V )}, (6a)
where for each (a, a′) ∈ {(0, 1), (1, 0)}, the function Ja is defined as
Ja(V ) = κa + θa √〈 Id − V V >, M̂a 〉 + ϑa′ √〈 Id − V V >, M̂a′ 〉 + 〈 Id − V V >, Ca 〉 , (6b)
and the parameters κ ∈ R, θ ∈ R, ϑ ∈ R and C ∈ Sd+ are defined as
κa = (p̂a + λ)εa + (p̂a′ − λ)εa′ , θa = 2|p̂a + λ| √ εa, ϑa′ = 2|p̂a′ − λ| √ εa′ , Ca = (p̂a + λ)M̂a + (p̂a′ − λ)M̂a′ . (6c)
We now briefly explain the steps that lead to the results in Theorem 3.2. Letting
J0(V ) = sup Q∈B(P̂)
(p̂0 + λ)EQ[`(V,X)|A = 0] + (p̂1 − λ)EQ[`(V,X)|A = 1],
J1(V ) = sup Q∈B(P̂)
(p̂0 − λ)EQ[`(V,X)|A = 0] + (p̂1 + λ)EQ[`(V,X)|A = 1],
then by expanding the term U(V,Q) using its definition, problem (4) becomes min
V ∈Rd×k,V >V=Ik max{J0(V ), J1(V )}.
Leveraging the definition the ambiguity set B(P̂), for any pair (a, a′) ∈ {(0, 1), (1, 0)}, we can decompose Ja into two separate supremum problems as follows
Ja(V ) = sup Qa:W(Qa,P̂a)≤εa (p̂a + λ)EQa [`(V,X)] + sup Qa′ :W(Qa′ ,P̂a′ )≤εa′ (p̂a′ − λ)EQa′ [`(V,X)].
The next proposition asserts that each individual supremum in the above expression admits an analytical expression. Proposition 3.3 (Reformulation). Fix a ∈ A. For any υ ∈ R, εa ∈ R+, it holds that
sup Qa:W(Qa,P̂a)≤εa υEQa [`(V,X)]
=  υ (√〈 Id − V V >, M̂a 〉 + √ εa )2 if υ ≥ 0, υ (√〈 Id − V V >, M̂a 〉 −√εa )2 if υ < 0 and 〈 Id − V V >, M̂a 〉 ≥ εa,
0 if υ < 0 and 〈 Id − V V >, M̂a 〉 < εa.
The proof of Theorem 3.2 now follows by applying Proposition 3.3 to each term in Ja, and balance the parameters to obtain (6c). A detailed proof is relegated to the appendix. In the next section, we study an efficient algorithm to solve (6a). Remark 3.4 (Recovery of the nominal PCA). If λ = 0 and εa = 0 ∀a ∈ A, our formulation (4) becomes the standard PCA problem (2). In this case, our robust fair principal components reduce to the standard principal components. On the contrary, existing fair PCA methods such as Samadi et al. (2018) and Olfat & Aswani (2019) cannot recover the standard principal components.

heading: 4 RIEMANNIAN GRADIENT DESCENT ALGORITHM
text: The distributionally robust fairness-aware PCA problem (4) is originally an infinite-dimensional min-max problem. Indeed, the inner maximization problem in (4) optimizes over the space of probability measures. Thanks to Theorem 3.2, it is reduced to the simpler finite-dimensional minimax problem (6a), where the inner problem is only a maximization over two points. Problem (6a) is, however, still challenging as it is a non-convex optimization problem over a non-convex feasible region defined by the orthogonality constraint V >V = Id. The purpose of this section is to devise an efficient algorithm for solving problem (6a) to local optimality based on Riemannian optimization.

heading: 4.1 REPARAMETRIZATION
text: As mentioned above, the non-convexity of problem (6a) comes from both the objective function and the feasible region. It turns out that we can get rid of the non-convexity of the objective function via a simple change of variables. To see that, we let U ∈ Rd×(d−k) be an orthonormal matrix complement to V , that is, U and V satisfy UU> + V V > = Id. Thus, we can express the objective function J via J(V ) = F (U) , max{F0(U), F1(U)}, where for (a, a′) ∈ {(0, 1), (1, 0)}, the function Fa is defined as
Fa(U) , κa + θa √〈 UU>, M̂a 〉 + ϑa′ √〈 UU>, M̂a′ 〉 + 〈 UU>, Ca 〉 .
Moreover, lettingM , {U ∈ Rd×(d−k) : U>U = Id−k}, we can re-express problem (6a) as min U∈M F (U). (7)
The feasible regionM of problem (7) is a Riemannian manifold, called the Stiefel manifold (Absil et al., 2007, Section 3.3.2). It is then natural to solve problem (7) by using Riemannian optimization algorithms (Absil et al., 2007). In fact, problem (6a) itself (before the change of variables) can also be seen as a Riemannian optimization problem over another Stiefel manifold. The change of variables above might seem unnecessary. Nonetheless, the upshot of problem (7) is that the objective function F is convex (in the traditional sense). This faciliates the application of the theoretical and algorithmic framework developed in Li et al. (2019) for (weakly) convex optimization over Stiefel manifolds.

heading: 4.2 THE RIEMANNIAN SUBGRADIENT
text: Note that the objective function F is non-smooth since it is defined as the maximum of two functions F0 and F1. To apply the framework in Li et al. (2019), we need to compute the Riemannian subgradient of the objective function F . Since the Stiefel manifoldM is an embedded manifold in Euclidean space, the Riemannian subgradient of F at any point U ∈ M is given by the orthogonal projection of the usual Euclidean subgradient onto the tangent space of the manifoldM at the point U , see Absil et al. (2007, Section 3.6.1) for example. Lemma 4.1. For any point U ∈ M, let3 aU ∈ arg maxa∈{0,1} Fa(U) and a′U = 1 − aU . Then, a Riemannian subgradient of the objective function F at the point U is given by
gradF (U) = (Id − UU>)  θaU√〈 UU>, M̂aU 〉M̂aUU + ϑa′U√〈 UU>, M̂a′U 〉M̂a′UU + 2CaUU  .
3 It is possible that the maximizer is not unique. In that case, choosing aU to be either 0 or 1 would work.

heading: 4.3 RETRACTIONS
text: Another important instrument required by the framework in Li et al. (2019) is a retraction of the Stiefel manifold M. At each iteration, the point U − γ∆ obtained by moving from the current iterate U in the opposite direction of the Riemannian gradient ∆ may not lie on the manifold in general, where γ > 0 is the stepsize. In Riemannian optimization, this is circumvented by the concept of retraction. Given a point U ∈ M on the manifold, the Riemannian gradient ∆ ∈ TUM (which must lie in the tangent space TUM) and a stepsize γ, the retraction map Rtr defines a point RtrU (−γ∆) which is guaranteed to lie on the manifoldM. Roughly speaking, the retraction RtrU ( · ) approximates the geodesic curve through U along the input tangential direction. For a formal definition of retractions, we refer the readers to (Absil et al., 2007, Section 4.1). In this paper, we focus on the following two commonly used retractions for Stiefel manifolds. The first one is the QR decomposition-based retraction
RtrqfU (∆) = qf(U + ∆), U ∈M,∆ ∈ TUM, where qf( · ) is the Q-factor in the QR decomposition. The second one is the polar decompositionbased retraction
RtrpolarU (∆) = (U + ∆)(Id−k + ∆ >∆)− 1 2 , U ∈M,∆ ∈ TUM. (8)

heading: 4.4 ALGORITHM AND CONVERGENCE GUARANTEES
text: Associated with any choice of retraction Rtr is a concrete instantiation of the Riemannian subgradient descent algorithm for our problem (7), which is presented in Algorithm 1.
Algorithm 1 Riemannian Subgradient Descent for (7) 1: Input: An initial point U0, a number of iterations τ and a retraction Rtr : (U,∆) 7→ RtrU (∆). 2: for t = 0, 1, . . . , τ − 1, do 3: Find at , arg maxa∈{0,1}{Fa(Ut)}. 4: Compute the Riemannian subgradient ∆t = gradF (Ut) using the formula
∆t = (I − UtUt>)  θat√〈 UtU>t , M̂at 〉M̂atUt + ϑa′t√〈 UtU>t , M̂a′t 〉M̂a′tUt + 2CatUt  .
5: Set Ut+1 = RtrUt(−γt∆t), where the step-size γt ≡ 1√τ+1 is constant. 6: end for 7: Output: Uτ .
The specific choice of the stepsizes γt is motivated by the theoretical results of (Li et al., 2019).
We now study the convergence guarantee of Algorithm 1. The following lemma shows that the objective function F is Lipschitz continuous (with respect to the Riemannian metric on the Stiefel manifoldM) with an explicit Lipschitz constant L. Lemma 4.2 (Lipschitz continuity). The function F is L-Lipschitz continuous onM, where L > 0 is given by
L , max { θ0
σmax(M̂0)√ σmin(M̂0) , θ1 σmax(M̂1)√ σmin(M̂1) , ϑ0 σmax(M̂0)√ σmin(M̂0) , ϑ1 σmax(M̂1)√ σmin(M̂1) ,
2 √ d− kσmax(C0), 2 √ d− kσmax(C1) } .
(9)
We now proceed to show that Algorithm 1 enjoys a sub-linear convergence rate. To state the result, we define the Moreau envelope
Fµ(U) , min U ′∈M
{ F (U ′) + 1
2µ ‖U ′ − U‖2F
} ,
where ‖ · ‖F denotes the Frobenius norm of a matrix. Also, to measure the progress of the algorithm, we need to introduce the proximal mapping on the Stiefel manifold (Li et al., 2019):
proxµF (U) ∈ arg min U ′∈M
{ F (U ′) + 1
2µ ‖U ′ − U‖2F
} .
From Li et al. (2019, Equation (22)), we have that ‖gradF (U)‖F ≤ ∥∥proxµF (U)− U∥∥F
µ , gapµ(U).
Therefore, the number gapµ(U) is a good candidate to quantify the progress of optimization algorithms for solving problem (7). Theorem 4.3 (Convergence guarantee). Let {Ut}t=1,...,τ be the sequence of iterates generated by Algorithm 1. Suppose that µ = 1/4L, where L is the Lipschitz constant of F in (9). Then, we have
min t=0,...,τ
gapµ(Ut) ≤ 2 √ Fµ(U0)−minU Fµ(U) + 2L3(L+ 1)
(τ + 1)1/4 .

heading: 5 NUMERICAL EXPERIMENTS
text: We compare our proposed method, denoted RFPCA, against two state-of-the-art methods for fair PCA: 1) FairPCA from Samadi et al. (2018)4, and 2) CFPCA from Olfat & Aswani (2019)5 with both cases: only mean constraint, and both mean and covariance constraints. We consider a wide variety of datasets with ranging sample sizes and number of features. Further details about the datatasets can be found in Appendix C. The code for all experiments is available in supplementary materials.
We include here some details about the hyper-parameters that we search in the cross-validation steps.
• RFPCA. We notice that the neighborhood size εa should be inversely proportional to the size of subgroup a. Indeed, a subgroup with large sample size is likely to have more reliable estimate of the moment information. Then we parameterize the neighborhood size εa by a common scalar α, and we have εa = α/ √ Na, where Na is the number of samples in group a. We search
α ∈ {0.05, 0.1, 0.15} and λ ∈ {0., 0.5, 1., 1.5, 2.0, 2.5}. For better convergence quality, we set the number of iteration for our subgradient descent algorithm to τ = 1000 and also repeat the Riemannian descent for 20 randomly generated initial point U0.
• FairPCA. According to Samadi et al. (2018), we only need tens of iterations for the multiplicative weight algorithm to provide good-quality solution; however, to ensure a fair comparison, we set the number of iterations to be 1000 for the convergence guarantee. We search the learning rate η of the algorithm from set of 17 values evenly spaced in [0.25, 4.25] and {0.1}.
• CFPCA. Followed Olfat & Aswani (2019), for the mean-constrained version of CFPCA, we search δ from {0., 0.1, 0.3, 0.5, 0.6, 0.7, 0.8, 0.9}, and for both mean and covariance constrained version, we fix δ = 0 while searching µ in {0.0001, 0.001, 0.01, 0.05, 0.5}.
Trade-offs. First, we examine the trade-off between the total reconstruction error and the gap between the subgroup error. In this experiment, we only compare our model with FairPCA and CFPCA mean-constraint version. We plot a pareto curve for each of them over the two criteria with different hyper-parameters (hyper-parameters test range are mentioned above). The whole datasets are used for training and evaluation. The results averaged over 5 runs are shown in Figure 2.
In testing methods with different principal components, we first split each dataset into training set and test set with equal size (50% each), the projection matrix of each method is learned from training set and tested over both sets. In this case, we only compare our method with traditional PCA and FairPCAmethod. We fix one set hyper-parameters for each method. For FairPCA, we set η = 0.1 and for RFPCA we set α = 0.15, λ = 0.5, others hyper-parameters are kept as discussed before. The results are averaged over 5 different splits. Figure 3 shows the consistence of our method performing fair projections over different values of k. Our method (cross) exhibits smaller gap of subgroup errors. More results and discussions on the effect of ε can be found in Appendix D.2. 4 https://github.com/samirasamadi/Fair-PCA 5 https://github.com/molfat66/FairML
Cross-validations. Next, we report the performance of all methods based on three criteria: absolute difference between average reconstruction error between groups (ABDiff.), average reconstruction error of all data (ARE.), and the fairness criterion defined by Olfat & Aswani (2019) with respect to a linear SVM’s classifier family (4FLin) 6. Due to the space constraint, we only include the first two criteria in the main text, see Appendix 4 for full results. To emphasize the generalization capacity of each algorithm, we split each dataset into a training set and a test set with ratio of 30% − 70% respectively, and only extract top three principal components from the training set. We find the best hyper-parameters by 3-fold cross validation, and prioritize the one giving minimum value of the summation (ABDiff.+ARE.). The results are averaged over 10 different training-testing splits. We report the performance on both training set (In-sample data) and test set (Out-of-sample data). The details results for Out-of-sample data is given in Table 1, more details about settings and performance can be found at Appendix D.
Results. Our proposed RFPCAmethod outperforms on 11 out of 15 datasets in terms of the subgroup error gap ABDiff, and 9 out of 15 with the totall error ARE. criterion. There are 5 datasets that RFPCA gives the best results for both criteria, and for the remaining datasets, RFPCA has small performance gaps compared with the best method.

heading: A PROOFS
text: A.1 PROOFS OF SECTION 2
Proof of Proposition 2.2. We first prove the “only if” direction. Suppose that there exists a fair projection matrix V ∈ Mk relative to Q. Let U ∈ Md−k be a complement matrix of V . Then, Definition 2.1 can be rewritten as 〈UU>, S〉 = 0, which implies that the null space of S has a dimension at least d − k. By the rank-nullity duality, we have rank(S) ≤ k. Next, we prove the “if” direction. Suppose that rank(S) ≤ k. Then, the matrix S has at least d− k (repeated) zero eigenvalues. Let U ∈Md−k be an orthonormal matrix whose columns are any d−k eigenvectors corresponding to the zero eigenvalues of S and V ∈ Mk be a complement matrix of U . Then, 〈Id − V V >, S〉 = 〈UU>, S〉 = 0. Therefore, V is a fair projection matrix relative to Q. This completes the proof.
A.2 PROOF OF SECTION 3
Proofs of Proposition 3.3. By exploiting the definition of the loss function `, we find
sup Qa:W(Qa,P̂a)≤εa υEQa [`(V,X)]
=  sup µa,Σa tr ( υ(I − V V >)(Σa + µaµ>a ) ) s.t. ‖µa − µ̂a‖22 + tr ( Σa + Σ̂a − 2 ( Σ̂ 1 2 aΣaΣ̂ 1 2 a ) 1 2 ) ≤ εa
=  inf γ(εa − tr ( Σ̂a ) ) + γ2 tr ( (γI − υ(I − V V >))−1Σ̂a ) + τ s.t. [ γI − υ(I − V V >) γµ̂a
γµ̂>a γ‖µ̂a‖22 + τ
] 0, γI υ(I − V V >), γ ≥ 0,
where the last equality follows from Nguyen (2019, Lemma 3.22). By the Woodbury matrix inversion, we have
(γI − υ(I − V V >))−1 = γ−1I − υ γ(υ − γ) (I − V V >).
Moreover, using the Schur complement, the semidefinite constraint is equivalent to
γ‖µ̂a‖22 + τ ≥ γ2µ̂>a (γI − υ(I − V V >))−1µ̂a,
which implies that at optimality, we have
τ = υγ
γ − υ µ̂>a (I − V V >)µ̂a.
At the same time, the constraint γI υ(I − V V >) is equivalent to γ > υ. Combining all previous equations, we have
sup Qa:W(Qa,P̂a)≤εa υEQa [`(V,X)] = inf γ>max{0,υ}
γεa + γυ γ − υ 〈 Id − V V >, M̂a 〉 .
The dual optimal solution γ? is given by
γ? =  υ ( 1 + √〈 Id−V V >,M̂a 〉 εa ) if υ ≥ 0, υ ( 1− √〈 Id−V V >,M̂a 〉 εa ) if υ < 0 and 〈 Id − V V >, M̂a 〉 ≥ εa,
0 if υ < 0 and 〈 Id − V V >, M̂a 〉 < εa.
Note that γ? ≥ max{0, υ} in all the cases. Therefore, we have
sup Qa:W(Qa,P̂a)≤εa υEQa [`(V,X)]
=  υ ( √ εa + √〈 Id − V V >, M̂a 〉)2 if υ ≥ 0, υ ( √ εa − √〈 Id − V V >, M̂a 〉)2 if υ < 0 and 〈 Id − V V >, M̂a 〉 ≥ εa,
0 if υ < 0 and 〈 Id − V V >, M̂a 〉 < εa.
This completes the proof.
We are now ready to prove Theorem 3.2.
Proof of Theorem 3.2. By expanding the absolute value, problem (4) is equivalent to
min V ∈Rd×k,V >V=Ik
max{J0(V ), J1(V )},
where for each (a, a′) ∈ {(0, 1), (1, 0)}, we can re-express Ja as
Ja(V ) = sup Qa:W(Qa,P̂a)≤εa (p̂a + λ)EQa [`(V,X)] + sup Qa′ :W(Qa′ ,P̂a′ )≤εa′ (p̂a′ − λ)EQa′ [`(V,X)]
Using Proposition 3.3 to reformulate the two individual supremum problems, we have Ja(V ) = (p̂a + λ)εa + 2|p̂a + λ| √ εa 〈 Id − V V >, M̂a 〉 + (p̂a + λ) 〈 Id − V V >, M̂a 〉 + (p̂a′ − λ)εa′ + 2|p̂a′ − λ| √ εa′ 〈 Id − V V >, M̂a′ 〉 + (p̂a′ − λ) 〈 Id − V V >, M̂a′ 〉 .
By defining the necessary parameters κ, θ, ϑ and C as in the statement of the theorem, we arrive at the postulated result.
A.3 PROOFS OF SECTION 4
Proof of Lemma 4.1. Let aU ∈ arg maxa∈{0,1} Fa(U) and a′U = 1 − aU . Then, an Euclidean subgradient of F is given by
∇F (U) = θaU√〈 UU>, M̂aU 〉M̂aUU + ϑa′U√〈 UU>, M̂a′U 〉M̂a′UU + 2CaUU ∈ Rd×(d−k). The tangent space of the Stiefel manifoldM at U is given by
TUM = {∆ ∈ Rd×(d−k) : ∆>U + U>∆ = 0},
whose orthogonal projection (Absil et al., 2007, Example 3.6.2) can be computed explicitly via
ProjTUM(D) = (Id − UU >)D +
1 2 U(U>D −D>U), D ∈ Rd×(d−k).
Therefore, a Riemannian subgradient of F at any point U ∈M is given by
gradF (U) = ProjTUM(∇F (U))
= (Id − UU>)  θaU√〈 UU>, M̂aU 〉M̂aUU + ϑa′U√〈 UU>, M̂a′U 〉M̂a′UU + 2CaUU  .
In the last line, we have used the fact that, if D = SU for some symmetric matrix S, then
U>D −D>U = U>SU − U>S>U = 0.
This completes the proof.
The proof of Lemma 4.2 relies on the following preliminary result.
Lemma A.1. Let M ∈ R(d−k)×(d−k) be a positive definite matrix. Then,∣∣〈UU>,M〉− 〈U ′U ′>,M〉∣∣ ≤ 2√d− kσmax(M)‖U − U ′‖F ∀U,U ′ ∈M, (10) and ∣∣∣∣√〈UU>,M〉−√〈U ′U ′>,M〉∣∣∣∣ ≤ σmax(M)√σmin(M)‖U − U ′‖F ∀U,U ′ ∈M, (11) where σmax(M) and σmin(M) denote the maximum and minimum eigenvalues of the matrix M .
Proof of Lemma A.1. For inequality (10),∣∣〈UU>,M〉− 〈U ′U ′>,M〉∣∣ ≤ ∣∣〈UU>,M〉− 〈UU ′>,M〉∣∣+ ∣∣〈UU ′>,M〉− 〈U ′U ′>,M〉∣∣ ≤ ∣∣〈U,M(U − U ′)〉∣∣+ ∣∣〈U ′,M(U − U ′)〉∣∣
≤ ‖U‖F ‖M(U − U ′)‖F + ‖U ′‖F ‖M(U − U ′)‖F = 2 √ d− kσmax‖U − U ′‖F .
For inequality (11), we first note that the function x 7→ √ x is 1/(2 √ xmin)-Lipschitz on [xmin,+∞)
and that 〈 UU>,M 〉 ≥ (d− k)σmin(M) ∀U ∈M. Therefore,∣∣∣∣√〈UU>,M〉−√〈U ′U ′>,M〉∣∣∣∣ ≤ 12√(d− k)σmin(M) ∣∣〈UU>,M〉− 〈U ′U ′>,M〉∣∣ ≤ σmax(M)√
σmin(M) ‖U − U ′‖F ,
where the last inequality follows from (10). This completes the proof.
We are now ready to prove Lemma 4.2.
Proof of Lemma 4.2. Let U , U ′ ∈M be two arbitrary points. We have
|F (U)− F (U ′)| = |max {F0(U), F1(U)} −max {F0(U ′), F1(U ′)}| ≤ max a∈{0,1} |Fa(U)− Fa(U ′)|
≤ max a∈{0,1} max θa σmax(M̂a)√ σmin(M̂a) , ϑ1−a σmax(M̂1−a)√ σmin(M̂1−a) , 2 √ d− kσmax(Ca)  ‖U − U ′‖F , where the last inequality follows from the definition of Fa and Lemma A.1. This completes the proof.
Proof of Theorem 4.3. The proof follows from the fact that F is convex on the Euclidean space Rd×(d−k), Lemma 4.2 and Li et al. (2019, Theorem 2) (and the remarks following it).

heading: B EXTENSION TO NON-BINARY SENSITIVE ATTRIBUTES
text: The main paper focuses on the case of a binary sensitive attribute withA = {0, 1}. In this appendix, we extend our approach to the case when the sensitive attribute is non-binary. Concretely, we suppose that the sensitive attribute A can take on any of the m possible values from 1 to m. In other words, the attribute space now becomes A = {1, . . . ,m}. Definition B.1 (Generalized unfairness measure). The generalized unfairness measure is defined as the maximum pairwise unfairness measure, that is,
Umax(V,Q) , max (a,a′)∈A×A
|EQ[`(V,X)|A = a]− EQ[`(V,X)|A = a′]|.
Notice that if A = {0, 1}, then Umax ≡ U recovers the unfairness measure for binary sensitive attribute defined in Section 2.2. We now consider the following generalized fairness-aware PCA problem
min V ∈Rd×k,V >V=Ik sup Q∈B(P̂)
EQ[`(V,X)] + λUmax(V,Q). (12)
Here we recall that the ambiguity set B(P̂) is defined in (5). The next theorem provides the reformulation of (12). Theorem B.2 (Reformulation of non-binary fairness-aware PCA). Suppose that for any a ∈ A, either of the following two conditions holds:
(i) 0 ≤ λ ≤ p̂a, (ii) the empirical second moment matrix M̂a = 1Na ∑ i∈Ia x̂ix̂ > i satisfies ∑d−k j=1 σj(M̂a) ≥ εa,
where σj(M̂a) is the j-th smallest eigenvalues of M̂a.
Then problem (12) is equivalent to
min V ∈Rd×k,V >V=Ik max a 6=a′ {∑ b∈A 2ca,a′,b √ εb〈Id − V V >, M̂b〉+ λ〈Id − V V >, M̂a − M̂a′〉+ λ(εa − εa′) } ,
where the parameter ca,a′,b admits values
ca,a′,b =  p̂a + λ if b = a, |p̂a′ − λ| if b = a′, p̂b otherwise.
Proof of Theorem B.2. For simplicity, we let E(V,Q, b) = EQ[`(V,X)|A = b]. Then, the objective function of problem (12) can be re-written as
sup Q∈B(P̂)
EQ[`(V,X)] + λUmax(V,Q)
= sup Q∈B(P̂) ∑ b∈A p̂bE(V,Q, b) + λmax a6=a′ {E(V,Q, a)− E(V,Q, a′)}
= max a 6=a′  ∑ b6=a,a′ sup W(Qb,P̂b)≤εb p̂bE(V,Qb, b) + sup W(Qa,P̂a)≤εa (p̂a + λ)E(V,Qa, a) + sup W(Qa′ ,P̂a′ )≤εa′ (p̂a′ − λ)E(V,Qa′ , a′)  = max a 6=a′ { ∑ b 6=a,a′ p̂b (√ 〈Id − V V >, M̂b〉+ √ εb )2 + (p̂a + λ) (√ 〈Id − V V >, M̂a〉+ √ εa )2
+ (p̂a′ − λ) (√ 〈Id − V V >, M̂a′〉+ sgn(p̂a′ − λ) √ εa′ )2}
= max a 6=a′ {∑ b∈A p̂b ( 〈Id − V V >, M̂b〉+ εb ) + ∑ b∈A 2ca,a′,b √ εb〈Id − V V >, M̂b〉
+ λ ( 〈Id − V V >, M̂a − M̂a′〉+ εa − εa′ )} ,
where the first equality follows from the definition of Umax(V,Q) and E(V,Q, b), the second from the definition (5) of the ambiguity set B(P̂), the third from Proposition 3.3 and the fourth from the definition of ca,a′,b. Noting that the first sum in the above maximization is independent of a and a′, the proof is completed.
Theorem 12 indicates that if the sensitive attribute admits finite values, then the distributionally robust fairness-aware PCA problem using an Umax unfairness measure can be reformulated as an optimization problem over the Stiefel manifold, where the objective function is a pointwise maximization of finite number of individual functions. It is also easy to see that each individual function can be reparametrized using U , and the Riemannian gradient descent algorithm in Section 4 can be adapted to solve for the optimal solution. The details on the algorithm are omitted.
C INFORMATION ON DATASETS
We summarize here the number of observations, dimensions, and the sensitive attribute of the data sets. For further information about the data sets and pre-processing steps, please refer to Samadi et al. (2018) for Default Credit and Labeled Faces in the Wild (LFW) data sets, and Olfat & Aswani (2019) for others. For each data set, we further remove columns with too small standard deviation
(≤ 1e−5) as they do not significantly affect the results, and ones with too large standard deviation (≥ 1000) which we consider as unreliable features.

heading: D ADDITIONAL RESULTS
text: D.1 DETAIL PERFORMANCES
Table 3 shows the performances of four examined methods with two criteria ABDiff. and ARE.. It is clear that our method achieves the best results over all 14 datasets w.r.t. ABDiff., and 7 datasets on ARE., which is equal to the number of datasets FairPCA out-perform others.
Table 4 complements Table 1 from the main text, from which we can see that two versions of CFPCA out-perform others over all datasets w.r.t.4FLin, which is the criteria they optimize for. Adjustment for LFW dataset. To demonstrate the efficacy of our method on high-dimensional data sets, we also do experiments on a subset of 2000 faces for each of male and female group (4000 in total) from LFW dataset 22, all images are are rescaled to have the size of 24 × 24 (dimensions d = 576). The experiment follows the same procedure in Section 5, with reducing the number of iterations to 500 for both RFPCA and FairPCA and 2-fold cross validation, results are averaged over 10 train-test simulations. Due to the high dimension of the input, the implementation of Olfat & Aswani (2019) fails to return any result.
7 Yeh & Lien (2009) 8 Mansouri (2013) 9 Nakai (1996) 10 Tsanas (2012) 11 Hofmann (1994) 12 mis (1990) 13 mis (1991) 14 Bock (2007) 15 mis (2009) 16 Blair (2013) 17 Srinivasan (1993) 18 Buscema (2010) 19 Yeh & Lien (2009) 20 Cortez (2009) 21 Afifi & Abdelhamed (2019) 22 https://github.com/samirasamadi/Fair-PCA
D.2.1 EFFECTS OF THE AMBIGUITY SET RADIUS
We examine the change of model’s performance with respect to the change of the radius of ambiguity sets. To generate toy data (also used for Figure 1), we use two 2-dimensional Gaussian distributions to represent two groups of a sensitive attribute, A = 0 and A = 1, or groups 0 and 1 for simplicity. The two distributions both have the mean at (0, 0) and covariance matrices for group 0, 1 are(
4.0 0 0 0.2
) and ( 0.2 0.4 0.4 3.0 ) ,
respectively. For the test set, the number of samples is 8000 for group 0 and 4000 for group 1, while for the training set, we have 200 for group 0 and 100 for group 1. We average the results over 100 simulations, for each simulation, test data is fixed, training data is randomly generated with the number of samples mentioned above. The projections are learned on training data and measured on test data by the summation of ARE. and ABDiff. We fixed λ = 0.1, which is not too small for achieving fair projections, and not too large to clearly observe the effects of ε, and we also fixed ε0 for better visualization. Note that we still compute ε1 = α/ √ N1 in which, α is tested with 100 values evenly spaced in [0, 10].
The experiment results are visualized in Figure 4. The result suggests that increasing the ambiguity set radius can improve the overall model’s performance. This justifies the benefit of adding distribu-
tional robustness to the fairness-aware PCA model. After a saturation point, a too large radius can lessen the role of empirical data, and the model prioritizes a more extreme distribution that is far from the target distribution, which causes the reduction in the model’s performance on target data.
D.2.2 PARETO CURVES
Figure 5: Pareto curves on Biodeg dataset (all data) with 3 principal components
Figure 6: Pareto curves on German Credit dataset (all data) with 3 principal components
D.2.3 PERFORMANCE WITH DIFFERENT PRINCIPAL COMPONENTS
Figure 9: Subgroup average error with different k on E. Coli dataset
Figure 10: Subgroup average error with different k on E. Coli dataset
Figure 11: Subgroup average error with different k on Magic dataset
Figure 12: Subgroup average error with different k on Magic dataset
Figure 13: Subgroup average error with different k on Steel dataset
Figure 14: Subgroup average error with different k on Steel dataset

Review summaries:
This paper studies the formulation, reformulation and algorithm for distributionally robust fairness-aware PCA. The reformulation exploits techniques from distributionally robust optimization, and the algorithm is based on reimannian sub-gradient descent. The theory is tested on UCI datasets.

This paper aims to improve fairness and distributional robustness in dimensionality reduction techniques. More specifically, it proposes a regularization term to principal compoenent analysis (PCA) such that the expected reconstruction errors for each groups (conditioned on a binary sensitive variable) are similar, and relaxes constraint on the empirical probability distribution by maximizing over a Wasserstein-type ambiguity set. This is optimized with subgradient descent on the Stiefel manifold, and numerical experiments comparing the proposed method to previous fair PCA are conducted.

This paper considers the fair PCA with distribution ally robust optimization algorithm. The problem is of great interest. The key idea of the proposed method is to consider a penalized loss function with the penalty on the fairness criterion. The distributionally robust optimization and its reformulation is developed for parameter estimation.