Fairness of linear regression in decision making

Abstract

Ranking systems conceived on historical data are central to our societies. Given a set of applicants and the information as to whether a past-applicant should have been selected or not, the task of fairly ranking the applicants (either by humans or by computers) is critical to the success of any institution. These tasks are typically carried out using regression methods, and considering the impact of these selection processes on our lives, it is natural to expect various fairness guarantees. In this article, we assume that affirmative action is enforced and that the number of candidates to admit from each protected group is predetermined. We demonstrate that even with this safety-net, classical linear regression methods may increase discrimination in the selection process, reinforcing implicit biases against minorities, in particular by poorly ranking the top minority applicants. We show that this phenomenon is intrinsic to linear regression methods and may happen even if the sensitive attribute is explicitly part of the input, or if a linear regression is computed on each minority group individually. We show that to better rank applicants it might be needed to adapt the choice of the regression methods (linear, polynomial, etc.) to each minority group individually.

Proof of theorem 1

Suppose the linear regression predictor tries to fit the population to the equation $A\mathbf {\beta } = \textbf{y}$, where each row of A corresponds to the non-sensitive attributes of an applicant (i.e., each row of A is a uniformly random vector in A) and the $i^{\text {th}}$ coordinate of y is given by the p th norm of the $i^{\text {th}}$ row of A. Note that $\mathbf {\beta }$ is the vector that minimizes the least square error, $\Vert A\mathbf {\beta } - \textbf{y}\Vert _2$.

First, we show below that $\mathbf {\beta }$ must have (almost) the same entry in all coordinates.

Theorem 2

As $n\rightarrow \infty$, with high probability $1-o(1)$, the optimal vector $\mathbf {\beta }=(b,\ldots ,b)+o(1)$ for some $b\in {\mathbb {R}}$.

Proof

The error $\Vert A\mathbf {\beta } -\textbf{y}\Vert _2^2$ is equal to $\underset{i=1}{\overset{n}{\sum }} |A_i \cdot \mathbf {\beta } -y_i |^2$, where $y_i$ is equal to the (scaled) p th norm of the row vector $A_i$. So we have a sum of N independent identically distributed copies which by large of large numbers tends to $\sim N \cdot \mathbb {E}(|A_i \cdot \mathbf {\beta } -y_i |^2)$. The expected value $\mathbb {E}(|A_i \cdot \mathbf {\beta } -y_i |^2)$ is a quadratic form in $\mathbf {\beta }$ given by:

$$\begin{aligned}{} & {} \mathbb {E}(|A_i \cdot \mathbf {\beta } -y_i |^2) = \mathbb {E}( (\left( \mathbf {\beta }^{T}\cdot A_i^T -y_i\right) \left( A_i\cdot \mathbf {\beta }-y_i\right) ) \\{} & {} = \mathbb {E}(\mathbf {\beta }^{T} A_i^T A_i \mathbf {\beta }) -2\mathbb {E}(y_iA_i\cdot \mathbf {\beta })+\mathbb {E}(y_i^2). \end{aligned}$$

Because the form $Q(\mathbf {\beta }) =\mathbb {E}(|A_i \cdot \mathbf {\beta } -y_i |^2)= \mathbb {E}(\mathbf {\beta }^{T} A_i^T A_i \mathbf {\beta }) -2\mathbb {E}(y_iA_i\cdot \mathbf {\beta })+\mathbb {E}(y_i^2)$ is non-negative, it has a unique minimum $\mathbf {\beta }$ which satisfies gradient condition $\frac{dQ(\mathbf {\beta } +t \textbf{v})}{dt}=0$ for every vector $\textbf{v}$.

$$\begin{aligned}&Q(\mathbf {\beta } + t\textbf{v}) = \mathbb {E}((\mathbf {\beta }+t \textbf{v})^{T} A_i^T A_i ( \mathbf {\beta }+ t\textbf{v}) \\&-2\mathbb {E}(y_iA_i\cdot (\mathbf {\beta } + t\textbf{v}))+\mathbb {E}(y_i^2)\\ =&Q(\mathbf {\beta }) + t \left( \textbf{v}^{T} \mathbb {E}(A_i^T A_i)\mathbf {\beta } + \mathbf {\beta }^T \mathbb {E}(A_i^T A_i) \cdot \textbf{v}- 2\mathbb {E}(y_iA_i)\cdot \textbf{v}\right) \\&+O(t^2). \end{aligned}$$

Therefore, for every $\textbf{v}$, we have

$$\begin{aligned} \textbf{v}^{T} \mathbb {E}(A_i^T A_i)\mathbf {\beta } + \mathbf {\beta }^T \mathbb {E}(A_i^T A_i) \textbf{v} -2\mathbb {E}(y_iA_i)\cdot \textbf{v}=0\\ \Rightarrow 2 \mathbb {E}(A_i^T A_i) \mathbf {\beta } - 2\mathbb {E}(y_iA_i)=\textbf{0}. \end{aligned}$$

Hence, the minimizer $\mathbf {\beta }$ satisfies

$$\begin{aligned} \mathbb {E}(A_i^T A_i) \mathbf {\beta } =\mathbb {E}(y_iA_i). \end{aligned}$$

If $y_i$ is any symmetric function of the coordinates of $A_i$ (like pth norm in this case), we have that $y_iA_i$ is vector with identically distributed coordinates, so the vector $\mathbb {E}(y_iA_i)$ has all equal coordinates.

If we use that $A_i$ is vector with iid coordinates from [0, 1], then $\mathbb {E}(A_i^T A_i)$ is given by a matrix M with the (r, s) entry given by $M_{r,s} = \mathbb {E}(X_rX_s)$, where $X_i$ are iid uniform variables on [0, 1].

So we have $M = \frac{1}{3} I + (\frac{1}{4} J - \frac{1}{4} I) = \frac{1}{12} I + \frac{1}{4}J$, where J is the all ones matrix, and thus, we see that $\mathbf {\beta }$ satisfies $\mathbf {\beta } = (b, b, \cdots , b)$, where $\frac{b}{12} + \frac{db}{4} = \mathbb {E}(y_i X_r)$.

Thus, we have that the minimizer of this quadratic form is of the form $\mathbf {\beta }= (b, b, \cdots , b)$ for some $b\in {\mathbb {R}}$. Therefore, with probability $1-o(1)$, the minimizer $\mathbf {\beta } = (b, b, \cdots , b) + o(1)$ (Because the quadratic form is very close to this expectation quadratic form with high probability, and the minimizer does not change under slight perturbations to the form) $\square$

Therefore, informally, we may conclude that linear regression simply selects the top $50\%$ of the applicants based on their $\ell _1$ norm.

Thus, ranking applicants using regression is equivalent to ranking according to $A_i \cdot \mathbf {\beta }$ which is proportional to the $\sum _{j=1}^d \mathbf {\beta } A_{i}(j)$ which essentially the $\ell _1$ norm of $A_i$. (We are going to rank according to $\mathbf {\beta } =(b, b, \cdots , b) + o(1)$, so the ranking which is determined by the volumes of the region $\mathbf {\beta } \cdot X > \tau$ is essentially equal to the volumes $(b, b, \cdots , b) \cdot X >\tau$- which corresponds to the rank by $\ell _1$ because in this case (non-negative entries) $\Vert X\Vert _1 = (1, 1, \cdots , 1) \cdot X$.)

Let $S_p$ be a ranking of vectors in $[0,1]^{d}$ based on their $\ell _p$-norm. Let $S_p^\tau$ be the restriction of the ranking to the top $\tau$ fraction of applicants. The value of $|S_1^\tau {\setminus } S_p^\tau |$ gives us the recall of the theorem statement, and this is calculated below.

Theorem 3

Let S be a uniformly random sample of N points in $[0,1]^d$. Let $p\in {\mathbb {R}}_{\ge 1}\cup \{\infty \}$. After ranking the points in S by their $\ell _1$-norm (resp. $\ell _p$-norm), let $S_1\subset S$ (resp. $S_p\subset S$) be all points in S ranked in the top half (breaking ties randomly). Then, for large enough d, n we have that

$$\begin{aligned} \frac{ |S_p\cap S_1|}{|S_1|}\sim 1-\frac{\tan ^{-1}\left( \sqrt{\frac{(p+2)^2}{6p+3}-1}\right) }{\pi }. \end{aligned}$$

Proof

By taking n large enough, it is enough to consider the case where you pick a point ${{\textbf {x}}}$ randomly from $[0,1]^d$ and compute the following probability:

$$\begin{aligned} \Pr _{{{\textbf {x}}}\sim [0,1]^d}\left[ \Vert {{\textbf {x}}}\Vert _1\ge m_1\text { and } \Vert {{\textbf {x}}}\Vert _p\ge m_p\right] , \end{aligned}$$

where $m_1$ (resp. $m_p$) is the median of the distribution of $\Vert {{\textbf {x}}}\Vert _1$ (resp. $\Vert x\Vert _p$). Asymptotically, for large d, we have $m_p\sim \root p \of {\frac{d}{p+1}}$. If ${{\textbf {x}}}:=(x_1,\ldots , x_d)$, then the above probability is essentially the following:

$$\begin{aligned} \Pr _{{{\textbf {x}}}\sim [0,1]^d}\left[ x_1+\cdots +x_d\ge m_1\text { and } x_1^p+\cdots x_d^p\ge m_p^p\right] . \end{aligned}$$

For every $i\in [d]$, let $\mathbf {y_i}:=(x_i,x_i^p)$. Then, the probability can be seen as:

$$\begin{aligned} \Pr _{{{\textbf {x}}}\sim [0,1]^d}\left[ \mathbf {y_1}+\cdots +\mathbf {y_d}\in {\mathcal {R}}\right] , \end{aligned}$$

(1)

where ${\mathcal {R}}:=[m_1,\infty )\times [m_p^p,\infty )$. Also, note that for every $i\in [d]$, we have

$$\begin{aligned} \underset{{{{\textbf {x}}}\sim [0,1]^d}}{{\mathbb {E}}}\left[ \mathbf {y_i}\right] = \left( \frac{d}{2},\frac{d}{p+1}\right) \sim (m_1,m_p^p). \end{aligned}$$

Thus, applying central limit theorem to all the $\mathbf {y_i}$s, as $d\rightarrow \infty$, we have:

$$\begin{aligned} \frac{1}{\sqrt{d}}\cdot \underset{i\in [d]}{\sum }\left( \mathbf {y_i}-{\mathbb E}\left[ \mathbf {y_i}\right] \right) \rightarrow {\mathcal {N}}(0,\Sigma ), \end{aligned}$$

where $\Sigma$ is the covariance matrix of $\mathbf {y_i}$s given by ${\mathbb {E}}[\mathbf {y_i}^T\mathbf {y_i}]$. We can thus compute $\Sigma$ to be:

$$\begin{aligned} \Sigma ={\mathbb {E}}\left[ \begin{array}{cc} x_i^2&{}x_i^{p+1}\\ x_i^{p+1}&{}x_i^{2p} \end{array}\right] =\left[ \begin{array}{cc} \nicefrac {1}{3}&{}\nicefrac {1}{p+2}\\ \nicefrac {1}{p+2}&{}\nicefrac {1}{2p+1} \end{array}\right] . \end{aligned}$$

So the probability in (1) converges to

$$\begin{aligned} \Pr _{{\textbf{Y}}\sim {\mathcal {N}}(0,\Sigma )}[{\textbf{Y}}\in [0,\infty )\times [0,\infty )]. \end{aligned}$$

Note that the distribution of ${\mathcal {N}}(0, \varvec{\Sigma })$ is given by

$$\begin{aligned} \displaystyle {\frac{\exp \left( -{\frac{1}{2}}({{\textbf{Y}}^T}{ \varvec{\Sigma } }^{-1}{{\textbf{Y}}})\right) }{ {2\pi \cdot |{ \det {\varvec{\Sigma }}}|}}} \end{aligned}$$

Moreover, we have the inverse of $\Sigma$ is:

$$\begin{aligned} \varvec{\Sigma } ^{-1} = \frac{1}{\frac{1}{3(2p+1)} -\frac{1}{(p+2)^2}}\left[ \begin{array}{cc} \frac{1}{2p+1}&{} -\frac{1}{p+2}\\ -\frac{1}{p+2}&{} \frac{1}{3}\end{array}\right] . \end{aligned}$$

Thus, using the integral

$$\begin{aligned}{} & {} \int _0^{\infty }\int _0^{\infty } \exp (a x^2 + b x y + c y^2) \, dx \, dy = \\{} & {} \frac{1}{2\sqrt{4ac-b^2}} \left( \pi + 2 \arctan \left( \frac{b}{\sqrt{4ac-b^2}}\right) \right) , \end{aligned}$$

we can compute probability to be the expression given in the theorem statement. $\square$