1 Introduction↩︎

In recent years, the ML community has witnessed an onslaught of charges that real-world machine learning algorithms have produced “biased” outcomes. The examples come from diverse and impactful domains. Google Photos labeled African Americans as gorillas [1], [2] and returned queries for CEOs with images overwhelmingly male and white [3], searches for African American names caused the display of arrest record advertisements with higher frequency than searches for white names [4], facial recognition has wildly different accuracy for white men than dark-skinned women [5], and recidivism prediction software has labeled low-risk African Americans as high-risk at higher rates than low-risk white people [6].

The community’s work to explain these observations has roughly fallen into either “biased data” or “biased algorithm” bins. In some cases, the training data might under-represent (or over-represent) some group, or have noisier labels for one population than another, or use an imperfect proxy for the prediction label (e.g., using arrest records in lieu of whether a crime was committed). Separately, issues of imbalance and bias might occur due to an algorithm’s behavior, such as focusing on accuracy across the entire distribution rather than guaranteeing similar false positive rates across populations, or by improperly accounting for confirmation bias and feedback loops in data collection. If an algorithm fails to distribute loans or bail to a deserving population, the algorithm won’t receive additional data showing those people would have paid back the loan, but it will continue to receive more data about the populations it (correctly) believed should receive loans or bail.

Many of the proposed solutions to “biased data” problems amount to re-weighting the training set or adding noise to some of the labels; for “biased algorithms”, most work has focused on maximizing accuracy subject to a constraint forbidding (or penalizing) an unfair model. Both of these concerns and approaches have significant merit, but form an incomplete picture of the ML pipeline and where unfairness might be introduced therein. Our work takes another step in fleshing out this picture by analyzing when dimensionality reduction might inadvertently introduce bias. We focus on principal component analysis (henceforth PCA), perhaps the most fundamental dimensionality reduction technique in the sciences [7]–[9]. We show several real-world data sets for which PCA incurs much higher average reconstruction error for one population than another, even when the populations are of similar sizes. Figure 1 shows that PCA on labeled faces in the wild data set (LFW) has higher reconstruction error for women than men even if male and female faces are sampled with equal weight.

This work underlines the importance of considering fairness and bias at every stage of data science, not only in gathering and documenting a data set [10] and in training a model, but also in any interim data processing steps. Many scientific disciplines have adopted PCA as a default preprocessing step, both to avoid the curse of dimensionality and also to do exploratory/explanatory data analysis (projecting the data into a number of dimensions that humans can more easily visualize). The study of human biology, disease, and the development of health interventions all face both aforementioned difficulties, as do numerous economic and financial analysis. In such high-stakes settings, where statistical tools will help in making decisions that affect a diverse set of people, we must take particular care to ensure that we share the benefits of data science with a diverse community.

We also emphasize this work has implications for representational rather than just allocative harms, a distinction drawn by [11] between how people are represented and what goods or opportunities they receive. Showing primates in search results for African Americans is repugnant primarily due to its representing and reaffirming a racist painting of African Americans, not because it directly reduces any one person’s access to a resource. If the default template for a data set begins with running PCA, and PCA does a better job representing men than women, or white people over minorities, the new representation of the data set itself may rightly be considered an unacceptable sketch of the world it aims to describe.

Our work proposes a different linear dimensionality reduction which aims to represent two populations \(A\) and \(B\) with similar fidelity—which we formalize in terms of reconstruction error. Given an \(n\)-dimensional data set and its \(d\)-dimensional approximation, the reconstruction error of the data with respect to its low-dimensional approximation is the sum of squares of distances between the original data points and their approximated points in the \(d\)-dimensional subspace. To eliminate the effect of size of a population, we focus on average reconstruction error over a population. One possible objective for our goal would find a \(d\)-dimensional approximation of the data which minimizes the maximum reconstruction error over the two populations. However, this objective doesn’t avoid grappling with the fact that population \(A\) may perfectly embed into \(d\) dimensions, whereas \(B\) might require many more dimensions to have low reconstruction error. In such cases, this objective would not necessarily favor a solution with average reconstruction error of \(\epsilon\) for \(A\) and \(y\gg\epsilon\) for \(B\) over one with \(y\) error for \(A\) and \(y\) error for \(B\). This holds even if \(B\) requires \(y\) reconstruction error to be embedded into \(d\) dimensions and thus the first solution is nearly optimal for both populations in \(d\) dimensions.

This motivates our focus on finding a projection which minimizes the maximum additional or marginal reconstruction error for each population above the optimal \(n\) into \(d\) projection for that population alone. This quantity captures how much a population’s reconstruction error increases by including another population in the dimensionality reduction optimization. Despite this computational problem appearing more difficult than solving “vanilla” PCA, we introduce a polynomial-time algorithm which finds an \(n\) into \((d+1)\)-dimensional embedding with objective value better than any \(d\)-dimensional embedding. Furthermore, we show that optimal solutions have equal additional average error for populations \(A\) and \(B\).

1.0.0.1 Summary of our results

We show PCA can overemphasize the reconstruction error for one population over another (equally sized) population, and we should therefore think carefully about dimensionality reduction in domains where we care about fair treatment of different populations. We propose a new dimensionality reduction problem which focuses on representing \(A\) and \(B\) with similar additional error over projecting \(A\) or \(B\) individually. We give a polynomial-time algorithm which finds near-optimal solutions to this problem. Our algorithm relies on solving a semidefinite program (SDP), which can be prohibitively slow for practical applications. We note that it is possible to (approximately) solve an SDP with a much faster multiplicative-weights style algorithm, whose running time in practice is equivalent to solving standard PCA at most 10-15 times. The details of the algorithm are given in the full version of this work. We then evaluate the empirical performance of this algorithm on several human-centric data sets.

2 Related work↩︎

This work contributes to the area of fairness for machine learning models, algorithms, and data representations. One interpretation of our work is that we suggest using Fair PCA, rather than PCA, when creating a lower-dimensional representation of a data set for further analysis. Both pieces of work which are most relevant to our work take the posture of explicitly trying to reduce the correlation between a sensitive attribute (such as race or gender) and the new representation of the data. The first piece is a broad line of work  [12]–[16] that aims to design representations which will be conditionally independent of the protected attribute, while retaining as much information as possible (and particularly task-relevant information for some fixed classification task). The second piece is the work by [17], who also look to design PCA-like maps which reduce the projected data’s dependence on a sensitive attribute. Our work has a qualitatively different goal: we aim not to hide a sensitive attribute, but instead to maintain as much information about each population after projecting the data. In other words, we look for representation with similar richness for population \(A\) as \(B\), rather than making \(A\) and \(B\) indistinguishable.

Other work has developed techniques to obfuscate a sensitive attribute directly [18]–[29]. This line of work diverges from ours in two ways. First, these works focus on representations which obfuscate the sensitive attribute rather than a representation with high fidelity regardless of the sensitive attribute. Second, most of these works do not give formal guarantees on how much an objective will degrade after their transformations. Our work directly minimizes the amount by which each group’s marginal reconstruction error increases.

Much of the other work on fairness for learning algorithms focuses on fairness in classification or scoring [30]–[33], or online learning settings [34]–[37]. These works focus on either statistical parity of the decision rule, or equality of false positives or negatives, or an algorithm with a fair decision rule. All of these notions are driven by a single learning task rather than a generic transformation of a data set, while our work focuses on a ubiquitous, task-agnostic preprocessing step.

3 Notation and vanilla PCA↩︎

We are given \(n\)-dimensional data points represented as rows of matrix \(M \in \mathbb{R}^{m \times n}\). We will refer to the set and matrix representation interchangeably. The data consists of two subpopulations \(A\) and \(B\) corresponding to two groups with different value of a binary sensitive attribute (e.g., males and females). We denote by \(\left[ \begin{array}{c} A \\ B \end{array} \right]\) the concatenation of two matrices \(A,B\) by row. We refer to the \(i^{th}\) row of \(M\) as \(M_i\), the \(j^{th}\) column of \(M\) as \(M^j\) and the \((i,j)^{th}\) element of \(M\) as \(M_{ij}\). We denote the Frobenius norm of matrix \(M\) by \(\|M\|_F\) and the \(2\)-norm of the vector \(M_i\) by \(\|M_i\|\). For \(k\in \mathbb{N}\), we write \([k]:=\{1, \ldots, k\}\). \(|A|\) denotes the size of a set \(A\). Given two matrices \(M\) and \(N\) of the same size, the Frobenius inner product of these matrices is defined as \(\langle M,N \rangle= \sum_{ij} M_{ij}N_{ij} = \mathop{\mathrm{Tr}}(M^TN)\).

3.1 PCA↩︎

This section recalls useful facts about PCA that we use in later sections. We begin with a reminder of the definition of the PCA problem in terms of minimizing the reconstruction error of a data set.

We will refer to \(\widehat{M}\) as an optimal rank-\(d\) approximation of \(M\). The following well-known fact characterizes the solutions to this classic problem [e.g., [38]].

The matrix \(WW^T\in \mathbb{R}^{n\times n}\) is called a projection matrix.

4 Fair PCA↩︎

Given the \(n\)-dimensional data with two subgroups \(A\) and \(B\), let \(\widehat{M}, \widehat{A}, \widehat{B}\) be optimal rank-\(d\) PCA approximations for \(M, A,\) and \(B\), respectively. We introduce our approach to fair dimensionality reduction by giving two compelling examples of settings where dimensionality reduction inherently makes a tradeoff between groups \(A\) and \(B\). Figure ¿fig:fig:ex1? shows a setting where projecting onto any single dimension either favors \(A\) or \(B\) (or incurs significant reconstruction error for both), while either group separately would have a high-fidelity embedding into a single dimension. This example suggests any projection will necessarily make a trade off between error on \(A\) and error on \(B\).

Our second example (shown in Figure ¿fig:fig:ex2?) exhibits a setting where \(A\) and \(B\) suffer very different reconstruction error when projected onto one dimension: \(A\) has high reconstruction error for every projection while \(B\) has a perfect representation in the horizontal direction. Thus, asking for a projection which minimizes the maximum reconstruction error for groups \(A\) and \(B\) might require incurring additional error for \(B\) while not improving the error for \(A\). So, minimizing the maximum reconstruction error over \(A\) and \(B\) fails to account for the fact that two populations might have wildly different representation error when embedded into \(d\) dimensions. Optimal solutions to such objective might behave in a counterintuitive way, preferring to exactly optimize for the group with larger inherent representation error rather than approximately optimizing for both groups simultaneously. We find this behaviour undesirable—it requires sacrifice in quality for one group for no improvement for the other group.

We therefore turn to finding a projection which minimizes the maximum deviation of each group from its optimal projection. This optimization asks that \(A\) and \(B\) suffer a similar loss for being projected together into \(d\) dimensions compared to their individually optimal projections. We now introduce our notation for measuring a group’s loss when being projected to \(Z\) rather than to its optimal \(d\)-dimensional representation:

Then, the optimization that we study asks to minimize the maximum loss suffered by any group. This captures the idea that, fixing a feasible solution, the objective will only improve if it improves the loss for the group whose current representation is worse. Furthermore, considering the reconstruction loss and not the reconstruction error prevents the optimization from incurring error for one subpopulation without improving the error for the other one as described in Figure ¿fig:fig:ex2?.

This definition does not appear to have a closed-form solution (unlike vanilla PCA—see Fact 1). To take a step in characterizing solutions to this optimization, Theorem 2 states that a fair PCA low dimensional approximation of the data results in the same loss for both groups.

Before proving Theorem 2, we need to state some building blocks of the proof, Lemmas 1, 2, and 3. For the proofs of the lemmas please refer to the appendix 9.

The next lemma presents some equalities that we will use frequently in the proofs.

Let the function \(g_A=g_A(U)\) measure the reconstruction error of a fixed matrix \(A\) with respect to its orthogonal projection to the input subspace \(U\). The next lemma shows that the value of the function \(g_A\) at any local minimum is the same.

:

Consider the functions \(g_A\) and \(g_B\) defined in Lemma 3. It follows from Lemma 1 and Lemma 2 that for \(V\in \mathbb{R}^{n\times d}\) with \(V^TV=I\) we have \[\begin{align} \label{eq:lossEq} loss(A, AVV^T) = \|\widehat{A}\|_F^2 - \|A\|_F^2 + g_A(V), \\ loss(B, BVV^T) = \|\widehat{B}\|_F^2 - \|B\|_F^2 + g_B(V). \nonumber \end{align}\tag{1}\]

Therefore, the Fair PCA problem is equivalent to \[\begin{align} \min_{V \in \mathbb{R}^{n\times d}, V^TV=I} \, f(V):=\max \left\{ \frac{1}{|A|} loss(A, AVV^T) , \frac{1}{|B|} loss(B, BVV^T) \right\}. \end{align}\]

We proceed to prove the claim by contradiction. Let \(W\) be a global minimum of \(f\) and assume that \[\label{eq:bigger} \frac{1}{|A|}loss(A, AWW^T) > \frac{1}{|B|}loss(B, BWW^T).\tag{2}\]

Hence, since \(loss\) is continuous, for any matrix \(W_\epsilon\) with \(W_\epsilon ^T W_\epsilon =I\) in a small enough neighborhood of \(W\), \(f(W_\epsilon)= \frac{1}{|A|}loss(A, AW_\epsilon W_\epsilon ^T)\). Since \(W\) is a global minimum of \(f\), it is a local minimum of \(\frac{1}{|A|} loss(A, AWW^T)\) or equivalently a local minimum of \(g_A\) because of 1 .

Let \(\{v_1, \ldots, v_n\}\) be an orthonormal basis of the eigenvectors of \(A^TA\) corresponding to eigenvalues \(\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_n\). Let \(V^*\) be the subspace spanned by \(\{v_1,\ldots,v_d\}\). Note that \(loss(A, A{V^*}^TV^*)=0\). Since the loss is always non-negative for both \(A\) and \(B\), 2 implies that \(loss(A, AWW^T) > 0\). Therefore, \(W\neq V^*\) and \(g_A(V^*) < g_A(W)\). By Lemma 3, this is in contradiction with \(V^*\) being a global minimum and \(W\) being a local minimum of \(g_A\). \(\hfill \Box\)

5 Algorithm and analysis↩︎

In this section, we present a polynomial-time algorithm for solving the fair PCA problem. Our algorithm outputs a matrix of rank at most \(d+1\) and guarantees that it achieves the fair PCA objective value equal to the optimal \(d\)-dimensional fair PCA value. The algorithm has two steps: first, relax fair PCA to a semidefinite optimization problem and solve the SDP; second, solve an LP designed to reduce the rank of said solution. We argue using properties of extreme point solutions that the solution must satisfy a number of constraints of the LP with equality, and argue directly that this implies the solution must lie in \(d+1\) or fewer dimensions. We refer the reader to  [40] for basics and applications of this technique in approximation algorithms.

The algorithm to prove Theorem 3 is presented in Algorithm ¿fig:alg:fairPCA?. Using Lemma 2, we can write the semi-definite relaxation of the fair PCA objective (Def. 4) as SDP [sdp]. This semi-definite program can be solved in polynomial time. The system of constraints [eq:two]-[eq:four] is a linear program in the variables \(\lambda_i\) (with the \(u_i\)’s fixed). Therefore, an extreme point solution \((\bar{\lambda},z^*)\) is defined by \(n+1\) equalities, at most three of which can be constraints in [eq:con1]-[eq:three] and the rest (at least \(n-2\) of them) must be from the \(\bar{\lambda}_i = 0\) or \(\bar{\lambda}_i = 1\) for \(i\in[n]\). Given the upper bound of \(d\) on the sum of the \(\bar{\lambda}_i\)’s, this implies that at least \(d-1\) of them are equal to \(1\), i.e., at most two are fractional and add up to \(1\).

5.0.0.1 Case 1.

All the eigenvalues are integral. Therefore, there are \(d\) eigenvalues equal to \(1\). This results in orthogonal projection to \(d\)-dimension.

5.0.0.2 Case 2.

\(n-2\) of eigenvalues are in \(\{0,1\}\) and two eigenvalues \(0 < \bar{\lambda}_d, \bar{\lambda}_{d+1} < 1\). Since we have \(n+1\) tight constraints, this means that both of the first two constraints are tight. Therefore \[\frac{1}{|A|}(\|\widehat{A}\|_F^2 - \sum_{i=1}^{n}\bar{\lambda}_i \langle {A^TA,u_i u_i^T} \rangle ) = \frac{1}{|B|}(\|\widehat{B}\|_F^2 - \sum_{i=1}^{n}\bar{\lambda}_i \langle {B^TB,u_i u_i^T} \rangle)=z^*\leq \hat{z},\] where the inequality is by observing that \((\hat{\lambda},\hat{z})\) is a feasible solution. Note that the loss of group \(A\) given by an affine projection \(P^*= \sum_{j=1}^n \lambda^* u_j u_j^T\) is \[\begin{align} loss(A,AP^*) = \|{A-AP^*}\|_F^2 - \|A- \widehat{A}\|^2_F = \text{Tr}\left((A-AP^*)(A-AP^*)^\top \right) - \|{A}\|_F^2 +\| \widehat{A}\|_F^2 \\ = \text{Tr}\left((A-AP^*)(A-AP^*)^\top \right) - \|{A}\|_F^2 +\| \widehat{A}\|_F^2 = \| \widehat{A}\|_F^2 - 2\text{Tr}(AP^*A^\top)+\text{Tr}(A{P^*}^2A^\top)\\ =\|\widehat{A}\|_F^2 - \sum_{i=1}^{n}(2\lambda_i^* -{\lambda_i^*}^2)\langle {A^TA,u_i u_i^T} \rangle =\|\widehat{A}\|_F^2 - \sum_{i=1}^{n}\bar{\lambda}\langle {A^TA,u_i u_i^T} \rangle , \\ \end{align}\] where the last inequality is by the choice of \(\lambda_j^* = 1-\sqrt{1-\bar{\lambda_j}}\). The same equality holds true for group \(B\). Therefore, \(P^*\) gives the equal loss of \(z^*\leq\hat{z}\) for two groups. The embedding \(x \rightarrow (x\cdot u_1, \ldots, x\cdot u_{d-1}, \sqrt{\lambda_d^*}\, x\cdot u_d, \sqrt{\lambda_{d+1}^*} \, x \cdot u_{d+1})\) corresponds to the affine projection of any point (row) of \(A,B\) defined by the solution \(P^*\).

In both cases, the objective value is at most that of the original fairness objective. \(\hfill \Box\)

The result of Theorem 3 in two groups generalizes to more than two groups as follows. Given \(m\) data points in \(\mathbb{R}^n\) with \(k\) subgroups \(A_1,A_2,\ldots,A_k\), and \(d\leq n\) the desired number of dimensions of projected space, we generalize Definition 4 of fair PCA problem as optimizing \[\begin{align} \label{minimization} \min_{U\in \mathbb{R}^{m\times n} \text{, } \mathsf{rank}(U) \leq d} \;\max_{i\in\{1,\ldots,k\}} & \left\{ \frac{1}{|A_i|} loss(A_i,U_{A_i}))\right\}, \end{align}\tag{3}\] where \(U_{A_i}\) are matrices with rows corresponding to rows of \(U\) for groups \(A_i\).

In contrast to the case of two groups, when there are more than two groups in the data, it is possible that all optimal solutions to fair PCA will not assign the same loss to all groups. However, with \(k-1\) extra dimensions, we can ensure that the loss of each group remains at most the optimal fairness objective in \(d\) dimension. The result of Theorem 4 follows by extending algorithm in Theorem 3 by adding linear constraints to SDP and LP for each extra group. An extreme solution \((\bar{\lambda},z^*)\) of the resulting LP contains at most \(k\) of \(\lambda_i\)’s that are strictly in between 0 and 1. Therefore, the final projection matrix \(P^*\) has rank at most \(d+k-1\).

5.0.0.3 Runtime

We now analyze the runtime of Algorithm ¿fig:alg:fairPCA?, which consists of solving SDP [sdp] and finding an extreme solution to an LP [eq:two]-[eq:four]. The SDP and LP can be solved up to additive error of \(\epsilon>0\) in the objective value in \(O(n^{6.5}\log(1/\epsilon))\) [41] and \(O(n^{3.5}\log(1/\epsilon))\) [42] time, respectively. The running time of SDP dominates the algorithm both in theory and practice, and is too slow for practicaluses for moderate size of \(n\).

We propose another algorithm of solving SDP using the multiplicative weight (MW) update method. In theory, our MW takes \(O(\frac{1}{\epsilon^2})\) iterations of solving standard PCA, giving a total of\(O(\frac{n^3}{\epsilon^2})\) runtime, which may or may not be faster than \(O(n^{6.5}\log(1/\epsilon))\) depending on \(n,\epsilon\). In practice, however, we observe that after appropriately tuning one parameter in MW, the MWalgorithm achieves accuracy \(\epsilon<10^{-5}\) within tens of iterations, and therefore is used to obtain experimental results in this paper. Our MW can handle data of dimension up to a thousand with running time in less than a minute. The details of implementation and analysis of MW method are in Appendix 8.

6 Experiments↩︎

We use two common human-centric data sets for our experiments. The first one is labeled faces in the wild (LFW) [43], the second is the Default Credit data set [44]. We preprocess all data to have its mean at the origin. For the LFW data, we normalized each pixel value by \(\tfrac{1}{255}\). The gender information for LFW was taken from [45], who manually verified the correctness of these labels. For the credit data, since different attributes are measurements of incomparable units, we normalized the variance of each attribute to be equal to 1. The code of all experiments is publicly available at https://github.com/samirasamadi/Fair-PCA.

6.0.0.1 Results

We focus on projections into relatively few dimensions, as those are used ubiquitously in early phases of data exploration. As we already saw in Figure 1 left, at lower dimensions, there is a noticeable gap between PCA’s average reconstruction error for men and women on the LFW data set. This gap is at the scale of up to 10% of the total reconstruction error when we project to 20 dimensions. This still holds when we subsample male and female faces with equal probability from the data set, and so men and women have equal magnitude in the objective function of PCA (Figure 1 right).

Figure 2 shows the average reconstruction error of each population (Male/Female, Higher/Lower education) as the result of running vanilla PCA and Fair PCA on LFW and Credit data. As we expect, as the number of dimensions increase, the average reconstruction error of every population decreases. For LFW, the original data is in 1764 dimensions (42\(\times\)​42 images), therefore, at 20 dimensions we still see a considerable reconstruction error. For the Credit data, we see that at 21 dimensions, the average reconstruction error of both populations reach 0, as this data originally lies in 21 dimensions. In order to see how fair are each of these methods, we need to zoom in further and look at the average loss of populations.

Figure 3 shows the average loss of each population as the result of applying vanilla PCA and Fair PCA on both data sets. Note that at the optimal solution of Fair PCA, the average loss of two populations are the same, therefore we have one line for “Fair loss”. We observe that PCA suffers much higher average loss for female faces than male faces. After running fair PCA, we observe that the average loss for fair PCA is relatively in the middle of the average loss for male and female. So, there is improvement in terms of the female average loss which comes with a cost in terms of male average loss. Similar observation holds for the Credit data set. In this context, it appears there is some cost to optimizing for the less well represented population in terms of the better-represented population.

7 Future work↩︎

This work is far from a complete study of when and how dimensionality reduction might help or hurt the fair treatment of different populations. Several concrete theoretical questions remain using our framework. What is the complexity of optimizing the fairness objective? Is it NP-hard, even for \(d=1\)? Our work naturally extends to \(k\) predefined subgroups rather than just \(2\), where the number of additional dimensions our algorithm uses is \(k-1\). Are these additional dimensions necessary for computational efficiency?

In a broader sense, this work aims to point out another way in which standard ML techniques might introduce unfair treatment of some subpopulation. Further work in this vein will likely prove very enlightening.

Acknowledgements↩︎

This work was supported in part by NSF awards CCF-1563838, CCF-1717349, and CCF-1717947.

8 Improved runtime of semi-definite relaxation by multiplicative weight update method↩︎

In this section, we show the multiplicative weight (MW)algorithm and runtime analysis to solve the fair PCA relaxation in two groups for \(n\times n\) matrix up to \(\epsilon\) additive error in \(O(\frac{1}{\epsilon^2})\) iterations of solving a standard PCA, such as Singular Value Decomposition (SVD). Because SVD takes \(O(n^3)\) time, the SDP relaxation [sdp] for two groups can be solved in \(O(\frac{n^3}{\epsilon^2})\). Comparing to \(O(n^{6.5}\log(1/\epsilon))\) runtime of an SDP solver that is commonly implemented with the interior point method [41], our algorithm may be faster or slower depending on \(n,\epsilon\). In practice, however, we tune the parameter of MW algorithm much more aggressively than in theory, and often take the last iterate solution of MW rather the average when the last iterate performs better, which gives a much faster convergence rate. Our runs of MW show that MW converges in at most 10-20 iterations. Therefore, we use MW to implement our fair PCA algorithm. We note at the conclusion of this section that the algorithm and analysis can be extended to solving fair PCA in \(k\) groups up to additive error \(\epsilon\) in \(O(\frac{\log k}{\epsilon^2})\) iterations.

Technically, the number of iterations for \(k\) groups is \(O(\frac{W^2\log k}{\epsilon^2})\), where \(W\) is the width of the problem, as defined in [46]. \(W\) can usually be bounded by the maximum number of input or the optimal objective value. For our purpose, if the total variance of input data over all dimension is \(L\), then the width \(W\) is at most \(L\). For simplicity, we assume \(L\leq 1\) (e.g. by normalization in prepossessing step), hence obtaining the \(O(\frac{\log k}{\epsilon^2})\) bound on number of iterations.

We first present an algorithmic framework and the corresponding analysis in the next two subsections, and later apply those results to our specific setting of solving the SDP [sdp] from fair PCA problem. The previous work by [46] shows how we may solve a feasibility problem of an LP using MW technique. Our main theoretical contribution is to propose and analyze the optimization counterpart of the feasibility problem, and the MW algorithm we need to solve such problem. The MW we develop fits more seamlessly into our fair PCA setting and simplifies the algorithm to be implemented for solving the SDP [sdp].

8.1 Problem setup and oracle access↩︎

We first formulate the feasibility problem and its optimization counterpart in this section. The previous and new MW algorithms and their analysis are presented in the following Section 8.2.

8.1.1 Previous work: multiplicative weight on feasibility problem↩︎

8.1.1.1 Problem

As in [46], we are given \(A\in \mathbb{R}^{m\times n}\) as an \(m\times n\) real matrix, \(x\in \mathbb{R}^n,b\in\mathbb{R}^m\), and \(\mathcal{P}\) as a convex set in \(\mathbb{R}^n\), and the goal is to check the feasibility problem \[\exists? x\in \mathcal{P}\;: \;Ax \geq b \label{feasibility}\tag{4}\] by giving a feasible \(x\in\mathcal{P}\) or correctly deciding that such \(x\) does not exist.

8.1.1.2 Oracle Access

We assume the existence if an oracle that, given any probability vector \(p\in \Delta_m\) over \(m\) constraints of 4 , correctly answers a single-constraint problem \[\exists? x\in \mathcal{P}\;: \;p^\top Ax \geq p^\top b \label{feasibility-oracle}\tag{5}\] by giving a feasible \(x\in\mathcal{P}\) or correctly deciding that such \(x\) does not exist. We may think of 5 as a weighted version of 4 , with weights on each constraint \(i\in[m]\) being \(p_i\).

As 5 consists of only one constraint, solving 5 is much easier than 4 in many problem settings. For example, in our PCA setting, solving [sdp] directly is non-trivial, but the weighted version 5 is a standard PCA problem: we weight each group \(A,B\) based on \(p\), and then apply a PCA algorithm (Singular Value Decomposition) on the sum of two weighted groups. The solution gives an optimal value of \(p^\top Ax-p^\top\;b\) in 5 . More details of application in fair PCA settings are in Section 8.3

8.1.2 New setting: multiplicative weight on optimization problem↩︎

8.1.2.1 Problem

The previous work gives an MW framework for the feasibility question. Here we propose an optimization framework, which asks for the best \(x\in\mathcal{P}\) rather than an existence of \(x\in\mathcal{P}\). The optimization framework can be formally stated as, given \(A\in \mathbb{R}^{m\times n}\) as an \(m\times n\) real matrix, \(x\in \mathbb{R}^n,b\in\mathbb{R}^m\), and \(\mathcal{P}\) as a convex set in \(\mathbb{R}^n\), we need to solve \[\min z \;: \;Ax-b+z \cdot \mathbf{1}\geq0 , \text{ s.t. }x\in \mathcal{P}\label{optimize}\tag{6}\] where \(\mathbf{1}\) denotes the \(m\times 1\) vector with entries 1. Denote \(z^*\) the optimum of 6 .

With the same type of oracle access, we may run 4 for \(O(\log \frac{n}{\epsilon})\) iterations to do binary search for the correct value of optimum \(z^*\) up to an additive error \(\epsilon\). However, our main contribution is to modify the previous multiplicative weight algorithm and the definition of the oracle to solve 6 without guessing the optimum \(z^*\). This improves the runtime slightly (reduce the \(\log(n/\epsilon)\) factor) and simplifies the algorithm.

8.1.2.2 Feasibility Oracle Access

We assume the existence of an oracle that, given any probability vector \(p\in \Delta_m\) over \(m\) constraints of 6 , correctly answers a single-constraint problem \[\text{Find } x\in \mathcal{P}\;: \;p^\top Ax - p^\top b+z^* \geq 0 \label{optimize-oracle2}\tag{7}\]

There is always such \(x\) because multiplying 6 on the left by \(p^\top\) shows that one of such \(x\) is the optimum \(x^*\) of 6 . However, finding one may not be as trivial as asserting problem’s feasibility. In general, 7 can be tricky to solve since we do not yet know the value of \(z^*\).

8.1.2.3 Optimization Oracle Access

We define the oracle that, given \(p\in \Delta_m\) over \(m\) constraints of 6 , correctly answers one maximizer of \[\min z \;: \;p^\top Ax-p^\top b+z \geq0 , \text{ s.t. }x\in \mathcal{P}\label{optimize-oracle}\tag{8}\] which is stronger than and is sufficient to solve 7 . This is because \(x^*\) of 6 is one feasible \(x\) to 8 , so the optimum \(\hat{z}\) of 8 is at most \(z^*\). Therefore, the optimum \(x\) by 8 can be a feasible solution to 7 . In many setting, because 7 is only one-constraint problem, it is possible to solve the optimization version 8 instead. For example, in our fair PCAon two groups setting, we can solve the 8 by standard PCA on the union of two groups after an appropriate weighting on each group. More details of application in fair PCA settings are in Section 8.3.

8.2 Algorithm and Analysis↩︎

The line of proof follows similarly from [46]. We first state the technical property that the oracle satisfies in our optimization framework, then show how to use that property to bound the number of iterations. We fix \(A\in \mathbb{R}^{m\times n}\) as an \(m\times n\) real matrix, \(x\in \mathbb{R}^n,b\in\mathbb{R}^m\), and \(\mathcal{P}\) is a convex set in \(\mathbb{R}^n\)

Note that even though we do not know \(z^*\), if we know the range of \(A_ix-b_i\) for all \(i\), we can bound the range of \(z^*\). Therefore, we can still find a useful \(\ell,\rho\) that an oracle satisfies.

Now we are ready to state the main result of this section: that we may solve the optimization version by multiplicative update as quickly as solving the feasibility version of the problem.

8.3 Application of multiplicative update method to the fair PCA problem↩︎

In this section, we apply MW results for solving LP to solve the SDP relaxation [sdp] of fair PCA.

8.3.0.1 LP formulation of fair PCA relaxation

The SDP relaxation [sdp] of fair PCA can be written in the form 6 as an LP with two constraints \[\begin{align} \tag{13} \min_{P\in\mathcal{P},z\in\mathbb{R}} &z \text{ s.t. } \\ z&\geq \alpha - \frac{1}{m_1}\langle A^\top A,P\rangle \tag{14} \\ z&\geq \beta -\frac{1}{m_2} \langle B^\top B,P\rangle \tag{15} \end{align}\] for some constants \(\alpha,\beta\), where the feasible region of variables is over a set of PSD matrices: \[\mathcal{P}=\{M\in\mathbb{R}^{n\times n}:0\preceq M \preceq I,\text{tr}(M)\leq d\}\] We will apply the multiplicative weight algorithm to solve 13 15 .

8.3.0.2 Oracle Access

First, we present an the oracle in Algorithm 4, which is in the form 8 and therefore can be used to solve 7 . As defined in 8 , the optimization oracle, given a weight vector \(p=(p_1,p_2)\in\Delta_2\), should be able to solve the LP with one weighted constraint obtained from weighting two constraints 14 and 15 by \(p\). However, because both constraints involve only dot products of same variable \(P\) with constant matrices \(A^\top A\) and \(B^\top B\), which are linear functions, the weighted constraint will involve the dot product of the same variable \(P\) with weighted sum of those constant matrices \(\frac{p_1}{m_1}A^\top A +\frac{p_2}{m_2} B^\top B\).

8.3.0.3 MW Algorithm

Our multiplicative weight update algorithm for solving fair PCA relaxation 13 15 is presented in Algorithm 5. The algorithm follows exactly from the construction in Theorem 5. The runtime analysis of our MW Algorithm 5 follows directly from the same theorem.

8.3.0.4 MW for More Than Two Groups

Algorithms 4 and 5 can be naturally extended to \(k\) groups. Theorem 5 states that we need \(O(\frac{\log k}{\epsilon^2})\) calls to the oracle with additional \(O(k)\) time per call (to update the weight for each loss). In each call, we must compute the weighted sum of \(k\) matrices of dimension \(n \times n\), which takes \(O(kn^2)\) arithmetic operations and perform SVD. In natural settings, \(k\) is much smaller than \(n\), and hence the runtime \(O(n^3)\) of SVD in each oracle call will dominate.

9 Proofs↩︎

: Since \(\mathsf{rank}(U) \leq d\), \(\mathsf{rank}(V) \leq d\) and thus \(\mathsf{rank}(\left[ \begin{array}{c} A \\ B \end{array} \right] VV^T) \leq d\). We will first show that \(loss(A, AVV^T) \leq loss(A, U_A)\).

Step 1. Since \(\{v_1,\ldots,v_d\}\) is an orthonormal basis of row space of \(U\), for every row of \(U_A\), we have that \((U_A)_i = c_iV^T\) for some \(c_i \in \mathbb{R}^{1 \times d}\).

Step 2. We show if we \(c_i \rightarrow A_iV\) and consequently substitute the row \((U_A)_i \rightarrow A_i VV^T\), the value of \(\|A_i - (U_A)_i\|\) decreases. \[\|A_i - (U_A)_i \|^2 = \|A_i - c_i V^T\|^2 = A_i {A_i}^T - 2A_i V {c_i}^T + c_i {c_i}^T\] Here we used the fact that \(V^TV=I\). Minimizing the right hand side with respect to \(c_i\), we get that \(c_i = A_i V\).

Step 3. Step 2 proved that for every \(i\), \(\|A_i - A_i VV^T \|^2 \leq \|A_i - (U_A)_i \|^2\). Remember that

\[\begin{align} loss(A, U_A) = \|A-U_A\|_F^2 - \|A-\widehat{A}\|_F^2 &= \sum \|A_i - (U_A)_i \|^2 - \|A-\widehat{A}\|_F^2 \\ loss(A, AVV^T) = \|A-AVV^T\|_F^2 - \|A-\widehat{A}\|_F^2 &= \sum \|A_i - A_iVV^T \|^2 - \|A-\widehat{A}\|_F^2 \end{align}\] This finished the proof that \(loss(A, AVV^T) \leq loss(A, U_A)\). Similarly, we can see that \(loss(B, BVV^T) \leq loss(B, U_B)\). Therefore \[\begin{align} f(\left[ \begin{array}{c} A \\ B \end{array} \right] VV^T) & = \max \big( \frac{1}{|A|} loss(A,AVV^T), \frac{1}{|B|} loss (B, BVV^T)\big) \\ & \leq \max \big( \frac{1}{|A|} loss(A,U_A), \frac{1}{|B|} loss (B, U_B)\big)\\ &= f(U) \end{align}\] \(\hfill \Box\)

: From Lemma 1, we know that there exist a matrix \(W_A \in \mathbb{R}^{n\times d}\) such that \(W_A^T W_A = I\) and \(\widehat{A} = A W_A W_A^T\). Considering this and the fact that \(V^TV=I\) \[\begin{align} loss(A, AVV^T) &= \|A-AVV^T\|_F^2 - \|A-AW_AW_A^T\|_F^2 \\ &= \sum_{i} \|A_i - A_i VV^T\|^2 - \|A_i - A_i W_A W_A^T\|^2 \\ & = \sum_i A_i A_i^T - A_i VV^T A_i^T - (\sum_i A_i A_i^T - \sum_i A_i W_A W_A^T)\\ & =\sum_{i} A_i W_A W_A^T A_i^T - \sum_{i} A_i VV^T A_i^T \\ \sum_{i} A_i W_A W_A^T A_i^T &= \sum_i \|A_i W_A\|^2 = \|AW_A\|_F^2 = \|AW_A W_A^T\|_F^2 = \|\widehat{A}\|_F^2 \\ \sum_{i} A_i VV^T A_i^T &= \sum_i \|A_i V\|^2 = \|AV\|_F^2 = \sum_i \|Av_i\|^2 \\ \sum_{i} A_i VV^T A_i^T &= \sum_i \|A_i V\|^2 = \|AV\|_F^2 = \mathop{\mathrm{Tr}}(V^TA^TAV) = \mathop{\mathrm{Tr}}(VV^TA^TA) = \langle A^TA, VV^T\rangle \end{align}\] Therefore \(loss(A, AVV^T) = \|\widehat{A}\|_F^2 - \sum_{i=1}^{d} \|Av_i\|^2 = \|\widehat{A}\|_F^2 - \langle A^TA, VV^T\rangle\).

\[\begin{align} \|A - AVV^T\|_F^2 &= \sum_i \|A_i - A_iVV^T\|^2 = \sum_i A_i A_i^T - \sum_i A_iVV^TA_i^T\\ &= \|A\|_F^2 - \sum_i \|Av_i\|^2 = \|A\|_F^2 - \|AV\|_F^2 \\ \end{align}\]

\(\hfill \Box\)

:

We prove that the value of function \(g_A\) at its local minima is equal to its value at its global minimum, which we know is the subspace spanned by a top \(d\) eigenvectors of \(A^TA\). More precisely, we prove the following: Let \(\{v_1, \ldots, v_n\}\) be an orthonormal basis of eigenvectors of \(A^TA\) with corresponding eigenvalues \(\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_n\) where ties are broken arbitrarily. Let \(V^*\) be the subspace spanned by \(\{v_1,\ldots, v_d\}\) and let \(U\) be some \(d\)-dimensional subspace s.t. \(g_A(U) > g_A(V^*)\). There is a continuous path from \(U\) to \(V^*\) s.t. the value of \(g_A\) is monotonically decreasing for every \(d\)-dimensional subspace on the path.

Before starting the proof, we will make a couple of notes which would be used throughout the proof. First note that \(g_A(V)\) is well-defined i.e., the value of \(g_A(V)\) is only a function of the subspace \(V\). More precisely, \(g_A(V)\) is invariant with respect to different choices of orthonormal basis of the subspace \(V\). Second, given Lemma 2, \(g_A(V) = \|A\|_F^2 - \sum_{i}\|Av_i\|^2\). Therefore, proving that \(g_A(V)\) is decreasing is equivalent to proving that \(\sum_i \|Av_i\|^2\) is increasing as a function of any choice of orthonormal basis of the subspaces on the path.

\(g_A(U) > g_A(V^*)\) therefore \(U\neq V^*\). Let \(k\) be the smallest index such that \(v_k \notin U\). Extend \(\{v_1,\ldots, v_{k-1}\}\) to an orthonormal basis of \(U\): \(\{v_1, \ldots, v_{k-1}, v'_k, \ldots, v'_d\}\). Let \(q \geq k\) be the smallest index such that \(\|Av_q\|^2 > \|Av'_q\|^2\). Such an index \(q\) must exist given that \(g_A(U)> g_A(V^*)\). Without loss of generality we can assume that \(q=1\) (this will be clear throughout the proof). Therefore, we assume that \(v_1\), the top eigenvector of \(A^TA\), is not in \(U\) and that it strictly maximizes the function \(\|Au\|^2\) over the space of unit vectors \(u\). Specifically, for any unit vector \(u\in U\), \(\|Au\|^2< \|Av_1\|^2 = \lambda_1\).

Let \(v_1 = \sqrt{1-a^2}z_1 + az_2\) where \(z_1 \in U\) and \(z_2 \perp U\), \(\|z_1\|=\|z_2\|=1\) i.e., the projection of \(v_1\) to \(U\) is \(\sqrt{1-a^2}z_1\). We distinguish two cases:

9.0.0.1 Case \(z_1 = 0\).

\(v_1\perp U\). Let \(w = \sqrt{1-\epsilon^2} u_1 + \epsilon v_1\). \(\|w\|=1\). Note that \(\{w, u_2, \ldots, u_d\}\) is an orthonormal set of vectors. We set \(U_{\epsilon}= span\{w, u_2, \ldots, u_d\}\). We show that \(g_A(U_\epsilon) < g_A (U)\). Using the formulation of \(g\) from Lemma 2, we need to show that \(\|Aw\|^2 + \|Au_2\|^2+\ldots+\|Au_d\|^2 > \|Au_1\|^2+\|Au_2\|^2+\ldots + \|Au_d\|^2\) or equivalently that \(\|Aw\|^2 > \|Au_1\|^2\).

\[\begin{align} \|Aw\|^2 - \|Au_1\|^2 &= \|A(\sqrt{1-\epsilon^2}u_1+\epsilon v_1)\|^2 - \|Au_1\|^2\\ &=(\sqrt{1-\epsilon^2}u_1^T+\epsilon v_1^T)A^TA(\sqrt{1-\epsilon^2}u_1+\epsilon v_1) - \|Au_1\|^2\\ &=(1-\epsilon^2)u_1^TA^TAu_1+\epsilon^2v_1^TA^TAv_1+2\sqrt{1-\epsilon^2}\epsilon u_1^TA^TAv_1 - \|Au_1\|^2\\ &=(1-\epsilon^2)\|Au_1\|^2+\epsilon^2\lambda_1+2\epsilon \sqrt{1-\epsilon^2} u_1^TA^TAv_1 - \|Au_1\|^2 \\ &= \epsilon^2 (\lambda_1 - \|Au_1\|^2) + 2\epsilon \sqrt{1-\epsilon^2} u_1^TA^TAv_1 \\ \end{align}\] where \(u_1^TA^TAv_1 = u_1^T (\lambda_1 v_1)= \lambda_1 u_1^Tv_1 =0\) since \(v_1\) is an eigenvector of \(A^TA\) and \(v_1 \perp u_1\). This, and considering the fact that \(\|Au_1\|^2 < \lambda_1\) \[\begin{align} \|Aw\|^2 - \|Au_1\|^2 &= \epsilon^2(\lambda_1 - \|Au_1\|^2) > 0 \end{align}\] Therefore, \(\|Aw\|^2 > \|Au_1\|^2\) and thus \(g_A(U_\epsilon) < g_A (U)\).

9.0.0.2 Case \(z_1 \neq 0\)

. Note that \(z_2 \neq 0\) either since we picked \(v_1 \notin U\). Let’s extend \(\{z_1\}\) to an orthonormal basis of \(U\): \(\{z_1,u_2, \ldots, u_k\}\). We will transform \(U\) s.t. the resulting subspace \(U_1\) is the span of \(v_1, u_2,\ldots, u_k\). This can then be repeated orthogonal to \(v_1\) till the subspace becomes \(V^*\).

For small enough \(\epsilon > 0\), consider the unit vector \(w = \sqrt{1-\epsilon^2}z_1 + \epsilon z_2\). We will move \(U\) to \(U_\epsilon:=span \{w,u_2,\ldots,u_d\}\). The latter is an orthonormal representation since both \(z_1\) and \(z_2\) are orthogonal to all of \(u_2, \ldots, u_d\) and \(w\) is in the span of \(z_1, z_2\). We will prove that \(g_A(U_\epsilon) < g_A(U)\). Given Lemma 2, since the chosen orthonormal basis of these two subspaces differ only in \(w\) and \(z_1\), it suffices to show that \(\|Aw\|^2 > \|Az_1\|^2\). We can write \[\begin{align} w &= \left(\sqrt{1-\epsilon^2} - \frac{\epsilon\sqrt{1-a^2}}{a}\right)z_1 + \frac{\epsilon}{a}(\sqrt{1-a^2}z_1 + az_2 )=\left(\sqrt{1-\epsilon^2} - \frac{\epsilon\sqrt{1-a^2}}{a}\right)z_1 + \frac{\epsilon}{a}v_1. \end{align}\] Therefore, noting that \(A^TAv_1 = \lambda_1 v_1\) since \(v_1\) is an eigenvector with eigenvalue \(\lambda_1\) and \(z_1^Tv_1=\sqrt{1-a^2}\), \[\begin{align} &\|Aw\|^2 \\ &= \left( \sqrt{1-\epsilon^2} - \frac{\epsilon\sqrt{1-a^2}}{a} \right)^2\|Az_1\|^2 + \frac{\epsilon^2}{a^2} \|Av_1\|^2 + 2\frac{\epsilon}{a}\left(\sqrt{1-\epsilon^2} - \frac{\epsilon\sqrt{1-a^2}}{a}\right)z_1^TA^TAv_1\\ &=\left(1-\epsilon^2 + \frac{\epsilon^2(1-a^2)}{a^2}-2\frac{\epsilon\sqrt{(1-\epsilon^2)(1-a^2)}}{a}\right)\|Az_1\|^2 + \frac{\epsilon^2}{a^2}\lambda_1 \\ & + 2\frac{\epsilon}{a}\left(\sqrt{1-\epsilon^2} - \frac{\epsilon\sqrt{1-a^2}}{a}\right)\lambda_1 z_1^Tv_1 =\left(1-2\epsilon^2 + \frac{\epsilon^2}{a^2}-2\frac{\epsilon\sqrt{(1-\epsilon^2)(1-a^2)}}{a}\right)\|Az_1\|^2 \\ &+ \left(\frac{\epsilon^2}{a^2} + 2\frac{\epsilon\sqrt{(1-\epsilon^2)(1-a^2)}}{a} - 2\frac{\epsilon^2(1-a^2)}{a^2}\right)\lambda_1 \\ &=\|Az_1\|^2 + (\lambda_1-\|Az_1\|^2)\left(2\frac{\epsilon\sqrt{(1-\epsilon^2)(1-a^2)}}{a}+2\epsilon^2 - \frac{\epsilon^2}{a^2}\right)\\ &>\|Az_1\|^2 \end{align}\] where the last inequality follows since \(\lambda_1 > \|Az_1\|^2\) and we can choose \(0 < \epsilon < \frac{1}{1+C}\) for \(C= 4a^2(1-a^2)\) so that \(2\frac{\epsilon\sqrt{(1-\epsilon^2)(1-a^2)}}{a} > \frac{\epsilon^2}{a^2}\) . Thus, \(\|Aw\|^2 > \|Az_1\|^2\) and therefore \(g_A(U_\epsilon) < g_A(U)\). proving the claim.

\(\hfill \Box\)

References↩︎