\section{Introduction}

Fix a closed, convex set $\Omega \subseteq \mathbf{R}^d$. For $p \ge 1$, let $\mathcal{P}_p(\Omega)$ denote the set of nonnegative Borel probability measures $\mu$ supported in $\Omega$ with finite $p^{th}$ moment $\int_\Omega |x|^p \d\mu < \infty$. Equip this space with the $p$-Wasserstein distance $W_p$, i.e.
$$
W_p^p(\mu,\nu):=\inf\left\{\int_{\Om\times\Om}|x-y|^p\d\g(x,y):\ \g\in\Pi(\mu,\nu) \right\},
$$
where $\Pi(\mu,\nu):=\left\{\g\in\cP_p(\Om\times\Om):\ (\pi^x)_\sharp\g=\mu,\ (\pi^y)_\sharp\g=\nu\right\}$ denotes the set of transportation plans between $\mu$ and $\nu$ and $\pi^x,\pi^y:\Om\times\Om\to\Om$ stand for the canonical projections $\pi^x(a,b)=a$, $\pi^y(a,b)=b$. We denote by $\Pi_o(\mu,\nu)\subseteq\Pi(\mu,\nu)$ the set of optimal plans that realize the value $W_p(\mu,\nu)$. It is well-known (see for instance \cite{AmbGigSav}) that $(\cP_p(\Om),W_p)$ defines a geodesic metric space. Moreover, if $\Om$ is compact then $W_p$ metrizes the weak-$\star$ convergence of probability measures in $\cP_p(\Om).$

In this paper, we are interested in properties of projection operators ${\rm{P}}_{\Omega}^p:\cP_p(\Om)\to \cK$, where $\cK\subseteq\cP_p(\Om)$ is a given closed and geodesically convex proper subset of $\cP_p(\Om)$. In particular, the main question we are interested in is the so-called {\it nonexpansiveness property} that reads as

\begin{align}\label{question:nonexpansive?}
	{\rm Is\ it\ true\ that\ \ }W_p\left({\rm{P}}_{\Omega}^p[\mu], {\rm{P}}_\Om^p[\nu]\right)\le W_p(\mu,\nu),\ \forall\ \mu,\nu\in\cP_p(\Om)\ ?\tag{Q}
\end{align}

For $\mu\in\cP_p(\Om)$, the projection ${\rm{P}}_\Om^p[\mu]$  is defined as the solution of the variational problem
\begin{equation}\label{def:proj}
{\rm{P}}_{\Omega}^p[\mu]:={\rm{argmin}}\left\{\frac{1}{p}W_p^p(\rho,\mu):\ \rho\in\cK\right\}.
\end{equation}
A few comments on the definition of this operator are necessary. The existence of a solution in this minimization problem is an easy consequence of the direct method of calculus of variations. Indeed, for $\mu\in\cP_p(\Om)$ and $C>0$, the set $\{\rho\in\cP_p(\Om): \ W_p(\rho,\mu)\le C\}$ is tight and the objective functional is weakly lower semicontinuous with respect to the narrow convergence of probability measures. However, for ${\rm{P}}_{\Omega}^p[\mu]$ to be well-defined, we would need to have the uniqueness of a minimizer in \eqref{def:proj}. This turns out to be a subtle question and it is linked to the strict convexity of $\rho\mapsto W_p^p(\rho,\mu)$ and/or the curvature properties of $(\cP_p(\Om),W_p)$.

While $\rho\mapsto W_p^p(\rho,\mu)$ is known to be convex with respect to the `flat' convex combination of probability measures, i.e. along $[0,1]\ni t\mapsto (1-t)\rho_0+ t\rho_1$, its strict convexity typically fails, unless additional conditions are imposed on $\mu$ (for instance absolute continuity with respect to $\cL^d\mres\Om$; see \cite[Proposition 7.17-7.19]{San}). From the geometric viewpoint however, when studying properties of projection operators, it is more natural to consider the notion of geodesic convexity (which is also referred to as {\it displacement convexity} in the case of $(\cP_p(\Om),W_p)$; see \cite{McC, AmbGigSav}). This notion is intimately linked to the curvature properties of the space. By \cite[Section 7.3]{AmbGigSav} we know that when $d\ge 2$, $(\cP_2(\Om),W_2)$ is a positively curved space in the sense of Alexandrov, and so the mapping $\rho\mapsto W_2^2(\rho,\mu)$ in general is not geodesically  $\lambda$-convex, for any $\lambda\in\R$. Similarly, it could be expected that  $(\cP_p(\Om),W_p)$ is also non-negatively curved for $p\neq 2$. However, to the best of our knowledge, a precise result in this direction is not available in the literature at this point.

These considerations let us conclude that the uniqueness of the projection onto closed and geodesically convex sets $\cK$ fails in general. To illustrate this fact, let us consider the following example. Let $\Om=\R^2$ and let $\cK:=\left\{\frac12\delta_{(-x,1)}+\frac12\delta_{(x,-1)}:\ x\in[-1,1]\right\}$. Then, $\cK$ is a closed geodesically convex set in $(\cP_2(\Om),W_2)$. Let $\mu:=\frac12\delta_{(-1,0)}+\frac12\delta_{(1,0)}$. Clearly, both measures $\rho_0:=\frac12\delta_{(-1,1)}+\frac12\delta_{(1,-1)}$ and $\rho_1:=\frac12\delta_{(1,1)}+\frac12\delta_{(-1,-1)}$ belong to $\cK$ and have the same minimal $W_2$ distance from $\mu$. So the projection of $\mu$ onto $\cK$ cannot be defined in a unique way in this case.

Because of this reason, in our study we will focus on some particular geodesically convex closed subsets $\cK\subset\cP_p(\Om)$ onto which we can guarantee the uniqueness of the projected measure in \eqref{def:proj}.

For a given $\lambda>0$, we consider
\begin{align}\label{def:Kp}
	\cK^\lambda_p(\Omega) := \left\{\rho \in \mathcal{P}_p(\Omega) \cap L^1(\Omega) \text{ and } 0 \leq \rho \leq \lambda {\rm{\ a.e.}}\right\},
\end{align}
that is the subset of absolutely continuous probability measures having densities uniformly bounded above by $\lambda$. Since the value of $\lambda$ will not play any role in our analysis, in the rest of the paper for simplicity of the exposition we will simply set $\lambda=1$ and we use the notation $\cK_p(\Omega)$ for $\cK^1_p(\Omega)$.

As we show in Lemma \ref{lem:admissibles_convex}, $\cK_p(\Om)$ is a closed geodesically convex subset of $(\cP_p(\Om),W_p)$. More importantly, arguments verbatim to the ones in \cite[Proposition 5.2]{DePMesSanVel} (that considered only the case $p=2$) let us conclude that the projection problem \eqref{def:proj} onto $\cK_p(\Om)$ has a unique solution for any $p>1$. A secondary motivation behind the consideration of these particular subsets is the following: in recent years the set $\cK_2(\Om)$ received some special attention in applications of optimal transport techniques to study the well-posedness and further properties of PDEs arising in crowd motion models under {\it density constraints}. For a non-exhaustive list of references on this subject we refer to \cite{MauRouSan,MesSan, DePMesSanVel,Mes}.

Coming back to the original motivation of our study, on the next pages of this note we investigate the question of nonexpansiveness of the projection operator ${\rm{P}}_\Om^p$ onto the set $\cK_p(\Om)$. When $d=1$, it is well-known that $\cP_2(\R)$ is isometrically isomorphic to a closed convex subset of a Hilbert space (the space of nondecreasing functions belonging to $L^2([0,1];\R)$, see \cite[Section 9.1]{AmbGigSav}). Therefore, $(\cP_2(\R),W_2)$ can be regarded as a flat space and so it is expected that ${\rm{P}}_\R^2$ is nonexpansive onto closed geodesically convex subsets $\cK\subset\cP_2(\R)$. Indeed, every closed geodesically convex subset $\cK\subset\cP_2(\R)$ corresponds to a closed convex subset of $L^2([0,1];\R)$. For instance, the space $\cK_2(\R)$ defined in \eqref{def:Kp} corresponds to $\{X\in L^2([0,1];\R):\ X'\ge 1\ {\rm{a.e.}}\}$. Therefore, the projection problem from $(\cP_2(\R),W_2)$ onto $\cK_2(\R)$ can be transferred to a projection problem in a Hilbertian setting, which has the nonexpansive property. Returning to the original setting via the isometric isomorphism, it follows that $\textnormal{P}_{\mathbf{R}}^2$ is nonexpansive. When $p\neq 2$, the nonexpansiveness property of projection operators on $L^p$ spaces is a more subtle question (see for instance \cite{Bru}) and therefore a conclusion similar to the one when $p=2$ seems to be nontrivial. In the case when $d>1$, even for $p=2$, it is not possible to identify $(\cP_2(\R^d),W_2)$ with the subset of a Hilbert space (and in particular, as discussed before, this space will not be flat). Therefore in those cases, the questions of nonexpansiveness seems to be highly nontrivial.

When $p=2$, Theorem \ref{thm:almost_lipschitz} presents a sort of {\it weak nonexpansiveness property} of the projection in arbitrary dimensions. Here we show that the left hand side of the inequality in \eqref{question:nonexpansive?} is always bounded above by the transportation cost of a certain suboptimal plan between the original measures. This suboptimal plan becomes optimal in two borderline scenarios: either when $d=1$ or when one of the original measures $\mu,\nu$ is a Dirac mass. So, this yields the nonexpansiveness of the projection operator when $d=1$ (see Corollary \ref{cor:1d}) or when one of the measures is a Dirac mass (see Corollary \ref{cor:dirac_nonexpansive}).

By \cite[Theorem 2.3]{KriRep} and \cite[Proposition 2.4 of Chapter II.2]{BriHae} we know that in the case of smooth Riemannian manifolds with non-positive sectional curvature or in the case of Alexandrov spaces having non-positive curvature the projection operator onto closed geodesically convex sets is nonexpansive. To the best of our knowledge, it is unclear whether the non-positive curvature condition of these spaces in general (beyond Riemannian manifolds) is also a necessary condition to ensure the nonexpansiveness of the projection operator in general. When $d\ge 2$, as we previously discussed, $(\cP_p(\Om),W_p)$ is expected to be non-negatively curved (even for $p\neq 2$). Therefore, there is a good reason to anticipate that there exist closed geodesically convex subsets of $\cP_p(\Om)$ such that the projection operator ${\rm{P}}_\Om^p$ onto these sets (whenever it is well-defined) fails to be nonexpansive.
This is precisely what we show in Proposition \ref{prop:small_p_counterexample}. Here, we will show in particular that there exists $p(d)>1$ small such that for any $p\in (1,p(d))$, the projection operator ${\rm{P}}_{\R^d}^p$ onto $\cK_p(\R^d)$ fails to be nonexpansive. Our proof is constructive, i.e. we construct a counterexample to the nonexpansiveness property. In our construction, we provide a quantitative asymptotic description of $p(d)$ as a function of $d$, for $d$ large.  Interestingly, relying again on Corollary \ref{cor:dirac_nonexpansive}, our construction does not provide a counterexample for $p=2$. Heuristically, our result would provide an argument  (in combination with \cite[Proposition 2.4 of Chapter II.2]{BriHae}) for the fact that $(\cP_p(\Om),W_p)$ is positively curved, when $p>1$ is close to $1$.

The structure of the rest of the paper is simple: in Section \ref{sec:prelim} we recall some preliminary results from optimal transport and we study some geometric properties of the projection operator ${\rm{P}}_\Om^2$. These properties seem to be interesting on their own right: we show that ${\rm{P}}_{\R^d}^2$ preserves barycenters of measures (see Proposition \ref{prop:barycenters}), and it satisfies a certain translation invariance with respect to distances between measures (see Proposition \ref{prop:translation_invariance}). Section \ref{sec:results} contains the proofs of our main results: in Theorem \ref{thm:almost_lipschitz} we show the `weak nonexpansiveness' property of ${\rm{P}}_\Om^2$, and deduce the full nonexpansiveness in the two cases mentioned above (see Corollaries \ref{cor:1d} and \ref{cor:dirac_nonexpansive}) Finally, Proposition \ref{prop:small_p_counterexample} constructs the counterexample to the nonexpansiveness of ${\rm{P}}_{\R^d}^p$ onto $\cK_p(\R^d)$ when $d\ge 2$ and $p\in (1,p(d))$, and studies the asymptotic behavior of $p(d)$ as the dimension becomes large.

\section{Preliminary results and some geometric properties of ${\rm{P}}_\Om^2$}\label{sec:prelim}

\subsection{Preliminary results from optimal transport} Some properties of the projection operator ${\rm{P}}_\Om^2$ onto the set $\cK_2(\Om)$ were studied in \cite{DePMesSanVel}. In particular, arguments verbatim to the ones presented there (see in particular Proposition 5.2) yield the following lemma.

\begin{lem}
Let $p\in(1,+\infty)$ and let $\Om\subseteq\R^d$ be closed and convex. Then for $\cK=\cK_p(\Om)$ with $
\cL^d(\Om)\ge 1$ and for any $\mu\in\cP_p(\Om)$, the problem \eqref{def:proj} has a unique solution ${\rm{P}}_\Om^p[\mu]$. Moreover, there exists $B\subseteq\Om$ Borel measurable such that
$${\rm{P}}_\Om^p[\mu]=\mu^{\rm{ac}}{\mathbbm{1}}_B+ \mathbbm{1}_{\Om\setminus B}.$$
\end{lem}

\begin{rem}
In the case of $p=2$, the projection ${\rm{P}}_{\Omega}^2$ behaves well with respect to interpolation along {\it generalized geodesics}. Let $\mu,\nu_0,\nu_1 \in \mathcal{P}_2(\Omega)$ with $\mu$ absolutely continuous. Then there are optimal maps $T_0,T_1$ which send $\mu$ to $\nu_0,\nu_1$ respectively. For $t \in [0,1]$, define the generalized geodesic connecting $\nu_0$ and $\nu_1$ with respect to $\mu$ by $\nu_t = (T_t)_{\sharp}\mu$, where $T_t = (1-t)T_0 + tT_1$. In \cite{MauRouSan}, using the displacement 1-convexity of $W_2^2(\mu,\cdot)$ along generalized geodesics, i.e. for all $t\in[0,1]$
\begin{align*}
		W_2^2(\mu,\nu_t) \leq (1-t)W_2^2(\mu,\nu_0) + tW_2^2(\mu,\nu_1) - t(1-t)W_2^2(\nu_0,\nu_1),
	\end{align*}
it was shown that ${\rm{P}}_\Om^2$ is locally $\frac12$-H\"older continuous. Since this argument is relying on the ``Hilbertian like'' behavior of $W_2$, it is unclear to us whether such reasoning could be carried through for $p\neq 2$.
\end{rem}

\begin{lem}\label{lem:admissibles_convex}
	Let $\Om\subseteq\R^d$ be closed and convex such that $\cL^d(\Om)\ge 1$ and $p\in(1,+\infty)$. The subspace $\cK_p(\Omega)$ defined in \eqref{def:Kp} is closed and geodesically convex in $(\cP_p(\Om),W_p)$.
\end{lem}

\begin{proof}
The closedness of $\cK_p(\Om)$ is straight forward.

Let $\mu_0,\mu_1\in\cK_p(\Om)$. To show that $\cK_p(\Om)$ is geodesically convex, we will show that the constant speed geodesic $[0,1]\ni t\mapsto\mu_t$ connecting $\mu_0$ to $\mu_1$ is absolutely continuous with a density bounded above by $1$. This is a consequence of \cite[Theorem 7.28]{San}, however for completeness we provide a direct proof of this result.

First, by \cite[Lemma 3.14]{Kel}, we have that $\mu_t\ll\cL^d\mres\Om$ for all $t\in[0,1]$.

To show the upper bound on $\mu_t$, we rely on the interpolation inequality for the Jacobian determinants of optimal transport maps provided in \cite[Theorem 3.13]{Kel} (see also \cite{CorMcCSch} for $p=2$). Let $\phi:\Om\to\R$ be the unique $c$-concave Kantorovich potential in the transport of $\mu_0$ onto $\mu_1$. Then, by \cite[Theorem 3.4]{Kel}, $T(x):=x-|\nabla\phi(x)|^{q-2}\nabla\phi(x)$ (where $1/p+1/q=1$) is the unique optimal transport map between $\mu_0$ and $\mu_1$. Moreover, $T_t(x):=x-t|\nabla\phi(x)|^{q-2}\nabla\phi(x)$ is the unique optimal transport map between $\mu_0$ and $\mu_t$.

Let us denote by $\Om_{\rm{id}}\subset\Om$ the set where $\phi$ is differentiable and $\nabla\phi=0$. Then, reasoning as in \cite{Kel}, we know that there exists a set $B\subseteq \Om\setminus\Om_{\rm{id}}$ of full measure such that $\phi $ is twice differentiable on $B$ with $\det(DT(x))>0$ if $x\in B$. Then by \cite[Theorem 3.13]{Kel} we have
\begin{equation}\label{eq:Jac}
\det(DT_t(x))^{\frac{1}{d}}\ge (1-t) + t \det(DT(x))^{\frac1d}.
\end{equation}
We remark that because our underlying space $\Om$ is flat, the volume distortion coefficients present in the previous inequality (stated in \cite{Kel} for general Finslerian manifolds) become 1.

By \eqref{eq:Jac}, if $\det(DT(x))\ge 1$, we conclude that  $\det(DT_t(x))\ge 1$, while if $\det(DT(x))\le 1$, then $\det(DT_t(x))\ge \det(DT(x))$. In conclusion,
$$
\det(DT_t(x))\ge \min\{1,\det(DT(x))\}.
$$

Now, since $\mu_t=(T_t)_\sharp\mu_0$, when restricted to the set $B$, the change of variable formula yields
\begin{align*}
\mu_t(T_t(x)) = \frac{\mu_0(x)}{\det(DT_t(x))}\le \frac{\mu_0(x)}{\min\{1,\det(DT(x))\}}\le\max\{\mu_0(x),\mu_1(T(x))\}\le 1.
\end{align*}
When restricted to the relative complement of $B$, $T_t$ essentially is the identity map, where the upper bound is also clearly preserved. The result follows.
\end{proof}

\subsection{Barycenters and translation invariance of ${\rm{P}}_\Om^2[\mu]$}

Suppose for now that $\Omega = \mathbf{R}^d$, so we do not have to worry about boundary issues. Then there are a few symmetries which one can exploit in Question \ref{question:nonexpansive?}. First, the projection operator ${\rm{P}}_{\Omega}^p$ commutes with translations. When $p=2$, the projection also preserves barycenters:

\begin{prop}\label{prop:barycenters}
Let $\mu \in \mathcal{P}_2( \mathbf{R}^d)$ and $\rho:= {\rm{P}}_{ \mathbf{R}^d}^2\mu$. Then
	\begin{align*}
		\int_ {\mathbf{R}^d} x \d\rho = \int_{ \mathbf{R}^d} x \d\mu.
	\end{align*}
\end{prop}

\begin{proof}
	For $h \in \mathbf{R}^d$, let $\tau \colon x \mapsto x+h$ denote the translation map by $h$. Then $\tau_{\sharp}\rho \in \cK_2( \mathbf{R}^d)$. Thus if $\gamma$ is an optimal plan between $\mu$ and $\rho$, then by Lemma \ref{lem:plan_translate} $(\id,\tau)_{\sharp}\gamma$ is optimal for $W_2(\mu,\tau_{\sharp}\rho)$. Thus, by the optimality of both $\rho$ and $\gamma$,
	\begin{align*}
		W^2_2(\mu,\rho)&=\int_{ \mathbf{R}^d\times \mathbf{R}^d} |x-y|^2 \d\gamma \leq W_2^2(\mu,\tau_{\sharp}\rho)= \int_{ \mathbf{R}^d\times\mathbf{R}^d} |x-y|^2 \d[(\id,\tau)_{\sharp}\gamma]\\
		& = \int_{ \mathbf{R}^d\times \mathbf{R}^d} |x-y-h|^2 \d\gamma = \int_{ \mathbf{R}^d\times \mathbf{R}^d}|x-y|^2\d\gamma - 2h \cdot \int_{ \mathbf{R}^d\times  \mathbf{R}^d} (x-y)\d\gamma + |h|^2.
	\end{align*}
	We conclude that $\displaystyle\int_{ \mathbf{R}^d\times \mathbf{R}^d}(x-y)\d\gamma = 0$, since otherwise one could set $\displaystyle h:=\lambda\int_{ \mathbf{R}^d\times \mathbf{R}^d}(x-y)\d\gamma$ and by sending $\lambda\downarrow 0$, the previous inequality would yield a contradiction. The result follows.
\end{proof}

\begin{lem}\label{lem:plan_translate}
Let $\mu,\nu\in\mathcal{P}_2(\mathbf{R}^d)$ and let $\nu'\in\mathcal{P}_2(\mathbf{R}^d)$ a translation of $\nu$, i.e. $\nu'=\tau_\sharp\nu$, where $\tau:x\mapsto x+h$ (for some given $h\in\mathbf{R}^d$). If $\gamma\in\mathcal{P}_2(\mathbf{R}^d\times\mathbf{R}^d)$ is optimal for $W_2^2(\mu,\nu)$, then $(\id,\tau)_\sharp\gamma$ is optimal for $W_2(\mu,\nu').$
\end{lem}

\begin{proof}
It is immediate to check that $\tilde\gamma:=(\id,\tau)_\sharp\gamma$ is an admissible plan for $W_2(\mu,\nu').$

By \cite[Theorem 5.10]{Vil} (see also \cite[Section 1.6.2]{San}) it is enough to show that $\tilde\gamma$ has cyclic monotone support. Let $n\in\mathbb{N}$. We notice that a collection of $n$ points from $\spt(\tilde\gamma)$ has the form $(x_i,y_i+h)_{i=1}^n$, where $(x_i,y_i)\in\spt(\gamma)$, $i\in\{1,\dots,n\}$. Let $\sigma:\{1,\dots,n\}\to\{1,\dots,n\}$ be a permutation of $n$ letters. Then we have
\begin{align*}
\sum_{i=1}^n |x_i-y_i-h|^2 & = \sum_{i=1}^n |x_i-y_i|^2-2\sum_{i=1}^n(x_i-y_i)\cdot h + n|h|^2\\
&\le \sum_{i=1}^n |x_i-y_{\sigma(i)}|^2-2\sum_{i=1}^n(x_i-y_i)\cdot h + n|h|^2 = \sum_{i=1}^n |x_i-y_{\sigma(i)}-h|^2,
\end{align*}
where in the inequality we have used the cyclic monotonicity of $\spt(\gamma)$. The result follows.
\end{proof}

From these observations one obtains the following ``translation invariance" when $p = 2$ and the ambient space is $ \mathbf{R}^d$.

\begin{prop}\label{prop:translation_invariance}
Let $\mu,\nu \in \mathcal{P}_2(\mathbf{R}^d)$ and $\nu'$ a translate of $\nu$. Then
	\begin{align*}
		W_2^2(\mu,\nu) - W_2^2(\textnormal{P}_{ \mathbf{R}^d}^2[\mu],\textnormal{P}_{ \mathbf{R}^d}^2[\nu]) = W_2^2(\mu,\nu') - W_2^2(\textnormal{P}_{ \mathbf{R}^d}^2[\mu],\textnormal{P}_{ \mathbf{R}^d}^2[\nu']).
	\end{align*}
\end{prop}

\begin{proof}
	Denote $\rho := \mathcal{P}_{ \mathbf{R}^d}^2[\mu]$, $\sigma := \mathcal{P}_{ \mathbf{R}^d}^2[\nu]$, and $\sigma' := \mathcal{P}_{ \mathbf{R}^d}^2[\nu']$. Let $\gamma\in\Pi_o(\mu,\nu)$ and $\eta\in\Pi_o(\rho,\sigma)$. If $\tau \colon x \mapsto x+h$ is the translation map which pushes forward $\nu$ onto $\nu'$, then we can construct optimal plans $\gamma',\eta'$ from $\mu,\rho$ to $\nu',\sigma'$, respectively, by $\gamma' = (\id,\tau)_{\sharp}\gamma$ and $\eta' = (\id,\tau)_{\sharp}\eta$ (see Lemma \ref{lem:plan_translate}; here we have also used the fact that the projection of the translate of a measure is the translate of the projection).
Thus
	\begin{align*}
		W_2^2(\mu,\nu') - W_2^2(\rho,\sigma')
		&= \int_{ \mathbf{R}^d\times \mathbf{R}^d} |x-y-h|^2 \d\gamma - \int_{ \mathbf{R}^d\times \mathbf{R}^d} |x-y-h|^2 \d\eta
		\\&= \int_{ \mathbf{R}^d\times  \mathbf{R}^d} |x-y|^2 \d\gamma - \int_{ \mathbf{R}^d\times \mathbf{R}^d} |x-y|^2 \d\eta\\
		& - 2h \cdot \int_{ \mathbf{R}^d\times \mathbf{R}^d} (x-y) \d\gamma + 2h \cdot \int_{ \mathbf{R}^d\times \mathbf{R}^d} (x-y) \d\eta
		\\&= W_2^2(\mu,\nu) - W_2^2(\rho,\sigma) + 2h \cdot \paren{\int_{ \mathbf{R}^d} x \d\rho - \int_{ \mathbf{R}^d} x \d\mu} + 2h \cdot \paren{\int_{ \mathbf{R}^d} y \d\nu - \int_{ \mathbf{R}^d} y \d\sigma},
	\end{align*}
	and the last two terms vanish by Proposition \ref{prop:barycenters}.
\end{proof}

In particular, any counterexample $\mu,\nu$ to nonexpansiveness must remain a counterexample when $\mu,\nu$ are replaced by translates of themselves. This already eliminates several candidates $\mu,\nu$ that may seem like potential counterexamples at first sight.

\section{Main Results}\label{sec:results}

Throughout this section, let $\mu,\nu \in \mathcal{P}_p(\Omega)$, and set $\rho = {\rm{P}}_{\Omega}^p(\mu)$ and $\sigma = {\rm{P}}_{\Omega}^p(\nu)$. Denote the optimal transport plan from $\rho$ to $\sigma$ by $\eta$. Note that since $\rho,\sigma$ are absolutely continuous, $\eta$ is induced by a map.

\subsection{Weak nonexpansiveness of the projection when $p = 2$.}

In the theorem below one bounds the distance squared between $\mu$ and $\nu$ by the transportation cost of a slightly suboptimal transport plan. This is a sort of ``weak nonexpansiveness.''

\begin{thm}\label{thm:almost_lipschitz}
	Let $\Omega\subseteq \mathbf{R}^d$ be a closed convex set. Let $T,U:\Omega\to\Omega$ stand for the  optimal maps from $\rho,\sigma$ to $\mu,\nu$, respectively. Take $p = 2$ and $\gamma := (T,U)_{\sharp}\eta\in\Pi(\mu,\nu)$. Then
	\begin{align*}
		W_2^2(\rho,\sigma) \leq \int_{\Omega\times\Omega} |x-y|^2 \d\gamma(x,y).
	\end{align*}
\end{thm}

\begin{proof}
	One can write
	\begin{align*}
	\int_{\Omega\times\Omega} |x-y|^2 \d\gamma
	&= \int_{\Omega\times\Omega} |T(x)-U(y)|^2 \d\eta
	= \int_{\Omega\times\Omega} |x-y + T(x)-x + y-U(y)|^2 \d\eta
	\\&= \int_{\Omega\times\Omega} |x-y|^2 \d\eta + 2\int_{\Omega\times\Omega} (x-y) \cdot (T(x)-x+y-U(y)) \d\eta\\
	& + \int_{\Omega\times\Omega} |T(x)-x+y-U(y)|^2 \d\eta
	\\&\geq \int_{\Omega\times\Omega} |x-y|^2 \d\eta + 2\int_{\Omega\times\Omega} (x-y) \cdot (T(x)-x) \d\eta + 2\int_{\Omega\times\Omega} (y-x) \cdot (U(y)-y) \d\eta.
	\end{align*}
	Thus it suffices to show that
	\begin{align*}
	\int_{\Omega\times\Omega} (x-y) \cdot (T(x)-x) \d\eta \geq 0
	\qquad \text{and (by symmetry)} \qquad
	\int_{\Omega\times\Omega} (y-x) \cdot (U(y)-y) \d\eta \geq 0.
	\end{align*}
	For $t \in (0,1)$, let $\pi_t(x,y) = (1-t)x+ty$. Then $\rho_t:=(\pi_t)_{\sharp}\eta \in \cK_2(\Omega)$ by the geodesic convexity of $\mathcal{K}_2(\Omega)$. The optimality of $\rho$ in the definition of $\textnormal{P}_{\Omega}^2(\mu)$, together with the fact that $\tilde\eta:=(T,\pi_t)_\sharp\eta\in \Pi(\mu,\rho_t)$, implies
	\begin{align*}
	W_2^2(\mu,\rho)
	&\le W_2^2(\mu,\rho_t)\leq \int_{\Omega\times\Omega}|x-y|^2\d\tilde\eta= \int_{\Omega\times\Omega} |T(x)-\pi_t(x,y)|^2 \d\eta
	= \int_{\Omega\times\Omega} |T(x)-x + t(x-y)|^2 \d\eta\\
	&= \int_{\Omega\times\Omega} |T(x)-x|^2 \d\eta + 2t \int_{\Omega\times\Omega} (x-y) \cdot (T(x)-x) \d\eta + t^2 \int_{\Omega\times\Omega} |x-y|^2 \d\eta\\
	&=W_2^2(\mu,\rho)+ 2t \int_{\Omega\times\Omega} (x-y) \cdot (T(x)-x) \d\eta+ t^2W_2^2(\rho,\sigma).
	\end{align*}
Thus, we have obtained
\begin{align*}
-tW_2^2(\rho,\sigma)\le 2 \int_{\Omega\times\Omega} (x-y) \cdot (T(x)-x) \d\eta.
\end{align*}

Letting $t \rightarrow 0$, we conclude that
\begin{align*}
\int_{\Omega\times\Omega} (x-y) \cdot (T(x)-x) \d\eta \geq 0,
\end{align*}
as desired.
\end{proof}

When $\Omega \subseteq \mathbf{R}$, this theorem is enough to deduce that the answer to Question \ref{question:nonexpansive?} is \emph{yes}.

\begin{cor}\label{cor:1d}
	Suppose $\Omega \subseteq \mathbf{R}$. Then ${\rm{P}}_{\Omega}^2$ is nonexapansive.
\end{cor}

\begin{proof}
	The plan $\gamma$ defined in Theorem \ref{thm:almost_lipschitz} is monotonically increasing, hence optimal.
\end{proof}

Theorem \ref{thm:almost_lipschitz} also implies nonexpansiveness when $\Pi(\mu,\nu)$ is a singleton. This is for instance the case when one of the measures $\mu,\nu$ is a Dirac mass.

\begin{cor}\label{cor:dirac_nonexpansive}
    Suppose that $\mu,\nu$ are such that $\Pi(\mu,\nu)$ is a singleton. Then
    \begin{align*}
        W_2(\rho,\sigma) \leq W_2(\mu,\nu).
    \end{align*}
\end{cor}

\begin{proof}
    There is only one transport plan between $\mu$ and $\nu$, so using the notation of Theorem \ref{thm:almost_lipschitz}, $\gamma\in\Pi(\mu,\nu)$ must be this plan. The result follows.
\end{proof}

\begin{rem}
	In general, $\gamma\in\Pi(\mu,\nu)$ in the statement of Theorem \ref{thm:almost_lipschitz} does not need to be optimal. In the case of $\Omega = \mathbf{R}^2$, consider for instance $\mu = \frac{1}{2}\delta_{(R,0)} + \frac{1}{2}\delta_{(-R,0)}$ and $\nu = \frac{1}{2}\delta_{(t,1)} + \frac{1}{2}\delta_{(-t,-1)}$, where $R$ is large and $t$ is small. For $t > 0$, the optimal map from $\mu$ to $\nu$ sends all the mass from $(R,0)$ to $(t,1)$, and all the mass from $(-R,0)$ to $(-t,-1)$. On the other hand, for $t < 0$, the optimal map sends all the mass from $(R,0)$ to $(-t,-1)$, and all the mass from $(-R,0)$ to $(t,1)$. This means that the optimal plan from $\mu$ to $\nu$ does not vary continuously with $t$ (it is discontinuous at $t = 0$). However, one can see that the plan $\gamma$ does depend continuously on $t$, so it cannot be optimal.
\end{rem}

\subsection{Failure of nonexpansiveness of the projection for $p=1 + o(1)$}

A restriction on $p$ is necessary in the statement of Proposition \ref{prop:small_p_counterexample}. For $p$ very close to $1$, the proof of the following proposition illustrates a counterexample to the nonexpansive property of ${\rm{P}}_{\mathbf{R}^d}^p$ onto $\cK_p(\mathbf{R}^d)$.

\begin{prop}\label{prop:small_p_counterexample}
	Let $\Omega = \mathbf{R}^d$ with $d > 1$. Then there exists $p(d) > 1$ such that  ${\rm{P}}_{\Omega}^p$ is not nonexpansive with respect to the $p$-Wasserstein distance for $1 < p < p(d)$. In fact, one can take
	\begin{align*}
	    p(d) = 1 + \frac{1}{O(d^2 \log d)}.
	\end{align*}
\end{prop}
In the proof below we use the following conventions: for positive quantities $A,B$ possibly depending on various parameters, we write $A \lesssim B$ if $A \leq CB$ with $C$ an absolute constant, $A \gtrsim B$ if $B \lesssim A$, and $A \sim B$ if $A \lesssim B \lesssim A$.

\begin{proof}[Proof of Proposition \ref{prop:small_p_counterexample}]
	Let $R>0$ be the radius of the ball of volume $\frac{1}{2}$. Let $\mu = \frac{1}{2}\delta_{(0,\dots,0)} + \frac{1}{2}\delta_{(2R,0,\dots,0)}$ and $\nu = \delta_{(0,\dots,0)}$. Then $\rho = {\rm{P}}_{\Omega}^p(\mu)$ and $\sigma = {\rm{P}}_{\Omega}^p(\nu)$ are the restriction of Lebesgue measure to
	\begin{align*}
		\{x \in \mathbf{R}^d : |x| \leq R \text{ or } |x-(2R,0,\dots,0)| \leq R\}
		\qquad \text{and} \qquad
		\{x \in \mathbf{R}^d : |x| \leq 2^{1/d}R\}
	\end{align*}
	respectively. Let $\gamma\in\Pi(\mu,\nu)$ and $\eta\in\Pi(\rho,\sigma)$ be arbitrary transport plans.
	We will show that
	\begin{align}\label{eqn:W1_ctrexample}
		\int_{\mathbf{R}^d\times\mathbf{R}^d} |x-y| \d\eta > \int_{\mathbf{R}^d\times\mathbf{R}^d} |x-y| \d\gamma
	\end{align}
	with an explicit lower bound on the difference; then we obtain the desired inequality
	\begin{align*}
		\int_{\mathbf{R}^d\times\mathbf{R}^d} |x-y|^p \d\eta > \int_{\mathbf{R}^d\times\mathbf{R}^d} |x-y|^p \d\gamma
	\end{align*}
	for $p \in [1,p(d))$ by continuity in $p$.

	The right hand side of \eqref{eqn:W1_ctrexample} is necessarily equal to
	\begin{align}\label{eqn:W1_mu_nu}
		\int_{\mathbf{R}^d\times\mathbf{R}^d}|x-y|\d\gamma = R = \int_{\R^d} x_1 \d\rho - \int_{\R^d} y_1 \d\sigma = \int_{\mathbf{R}^d\times\mathbf{R}^d} (x_1-y_1) \d\eta.
	\end{align}
	Thus to get a quantitative form of \eqref{eqn:W1_ctrexample}, it is enough to estimate
	\begin{align}\label{eqn:to_lbd}
		\int_{\mathbf{R}^d\times\mathbf{R}^d} (|x-y|-|x_1-y_1|) \d\eta
	\end{align}
	from below. Given $x \in \mathbf{R}^d$, denote $x' = (x_2,\dots,x_d) \in \mathbf{R}^{d-1}$. Let
	\begin{align*}
		E = \{y \in \mathbf{R}^d \colon 1.1^{1/d}R \leq |y'| \leq 1.9^{1/d}R \text{ and } |y_1| \leq \sqrt{2^{2/d}-1.9^{2/d}}R\} \subseteq \{x \in \mathbf{R}^d \colon |x| \leq 2^{1/d}R\} = \spt\sigma.
	\end{align*}
	Suppose $(x,y) \in \spt\eta$ with $y \in E$. Then $|x'| \leq R$, so
	\begin{align*}
		|x'-y'| \geq (1.1^{1/d}-1)R
		= R \int_{0}^{1/d} 1.1^t \log 1.1 \d t
		\gtrsim R/d,
	\end{align*}
	On the other hand
	\begin{align*}
		|x-y| \leq \diam\{\spt\rho \cup \spt\sigma\} \lesssim R.
	\end{align*}
	Combining these two facts yields
	\begin{align*}
		|x-y| - |x_1-y_1| \gtrsim \Big(\sqrt{1+\frac{1}{d^2}} - 1\Big)R \gtrsim \frac{R}{d^2}.
	\end{align*}
	Plugging this into \eqref{eqn:to_lbd} and recalling \eqref{eqn:W1_mu_nu}, we deduce that
	\begin{align}\label{eqn:abstract_lbd}
		\int_{\R^d\times\R^d} |x-y| \d\eta - \int_{\R^d\times\R^d} |x-y| \d\gamma &\geq \int_{\R^d\times\R^d} (|x-y| - |x_1-y_1|) \d\eta\\
		&\nonumber \geq \int_{\R^d \times E} (|x-y|-|x_1-y_1|) \d\eta
		\gtrsim \sigma(E) \frac{R}{d^2}.
	\end{align}
	In computing $\sigma(E)$, it will be convenient to write $\Vol_n(r)$ for the volume of the $n$-dimensional ball of radius $r$, and $\Rad_n(v)$ for the radius of the $n$-dimensional ball of volume $v$. Then $R = \Rad_d(1/2)$ by definition, and by classical formulas,
	\begin{align*}
		\Vol_n(r) = \frac{\pi^{n/2}}{\Gamma(n/2+1)} r^n
		\qquad \text{ and } \qquad
		\Rad_n(v) = \frac{\Gamma(n/2+1)^{1/n}}{\sqrt{\pi}} v^{1/n} \sim \sqrt{n} v^{1/n}.
	\end{align*}
	Thus $R \sim \sqrt{d}$, and
	\begin{align*}
		\sigma(E) &= 2\sqrt{2^{2/d} - 1.9^{2/d}} R [\Vol_{d-1}(1.9^{1/d}R) - \Vol_{d-1}(1.1^{1/d}R)]
		\\&\sim \frac{R}{\sqrt{d}} \Vol_{d-1}(R)
		= \frac{1}{\sqrt{d}} \frac{R\Vol_{d-1}(R)}{\Vol_d(R)} \Vol_d(R)
		\sim \frac{1}{\sqrt{d}} \frac{\Gamma(d/2+1)}{\Gamma((d-1)/2+1)}
		\sim 1,
	\end{align*}
	where the final estimate follows from Stirling's asymptotic for the Gamma function and from the fact that by the choice of $R$, $\Vol_d(R)=\frac12$.
	From \eqref{eqn:abstract_lbd} we therefore conclude
	\begin{align}\label{eqn:W1_case}
		\int_{\R^d\times\R^d} |x-y| \d\eta - \int_{\R^d\times\R^d} |x-y| \d\gamma \gtrsim \frac{1}{d^{3/2}}.
	\end{align}

	The inequality we need to prove is
	\begin{align}\label{eqn:goal}
		\int_{\R^d\times\R^d} |x-y|^p \d\gamma < \int_{\R^d\times\R^d} |x-y|^p \d\eta
	\end{align}
	for $p \in (1,p(d))$. If $c > 0$ is the implied constant in \eqref{eqn:W1_case}, then \eqref{eqn:goal} will follow from \eqref{eqn:W1_case} as long as $p$ is small enough that
	\begin{align*}
		\bracket{\int_{\R^d\times\R^d}|x-y|^p\d\gamma - \int_{\R^d\times\R^d}|x-y|\d\gamma} + \bracket{\int_{\R^d\times\R^d}|x-y|\d\eta - \int_{\R^d\times\R^d}|x-y|^p\d\eta} < \frac{c}{d^{3/2}}.
	\end{align*}
	The first term in brackets is simply
	\begin{align*}
		\int_{\R^d\times\R^d} |x-y|^p \d\gamma - \int_{\R^d\times\R^d} |x-y| \d\gamma  = 2^{p-1}R^p - R.
	\end{align*}
	Because of the general inequality $t-t^p \leq p-1$ for all $t \geq 0$ and $p\ge 1$, the second term in brackets must be at most $p-1$. Thus \eqref{eqn:goal} holds whenever
	\begin{align*}
		2^{p-1}R^p - R + p - 1 < \frac{c}{d^{3/2}}.
	\end{align*}
	One can estimate
	\begin{align*}
		2^{p-1}R^p - R \lesssim (p-1) R^p \log R
		\lesssim (p-1) d^{p/2} \log d
	\end{align*}
	for $p$ bounded, so it is enough if
	\begin{align*}
		(p-1)(1+d^{p/2}\log d) < \frac{c'}{d^{3/2}}
	\end{align*}
	for some smaller absolute constant $c' > 0$. This is true for
	\begin{align*}
		p < 1 + \frac{1}{O(d^2 \log d)},
	\end{align*}
	as long as the implied constant is sufficiently large.
\end{proof}
\begin{rem}
Let us note as a consequence of Corollary \ref{cor:dirac_nonexpansive} that this construction is not a counterexample when $p =2$.
\end{rem}

\medskip

\noindent
{\sc Acknowledgements}\ \  We thank Hugo Lavenant for his feedback on an earlier version of the manuscript.


\section{introduction}

Graph neural networks (GNNs)~\cite{bruna2013spectral,kipf2016semi,velivckovic2017graph,hamilton2017inductive} have achieved remarkable performance on various domains such as knowledge graph~\cite{hamaguchi2017knowledge,wang2018cross}, social media mining~\cite{hamilton2017inductive}, nature language processing~\cite{kipf2016semi,yao2019graph}, and recommendation system~\cite{ying2018graph,berg2017graph}. Generally, message-passing process
is adopted in GNNs~\cite{kipf2016semi,hamilton2017inductive}, where information from neighbors is aggregated for every node in each layer. This process enriches node representations, and preserves both node feature characteristics and topological structures.

Despite the success in modeling graph data, GNNs trained on graphs may inherit the societal bias in data, which limits the adoption of GNNs in many real-world applications. \textit{First}, extensive studies~\cite{dwork2012fairness,beutel2017data,creager2019flexibly} have revealed that historical data may include patterns of previous discrimination and societal bias. Machine learning models trained on such data can inherit the bias on sensitive attributes such as ages, genders, skin color, and regions~\cite{dwork2012fairness,beutel2017data,creager2019flexibly,buolamwini2018gender}, which implies that GNNs could also exhibit the bias. \textit{Second}, the topology of graphs and the message-passing of GNNs could magnify the bias. Generally, in graphs such as social networks, nodes of similar sensitive attributes are more likely to connect to each other than nodes of different sensitive attributes~\cite{dong2016young,rahman2019fairwalk,mehrabi2019survey,olteanu2019social}. For example, young people tend to build friendship with people of similar age on the social network \cite{dong2016young}.
This makes the aggregation of neighbors' features in GNN have similar representations for nodes of similar sensitive information while different representations for nodes of different sensitive features, leading to severe bias in decision making, i.e., the predictions are highly correlated with the sensitive attributes of the nodes. Our preliminary experiments in Sec.~\ref{sec:dis} indicate that GNNs have larger bias due to the adoption of graph structure than models which only use node attributes, which verifies our intuition. The bias would largely limit the wide adoption of GNNs in domains such as ranking
of job applicants~\cite{mehrabi2019survey}, crime rate prediction~\cite{suresh2019framework}, and credit score estimation~\cite{zemel2013learning}. Thus, it is critically important to investigate fair GNNs.

However, developing fair GNNs is a non-trivial task. \textit{First}, to achieve fairness, we need to obtain abundant nodes with known sensitive attributes so that we can either revise the data or regularize the model; whereas people are unwilling to share their sensitive information in the real-world, and resulting in inadequate nodes with sensitive attributes known for fair model learning.
For example, only 14\% teen users public their complete profiles on Facebook~\cite{madden2013teens}.
The lacking of sensitive information challenges many existing work on fair models~\cite{beutel2017data,locatello2019fairness,louizos2015variational,creager2019flexibly,song2018learning,liao2019learning,madras2018learning,zhang2018mitigating}. \textit{Second}, though extensive efforts have been made to establish fair models by revising features \cite{zhang2017achieving,kamiran2009classifying,kamiran2012data}, disentanglement \cite{louizos2015variational,creager2019flexibly}, adversarial debiasing \cite{edwards2015censoring,beutel2017data,zhang2018mitigating,madras2018learning} and fairness constraints \cite{zafar2015fairness,zafar2017fairness},
they are overwhelmingly dedicated to independently and identically distributed (i.i.d) data, which cannot be directly applied on graph data for the absence of simultaneous consideration of the bias from node attributes and graph structures.
Recently, \cite{rahman2019fairwalk,bose2019compositional} aim to learn fair node representations from  graphs. These methods merely deal with plain graphs without any node attributes, and focus on fair node representations instead of fair node classifications.

Therefore, in this paper, we study a novel problem of learning fair graph neural networks with limited sensitive information. In essence, we need to solve two challenges: (i) how to overcome the shortage of sensitive attributes for eliminating discrimination; and (ii) how to ensure the fairness of the GNN classifier. In an attempt to address these challenges, we propose a novel framework named as \textbf{FairGNN} for fair node classification.
A GNN sensitive attribute estimator is adopted in FairGNN to predict plenty of sensitive attributes with noise for fair classification. Inspired by existing works of fair classification on i.i.d data with adversarial learning \cite{edwards2015censoring,beutel2017data,zhang2018mitigating,madras2018learning}, we deploy an adversary to ensure the GNN classifier make predictions independent with the estimated sensitive attributes.
To further stabilize the training process and performance in fairness, we introduce a fairness constraint to make the predictions invariant with the estimated sensitive attributes.
Our main contributions are:
\begin{itemize}
	\item We study a novel problem of fair graph neutral networks learning with limited sensitive information;
	\item A new framework, FairGNN, is proposed to settle the shortage of sensitive attributes for adversarial debiasing and fairness constraint by estimating users' sensitive attributes;
	\item We conduct theoretical analysis showing fairness achieves at the global minimum even with estimated sensitive attributes;
	\item Extensive experiments on different datasets demonstrate the effectiveness of our methods in eliminating discrimination while keeping high accuracy of GNNs.
\end{itemize}

The rest of the paper is organized as follows. In Sec.~\ref{sec:related_work}, we review related work. In Sec.~\ref{sec:preliminary_analysis}, we  conduct preliminary analysis to understand the bias issue of GNNs. In Sec.~\ref{sec:methdology}, we give the details of FairGNN. In Sec.~\ref{sec:experiments}, we conduct experiments to show the effectiveness of FairGNN. In Sec.~\ref{sec:conclusion}, we conclude with future work.

.

\section{related work} \label{sec:related_work}
In this section, we will review related work including graph neural networks and fairness in machine learning.

\subsection{Graph Neural Networks}
Graph neural networks (GNNs), which generalize neural networks for graph structured data, have shown great success for various applications~\cite{hamaguchi2017knowledge,wang2018cross,kipf2016semi,yao2019graph,ying2018graph,berg2017graph}.
Generally, GNNs can be categorized into two categories, i.e., spectral-based~\cite{bruna2013spectral,henaff2015deep,defferrard2016convolutional,kipf2016semi,levie2018cayleynets} and spatial-based~\cite{velivckovic2017graph,hamilton2017inductive,chen2018fastgcn,chiang2019cluster,ying2018graph}.
Spectral-based GNNs define graph convolution based on spectral graph theory, which is first explored by \citeauthor{bruna2013spectral} \cite{bruna2013spectral}. Since then, more spectral-based methods are developed for further improvements and extensions \cite{henaff2015deep,defferrard2016convolutional,kipf2016semi,levie2018cayleynets}. Graph Convolutional Network (GCN)~\cite{kipf2016semi} is a particularly popular method which simplifies the convolutional operation on the graph. Spatial-based graph convolution directly updates the node representation by aggregating its neighborhoods' representations \cite{niepert2016learning,gilmer2017neural,hamilton2017inductive,ying2018graph}.
\citeauthor{velivckovic2017graph} \cite{velivckovic2017graph} introduce the self-attention into the aggregation of spatial graph convolution by assigning higher weights to the more important nodes in graph attention network (GAT). Various spatial methods are proposed to solve the scalability issue of GCN \cite{hamilton2017inductive, chen2018fastgcn,chiang2019cluster}.  For example, a neighbor sampling method to make it possible to
train GNN with nodes in mini-batch instead of the whole graph is developed in GraphSAGE~\cite{hamilton2017inductive}.
Moreover, spatial-based methods have already been successfully deployed to deal with extremely large industrial datasets \cite{ying2018graph}.

The essential idea of GNNs is to propagate the information of nodes through the graph to get better representations. However, people tend to build relationship with those sharing the same sensitive attributes. Then, representations in GNNs are nearly propagated within the subgroup, which highly increases the risk of discrimination towards sensitive attributes. Despite the risk of discrimination in GNNs, there is no existing work to address this problem. Thus, we study the novel problem of learning fair GNNs
which will eliminate the potential discrimination.

\subsection{Fairness in Machine Learning}
Many works have been conducted to deal with the bias in the training data to achieve fairness in machine learning~\cite{zhang2017achieving,kamiran2009classifying,kamiran2012data,beutel2017data,locatello2019fairness,dwork2012fairness,hardt2016equality}. Based on which stage of the machine learning training process is revised, algorithms could be split into three categories: the pre-processing approaches, the in-processing approaches, and the post-processing approaches.
The pre-processing approaches are applied before training machine learning models. They could reduce the bias by modifying the training data through correcting labels \cite{zhang2017achieving,kamiran2009classifying}, revising attributes of data \cite{kamiran2012data,feldman2015certifying}, generating non-discriminatory labeled data \cite{xu2018fairgan,xu2019fairgan+,sattigeri2019fairness}, and obtaining fair data representations \cite{beutel2017data,locatello2019fairness,edwards2015censoring,zemel2013learning,louizos2015variational,creager2019flexibly,song2018learning,liao2019learning,madras2018learning}. The in-processing approaches are designed to revise the training of the state-of-the-art models.
Typically the machine learning models are trained with
additional regularization terms or a new objective function. \cite{dwork2012fairness,zafar2015fairness,kamishima2011fairness,zhang2018mitigating}. Finally, the post-processing approaches directly change the predictive labels to ensure fairness \cite{hardt2016equality,pleiss2017fairness}.
Recently, several works explore the learning of fair graph embeddings for recommendation \cite{rahman2019fairwalk,bose2019compositional}. Fairwalk \cite{rahman2019fairwalk} modifies the random walk procedure of node2vec \cite{grover2016node2vec} to obtain a more diverse network neighborhood representations. The sensitive attributes of all the nodes are required in the sampling procedure of FairWalk. \citeauthor{bose2019compositional} \cite{bose2019compositional} propose to add discriminators
to eliminate the sensitive information in the graph embeddings. Similar to Fairwalk, the training process of the discriminators is in need of the sensitive attributes of all the nodes.

Our work is inherently different from existing works: (i) we focus on learning fair GNNs for node classification instead of fair graph embeddings; (ii)  we address the problem that only a limited number of nodes are provided with sensitive attributes in practice.

\section{Preliminaries Analysis} \label{sec:preliminary_analysis}
In this section, we first conduct preliminary analysis on real-world datasets to show that GNNs could exhibit more serve bias due to the graph structure and the message-passing. Sequentially, We formally give the problem definition of fair node classification.

\subsection{Notations} \label{sec:notations}
We use $\mathcal{G}=(\mathcal{V},\mathcal{E}, \mathbf{X})$ to denote an attributed graph, where $\mathcal{V}=\{v_1,...,v_N\}$ is the set of $N$ nodes, $\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}$ is the set of edges, and $\mathbf{X}=\{\mathbf{x}_1,...,\mathbf{x}_N\}$ is the set of node features. $\mathbf{A} \in \mathbb{R}^{N \times N}$ is the adjacency matrix of the graph $\mathcal{G}$, where $\mathbf{A}_{ij}=1$ if nodes ${v}_i$ and ${v}_j$ are connected; otherwise, $\mathbf{A}_{ij}=0$. In the semi-supervised setting, part of nodes $v \in \mathcal{V}_L$ are provided with labels $y_v \in \mathcal{Y}$,
where $\mathcal{V}_L \subseteq \mathcal{V}$ denotes nodes with labels, and $\mathcal{Y}$ is the set of labels. Sensitive attributes of training nodes are required to achieve fairness of machine learning algorithms. In our setting, only a small set of nodes $\mathcal{V}_S \subset \mathcal{V}$ are provided with the sensitive attribute $s \in \{0,1\} $. The set of provided sensitive attributes is denoted by $\mathcal{S}$.

\subsection{Datasets} \label{sec:datasets}
\begin{table}[t]
	\small
	\centering
	\caption{The statistics of datasets.}
	\vskip -1.5em
	\begin{tabularx}{0.88\linewidth}{p{0.35\linewidth}XXX}
		\toprule
		Dataset   & Pokec-z & Pokec-n  &  NBA  \\
		\midrule
		\# of nodes  & 67,797 & 66,569 & 403\\
		\# of node attributes  & 59 & 59 & 39\\
		\# of edges  & 882,765 & 729,129 & 16,570 \\
		Size of $\mathcal{V}_L$  & 500 & 500 & 100\\
		Size of $\mathcal{V}_S$ & 200 & 200 & 50 \\

		Group ratio & 1.84 & 2.46  & 2.77\\
		\# of inter-group edges  & 39,804 & 31,515 & 4,401\\
		\# of intra-group edges  & 842,961 & 697,614 & 12,169\\
		\bottomrule
	\end{tabularx}
	\label{tab:stat}
	\vskip -1.5em
\end{table}
For the purpose of this study, we collect and sample datasets from Pokec and NBA. The details are described as below.

\textbf{Pokec \cite{takac2012data}}: It is the most popular social network in Slovakia, which is very similar to Facebook and Twitter. This dataset contains  anonymized data of the whole social network in 2012. User profiles of Pokec contain gender, age, hobbies, interest, education, working field and etc. The original Pokec dataset contains millions of users. Based on the provinces that users belong to, we sampled two datasets named as: \textbf{Pokec-z} and \textbf{Pokec-n}. Both Pokec-z and Pokec-n consist of users belonging to two major
regions of the corresponding provinces. We treat the region as the sensitive attribute. The classification task is to predict the working field of the users.

\textbf{NBA}: This extended from a Kaggle dataset~\footnote{https://www.kaggle.com/noahgift/social-power-nba} containing around 400 basketball players in NBA.
The performance statistics of players in the 2016-2017 season and other various information e.g., nationality, age, and salary are provided. To obtain the graph that links the NBA players together, we collect the relationships of the NBA basketball players on Twitter with its official crawling API~\footnote{https://developer.twitter.com/en}. We binarize the nationality to
two categories, i.e., U.S. players and oversea players, which is used as sensitive attribute.
The classification task is to predict whether the salary of the player is over median.

For all the datasets, we eliminate nodes without any links with others. We randomly sample labels and sensitive attributes separately to get $\mathcal{V}_L$ and $\mathcal{V}_S$.
We randomly sample 25\% and 50\%  of nodes containing both sensitive attributes and labels in Pokec-z, Pokec-n and NBA as validation sets and test sets. Note that the validation sets and test sets have no overlap with $\mathcal{V}_L$ and $\mathcal{V}_S$.
The key statistics of the datasets are given in Table \ref{tab:stat}. Apart from the basic statistics, we also report the ratio of the majority and minority group and the number of edges linking the same group and different groups. It is evident from the table that: (i) skew exists in sensitive attributes; (ii) most of relationships are between users who share the same sensitive attribute.

\subsection{Preliminaries of Graph Neural Networks}

Graph neural networks (GNNs) utilize the node attributes and edges to learn a representation $\mathbf{h}_v$ of the node $v \in \mathcal{V}$. The goal of learning representation in node classification is to predict the node $v$'s label as $y_v = f(\mathbf{h}_v)$ . Current GNNs are neighborhood aggregation approaches, which will update the representations of the nodes with the representations of the neighborhood nodes. The representations after $k$ layers'
aggregation would capture the
structural information of the $k$-hop network neighborhoods. The updating process of the $k$-th layer in GNN could be formulated as:
\begin{equation}
\begin{aligned}
\mathbf{a}^{(k)}_v & = \text{AGGREGATE}^{(k-1)}(\{\mathbf{h}^{(k-1)}_u: u \in \mathcal{N}(v)\}),
\label{eq:GNN_a} \\
\mathbf{h}^{(k)}_{v} & =\text{COMBINE}^{(k)}(\mathbf{h}^{(k-1)}_v, \mathbf{a}^{(k)}),
\end{aligned}
\end{equation}
where $\mathbf{h}^{(k)}_v$ is the representation vector of the node $v \in \mathcal{V}$ at $k$-th layer and $\mathcal{N}(v)$ is  a set of neighborhoods of $v$.

\subsection{Fairness Evaluation Metrics}
In this subsection, we will present two definitions of fairness for the binary label $y \in \{0,1\}$ and the sensitive attribute $s \in \{0,1\}$. $\hat{y} \in  \{0, 1\}$ denotes the prediction of the classifier $\eta$: $\mathbf{x} \rightarrow y$.
\begin{definition}
	(Statistical Parity \cite{dwork2012fairness}). Statistical parity requires the predictions to be independent with the sensitive attribute $s$, i.e., $\hat{y} \bot s$. It could be formally written as:
	\begin{equation}
	P(\hat{y}|s=0)=P(\hat{y}|s=1).
	\end{equation}
\end{definition}
\begin{definition} (Equal Opportunity \cite{hardt2016equality}).
	Equal opportunity requires the probability of an instance in a positive class being assigned to a positive outcome should be equal for both subgroup members. The property of equal opportunity is defined as:
	\begin{equation}
	P(\hat{y}=1|y=1,s=0) = P(\hat{y}=1|y=1,s=1).
	\end{equation}
	The equal opportunity expects the classifier to give equal true positive rates across the subgroups.
	According to~\cite{louizos2015variational,beutel2017data}, we apply the following metrics to quantitatively evaluate statistical parity and equal opportunity:
	\begin{equation}
	\Delta_{SP}  = |P(\hat{y}=1|s=0)-P(\hat{y}=1|s=1)|,
	\end{equation}
	\begin{equation}
	\Delta_{EO} = |P(\hat{y}=1|y=1,s=0) - P(\hat{y}=1|y=1,s=1)|,
	\end{equation}
	where the probabilities are evaluated on the test set.
\end{definition}
\label{sec:fairness}
\subsection{Discrimination in Graph Neural Networks}

\begin{table}[t]
	\small
	\centering
	\caption{Results of models w/ and w/o utilizing graph.}
	\vskip -1.5em
	\begin{tabularx}{\linewidth}{|p{0.1\linewidth}|p{0.12\linewidth}|X|X|X|X|}
		\hline
		Dataset & Metrics & MLP & MLP-e & GCN & GAT \\
		\hline
		\hline
		\multirow{4}{*}{Pokec-z}
		& ACC (\%) & 65.3 $\pm 0.5$ & 68.6 $\pm 0.3$ & 70.2  $\pm 0.1$ & 70.4 $\pm 0.1$ \\
		& AUC (\%)& 71.3 $\pm 0.3$ & 74.8 $\pm 0.3$ & 77.2  $\pm 0.1$ & 76.7 $\pm 0.1$\\
		& $\Delta_{SP}$ (\%)& 3.8 $\pm 1.3$ & 6.9 $\pm 1.0$ & 9.9  $\pm 1.1$ & 9.1 $\pm 0.9$  \\
		& $\Delta_{EO}$ (\%)& 2.2 $\pm 0.7$ & 4.0 $\pm 1.5$ & 9.1  $\pm 0.6$ & 8.4  $\pm 0.6$ \\
		\hline
		\hline
		\multirow{4}{*}{Pokec-n}
		& ACC (\%) & 63.1 $\pm 0.4$ & 66.3 $\pm 0.6$ & 70.5 $\pm 0.2$ & 70.3 $\pm 0.1$\\
		& AUC (\%)& 68.2 $\pm 0.3$ & 72.4 $\pm 0.6$  & 75.1 $\pm 0.2$ & 75.1 $\pm 0.2$\\
		& $\Delta_{SP}$ (\%)& 3.3 $\pm 0.6$ & 8.7 $\pm 1.0$ & 9.6 $\pm 0.9$ & 9.4 $\pm 0.7$ \\
		& $\Delta_{EO}$ (\%)& 7.1 $\pm 0.9$ & 9.9 $\pm 0.6$ & 12.8 $\pm 1.3$ & 12.0 $\pm 1.5$ \\
		\hline
		\hline
		\multirow{4}{*}{NBA}     & ACC (\%)& 63.6 $\pm 0.9$ & 66.1 $\pm 1.1$ & 71.2 $\pm 0.5$ & 71.9 $\pm 1.1$\\
		& AUC (\%)& 73.5 $\pm 0.3$ & 74.4 $\pm 1.2$ & 78.3 $\pm 0.3$ & 78.2 $\pm 0.6$\\
		& $\Delta_{SP}$ (\%)& 6.0$\pm 1.5$ & 10.9 $\pm 1.9$ & 7.9 $\pm 1.3$ & 10.2 $\pm 2.5$\\
		& $\Delta_{EO}$ (\%)& 6.1 $\pm 1.8$ & 8.8 $\pm 3.0 $ & 17.8 $\pm 2.6$ & 15.9 $\pm 4.0$ \\
		\hline
	\end{tabularx}

	\label{tab:pre}
	\vspace{-3mm}
\end{table}
Various machine learning algorithms such as logistic regression \cite{zafar2015fairness}, SVM \cite{zafar2015fairness}, and MLP \cite{edwards2015censoring} have been reported to have discrimination. The features of the instances may contain proxy variables of the sensitive attribute. It could result in biased predictions. For GNNs, edges in graph can bring linking bias, i.e., the misrepresentation due to the connections of users \cite{mehrabi2019survey}.
It has been proven that the embeddings of nodes within the connected component will be closer after one aggregation in GCN \cite{li2018deeper,wang2020unifying}. Since most of edges are intra-group as Table \ref{tab:stat} shows, embeddings of nodes sharing the same sensitive attribute will be closer after $k$-layer information aggregation. As a result, representations of the nodes may exhibit bias.
Intuitively, similar discrimination also exists in other GNNs that aggregate information of neighborhoods.

To empirically demonstrate the existence of discrimination in GNNs, we make comparisons between the following models:
\begin{itemize}
	\item \textbf{MLP}: A multi-layer perception model trained on $\mathcal{V}_L$.
	\item \textbf{MLP-e}: A MLP model utilizes graph structure by adding embeddings learned by deepwalk to the features.
	\item \textbf{GCN} \cite{kipf2016semi}: A state-of-the-art spectral graph neural network.
	\item \textbf{GAT} \cite{velivckovic2017graph}: A spatial graph neural network which utilizes attention to assign higher weights to more important edges.
\end{itemize}
For each model, we run the experiment 5 times. The classification results and discrimination scores on the test set are reported in Table \ref{tab:pre}.
From the table, we observe that (i) the performance of GCN and GAT are much better than that of MLP, which is as expected because GCN and GAT adopt both node attributes and the graph structure for classification; (ii) Compared with MLP, models utilizing graph structure, i.e., GCN and GAT, perform significantly worse in terms of fairness, which verifies that \textit{bias exists in GNNs and the graph structure could further aggravate the discrimination.}

\label{sec:dis}

\subsection{Problem definition}

Our preliminary analysis verifies that GNNs have severe bias issue. Thus, it is important to develop fair GNNs. Following existing work on fair models \cite{louizos2015variational,feldman2015certifying,beutel2017data,xu2018fairgan}, we focus on the binary class and binary sensitive attribute setting, i.e., both $y$ and $s$ can only be 0 or 1. We leave the extension to multi-class and multi-sensitive attribute setting as a future work.
With the notations given in Section \ref{sec:notations},
the fair GNN problem is formally defined as:

\vspace*{-0.3em}

\newtheorem{problem}{Problem}
\begin{problem}
	Given a graph $\mathcal{G}=(\mathcal{V},\mathcal{E}, \mathbf{X})$, small labeled node set $\mathcal{V}_L \in \mathcal{V}$ with the corresponding labels in $\mathcal{Y}$, and a small set of nodes $\mathcal{V}_S \in \mathcal{V}$ with corresponding sensitive attributes in $\mathcal{S}$, learn a fair GNN for fair node classification, i.e.,
	\begin{equation}
	f(\mathcal{G}, \mathcal{Y}, \mathcal{S}) \rightarrow \hat{\mathcal{Y}}
	\end{equation}
	where $f$ is the function we aim to learn and $\hat{\mathcal{Y}}$ is the set of predicted labels for unlabeled nodes.  $\hat{\mathcal{Y}}$ should maintain high accuracy whilst satisfying the
	fairness criteria such as statistical parity.
\end{problem}

\section{methodology} \label{sec:methdology}

In this section, we give the details of FairGNN. An illustration of the proposed framework is shown in Figure \ref{fig:framework}, which is composed of a GNN classifier $f_{\mathcal{G}}$, a GCN based sensitive attribute estimator $f_E$ and an adversary $f_A$. The classifier $f_{\mathcal{G}}$ takes $\mathcal{G}$ as input for node classification.  The sensitive attribute estimator $f_E$ is to predict the sensitive attributes for nodes whose sensitive attributes are unknown, which paves us a way to adopt adversarial learning to learn fair node representations and to regularize the predictions of $f_{\mathcal{G}}$. Specifically, the adversary $f_A$ aims to predict the known or estimated sensitive attributes by $f_E$ from the node representation learned by $f_\mathcal{G}$; while $f_\mathcal{G}$ aims to learn fair node representations that can fool the adversary $f_A$ to make wrong predictions. We theoretically prove that under mild conditions, such minmax game can guarantee that learned representations are fair.
In addition to make the representations fair, we directly add a regularizer on the predictions of $f_{\mathcal{G}}$ to guarantee that $f_{\mathcal{G}}$ gives fair predictions. Next, we introduce each component in detail along with theoretical proof

\begin{figure}
	\centering
	\includegraphics[width=0.90\linewidth]{framework.png}
	\vskip -1.5em
	\caption{The overall framework of FairGNN.}
	\vskip -1.5em
	\label{fig:framework}
\end{figure}

\subsection{The GNN Classifier $f_{\mathcal{G}}$}
The GNN classifier $f_{\mathcal{G}}$ takes $\mathcal{G}$ as input and predicts node labels. The proposed framework FairGNN is flexible. Any GNNs that follow the structure of Eq.(\ref{eq:GNN_a}) can be used such as GCN~\cite{kipf2016semi} and GAT~\cite{velivckovic2017graph}. Let $f_{\mathcal{G}}^{(k)}$ denote the operation of aggregating and combining the information of node $v$ and its $k$-hop neighborhoods through $k$ layers' iterations in GNN classifier $f_{\mathcal{G}}$. For a GNN with $K$ layers, the representation of node $v$ of the final layer could be written as:
\begin{equation}
\mathbf{h}_v = f_{\mathcal{G}}^{(K)}(\mathbf{x}_v,\mathcal{N}^{(K)}_v),
\end{equation}
where $\mathcal{N}^{(K)}_v$ represents the $K$-hop neighborhoods of $v$. To get the $\hat{y}_v$, i.e., the prediction of node $v$, a linear classification layer is applied to $\mathbf{h}_v$ as:
\begin{equation}
\hat{y}_v = \sigma(\mathbf{h}_v \cdot \mathbf{w}),
\label{eq:y_v}
\end{equation}
where $\mathbf{w} \in \mathbb{R}^{d}$ is the weights of the linear classification layer and $\sigma$ is the sigmoid function. The loss function for training $f_{\mathcal{G}}$ is
\begin{equation} \label{eq:obj_GNN_classifier}
\min_{\theta_{\mathcal{G}}}\mathcal{L}_C = -\frac{1}{|\mathcal{V}_L|}\sum_{v\in \mathcal{V}_L} [y_v \log{\hat{y}_v} + (1-y_v) \log{(1 - \hat{y}_v})],
\end{equation}
where $|\mathcal{V}_L|$ denotes the size of $\mathcal{V}_L$, $\theta_{f_{\mathcal{G}}}$ represents the parameters of $f_{\mathcal{G}}$ and $y_v$ is the groundtruth label of node $v$.

\subsection{Adversarial Debiasing with Estimator $f_E$}
The GNN classifier $f_\mathcal{G}$ can make biased predictions because the learned representations of $f_\mathcal{G}$  exhibit bias due to the node features, graph structure and aggregation mechanism of GNN. One way to make $f_\mathcal{G}$ fair is to eliminate the bias of the final layer representations $\mathbf{h}_{v}$. Recently, adversarial debiasing has been proven to be effective in alleviating the bias of representations  \cite{beutel2017data,edwards2015censoring,liao2019learning,madras2018learning}.
In the general process of adversarial debiasing, an adversary is used to predict sensitive attributes from the representations of the classifier; while the classifier is trained to learn representations to make the adversary unable to predict the sensitive attributes while keep high accuracy in the classification task. Such process requires \textit{abundant data samples with known sensitive attributes} so that we can judge if the adversary can make accurate predictions or not.

However, in practice people are reluctant to share their sensitive attributes, which leads to small size $\mathcal{V}_S$. Lacking of data with labeled sensitive attributes would result in poor improvement in fairness even with adversarial debiasing. Though we have limited nodes with sensitive attributes, i.e., small $\mathcal{V}_S$, generally, nodes with similar sensitive attributes are more likely connected to each other, which makes it possible to accurately predict the sensitive attributes for nodes in $\mathcal{V} - \mathcal{V}_S$ using the graph $\mathcal{G}$ and $\mathcal{V}_S$. Thus, we deploy a graph convolutional network $f_E: \mathcal{G} \rightarrow \mathcal{S}$ to estimate the sensitive attribute of node whose sensitive attribute is unavailable. The large amount of estimated sensitive attributes would greatly
benefit the adversarial debiasing. Note that it is important to use two separate GNNs for node label prediction and sensitive attribute prediction because we aim to learn fair representations $\mathbf{h}_v$ for $f_{\mathcal{G}}$, i.e., $\mathbf{h}_v$ does not contain the sensitive information. The objective function of training $f_E$ is
\begin{equation}
\min_{\theta_E} \mathcal{L}_E = -\frac{1}{|\mathcal{V}_S|} \sum_{v\in \mathcal{V}_S} [s_v \log{\hat{s}_v} + (1-s_v) \log{(1 - \hat{s}_v})],
\label{eq:L_E}
\end{equation}
where $\hat{s}_v$ is the predicted sensitive attribute of node $v \in \mathcal{V}_S$ by  $f_E$ and $\theta_E$ is the set of parameters of $f_E$.

With $f_E$, we could get the estimation of the sensitive attributes $\hat{\mathcal{S}}_u$ of the  nodes $u \in (\mathcal{V} - \mathcal{V}_S)$.
We use $\hat{\mathcal{S}}$ to denote the set of sensitive attributes by combining $\mathcal{S}$ and $\hat{\mathcal{S}}_u$, i.e., $\hat{\mathcal{S}} = \mathcal{S} \cup \hat{\mathcal{S}}_u$.
During the training process, for each node $v \in \mathcal{V}$, the adversary
$f_A$ tries to predict $v$'s sensitive attribute $\hat{s}_v$ given the representation $\mathbf{h}_v$ as $f_A(\mathbf{h}_v)$;
while  $f_{\mathcal{G}}$ aims to learn node representation $\mathbf{h}_v$ that makes the adversary $f_A$ unable to distinguish which sensitive group the node $v$ belong to. This min max game can be written as
\begin{equation}
\begin{aligned}
\min_{\theta_{\mathcal{G}}} \max_{\theta_A} \mathcal{L}_A & = \mathbb{E}_{\mathbf{h} \sim p(\mathbf{h}|\hat{s}=1)} [\log(f_A(\mathbf{h}))] \\
& + \mathbb{E}_{\mathbf{h} \sim p(\mathbf{h}|\hat{s}=0)} [ \log(1-f_A(\mathbf{h}))],
\end{aligned}
\label{eq:A_S}
\end{equation}
where $\mathbf{h} \sim p(\mathbf{h}|\hat{s}=1)$ means sampling a node with sensitive attribute as 1 from $\mathcal{G}$. $\theta_A$ is the parameters of $f_A$.

\textbf{\textbf{Theoretical Analysis.}}
Since the size of $\mathcal{V}_S$ is small, the estimation of sensitive attributes will introduce nonnegligible noise. The noise of the sensitive attributes may influence the adversarial debiasing. Thus, we conduct theoretical analysis to show that sensitive attributes containing noise could help to achieve statistical parity under mild conditions. Next, we give the details of the proof.
\begin{proposition}
	The global minimum of Eq.(\ref{eq:A_S}) is achieved if and only if $p(\mathbf{h}|\hat{s}=1)=p(\mathbf{h}|\hat{s}=0)$, where $\hat{s} \in \hat{\mathcal{S}}$ and $\mathbf{h}$ is final layer representation learned by the $K$-layer GNN classifier $f_{\mathcal{G}}$.
	\label{prop:A_s}
\end{proposition}
\begin{proof}
	According to Proposition 1. in \cite{goodfellow2014generative}, the optimal adversary is
	$f_A^*(\mathbf{h})=\frac{p(\mathbf{h}|\hat{s}=1)}{p(\mathbf{h}|\hat{s}=1) + p(\mathbf{h}|\hat{s}=0)}$. Then the min max game in Eq.(\ref{eq:A_S}) could be reformulated as minimizing this function:
	\begin{equation}
	\begin{aligned}
	C^s & = \mathbb{E}_{\mathbf{h} \sim p(\mathbf{h}|\hat{s}=1)} \big[\log{\frac{p(\mathbf{h}|\hat{s}=1)}{p(\mathbf{h}|\hat{s}=1) + p(\mathbf{h}|\hat{s}=0)}}\big] \\
	& + \mathbb{E}_{\mathbf{h} \sim p(\mathbf{h}|\hat{s}=0)} \big[ \log{\frac{p(\mathbf{h}|\hat{s}=0)}{p(\mathbf{h}|\hat{s}=1) + p(\mathbf{h}|\hat{s}=0)}} \big] \\
	& = -\log(4) + 2\cdot JSD(p(\mathbf{h}|\hat{s}=1)||p(\mathbf{h}|\hat{s}=0).
	\end{aligned}
	\end{equation}
	The Jensen-Shannon divergence between two distributions is non-negative, and become zero if the two distributions are equal. Thus, only if $p(\mathbf{h}|\hat{s}=1)=p(\mathbf{h}|\hat{s}=0)$, the objective function $C^s$ will reach the minimum, which completes our proof.
\end{proof}
\begin{theorem}
	Let $\hat{y}$ denote the prediction of $f_{\mathcal{G}}$.
	Suppose:
	\begin{enumerate}
		\item The estimated sensitive attribute $\hat{s} $ and $\mathbf{h}$ are independent conditioned on true sensitive attribute  $s$, i.e., $p(\hat{s},\mathbf{h}|s) = p(\hat{s}|s)p(\mathbf{h}|s)$;
		\item $p(s=1|\hat{s}=1) \neq p(s=1|\hat{s}=0)$.
	\end{enumerate}
	If Eq.(\ref{eq:A_S}) reaches the global minimum, the GNN classifier $f_{\mathcal{G}}$ will achieve statistical parity, i.e., $p(\hat{y}|s=0)=p(\hat{y}|s=1)$.
	\label{theorem:adv}
\end{theorem}
\begin{proof}
	Under the assumption that $p(\hat{s},\mathbf{h}|s) = p(\hat{s}|s)p(\mathbf{h}|s)$, we could obtain $p(\mathbf{h}|s,\hat{s}) = p(\mathbf{h}|s)$.
	From Proposition \ref{prop:A_s}, we have  $p(\mathbf{h}|\hat{s}=1) = p(\mathbf{h}|\hat{s}=0)$ when the algorithm converges, which is equivalent to $\sum_s p(\mathbf{h},s|\hat{s}=1) = \sum_s p(\mathbf{h},s|\hat{s}=0)$. Together with $p(\mathbf{h}|s,\hat{s}) = p(\mathbf{h}|s)$, we arrive at
	\begin{equation}
	\begin{aligned}

	\sum_s p(\mathbf{h}|s)p(s|\hat{s}=1) = \sum_s p(\mathbf{h}|s)p(s|\hat{s}=0)
	\end{aligned}
	\label{eq:13}
	\end{equation}
	Reordering the terms in Eq.(\ref{eq:13}), we can ge
	\begin{equation}
	\begin{aligned}
	\frac{p(\mathbf{h}|s=1)}{p(\mathbf{h}|s=0)} & = \frac{p(s=0|\hat{s}=1)-p(s=0|\hat{s}=0)}{p(s=1|\hat{s}=0)-p(s=1|\hat{s}=1)} \\
	& = \frac{(1-p(s=1|\hat{s}=1))-(1-p(s=1|\hat{s}=0))}{p(s=1|\hat{s}=0)-p(s=1|\hat{s}=1)} \\
	& = 1
	\end{aligned}
	\label{eq:ratio}
	\end{equation}
	Eq.(\ref{eq:ratio}) shows that at the global minimum $p(\mathbf{h}|s=1) = p(\mathbf{h}|s=1)$ under the assumption  $p(s=1|\hat{s}=1) \neq p(s=1|\hat{s}=0)$.
	Since $\hat{y} = \sigma(\mathbf{h} \cdot \mathbf{w})$, we could get $p(\hat{y}|s=1) = p(\hat{y}|s=0)$. Thus, the statistical parity is achieved when Eq.(\ref{eq:A_S}) converges.
\end{proof}

In our proof, two assumptions are made. For the first assumption, since we use $f_E$ to predict the sensitive attributes $\hat{s}$ and $f_{\mathcal{G}}$ to get the latent representation $\mathbf{h}$, and $f_E$ and $f_{\mathcal{G}}$ doesn't share any parameters, it is generally true that $\hat{s}$ is independent with the representation $\mathbf{h}$, i.e., $p(\hat{s},\mathbf{h}|s) = p(\hat{s}|s)p(\mathbf{h}|s)$.
As for the second assumption, it will be satisfied when we have a reasonable estimator $f_E$, i.e., $f_E$ doesn't give random predictions.

\subsection{Covariance Constraint}
The instability of the training process of adversarial learning is well known \cite{arjovsky2017towards}. In adversarial debiasing, failure to coverage may result in a classifier with discrimination.
To alleviate this issue, we add a covariance constraint~\cite{zafar2015fairness,zafar2017fairness} on the output of $f_{\mathcal{G}}$ to help the model achieve fairness.
The covariance constraint has been explored in~\cite{zafar2015fairness,zafar2017fairness} by minimizing the absolute covariance between users' sensitive attributes and the signed distance from the users' features to the decision boundary for fair linear classifiers .
In our problem, only a small portion of users' sensitive attributes are known and the decision boundary of GNN is hard to obtain. Thus, we propose to minimize the absolute covariance between the noisy sensitive attribute $\hat{s} \in \hat{\mathcal{S}}$ and prediction $\hat{y}$ as
\begin{equation}
\mathcal{L}_R = |\text{Cov}(\hat{s},\hat{y})| = |\mathbb{E}[(\hat{s}-\mathbb{E}(\hat{s}))(\hat{y}-\mathbb{E}(\hat{y}))]|,
\end{equation}
where $|\cdot|$ indicates the absolute value.

\textbf{Theoretical Analysis.}
Since $\mathcal{L}_R$ is the absolute value of covariance between $\hat{y}$ and $\hat{s}$, $\mathcal{L}_R=0$, i.e., the global minimum of $\mathcal{L}_R$, is the prerequisite that $\hat{y}$ and $\hat{s}$ are independent.
Thus, we will show that $\mathcal{L}_R=0$ is the prerequisite of the statistical parity under mild assumption with the following theorem.
\begin{theorem}
	Suppose that $p(\hat{s},\mathbf{h}|s) = p(\hat{s}|s)p(\mathbf{h}|s)$, when $f_{\mathcal{G}}$ satisfies statistical parity, i.e. $\hat{y} \bot s $,  $\hat{y}$ is independent with $\hat{s}$ and $\mathcal{L}_R = 0$.
	\label{Theorem:constraint}
\end{theorem}
\begin{proof}
	Through $p(\hat{s},\mathbf{h}|s) = p(\hat{s}|s)p(\mathbf{h}|s)$, we could get $p(\mathbf{h}|s,\hat{s})=p(\mathbf{h}|s)$. Then, $p(\hat{y}|s,\hat{s})=p(\hat{y}|s)$ could be derived.  When $\hat{y} \bot s$, the distribution $p(\hat{y},\hat{s})$ would be:
	\begin{equation}
	\begin{aligned}
	p(\hat{y},\hat{s}) & = \sum_s p(\hat{y}|s) p(\hat{s},s) =  \sum_s p(\hat{y}) p(\hat{s},s) = p(\hat{y})p(\hat{s}).
	\end{aligned}
	\end{equation}
	Thus, $\hat{y}$ is independent with $\hat{s}$ when the statistical parity is achieved. Then, we can get $\mathcal{L}_R = |\text{Cov}(\hat{s},\hat{y})| = |\mathbb{E}(\hat{s},\hat{y})-\mathbb{E}(\hat{s})\mathbb{E}(\hat{y})|=0$.
\end{proof}
In the proof, we use the first assumption in Theorem \ref{Theorem:constraint}, which is generally valid as discussed previously.

\subsection{Final Objective Function of FairGNN}
We now have $f_{\mathcal{G}}$ for label prediction, $f_E$ for sensitive attribute estimation, $f_A$ with adversarial debiasing to force the node representations learned by $f_{\mathcal{G}}$ are fair, and covariance constraint to further ensure that the prediction of $f_{\mathcal{G}}$ is fair. Combining all these together, the final objective function could be formulated as:
\begin{equation}
\min_{\theta_\mathcal{G},\theta_E} \max_{\theta_A} \mathcal{L}_C + \mathcal{L}_E + \alpha  \mathcal{L}_R - \beta \mathcal{L}_A,
\label{eq:total}
\end{equation}
where $\theta_\mathcal{G}$, $\theta_E$, and $\theta_A$ are the parameters of classifier, estimator, and adversary, respectively. $\alpha$ and $\beta$ are scalars to control the contributions of the covariance constraint and adversarial debiasing.

\begin{algorithm}[t]
	\caption{ Training Algorithm of FairGNN.}
	\label{alg:Framwork}
	\begin{algorithmic}[1]
		\REQUIRE
		$\mathcal{G}=(\mathcal{V},\mathcal{E}, \mathbf{X})$, $\mathcal{Y}$, $\mathcal{S}$, $\alpha$ and $\beta$.
		\ENSURE $f_{\mathcal{G}}$, $f_A$, and $f_E$
		\STATE Initialize $f_E$ by optimizing Eq.(\ref{eq:L_E}) w.r.t $\theta_E$
		\REPEAT
		\STATE Obtain the estimated sensitive attributes with $f_E$
		\STATE Optimize the GNN classifier parameters $\theta_{\mathcal{G}}$, the adversary parameters $\theta_A$, and the estimator parameters $\theta_E$ by Eq.(\ref{eq:total}).

		\UNTIL convergence
		\RETURN $f_{\mathcal{G}}$, $f_A$, and $f_E$
	\end{algorithmic}
\end{algorithm}

\subsection{An Training Algorithm of FairGNN}
The training algorithm of FairGNN is presented in Algorithm \ref{alg:Framwork}. Specially, we first pretrain $f_E$ to ensure it meets the second assumption in Theorem \ref{theorem:adv}. Sequentially, we optimize the whole model with Eq.(\ref{eq:total}) through the ADAM optimizer \cite{kingma2014adam}.
In the training process, we replace the hard labels in $\mathcal{L}_{A}$ with soft labels, i.e., the probability produced by $f_E$, to stabilize the adversarial learning \cite{salimans2016improved}.

\begin{table*}[t]
	\small
	\centering
	\caption{The comparisons of our proposed methods with the baselines.}
	\vskip -1.5em
	\begin{tabularx}{\textwidth}{|p{0.06\textwidth}|p{0.06\textwidth}|XX|XXXXXXX|}
		\hline
		Dataset & Metrics & GCN & GAT & ALFR & ALFR-e & Debias & Debias-e & FCGE & FairGCN & FairGAT \\
		\hline
		\hline

		\multirow{4}{*}{Pokec-z}
		& ACC (\%) & 70.2 $\pm 0.1$ & 70.4 $\pm 0.1$ & 65.4 $\pm 0.3$ & 68.0 $\pm 0.6$ & 65.2 $\pm 0.7$ & 67.5 $\pm 0.7$ & 65.9 $\pm 0.2$& \textbf{70.0} $\pm \mathbf{0.3}$ &  \textbf{70.1} $\pm \mathbf{0.1}$\\
		& AUC (\%)& 77.2 $\pm 0.1$ & 76.7 $\pm 0.1$ & 71.3 $\pm 0.3$ & 74.0 $\pm 0.7$ & 71.4 $\pm 0.6$ & 74.2 $\pm 0.7$ & 71.0 $\pm 0.2$ & \textbf{76.7} $\pm \mathbf{0.2}$ & \textbf{76.5} $\pm \mathbf{0.2}$\\
		& $\Delta_{SP}$ (\%)& 9.9 $\pm 1.1$ & 9.1 $\pm 0.9$ & 2.8 $\pm 0.5$ & 5.8 $\pm 0.4$ & 1.9 $\pm 0.6$ & 4.7 $\pm 1.0$ & 3.1 $\pm 0.5$ & \textbf{0.9} $\pm \mathbf{0.5}$ & \textbf{0.5} $\pm \mathbf{0.3}$\\
		& $\Delta_{EO}$ (\%)& 9.1 $\pm 0.6$ & 8.4 $\pm 0.6$ & 1.1 $\pm 0.4$ & 2.8 $\pm 0.8$ & 1.9 $\pm 0.4$ & 3.0 $\pm 1.4$ & 1.7 $\pm 0.6$ & \textbf{1.7} $\pm \mathbf{0.2}$ & \textbf{0.8} $\pm\mathbf{0.3}$\\
		\hline
		\hline
		\multirow{4}{*}{Pokec-n}
		& ACC (\%) & 70.5 $\pm 0.2$ & 70.3 $\pm 0.1$ & 63.1 $\pm 0.6$ & 66.2 $\pm 0.5$ & 62.6 $\pm 0.9$ & 65.6 $\pm 0.8$ & 64.8 $\pm 0.5$ & \textbf{70.1} $\pm \mathbf{0.2}$ & \textbf{70.0} $\pm \mathbf{0.2}$\\
		& AUC (\%)& 75.1 $\pm 0.2$ & 75.1 $\pm 0.2$ & 67.7 $\pm 0.5$ & 71.9 $\pm 0.3$ & 67.9 $\pm 0.7$ & 71.7 $\pm 0.7$ & 69.5 $\pm 0.4$ & \textbf{74.9} $\pm \mathbf{0.4}$& \textbf{74.9} $\pm \mathbf{0.4}$\\
		& $\Delta_{SP}$ (\%)& 9.6 $\pm 0.9$ & 9.4 $\pm 0.7$ & 3.05 $\pm 0.5$ & 4.1 $\pm 0.5$ & 2.4 $\pm 0.7$ & 3.6 $\pm 0.2$ & 4.1 $\pm 0.8$ & \textbf{0.8} $\pm \mathbf{0.2}$ & \textbf{0.6} $\pm \mathbf{0.3}$\\
		& $\Delta_{EO}$ (\%)& 12.8 $\pm 1.3$ & 12.0 $\pm 1.5$ & 3.9 $\pm 0.6$ & 4.6 $\pm 1.6 $ & 2.6 $\pm 1.0$ & 4.4 $\pm 1.2$ & 5.5 $\pm 0.9$& \textbf{1.1} $\pm \mathbf{0.5}$ & \textbf{0.8} $\pm \mathbf{0.2}$\\
		\hline
		\hline
		\multirow{4}{*}{NBA}
		& ACC (\%) & 71.2 $\pm 0.5$ & 71.9 $\pm 1.1$ & 64.3 $\pm 1.3$ & 66.0 $\pm 0.4$ & 63.1 $\pm 1.1$  & 65.6 $\pm 2.4$ & 66.0 $\pm 1.5$ & \textbf{71.1} $\pm \mathbf{1.0}$ & \textbf{71.5} $\pm \mathbf{0.8}$\\
		& AUC (\%)& 78.3 $\pm 0.3$ & 78.2 $\pm 0.6$ & 71.5 $\pm 0.3$ & 72.9 $\pm 1.0$ & 71.3 $\pm 0.7$ & 72.9 $\pm 1.2$ & 73.6 $\pm 1.5$ & \textbf{77.0} $\pm \mathbf{0.3}$ & \textbf{77.5} $\pm \mathbf{0.7}$\\
		& $\Delta_{SP}$ (\%)& 7.9 $\pm 1.3$ & 10.2 $\pm 2.5 $ & 2.3 $\pm 0.9$ & 4.7 $\pm 1.8$ & 2.5 $\pm 1.5$ & 5.3 $\pm 0.9$ & 2.9 $\pm 1.0$ & \textbf{1.0} $\pm \mathbf{0.5} $ & \textbf{0.7} $\pm \mathbf{0.5}$\\
		& $\Delta_{EO} (\%)$ & 17.8 $\pm 2.6$ & 15.9 $\pm 4.0$ & 3.2 $\pm 1.5$ & 4.7 $\pm 1.7$ & 3.1 $\pm 1.9$ & 3.1 $\pm 1.3$ & 3.0 $\pm 1.2$ & \textbf{1.2} $\pm \mathbf{0.4}$ & \textbf{0.7} $\pm \mathbf{0.3}$\\
		\hline
	\end{tabularx}

	\label{tab:results}
	\vskip -1em
\end{table*}

\section{experiments} \label{sec:experiments}
In this section, we conduct experiments to show the effectiveness of FairGNN for fair node classification. In particular, we aim to answer the following questions
\begin{itemize}
	\item \textbf{RQ1} Can the proposed FairGNN reduce the bias of GNNs while maintain high accuracy?

	\item \textbf{RQ2} How do the sensitive attribute estimator, adversarial loss, and covariance constraint affect FairGNN?
	\item \textbf{RQ3} Is FairGNN effective when different amount of sensitive attributes or labels are provided in the training set?

\end{itemize}
We use the same datasets introduced in Sec.~\ref{sec:datasets} for all the experiments. Next, we will begin by introducing compared methods.

\subsection{Compared Methods}
We compare our proposed framework with GCN, GAT, and the following representative and state-of-the-art methods for fair classification and fair graph embedding learning:
\begin{itemize}[leftmargin=*]

	\item \textbf{ALFR} \cite{edwards2015censoring}: This is a pre-processing method. A discriminator is applied to remove the sensitive information in the representations produced by a MLP-based autoencoder. Then, linear classifier is trained on the debiased representations.
	\item \textbf{ALFR-e}: To utilize the graph structure information, ALFR-e concatenates
	the graph embeddings learned by deepwalk \cite{perozzi2014deepwalk} with the user features in the ALFR.
	\item \textbf{Debias} \cite{zhang2018mitigating}: This is an in-processing fair classification method.  It directly applies an discriminator on the estimated probability of classifier $\eta: \mathbf{x} \rightarrow \mathbb{R}$. It would make the probability distribution $p(\eta(\mathbf{x})|s=0)$ closer to $p(\eta(\mathbf{x})|s=1)$.
	\item \textbf{Debias-e}: Similar to the ALFR-e, we also add the deepwalk embeddings to the features used in Debias.
	\item \textbf{FCGE} \cite{bose2019compositional}: FCGE is proposed to learn fair node embeddings  in graph without node features through edge prediction. The sensitive information in the embeddings is filtered by discriminators.
\end{itemize}
ALFR and ALFR-e are trained with features of all the users $\mathcal{V}$, labels of $\mathcal{V}_L$, and the sensitive attributes of $\mathcal{V}_S$ for fair classification.
Debis and Debias-e require the sensitive attributes of labeled nodes, which is on contrary with our setting that $\mathcal{V}_L$ could have no overlap with $\mathcal{V}_S$. Thus, we use the estimated labels of $\mathcal{V}_S$, features of $\mathcal{V}_L$, and labels of $\mathcal{V}_L$ to train Debias and Debias-e. FCGE utilizes $\mathcal{G}$, labels of $\mathcal{V}_L$, and sensitive attributes of $\mathcal{V}_S$.

For FairGNN, we deploy a one hidden layer GCN for $f_E$. The hidden dimension is set as 128. We use a linear classifier for $f_A$. To verify that our framework is useful for various GNNs, we adopt both GCN and GAT as the backbone of the FairGNN classifier $f_{\mathcal{G}}$, which are named as \textbf{FairGCN} and \textbf{FairGAT}. In FairGCN, the GCN classifier contains one hidden layer with dimension 128. The GAT
classifier in FairGAT also contains two layers in total. We set the number of heads as 1. The dimensions of the GAT classifiers' hidden layer for Pokec-z, Pokec-n and NBA are 64, 64 and 32, respectively.

\subsection{Fair Classification on Graph}

To answer \textbf{RQ1}, we evaluate our proposed FairGNN in terms of fairness and classification performance.
$\Delta_{SP}$ and $\Delta_{EO}$ are used to show the discrimination level, which are introduced in Section~\ref{sec:fairness}. The smaller $\Delta_{SP}$ and $\Delta_{EO}$ are, the more fair the classifier is. Accuracy (ACC) and ROC AUC score are used to evaluate the classification performance. For all the models, we tune the hyperparameters on the training set via cross validation.
For FairGCN, we set $\alpha$ to 100 and $\beta$ to 1.
For FairGAT, $\alpha$ is 2 and $\beta$ is 0.1. More details about hyperparameter selection will be discussed in Sec~\ref{sec:parameter_sensitivity}. All the experiments are conducted 5 times.
The mean and standard deviations for all the models on the three datasets are reported in Table \ref{tab:results}. From the table, we make the following observations:
\begin{itemize}[leftmargin=*]
	\item Compared with GCN and GAT, the general fair classification methods and graph embeddings learning method show poor performance in classification even with the help of graph information, while FairGCN and FairGAT perform very close to the based GNNs. This suggests the necessity of investigating fair classification algorithms on GNNs for accurate predictions;
	\item Under the condition of limited sensitive information, baselines show obvious bias and the ones utilizing graph information are even worse. On the contrary, our proposed models obtain $\Delta_{SP}$ and $\Delta_{EO}$ that are close to 0, which indicates that the discrimination is basically eliminated; and
	\item FairGAT is slightly better than FairGCN in Fairness. This is reasonable because the learnable edge coefficients in GAT could be helpful to reduce the weights of the edges that bring bias.
\end{itemize}
These observations demonstrate the effectiveness of our proposed framework in making fair and accurate predictions.

\subsection{Ablation Study}
To answer \textbf{RQ2}, we conduct ablation studies to understand the impacts of $f_E$, adversarial loss, and covariance constraint.

\subsubsection{Impact of $f_E$}  In our proposed framework, a GCN estimator is deployed to predict  sensitive attributes for adversarial debiasing.
To show the importance of the GCN estimator, we analyze it from two aspects.
Firstly, to demonstrate the effectiveness of the noisy sensitive attributes, we eliminate the estimator and only use the provided sensitive attributes $\mathcal{S}$ to get a variant denoted as FairGNN$\backslash$E. Secondly, to investigate how a weaker estimator would influence the fair classification, we train a variant $\text{FairGNN}_{MLP}$ by using MLP as the estimator. Hyperparameters of
these variants are determined by cross validation with gird search. Specifically, we vary $\alpha$ and $\beta$ among $\{0.0001,0.001,0.1,1\}$ and $\{1, 2, 5, 10, 20, 50, 100\}$, respectively.
For each variant, the experiments are conducted 5 times. The average performance of fairness in terms of $\Delta_{SP}$ and node classificaiton in terms of AUC on Pockec-z are presented in Fig. \ref{fig:abl}(a) and (b), respectively. We only show the results on Pockec-z as we have similar observations on the other datasets. From the figures, we make the following observations:
\begin{itemize}[leftmargin=*]
	\item The $\Delta_{SP}$ score of FairGNN$\backslash$E is much larger than that of FairGNN. which is because the provided sensitive attributes are inadequate. This shows that $f_E$ plays an important role in FairGNN; and
	\item The performance of sensitive attribute prediction in terms of AUC for MLP estimator is 0.69, which is much lower than that of GCN estimator, which is 0.80. Though $\text{FairGNN}_{MLP}$ adopts a much weaker estimator than FairGNN, the performance in terms of fairness is slightly worse than FairGNN. This aligns with our theoretical analysis that $f_E$ doesn't need to be very accurate. However, the marginal differences still indicate that too much noise in sensitive attributes may still slightly affect the fairness.
\end{itemize}

\begin{figure}[t]
	\centering
	\begin{subfigure}{0.49\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{SP_abl.pdf}
		\vskip -0.5em
		\caption{$\Delta_{SP}$}
		\label{fig:abl_sp}
	\end{subfigure}

	\begin{subfigure}{0.49\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{roc_abl.pdf}
		\vskip -0.5em
		\caption{ROC AUC}
		\label{fig:abl_roc}
	\end{subfigure}
	\vspace{-1.5em}
	\caption{Comparisons between FairGNN and its variants.}
	\vspace{-1em}
	\label{fig:abl}
\end{figure}

\subsubsection{Impacts of the adversarial debiasing and covariance constraint}
To demonstrate the effects of the adversarial loss and covariance constraint, we train two variants of FairGNN, i.e., FairGNN$\backslash$A and FairGNN$\backslash$C, where FairGNN$\backslash$A means FairGNN without the adversarial loss, and FairGNN$\backslash$C means FiarGNN without covariance constraint. Similarly, for each variant, we run the experiment 5 times on Pokec-z and the average performances are shown in Figure \ref{fig:abl}. From the figure, we observe:
\begin{itemize}[leftmargin=*]
	\item The $\Delta_{SP}$ scores for both FairGNN$\backslash$C and FairGNN$\backslash$A are much smaller than that of GNNs in Figure \ref{fig:abl}, which shows that covariance constraint and adversarial debiasing can improve fairness, respectively; and
	\item The $\Delta_{SP}$ scores for both FairGNN$\backslash$C and FairGNN$\backslash$A are much larger than that of FairGNN, which implies that using both covariance constraint and adversarial debiasing can achieve better fairness. This is because they regularize the GNN from two different perspectives, i.e., adversarial debiasing regularizes on the node representations while covariance cosntraint is directly on the predictions for fair classification.
\end{itemize}

\subsection{Impacts of Sizes of $\mathcal{V}_S$ and $\mathcal{V}_L$}
To answer \textbf{RQ3}, we study the impacts of the sizes of $\mathcal{V}_S$ and  $\mathcal{V}_L$ on FairGAT.
We set $\alpha=0.1$ and $\beta=2$ based on cross validation.
We vary $|\mathcal{V}_S|$ as $\{200,600,1000,1400,1800,2200,2600,3000\}$. Each experiment is conducted 5 times and the average results on Pokec-z with comparison to FairGAT$\backslash$E and ALFR-e are  shown in Fig.~\ref{fig:sens_compare}.
From the figure, we observe that:
(i) Generally, both FairGAT$\backslash$E and ALFR-e have high discrimination scores when  $|\mathcal{V}_S|$ is small. They need plenty of data with sensitive attributes to become effective. FairGAT could get very low $\Delta_{SP}$ even when  $|\mathcal{V}_S|$ is as small as 200. This implies that FairGAT is insensitive to the size of data with sensitive attributes, which is because we have $f_E$ to estimate the sensitive attributes. Though extremely small $|\mathcal{V}_S|$ would lead to a weak $f_E$, we still have similar $\Delta_{SP}$ score as that when $\mathcal{V}_S$ is large. This verifies  our theoretical analysis that we can achieve good fairness with a reasonable $f_E$;
(ii) FairGAT$\backslash$E and ALFR-e decrease slightly
in classification performance with the increasing of the size of $\mathcal{V}_S$, which is because more data with sensitive attribute would lead to a stricter regularization. In the contrary, FairGAT keeps high classification performance and even perform slightly better with more sensitive attributes. This is because the size of sensitive attributes $\hat{\mathcal{S}}$ used for training FairGAT are fixed to the size of $\mathcal{V}$, and less noise in  the estimation of the sensitive attributes is helpful to better learn representations for classification.

Similarly, we vary $|\mathcal{V}_L|$ as $\{500, 1000, 1500, 2000\}$ and each experiment is run for 5 times. The average results on Pokec-z are reported in Figure \ref{fig:sens_sp_l}.
We only report the results on Pokec-z as we have simialr observations on other datasets.
From the figure, we observe that: FairGAT consistently shows effectiveness in eliminating discrimination. The drop in classification performance is marginal.
This demonstrates that our proposed method could achieve fairness while keep high accuracy in general  scenarios which correspond to various sizes of $\mathcal{V}_S$ and $\mathcal{V}_L$.

\begin{figure}[t]
	\centering
	\begin{subfigure}{0.45\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{sens_sp.png}
		\vspace{-5mm}
		\caption{$\Delta_{SP}$}
		\label{fig:sens_sp}
	\end{subfigure}
	\begin{subfigure}{0.45\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{sens_roc.png}
		\vspace{-5mm}
		\caption{ROC AUC}
		\label{fig:sens_roc}
	\end{subfigure}
	\vskip -1.2em
	\caption{Impacts of the size of $\mathcal{V}_S$ to FairGAT.}
	\vskip -1em
	\vspace{-0.5em}
	\label{fig:sens_compare}
\end{figure}

\begin{figure}[t]

	\centering
	\begin{subfigure}{0.45\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{label_sp.png}
		\vspace{-5mm}
		\caption{$\Delta_{SP}$}
		\label{fig:sens_sp_l}
	\end{subfigure}
	\begin{subfigure}{0.45\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{label_roc.png}
		\vspace{-5mm}
		\caption{ROC AUC}
		\label{fig:sens_roc_l}
	\end{subfigure}
	\vskip -1.2em
	\caption{Impacts of the size of $\mathcal{V}_L$ to FairGAT.}
	\vspace{-1em}
	\label{fig:sens_compare_l}
\end{figure}

\subsection{Parameter Sensitivity} \label{sec:parameter_sensitivity}
There are two important hyperparameters in our proposed model, .i.e., $\alpha$ controlling the influence of the adversary to the GNN classifier, while $\beta$ controlling the contribution of the covariance constraint to ensure fairness.
To investigate the parameter sensitivity and find the ranges that achieve high accuracy with low discrimination score, we train FairGAT models on Pokec-z with various hyperparameters.
More specifically, we alter the values of $\alpha$ and $\beta$ among $\{0.0001,0.001,0.01,0.1,1\}$ and $\{1,2,5,10,20,50,100\}$.
The results are presented in Figure \ref{fig:para}.
From Figure \ref{fig:para} (\subref{fig:abl_roc}), we can find that when $\alpha \leq 0.01$ and $\beta \leq 20$ the classification performance is almostly unaffected. Once $\alpha$ and $\beta$ is too large, the classifier's performance will decay rapidly. The impacts of the hyperparameters to the discrimination score are presented in Figure \ref{fig:para} (\subref{fig:abl_sp}). When we increase the value of $\alpha$, $\Delta_{SP}$ will
firstly decrease as expected. Then, it would increase when the value of $\alpha$ is too large. Because it would be difficult to optimize the GNN classifier to the global minimum when the contribution of the adversary is extremely high. As for $\beta$, the discrimination score would consistently reduce when we increase its value. Combining the two figures, we could determine that when $\alpha \in [0.001, 0.01]$ and $\beta \in [5,20]$, the GNN classifier achieves fairness and maintains high node classification accuracy.
\begin{figure}[t]
	\centering
	\begin{subfigure}{0.45\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{parameter_sp.pdf}
		\vspace{-5mm}
		\caption{$\Delta_{SP}$ (\%)}
		\label{fig:abl_sp}
	\end{subfigure}
	\begin{subfigure}{0.45\columnwidth}
		\centering
		\includegraphics[width=\linewidth]{parameter_roc.pdf}
		\vspace{-5mm}
		\caption{ROC AUC (\%)}
		\label{fig:abl_roc}
	\end{subfigure} ~~~
	\vskip -1.3em
	\caption{Parameter sensitivity analysis.}

	\vskip -1.5em
	\label{fig:para}
\end{figure}

\section{conclusion and future work} \label{sec:conclusion}
In this paper, we study a novel problem of fair GNN learning with limited sensitive information. We empirically demonstrate that GNNs exhibit severe bias. We propose a novel and flexible framework FairGNN which is able to significantly alleviate the bias issue of GNNs meanwhile  maintain high performance on node classification. FairGNN adopts a sensitive attribute estimator to alleviate the issue of lacking sensitive attribute information. With the estimated sensitive attributes, FairGNN designs adversarial debiasing and covariance constraint to regularize the GNN to have fair node representations and predictions, respectively. We theoretically show that FairGNN can reduce the bias.
Experiment results on real-world datasets demonstrate the effectiveness of the proposed framework in terms of both fairness and classification performance.
There are several interesting directions which need further investigation. First, we assume the provided  sensitive attributes are clean. However, for some applications in social media, users might randomly input sensitive attributes such as gender due to privacy concern. Thus, we will extend FairGNN to deal with limited and inaccurate sensitive information. Second, the experiments show that the edges are possible to bring bias. Thus, we will also explore methods which add/delete links in graphs to improve the fairness and classification performance of FairGNN.

\bibliographystyle{ACM-Reference-Format}