# Debiased Graph Neural Networks with
Agnostic Label Selection Bias

Shaohua Fan, Xiao Wang, Chuan Shi2, Kun Kuang, Nian Liu, Bai Wang 2
Corresponding author.S. Fan, N. Liu and B. Wang are with the Department of
Computer Science, Beijing University of Posts and Telecommunications, Beijing,
China. E-mail: {fanshaohua, nianliu<EMAIL_ADDRESS>Wang and C. Shi
are with the Key Laboratory of Trustworthy Distributed Computing and Service
(MoE), Beijing University of Posts and Telecommunications, Beijing, China and
the Peng Cheng Laboratory, Shenzhen, China. E-mail:
<EMAIL_ADDRESS>Kuang is with the College of Computer
Science and Technology, Zhejiang University, Zhejiang, China. E-mail:
<EMAIL_ADDRESS>received 05 Feb. 2021; revised 02 Sep. 2021;
accepted 29 Dec 2021.

###### Abstract

Most existing Graph Neural Networks (GNNs) are proposed without considering
the selection bias in data, i.e., the inconsistent distribution between the
training set with test set. In reality, the test data is not even available
during the training process, making selection bias agnostic. Training GNNs
with biased selected nodes leads to significant parameter estimation bias and
greatly impacts the generalization ability on test nodes. In this paper, we
first present an experimental investigation, which clearly shows that the
selection bias drastically hinders the generalization ability of GNNs, and
theoretically prove that the selection bias will cause the biased estimation
on GNN parameters. Then to remove the bias in GNN estimation, we propose a
novel Debiased Graph Neural Networks (DGNN) with a differentiated
decorrelation regularizer. The differentiated decorrelation regularizer
estimates a sample weight for each labeled node such that the spurious
correlation of learned embeddings could be eliminated. We analyze the
regularizer in causal view and it motivates us to differentiate the weights of
the variables based on their contribution on the confounding bias. Then, these
sample weights are used for reweighting GNNs to eliminate the estimation bias,
thus help to improve the stability of prediction on unknown test nodes.
Comprehensive experiments are conducted on several challenging graph datasets
with two kinds of label selection biases. The results well verify that our
proposed model outperforms the state-of-the-art methods and DGNN is a flexible
framework to enhance existing GNNs.

###### Index Terms:

Graph Neural Networks, Casual Inference, Selection Bias.

##  I Introduction

Graph Neural Networks (GNNs) are powerful deep learning algorithms on graphs
with various applications [1, 2, 3, 4]. Existing GNNs mainly learn a node
embedding through aggregating the features from its neighbors, and such
message-passing framework is supervised by the node label in an end-to-end
manner. During this training procedure, GNNs will effectively learn the
correlation between the structure patterns and node features with the node
labels, so that GNNs are capable of learning the embeddings of new nodes and
inferring their labels.

One basic requirement of GNNs making precise predictions on unseen test nodes
is that the distribution of labeled training nodes and test nodes is the same,
i.e., the structure and feature of labeled training and test nodes follow the
similar pattern, so that the learned correlation between the current graph and
labels can be well generalized to the new nodes. However, in reality, there
are two inevitable issues. (1) Because it is difficult to control the graph
collection in an unbiased manner, the relationship between the collected real-
world graph and the labeled nodes is inevitably biased. Training on such a
graph will cause biased correlation with node labels. Taking a scientist
collaboration network as an example, if most scientists with “machine
learning” (ML) label collaborate with those with “computer vision” (CV) label,
existing GNNs may learn spurious correlation, i.e., scientists who cooperate
with CV scientists are ML scientists. If a new ML scientist only connects with
ML scientists or the scientists in other areas, it will be probably
misclassified. (2) The test nodes in the real scenario are usually not
available in the training phase, implying that the distribution of new nodes
is agnostic. Once the distribution is inconsistent with that in the training
nodes, the performance of all the current GNNs will be hindered. Even transfer
learning is able to solve the distribution shift problem, however, it still
needs the prior of test distribution, which actually cannot be obtained
beforehand. Therefore, the agnostic label selection bias greatly affects the
generalization ability of GNNs on unknown test data.

In order to observe selection bias in real graph data, we conduct an
experimental investigation to validate the effect of selection bias on GNNs
(see Section II-A). We select training nodes with different biased degrees for
each dataset, making the distribution of training nodes and test nodes
inconsistent. The results clearly show that selection bias drastically hinders
the performance of GNNs on unseen test nodes. Moreover, with heavier bias, the
performance drops more. Further, we theoretically analyze how the data
selection bias results in the estimation bias in GNN parameters (see Section
II-B). Based on the stable learning technique [5], we can assume that the
learned embeddings consist of two parts: stable variables and unstable
variables. The data selection bias will cause spurious correlation between
these two kinds of variables. Thereby we prove that with the inevitable model
misspecification, the spurious correlation will further cause the parameter
estimation bias. Once the weakness of the current GNNs with selection bias is
identified, one natural question is “how to remove the estimation bias in
GNNs?”

In this paper, we propose a novel Debiased Graph Neural Network (DGNN)
framework for stable graph learning by jointly optimizing a differentiated
decorrelation regularizer and a weighted GNN model. Specifically, the
differentiated decorrelation regularizer is able to learn a set of sample
weights under differentiated variable weights, so that the spurious
correlation between stable and unstable variables would be greatly eliminated.
Based on the causal view analysis of the decorrelation regularizer, we
theoretically prove that the weights of variables can be differentiated by the
regression coefficients. Compared with existing decorrelation methods [5, 6],
the proposed regularizer is able to remove the spurious correlation while
maintaining a higher effective sample size and requiring less prior knowledge.
Moreover, to better combine the decorrelation regularizer with existing GNN
architecture, the theoretical result shows that adding the regularizer to the
embeddings learned by the penultimate layer could be both theoretically sound
and flexible. Then the sample weights learned by the decorrelation regularizer
are used to reweight the GNN loss so that the parameter estimation could be
unbiased.

In summary, the contributions of this paper are three-fold: i) We investigate
a new problem of learning GNNs with agnostic label selection bias. The problem
setting is general and practical for real applications. ii) We bring the idea
of variable decorrelation into GNNs to relieve bias influence on model
learning and propose a general framework DGNN that could be adopted to various
GNNs. iii) We conduct experiments on real-world graph benchmarks with two
kinds of agnostic label selection biases, and the experimental results
demonstrate the effectiveness and flexibility of our model.

##  II Effect of Label Selection Bias on GNNs

In this section, we first summarize the main notations used in this paper in
Table I and then formulate our target problem:

######  Problem 1 (Semi-supervised Learning on Graph with Agnostic Label
Selection Bias).

Given a training graph
$\mathcal{G}_{train}=\\{\mathbf{A}_{train},\mathbf{X}_{train},\mathbf{Y}_{train}\\}$,
where $\mathbf{A}_{train}\in\mathbb{R}^{N\times N}$ ($N$ nodes) represents the
adjacency matrix, $\mathbf{X}_{train}\in\mathbb{R}^{N\times D}$ ($D$ features)
refers to the node feature vectors and
$\mathbf{Y}_{train}\in\mathbb{R}^{n\times C}$ ($n$ labeled nodes, $C$ classes)
refers to the available labels for training ($n\ll N$), the task is to learn a
GNN $g_{\theta}(\cdot)$ with parameter $\theta$ to precisely predict the label
of nodes on test graph
$\mathcal{G}_{test}=\\{\mathbf{A}_{test},\mathbf{X}_{test},\mathbf{Y}_{test}\\}$,
where distribution $\Psi(\mathcal{G}_{train})\neq\Psi(\mathcal{G}_{test})$.

TABLE I: Glossary of Notations. Notation | Description
---|---
$\mathcal{G}_{train}$ | Training graph
$\mathcal{G}_{test}$ | Test graph
$\mathbf{A}_{train/test}$ | The adjacency matrix of $\mathcal{G}_{train}$ or $\mathcal{G}_{test}$
$\mathbf{X}_{train/test}$ | The node feature vectors of $\mathcal{G}_{train}$ or $\mathcal{G}_{test}$
$\mathbf{Y}_{train/test}$ | The node label vectors of $\mathcal{G}_{train}$ or $\mathcal{G}_{test}$
$\mathbf{H}$/$\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})$ | Node embeddings matrix learned by GNNs
$\mathbf{S}$ | The stable variables in $\mathbf{H}$
$\mathbf{V}$ | The unstable variables in $\mathbf{H}$
$\bar{\mathbf{S}}$ | The latent stable variables to generate $\mathbf{Y}$
$\bar{\mathbf{V}}$ | The unstable variables to generate label $\mathbf{Y}$
$\tilde{\beta}_{S}$ | The linear coefficients for $\mathbf{S}$
$\tilde{\beta}_{V}$ | The linear coefficients for $\mathbf{V}$
$g(\cdot)$ | The non-linear transformation for stable variables $\mathbf{S}$
$\beta_{S}$ | The linear coefficients for $\bar{\mathbf{S}}$
$\beta_{V}$ | The linear coefficients for $\bar{\mathbf{V}}$
$T$ | Treatment variable
$\mathbf{X}$ | Confounders
$\mathbf{w}$ | Sample weights
$\alpha$ | Variable weights in DVD term
$\xi$ | The linear coefficients for confounders $\mathbf{X}$
$\gamma$ | The linear coefficient for treatment $T$
$Y_{i}(t)$ | The potential outcome of sample $i$ with treatment $t$

###  II-A Experimental Investigation

(a) Cora

(b) Citeseer

(c) Pubmed

Figure 1: Effect of selection bias on GCN and GAT.

We conduct an experimental investigation to examine whether the existing GNNs
are sensitive to the distribution shifts caused by the label selection bias.
One motivating example is that due to the research interests of the researcher
that they are more likely to label the interdisciplinary documents that cite
more papers from different subjects in a citation network, while testing may
be conducted on nodes with any neighborhood distribution. Based on this
example, the main idea of the experimental investigation is that we will
perform two representative GNNs: GCN [2] and GAT [3] on three widely used
graph datasets: Cora, Citeseer, Pubmed [7] with different bias degrees. If the
performance drops sharply comparing with the scenarios without selection bias,
this will demonstrate that GNNs cannot generalize well in selection bias
settings.

To simulate the agnostic selection bias scenario, we first follow the
inductive setting in [8] that masks the validation and test nodes as the
training graph $\mathcal{G}_{train}$ in the training phase, and then infer the
labels of validation and test nodes with whole graph $\mathcal{G}_{test}$. In
this way, the distribution of test node can be considered agnostic. Following
[9], we design a biased label selection method on training graph
$\mathcal{G}_{train}$. The selection variable $e$ is introduced to control
whether the node will be selected as labeled nodes, where $e=1$ means selected
and 0 otherwise. For node $i$, we calculate its neighbor distribution ratio:
$r_{i}=|\\{j|j\in\mathcal{N}_{i},y_{j}\neq y_{i}\\}|/|\mathcal{N}_{i}|$, where
$\mathcal{N}_{i}$ is neighborhood of node $i$ in $\mathcal{G}_{train}$ and
$y_{j}\neq y_{i}$ means the label of central node $i$ is not the same as the
label of its neighborhood node $j$. And $r_{i}$ measures the difference
between the label of central node $i$ with the labels of its neighborhood.
Then we average all the nodes’ $r$ to get a threshold $t$. For each node, the
probability to be selected is:
$P(e_{i}=1|r_{i})=\left\\{\begin{array}[]{rcl}\epsilon&{r_{i}\geq t}\\\
1-\epsilon&{r_{i}<t}\\\ \end{array}\right.$, where $\epsilon\in(0.5,1)$ is
used to control the degree of selection bias and a larger $\epsilon$ means a
heavier bias. We set $\epsilon$ as {0.7, 0.8, 0.9} to get three bias degrees
for each dataset, termed as Light, Medium, Heavy, respectively. We select 20
nodes for each class for training and the validation and test nodes are the
same as [10]. Furthermore, we take the unbiased datasets as baselines, where
the labeled nodes are selected randomly.

Figure 1 is the results of GCN, GAT and our proposed method, GCN/GAT-DVD, on
these datasets with four bias degrees. We can find that: i) Compared with the
unbiased scenario, when performing GCN/GAT on biased datasets, they suffer
from serious performance decrease, indicating that selection bias greatly
affects the GNNs’ performance. ii) All lines decrease monotonically with the
increase of bias degree, demonstrating that heavier biases will cause larger
performance reduction. iii) GCN/GAT-DVD outperforms the corresponding base
models (i.e., GCN/GAT) consistently and achieves larger improvements in
heavier bias degree scenarios, indicating that our proposed method could
relieve the effect of selection bias.

###  II-B Theoretical Analysis

The above experiment empirically verifies the effect of selection bias on
GNNs. Here we theoretically analyze the effect of selection bias on estimating
the parameters in GNNs. First, because biased labeled nodes have biased
neighborhood structure and features, GNNs will encode this biased information
into the node embeddings, which is validated by the experimental
investigation. Based on stable learning technique [5], we make the following
assumption:

######  Assumption 1.

The node embeddings learned by GNNs for each node can be decomposed as
$\mathbf{H}=\\{\mathbf{S},\mathbf{V}\\}$, where $\mathbf{S}$ represents the
stable variables and $\mathbf{V}$ represents the unstable variables.
Specifically, for both training and test environments,
$\mathbb{E}(\mathbf{Y}|\mathbf{S}=s,\mathbf{V}=v)=\mathbb{E}(\mathbf{Y}|\mathbf{S}=s)$.

Under Assumption 1, the distribution shift between training set and test set
is mainly induced by the variation in the joint distribution over
($\mathbf{S}$, $\mathbf{V}$), i.e.,
$\mathbb{P}(\mathbf{S}_{train},\mathbf{V}_{train})\neq\mathbb{P}(\mathbf{S}_{test},\mathbf{V}_{test})$.
However, there is an invariant relationship between stable variables
$\mathbf{S}$ and outcome $\mathbf{Y}$ in both training and test environments,
which can be expressed as
$\mathbb{P}(\mathbf{Y}_{train}|\mathbf{S}_{train})=\mathbb{P}(\mathbf{Y}_{test}|\mathbf{S}_{test})$.
Assumption 1 can be guaranteed by $\mathbf{Y}\bot\mathbf{V}|\mathbf{S}$. Thus,
one can solve the stable prediction problem by developing a function
$g(\cdot)$ based on $\mathbf{S}$. However, one can hardly identify such
variables in GNNs.

Without loss of generality, we take $\mathbf{Y}$ as a continuous variable for
analysis and have the following assumption:

######  Assumption 2.

The true generation process of target variable $\mathbf{Y}$ contains not only
the linear combination of stable variables $\mathbf{S}$, but also the
nonlinear transformation of stable variables.

Based on the above assumptions, we formalize the label generation process as
follows:

$\displaystyle\mathbf{Y}=f(\mathbf{X},\mathbf{A})+\varepsilon$
$\displaystyle=\bar{\mathbf{S}}\beta_{S}+\bar{\mathbf{V}}\beta_{V}+g(\bar{\mathbf{S}})+\varepsilon,$
(1)

where $\bar{\mathbf{S}}$ and $\bar{\mathbf{V}}$ are latent stable variables
and unstable variables to generate label $\mathbf{Y}$, which can be learned by
GNNs from raw graph data, $\beta_{S}$ and $\beta_{V}$ are the corresponding
linear coefficients and they represent the effect of each latent variable on
outcome $\mathbf{Y}$, $\varepsilon$ is the independent random noise, and
$g(\cdot)$ is the nonlinear transformation function of stable variables.
According to Assumption 1, we know that coefficients of unstable variables
$\bar{\mathbf{V}}$ are actually $\mathbf{0}$ (i.e., $\beta_{V}=\mathbf{0}$).

For a classical GNN model with a linear regression predictor, its prediction
function can be formulated as:

$\hat{\mathbf{Y}}=\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})_{S}\hat{\beta}_{S}+\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})_{V}\hat{\beta}_{V}+\varepsilon,$
(2)

where
$\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})\in\mathbb{R}^{N\times p}$
denotes the node embeddings learned by a GNN, such as GCN and GAT, the output
variables of $\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})$ can be
decomposed as stable variables
$\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})_{S}\in\mathbf{R}^{N\times
m}$ and unstable variables
$\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})_{V}\in\mathbf{R}^{N\times
q}$ ($m+q=p$) corresponding to $\bar{\mathbf{S}}$ and $\bar{\mathbf{V}}$ in
Eq. (1). Compared with Eq. (1), we can find that the parameters of GNN model
could be unbiasedly estimated if the nonlinear term $g(\bar{\mathbf{S}})=0$
(i.e., there does not exist any non-linear relationship in the label
generation process that cannot be learned by GNNs), because the GNN model in
Eq. (2) will have the same label generation mechanism as Eq. (1). However, as
common used GNNs only have several layers which may limit their nonlinear
power and the real-world graph data is far more complicated, it is reasonable
to assume that there is a nonlinear term $g(\bar{\mathbf{S}})\neq 0$ that
cannot be fitted by the GNNs. Under this assumption, next, we taking a vanilla
GCN [2] as an example to illustrate how the distribution shift will induce
parameter estimation bias. A two-layer GCN can be formulated as
$\hat{\mathbf{A}}\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})\mathbf{W}^{(1)}$,
where $\hat{\mathbf{A}}$ is the normalized adjacency matrix, $\mathbf{W}$ is
the transformation matrix at each layer and $\sigma(\cdot)$ is the Relu
activation function. We decompose GCN as two parts: one is embedding learning
part $\hat{\mathbf{A}}\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})$,
which can be decomposed as
$[\mathbf{S}^{\mathrm{T}},\mathbf{V}^{\mathrm{T}}]$, corresponding to
$\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})_{S}$ and
$\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})_{V}$ in Eq. (2), and the
other part is $\mathbf{W}^{(1)}$, where the learned parameters can be
decomposed as $[\tilde{\beta}_{S},\tilde{\beta}_{V}]$, corresponding to
$[\hat{\beta}_{S},\hat{\beta}_{V}]$ in Eq. (2). We aim at minimizing the
least-square loss:

$\mathcal{L}_{GCN}=\sum^{n}_{i=1}(\mathbf{S}_{i}^{\mathrm{T}}\tilde{\beta}_{S}+\mathbf{V}^{\mathrm{T}}_{i}\tilde{\beta}_{V}-\mathbf{Y}_{i})^{2}.$
(3)

According to the derivation rule of partitioned regression model [5, 11], with
$\mathbf{S}=\bar{\mathbf{S}}$ and $\mathbf{V}=\bar{\mathbf{V}}$, we have:

$\displaystyle\tilde{\beta}_{V}-\beta_{V}$
$\displaystyle=(\frac{1}{n}\sum^{n}_{i=1}\mathbf{V}_{i}^{\mathrm{T}}\mathbf{V}_{i})^{-1}(\frac{1}{n}\sum^{n}_{i=1}\mathbf{V}_{i}^{\mathrm{T}}g(\mathbf{S}_{i}))$
(4)
$\displaystyle+(\frac{1}{n}\sum^{n}_{i=1}\mathbf{V}_{i}^{\mathrm{T}}\mathbf{V}_{i})^{-1}(\frac{1}{n}\sum^{n}_{i=1}\mathbf{V}_{i}^{\mathrm{T}}\mathbf{S}_{i})(\beta_{S}-\tilde{\beta}_{S}),$

$\displaystyle\tilde{\beta}_{S}-\beta_{S}$
$\displaystyle=(\frac{1}{n}\sum^{n}_{i=1}\mathbf{S}_{i}^{\mathrm{T}}\mathbf{S}_{i})^{-1}(\frac{1}{n}\sum^{n}_{i=1}\mathbf{S}_{i}^{\mathrm{T}}g(\mathbf{S}_{i}))$
(5)
$\displaystyle+(\frac{1}{n}\sum^{n}_{i=1}\mathbf{S}_{i}^{\mathrm{T}}\mathbf{S}_{i})^{-1}(\frac{1}{n}\sum^{n}_{i=1}\mathbf{S}_{i}^{\mathrm{T}}\mathbf{V}_{i})(\beta_{V}-\tilde{\beta}_{V}),$

where $n$ is labeled node size, $\mathbf{S}_{i}$ is $i$-th sample of
$\mathbf{S}$,
$\frac{1}{n}\sum^{n}_{i=1}\mathbf{V}_{i}^{\mathrm{T}}g(\mathbf{S}_{i})=\mathbb{E}(\mathbf{V}^{\mathrm{T}}g(\mathbf{S}))+o_{p}(1)$,
$\frac{1}{n}\sum^{n}_{i=1}\mathbf{V}_{i}^{\mathrm{T}}\mathbf{S}_{i}=\mathbb{E}(\mathbf{V}^{\mathrm{T}}\mathbf{S})+o_{p}(1)$
and $o_{p}(1)$ is the error which is negligible. Ideally,
$\tilde{\beta}_{V}-\beta_{V}=0$ indicates that there is no bias between the
estimated and the real parameter. However, if
$\mathbb{E}(\mathbf{V}^{\mathrm{T}}\mathbf{S})\neq 0$ or
$\mathbb{E}(\mathbf{V}^{\mathrm{T}}g(\mathbf{S}))\neq 0$ in Eq. (4),
$\tilde{\beta}_{V}$ will be biased, leading to the biased estimation on
$\tilde{\beta}_{S}$ in Eq. (5) as well, i.e, the true effect of learned
embeddings on label $\mathbf{Y}$ can not be estimated precisely. Since the
correlation between $\mathbf{V}$ and $\mathbf{S}$ (or $g(\mathbf{S})$) 111We
assume all variables are centered with zero mean. This assumption could be
satisfied by adding a normalization layer after the learned embeddings. might
shift in test phase, the biased parameters learned in training set is not the
optimal parameters for predicting testing nodes. Therefore, to increase the
stability of prediction, we need to unbiasedly estimate the parameters of
$\tilde{\beta}_{V}$ by removing the correlation between $\mathbf{V}$ and
$\mathbf{S}$ (or $g(\mathbf{S})$) on training graph, making
$\mathbb{E}(\mathbf{V}^{\mathrm{T}}\mathbf{S})=0$ and
$\mathbb{E}(\mathbf{V}^{\mathrm{T}}g(\mathbf{S}))=0$. Note that
$\frac{1}{n}\sum^{n}_{i=1}\mathbf{S}_{i}^{\mathrm{T}}g(\mathbf{S}_{i})$ in Eq.
(5) can also cause estimation bias, but the relation between $\mathbf{S}$ and
$g(\mathbf{S})$ is stable across environments, which do not affect the
stability to some extent.

Figure 2: (a) Diagram of decorrelating node embeddings with confounding
balance. $\tilde{\mathbf{H}}^{(K-1)}$ is the node embedding matrix to be
decorrelated. $T$ is the treatment variable, corresponding to one target
variable in $\tilde{\mathbf{H}}^{(K-1)}$. $\mathbf{X}$ means the confounders,
corresponding to the remaining variables of the target variable in
$\tilde{\mathbf{H}}^{(K-1)}$. $Y$ is the outcome, corresponding to labels. (b)
The framework of GNN-DVD. The same color in the two figures represents the
same kind of variable.

##  III Proposed Model

###  III-A Revisiting on Variable Decorrelation in Causal View

To decorrelate $\mathbf{V}$ and $\mathbf{S}$ (or $g(\mathbf{S})$), we should
decorrelate the output variables of
$\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})$. [5] proposes a Variable
Decorrelation (VD) term with sample reweighting technique to eliminate the
correlation between each variable pair, in which the sample weights are
learned by jointly minimizing the moment discrepancy between each variable
pair:

$\mathcal{L}_{VD}(\mathbf{H})=\sum^{p}_{j=1}||\mathbf{H}_{.j}^{\mathrm{T}}\Lambda_{\mathbf{w}}\mathbf{H}_{.-j}/n-\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{w}/n\cdot\mathbf{H}_{.-j}^{\mathrm{T}}\mathbf{w}/n||^{2}_{2},$
(6)

where $\mathbf{H}\in\mathbb{R}^{n\times p}$ means the variables needed to be
decorrelated, i.e., $\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})$ of
GNNs, $\mathbf{H}_{.j}$ is $j$-th variable of $\mathbf{H}$,
$\mathbf{H}_{.-j}=\mathbf{H}\backslash{\mathbf{H}_{.j}}$ means all the
remaining variables by setting the value of $j$-th variable in $\mathbf{H}$ as
zero, $\mathbf{w}\in\mathbb{R}^{n\times 1}$ is the sample weights,
$\sum_{i=1}^{n}\mathbf{w}_{i}=n$ and
$\Lambda_{\mathbf{w}}=\text{diag}(\mathbf{w}_{1},\cdots,\mathbf{w}_{n})$ is
the corresponding diagonal matrix. As we can see,
$\mathcal{L}_{VD}(\mathbf{H})$ can be reformulated as $\sum_{j\neq
k}||\mathbf{H}_{.j}^{\mathrm{T}}\Lambda_{\mathbf{w}}\mathbf{H}_{.k}/n-\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{w}/n\cdot\mathbf{H}_{.k}^{\mathrm{T}}\mathbf{w}/n||^{2}_{2}$,
and it aims to let
$\mathbb{E}(\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{H}_{.k})=\mathbb{E}(\mathbf{H}_{.j}^{\mathrm{T}})\mathbb{E}(\mathbf{H}_{.k})$
for each variable pair $j$ and $k$. $\mathcal{L}_{VD}(\mathbf{H})$
decorrelates all the variable pairs equally. However, decorrelating all the
variables requires sufficient samples [5], i.e., $n\rightarrow\infty$, which
is hard to be satisfied, especially in the semi-supervised setting. In this
scenario, we cannot guarantee $\mathcal{L}_{VD}(\mathbf{H})=0$. Therefore, the
key challenge is the difficulties of removing the spurious correlation that
has the largest impact on the unbiased estimation when
$\mathcal{L}_{VD}(\mathbf{H})\neq 0$.

Inspired by confounding balancing technique in observational studies [12], we
revisit the VD regularizer in causal view and show how to differentiate each
variable pair. Confounding balancing techniques are often used for causal
effect estimation of treatment $T$, where the distributions of confounders
$\mathbf{X}$ are different between treated ($T=1$) and control ($T=0$) groups
because of non-random treatment assignment. One could balance the distribution
of confounders between treatment and control groups to unbiasedly estimate
causal treatment effects [13, 14, 15]. Most balancing approaches exploit
moments to characterize distributions, and balance them by adjusting sample
weights $\mathbf{w}$ as follows:
$\mathbf{w}=\arg\min\limits_{\mathbf{w}}||\sum_{i:T_{i}=1}\mathbf{X}_{i}-\sum_{i:T_{i}=0}\mathbf{w}_{i}\cdot\mathbf{X}_{i}||^{2}_{2}.$
After balancing, the treatment $T$ and the confounders $\mathbf{X}$ tend to be
independent.

Given one targeted variable $j$, its decorrelation term,
$\mathcal{L}_{VD_{j}}=||\mathbf{H}_{.j}^{\mathrm{T}}\Lambda_{\mathbf{w}}\mathbf{H}_{.-j}/n-\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{w}/n\cdot\mathbf{H}_{.-j}^{\mathrm{T}}\mathbf{w}/n||^{2}_{2}$,
is to make $\mathbf{H}_{.j}$ independent from $\mathbf{H}_{.-j}$ 222Nonlinear
relationships between variables can be incorporated by considering high-order
moments in Eq. (6), for example, a polynomial augmented function
$f(\mathbf{H})=(\mathbf{H},\mathbf{H}^{2},\mathbf{H}_{.i}\mathbf{H}_{.j},\mathbf{H}^{3},\mathbf{H}_{.i}\mathbf{H}_{.j}\mathbf{H}_{.k},\cdots)$.,
which is same as the confounding balancing term making treatment and
confounders independent. Thereby, $\mathcal{L}_{VD_{j}}$ can also be viewed as
a confounding balancing term, where $\mathbf{H}_{.j}$ is treatment and
$\mathbf{H}_{.-j}$ is confounders, illustrated in Fig. 2(a). Hence, our target
can be explained as unbiased estimation of causal effect of each variable
which is invariant across training and test set. As different variable may
contribute unequally to the confounding bias, it is necessary to differentiate
the confounders. The target of differentiating confounders exactly matches our
target that removing the correlation of variables that has the largest impact
on the unbiased estimation.

###  III-B Differetiated Variable Decorrelation

Considering a continuous treatment, the causal effect of treatment can be
measured by Marginal Treatment Effect Function (MTEF) [16], and defined as:
$MTEF=\frac{\mathbb{E}[Y_{i}(t)]-\mathbb{E}[Y_{i}(t-\Delta t)]}{\Delta t}$,
where $Y_{i}(t)$ represents the potential outcome of sample $i$ with treatment
status $T=t$, and $\Delta t$ denotes the increasing level of treatment. With
the sample weights $\mathbf{w}$ decorrelating the treatment and the
confounders, we can estimate the MTEF by:

$\widehat{MTEF}=\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}\cdot
Y_{i}(T_{i})-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta t}\mathbf{w}_{j}\cdot
Y_{j}(T_{j})}{\Delta t},$ (7)

where $n_{i}$ and $n_{j}$ are the number of samples for two groups,
respectively. Next, we theoretically analyze how to differentiate confounders’
weights with the following theorem.

######  Theorem 1.

In observational studies, different confounders make unequal confounding bias
on Marginal Treatment Effect Function (MTEF) with their own weights, and the
weights can be learned via regressing outcome $Y$ on confounders $\mathbf{X}$
and treatment variable $T$.

###### Proof.

Recalling the Assumption 1, we rewrite the label generation process Eq. (1)
under MTEF setting as:

$Y(t)=\sum_{k}\xi_{k}\mathbf{X}_{.k}+\gamma t+g_{t}(\mathbf{S})+\varepsilon,$
(8)

where $\alpha=[\xi,\gamma]$ are the linear coefficients, and
$g_{T=t}(\mathbf{S})$ is the output of nonlinear transformation of the stable
variables $\mathbf{S}$ when treatment $T$ is $t$. Note that if
$T\notin\mathbf{S}$, changing the value of $T$ will not change the value of
$g_{t}(\mathbf{S})$, $T\in\mathbf{S}$ otherwise.

Under above formulation, we write estimator of $\widehat{MTEF}$ as:

$\displaystyle\widehat{MTEF}=\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}Y_{i}(T_{i})-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}Y_{j}(T_{j})}{\Delta t}$ (9)
$\displaystyle=\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}(\sum_{k}\xi_{k}\mathbf{X}_{ik}+\gamma
t+g_{T=t}(\mathbf{S}_{i})+\epsilon)}{\Delta t}$
$\displaystyle-\frac{\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}(\sum_{k}\xi_{k}\mathbf{X}_{jk}+\gamma(t-\Delta t)+g_{t-\Delta
t}(\mathbf{S}_{i})+\epsilon)}{\Delta t}$
$\displaystyle=\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}\gamma
t-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta t}\mathbf{w}_{j}\gamma(t-\Delta
t)}{\Delta t}$
$\displaystyle+\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}\sum_{k}\xi_{k}\mathbf{X}_{ik}-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}\sum_{k}\xi_{k}\mathbf{X}_{ik}}{\Delta t}$
$\displaystyle+\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}g_{t}(\mathbf{S}_{i})-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}g_{t-\Delta t}(\mathbf{S}_{i})}{\Delta t}+\phi(\epsilon)$
$\displaystyle=MTEF$ $\displaystyle+\sum_{k\neq
t}\xi_{k}(\frac{\sum_{i:T_{i}=t}\frac{1}{n_{i}}\mathbf{w}_{i}\mathbf{X}_{ik}-\sum_{j:T_{j}=t-\Delta
t}\frac{1}{n_{j}}\mathbf{w}_{j}\mathbf{X}_{jk}}{\Delta t})$
$\displaystyle+\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}g_{t}(\mathbf{S}_{i})-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}g_{t-\Delta t}(\mathbf{S}_{i})}{\Delta t}+\phi(\epsilon),$

where $\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}\gamma
t-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta t}\mathbf{w}_{j}\gamma(t-\Delta
t)}{\Delta t}$ is the ground truth of $MTEF$, $\phi(\varepsilon)$ means the
noise term, and $\phi(\varepsilon)\simeq 0$ with Gaussian noise. According to
the last equation, to reduce the bias of $\widehat{MTEF}$, we need regulate
the term
$\sum\nolimits_{k}\xi_{k}(\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}\mathbf{X}_{ik}-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}\mathbf{X}_{jk}}{\Delta t})$ and
$\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}g_{t}(\mathbf{S}_{i})-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}g_{t-\Delta t}(\mathbf{S}_{i})}{\Delta t}$, where the second
term has the unknown term $g_{T}(\mathbf{S}_{i})$ so that we can only try to
reduce the first term.
$\frac{\frac{1}{n_{i}}\sum_{i:T_{i}=t}\mathbf{w}_{i}\mathbf{X}_{ik}-\frac{1}{n_{j}}\sum_{j:T_{j}=t-\Delta
t}\mathbf{w}_{j}\mathbf{X}_{jk}}{\Delta t}$ means the difference of the $k$-th
confounder between treated and control samples. The parameter $\xi_{k}$
represents the confounding bias weight of the $k$-th confounder, and it is the
coefficient of $\mathbf{X}_{.k}$. Moreover, because our target is to learn the
weight of each variable pair, i.e., between treatment and each confounder, we
also need to learn the weight $\gamma$ of treatment $T$. Hence, according to
Eq. (8), the confounder weights and treatment weight can be learned by
regressing observed outcome $Y$ on confounders $\mathbf{X}$ and treatment $T$.
∎

Due to the connection between treatment effect estimation with variable
decorrelation as analyzed in Section III-A, we utilize Theorem 1 to reweight
the variable weight in variable decorrelation term. When apply the Theorem 1
to GNNs, the confounders $\mathbf{X}$ should be $\mathbf{H}_{.-j}$ and
treatment $T$ is $\mathbf{H}_{.j}$, where the embedding $\mathbf{H}$ is
learned by $\hat{\mathscr{G}}(\mathbf{X},\mathbf{A};\theta_{g})$ in Eq. (2).
And the variable weights $\alpha$ is equal to the regression coefficients for
$\mathbf{H}$, i.e, $\hat{\beta}$ in Eq. (2). Then the Differentiated Variable
Decorrelation (DVD) term can be formulated as follows:

$\displaystyle\min_{\mathbf{w}}\mathcal{L}_{DVD}(\mathbf{H})=$
$\displaystyle\sum\nolimits_{j=1}^{p}(\alpha^{\mathrm{T}}\cdot\text{abs}(\mathbf{H}_{.j}^{\mathrm{T}}\Lambda_{\mathbf{w}}\mathbf{H}_{.-j}/n$
(10)
$\displaystyle-\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{w}/n\cdot\mathbf{H}_{.-j}^{\mathrm{T}}\mathbf{w}/n))^{2}$
$\displaystyle+\frac{\lambda_{1}}{n}\sum\nolimits_{i=1}^{n}\mathbf{w}_{i}^{2}+\lambda_{2}(\frac{1}{n}\sum\nolimits_{i=1}^{n}\mathbf{w}_{i}-1)^{2},$
$\displaystyle s.t.\mathbf{w}\succeq 0$

where $\text{abs}(\cdot)$ means the element-wise absolute value operation,
preventing positive and negative values from eliminating. Term
$\frac{\lambda_{1}}{n}\sum_{i=1}^{n}\mathbf{w}_{i}^{2}$ is added to reduce the
variance of sample weights to achieve stability, and the formula
$\lambda_{2}(\frac{1}{n}\sum_{i=1}^{n}\mathbf{w}_{i}-1)^{2}$ avoids all the
sample weights to be 0. The term $\mathbf{w}\succeq 0$ constrains each sample
weight to be non-negative. After variable reweighting, the weighted
decorrelation term in Eq. (10) can be rewritten as $\sum_{j\neq
k}\alpha_{j}^{2}\alpha_{k}^{2}||\mathbf{H}_{.j}^{\mathrm{T}}\Lambda_{\mathbf{w}}\mathbf{H}_{.k}/n-\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{w}/n\cdot\mathbf{H}_{.k}^{\mathrm{T}}\mathbf{w}/n||^{2}_{2}$,
and the weight for variable pair $j$ and $k$ would be
$\alpha_{j}^{2}\alpha_{k}^{2}$, hence it considers both the weights of
treatment and confounder. Then we derive the uniqueness property of
$\mathbf{w}$ as follows:

######  Theorem 2 (Uniqueness).

If $\lambda_{1}n\gg p^{2}+\lambda_{2}$,
$p^{2}\gg\max(\lambda_{1},\lambda_{2})$, $|\mathbf{H}_{i,j}|\leq c$ and
$|\alpha_{i}|\leq c$ for some constant $c$, the solution
$\hat{\mathbf{w}}\in\\{\mathbf{w}:|\mathbf{w}_{i}|\leq c\\}$ to minimize Eq.
(10) is unique.

###### Proof.

See Appendix A.∎

###  III-C Debiased GNN Framework

In this section, we describe the framework of Debiased GNN (DGNN) that
incorporates VD/DVD term with GNNs in a seamless way. As analyzed in Section
II-B, decorrelating
$\hat{\mathbf{A}}\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})$ could
make parameter estimation of GCN unbiased. However, most GNNs follow a layer-
by-layer stacking architecture, and the output embedding of each layer is more
easy to obtain in implementing. Since
$\hat{\mathbf{A}}\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})$ is the
aggregation of the first layer embedding
$\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})$, decorrelating
$\hat{\mathbf{A}}\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})$ may lack
the flexibility that incorporates VD/DVD term with other GNN architectures.
Fortunately, we have the following theorem to identify a more flexible way to
combine variable decorrelation with GNNs.

######  Theorem 3.

Given $p$ pairwise uncorrelated variables
$\mathbf{Z}=(\mathbf{Z}_{1},\mathbf{Z}_{2},\cdots,\mathbf{Z}_{p})$, with a
linear aggregation operator $\hat{\mathbf{A}}$, the variables of
$\mathbf{Y}=\hat{\mathbf{A}}\mathbf{Z}$ are still pairwise uncorrelated.

###### Proof.

Let $\mathbf{Z}=\\{\mathbf{Z}_{1},\mathbf{Z}_{2},\cdots,\mathbf{Z}_{p}\\}$ be
$p$ pairwise uncorrelated variables.
$\forall\mathbf{Z}_{i},\mathbf{Z}_{j}\in\mathbf{Z}$,
$(\mathbf{Z}_{i}^{(1)},\mathbf{Z}_{i}^{(2)},\cdots,\mathbf{Z}_{i}^{(n)})$ and
$(\mathbf{Z}_{j}^{(1)},\mathbf{Z}_{j}^{(2)},\cdots,\mathbf{Z}_{j}^{(n)})$ are
$n$ simple random samples drawn from $\mathbf{Z}_{i}$ and $\mathbf{Z}_{j}$
respectively, and have same distribution with $\mathbf{Z}_{i}$ and
$\mathbf{Z}_{j}$. Given a linear aggregation matrix
$\hat{\mathbf{A}}=(a_{ij})$, $\forall s,v\in(1,2,\cdots,n)$, let
$\mathbf{Y}_{i}^{(s)}=\sum_{k=1}^{n}a_{sk}\mathbf{Z}_{i}^{(k)}$ and
$\mathbf{Y}_{j}^{(v)}=\sum_{l=1}^{n}a_{vl}\mathbf{Z}_{j}^{(l)}$, and we have
following derivation:

$\displaystyle\text{Cov}(\mathbf{Y}_{i}^{(s)},\mathbf{Y}_{j}^{(v)})$
$\displaystyle=\text{Cov}(\sum_{k=1}^{n}a_{sk}\mathbf{Z}_{i}^{(k)},\sum_{l=1}^{n}a_{vl}\mathbf{Z}_{j}^{(l)})$
$\displaystyle=\sum_{k=1}^{n}\sum_{l=1}^{n}a_{sk}a_{vl}\text{Cov}(\mathbf{Z}_{i}^{(k)},\mathbf{Z}_{j}^{(l)})$
$\displaystyle=\sum_{k=1}^{n}\sum_{l=1}^{n}a_{sk}a_{vl}\delta_{ij},$

where $\delta_{ij}=0$ when $i\neq j$, otherwise $\delta_{ij}=1$. Therefore,
when $i\neq j$, we have
$\text{Cov}(\mathbf{Y}_{i}^{(s)},\mathbf{Y}_{j}^{(v)})=0$ and
$\text{Cov}(\mathbf{Y}_{i},\mathbf{Y}_{j})=0$. Extended the conclusion to
multiple variable,
$\mathbf{Y}=(\mathbf{Y}_{1},\mathbf{Y}_{2},\cdots,\mathbf{Y}_{n})$ are
pairwise uncorrelated. ∎

The theorem indicates that if the variables of embeddings $\mathbf{Z}$ are
uncorrelated, after any form of linear neighborhood aggregation
$\hat{\mathbf{A}}$, e.g., average, attention or sum, the variables of
transformed embeddings $\mathbf{Y}$ would be also uncorrelated. Therefore,
decorrelating $\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})$ can also
reduce the estimation bias. For a $K$ layers of GNN, we can directly
decorrelate the output of $(K-1)$-th layer, i.e.,
$\sigma(\hat{\mathbf{A}}\cdots\sigma(\hat{\mathbf{A}}\mathbf{X}\mathbf{W}^{(0)})\cdots\mathbf{W}^{(K-2)})$
for a $K$ layers of GCN.

The previous analysis finds a flexible way to incorporate VD/DVD term with
GNNs, however, recall that we analyze GNNs based on the least-squares loss in
Eq. (3), and most existing GNNs are designed for the classification task.
Therefore, in the following, we analyze that the previous conclusions are
still applicable in classification. We consider the cases that the softmax
layer is used as the output layer of GNNs and loss is the cross-entropy error
function. We use the Newton-Raphson update rule [17] to bridge the gap between
linear regression and multi-classification. According to the Newton-Raphson
update rule, the update formula for transformation matrix $\mathbf{W}^{(K-1)}$
of the last layer of GCN can be derived:

$\displaystyle\mathbf{W}_{.j}^{(\text{new})}$
$\displaystyle=\mathbf{W}_{.j}^{(\text{old})}-(\mathbf{H}^{\mathrm{T}}\mathbf{R}\mathbf{H})^{-1}\mathbf{H}^{\mathrm{T}}(\mathbf{H}\mathbf{W}_{.j}^{(\text{old})}-\mathbf{Y}_{.j})$
(11)
$\displaystyle=(\mathbf{H}^{\mathrm{T}}\mathbf{R}\mathbf{H})^{-1}\\{\mathbf{H}^{\mathrm{T}}\mathbf{R}\mathbf{H}\mathbf{W}_{.j}^{(\text{old})}-\mathbf{H}^{\mathrm{T}}(\mathbf{H}\mathbf{W}_{.j}^{(\text{old})}-\mathbf{Y}_{.j})\\}$
$\displaystyle=(\mathbf{H}^{\mathrm{T}}\mathbf{R}\mathbf{H})^{-1}\mathbf{H}^{\mathrm{T}}\mathbf{R}\mathbf{z},$

where
$\mathbf{R}_{kj}=-\sum_{n=1}^{N}\mathbf{H}_{n}\mathbf{W}_{.k}^{(\text{old})}(\mathbf{I}_{kj}-\mathbf{H}_{n}\mathbf{W}_{.j}^{(\text{old})})$
is a weighing matrix and $\mathbf{I}_{kj}$ is the element of the identity
matrix, and
$\mathbf{z}=\mathbf{H}\mathbf{W}_{.j}^{(\text{old})}-\mathbf{R}^{-1}(\mathbf{Y}_{.j}-\mathbf{W}_{.j}\mathbf{H})$
is an effective target value. Eq. (11) takes the form of a set of normal
equations for a weighted least-squares problem. As the weighing matrix
$\mathbf{R}$ is not constant but depends on the parameter vector
$\mathbf{W}^{(\text{old})}_{.j}$, we must apply the normal equations
iteratively. Each iteration uses the last iteration weight vector
$\mathbf{W}^{(\text{old})}_{.j}$ to calculate a revised weighing matrix
$\mathbf{R}$ and regresses the target value $\mathbf{z}$ with
$\mathbf{H}\mathbf{W}_{.j}^{(\text{new})}$. Therefore, the variable
decorrelation can also be applied to the GNNs with softmax classifier to
reduce the estimation bias in each iteration. Note that according to update
formula Eq. (11), we should calculate the inverse matrix
$(\mathbf{H}^{\mathrm{T}}\mathbf{R}\mathbf{H})^{-1}$ in each iteration, which
requires high computation. In practice, we use gradient descent methods to
approximate Newton-Raphson update rule and it works well in experiments.

Fig. 2(b) is the framework of GNN-DVD, where we input the labeled nodes’
embeddings $\tilde{\mathbf{H}}^{(K-1)}$ into the regularizer
$\mathcal{L}_{DVD}(\tilde{\mathbf{H}}^{(K-1)})$. As GCN has the formula
$\text{softmax}(\hat{\mathbf{A}}\mathbf{H}^{(K-1)}\mathbf{W}^{(K-1)})$, the
variable weights of $\tilde{\mathbf{H}}^{(K-1)}$ used for differentiating
$\mathcal{L}_{DVD}(\tilde{\mathbf{H}}^{(K-1)})$ can be computed from

$\alpha=\text{Var}(\mathbf{W}^{(K-1)},\text{axis}=1)$ (12)

where $\text{Var}(\cdot,\text{axis}=1)$ refers to calculating the variance of
each row of some matrix and it reflects each variable’s weight for
classification which is similar to the regression coefficients. Note that when
incorporating VD term with GNNs, we do not need compute the variable weights.
Then the sample weights $\mathbf{w}$ learned by DVD term have the ability to
remove the correlation in $\tilde{\mathbf{H}}^{(K-1)}$. Inspired by sample
weighting methods [18], we propose to use this sample weights to reweight
softmax loss:

$\min\limits_{\theta}\mathcal{L}_{G}=\sum\limits_{l\in\mathcal{Y}_{L}}\mathbf{w}_{l}\cdot\ln(q(\tilde{\mathbf{H}}^{(K)}_{l})\cdot\mathbf{Y}_{l}),\\\
$ (13)

where $q(\cdot)$ is the softmax function, $\mathcal{Y}_{L}$ is the set of
labeled node indices and $\theta$ is the set of parameters of GCN. After
reweighting nodes by these weights, we can create a pseudo-population where
the biases in node neighborhood are effectively reduced, with which the off-
the-shelf GNN models can achieve more accurate prediction under agnostic
environments. The whole algorithm is summarized in Algorithm 1.

Input : Training graph
$\mathcal{G}_{train}=\\{\mathbf{A},\mathbf{X},\mathbf{Y}\\}$, and indices of
labeled nodes $\mathcal{Y}_{L}$; Max iteration:$maxIter$

Output : GNN parameter $\theta$ and sample weights $\mathbf{w}$

Initialization : Let $\mathbf{w}=\omega\odot\omega$ and initialize sample
weights $\omega$ with $\mathbf{1}$; Initialize GNN’s parameters $\theta$ with
random uniform distribution; Iteration $t\leftarrow 0$

1 while  _not converged or $t<maxIter$_ do

2

3 Optimize $\theta^{(t)}$ to minimize $\mathcal{L}_{G}$ via Eq. (13);

4 Calculate variable weights $\alpha^{(t)}$ from $\mathbf{W}^{(K-1)}$ via Eq.
(12);

5 Optimize $\omega^{(t)}$ to minimize
$\mathcal{L}_{DVD}(\tilde{\mathbf{H}}^{(K-1)})$ via Eq. (10);

6 $t=t+1$;

7 end while

Return: $\theta$ and $\mathbf{w}=\omega\odot\omega$

Algorithm 1 GNN-DVD Algorithm

To optimize our GNN-DVD algorithm, we propose an iterative method. Firstly, we
let $\mathbf{w}=\omega\odot\omega$ to ensure non-negativity of $\mathbf{w}$
and initialize sample weight $\omega_{i}=1$ for each sample $i$ and GNN’s
parameters $\theta$ with random uniform distribution. Once the initial values
are given, in each iteration, we fix the sample weights $\omega$ and update
the GNN’s parameters $\theta$ by $\mathcal{L}_{G}$ with gradient descent, then
compute the confounder weights $\alpha$ from the linear transform matrix
$\mathbf{W}^{(K-1)}$. With $\alpha$ and fixing the GNN’s parameters $\theta$,
we update the sample weights $\omega$ with gradient descent to minimize
$\mathcal{L}_{DVD}(\mathbf{H}^{(K-1)})$. We iteratively update the sample
weights $\mathbf{w}$ and GNN’s parameters $\theta$ until $\mathcal{L}_{G}$
converges.

###  III-D Extension to GAT

We can easily incorporate the VD/DVD term with other GNNs. We combine them
with GAT and more extensions leave as future work. GAT utilizes an attention
mechanism to aggregate neighbor information. It also follows the linear
aggregation and transformation steps. Similar to GCN, the hidden embedding
$\tilde{\mathbf{H}}^{(K-1)}$ is the input of VD/DVD term, and the variable
weights $\alpha$ are calculated from the transformation matrix
$\mathbf{W}^{(K-1)}$ and the sample weights $\mathbf{w}$ are used to reweight
the softmax loss. Note that the original paper utilizes the same
transformation matrix $\mathbf{W}^{(K-1)}$ for transforming embedding and
learning attention values. Because $\alpha$ means the importance of each
variable for classification, and it should be calculated from transformation
matrix $\mathbf{W}^{(K-1)}$ for transforming embedding, hence we use separate
matrices for transforming embedding and learning attention values,
respectively. This modification does not change the performance of GAT in
experiments.

###  III-E Complexity Analysis

Compared with base models (e.g., GCN and GAT), the main incremental time cost
is the complexity from VD/DVD term. For a training graph with $n$ labeled
nodes, we analyze the time complexity of the VD/DVD term in each iteration.
For calculating the VD loss, its complexity is $\mathcal{O}(np^{2})$, where
$p$ is the dimension of embeddings. And for DVD loss, its complexity is the
same as VD, as the complexity of calculating variable weights $\alpha$ is
$\mathcal{O}(np)$, which is relatively small comparing with
$\mathcal{O}(np^{2})$. For updating $\mathbf{w}$, the complexity is dominated
by the step of calculating the partial gradients of the function
$\mathcal{L}_{DVD}(\mathbf{H})$ with respect to variable $\mathbf{w}$. The
complexity of
$\frac{\partial\mathcal{L}_{DVD}(\mathbf{H})}{\partial\mathbf{w}}$ is
$\mathcal{O}(np^{2})$. In total, the complexity of each iteration for VD/DVD
term in Algorithm 1 is $\mathcal{O}(np^{2})$. And it is quite smaller than the
base models (e.g., the complexity of GCN is
$\mathcal{O}({\mathcal{E}cp^{2}})$, where $\mathcal{E}$ is the number of edges
and $c$ is the dimension of input node features).

###  III-F Discussion

In our paper, we propose to integrate two decorrelation terms (i.e., VD/DVD
term) with GNN models to eliminate the estimation bias. Here we discuss the
advantages and disadvantages of these two terms. VD term aggressively
decorrelates all the variables learned by GNNs, however, it theoretically
requires a large number of samples to achieve this goal. To overcome this
dilemma, the DVD term is proposed to differentiate the variable weights in the
VD term, aiming to remove the most unexpected correlation. However, due to the
existence of unknown term $g(\mathbf{S})$ in Eq. (8), introducing more
parameters to optimize may increase the instability of the model. Hence, when
the number of labeled samples is large, performing GNN-VD may induce more
stable results.

##  IV Experiments

###  IV-A Datasets

Here, we validate the effectiveness of our methods on node classification with
two kinds of selection bias, i.e., label selection bias and small sample
selection bias. For label selection bias, we employ three widely used graph
datasets: Cora, Citeseer and Pubmed [7]. As in Section II-A, we get three
biased degrees as well as the original unbiased labeled nodes for each
dataset. For small sample selection bias, we conduct the experiments on NELL
dataset [19], where each class only has at most 1/5/10 labeled nodes for
training. Due to the large scale of this dataset, the test nodes are easily to
have distribution shifts from training nodes. The details of the datasets are
summarized in Table II.

TABLE II: Dataset statistics Dataset | Type | Nodes | Edges | Classes | Features | Bias degree ($\epsilon$) | Bias type
---|---|---|---|---|---|---|---
Cora | Citation network | 2,708 | 5,429 | 7 | 1,433 | 0.7/0.8/0.9 | Label selection bias
Citeseer | Citation network | 3,327 | 4,732 | 6 | 3,703 | 0.7/0.8/0.9 | Label selection bias
Pubmed | Citation network | 19,717 | 44,338 | 3 | 500 | 0.7/0.8/0.9 | Label selection bias
NELL | Knowledge graph | 65,755 | 266,144 | 210 | 5,414 | 1/5/10 labeled nodes per class | Small sample selection bias

###  IV-B Baselines

We compare our proposed framework with several related baselines:

  * •

Base models: GCN [2] and GAT [3] are classical GNN methods. We utilize them as
the base models in our framework, so they are the most related baselines to
validate the effectiveness of the proposed framework.

  * •

GNM-GCN/GAT [20]: A GNN method which considers unbalanced label selection bias
problem in transductive setting. They also utilize GCN/GAT as their base
models.

  * •

Chebyshev [2]: It is a GCN-based method utilizing third-order Chebyshev
filters.

  * •

SGC [8]: It is a simplified GCN-based method, which reduces the excess
complexity through successively removing nonlinearities and collapsing weight
matrices between consecutive layers.

  * •

APPNP [21]: It is one of the state-of-the-art GNN methods that combines
PageRank with GCN.

  * •

Planetoid [10]: It is a classical semi-supervised graph embedding method. We
use its inductive variant.

  * •

MLP: It is a two-layer multilayer perceptron trained on the labeled nodes with
only node features as input.

  * •

DGNN: It is the debiased GNN framework proposed in this paper. We incorporate
the VD/DVD term with GCN/GAT under our proposed framework called GCN/GAT-
VD/DVD.

###  IV-C Experimental Setup

As the Section II-A has described, for all datasets, to simulate the agnostic
selection bias scenario, we first follow the inductive setting in [8] that
masks the validation and test nodes in the training phase and validation and
test with the whole graph so that the test nodes will be agnostic. For GCN and
GAT, we utilize the same two-layer architecture as their original paper [2,
3]. We use the following sets of hyperparameters for GCN on Cora, Citeseer,
Pubmed: 0.5 (dropout rate), $5\cdot 10^{-4}$ (L2 regularization) and 32
(number of hidden units); and for NELL: 0.1 (dropout rate), $1\cdot 10^{-5}$
(L2 regularization) and 64 (number of hidden units). For GAT on Cora,
Citeseer, we use: 8 (first layer attention heads), 8 (features each head), 1
(second layer attention head), 0.6 (dropout), $5\cdot 10^{-4}$ (L2
regularization); and for Pubmed: 8 (second layer attention head), $1\cdot
10^{-3}$ (L2 regularization), other parameters are the same as Cora and
Citeseer. To fair comparison, the GNN part of our model uses the same
architecture and hyper-parameters with the base model and we grid search
$\lambda_{1}$ and $\lambda_{2}$ from $\\{0.01,0.1,1,10,100\\}$. For other
baselines, we use the optimal hyper-parameters in the original literatures on
each dataset. For all the experiments, we run each model 10 times with
different random seeds and report its average Accuracy results.

TABLE III: Performance of three citation networks. The ‘*’ indicates the best
results of the baselines. Best results of all methods are indicated in bold.
‘% gain over GCN/GAT’ means the improvement percent of GCN/GAT-DVD against
GCN/GAT.

Method | Cora | Citeseer | Pubmed
---|---|---|---
Unbiased | Light | Medium | Heavy | Unbiased | Light | Medium | Heavy | Unbiased | Light | Medium | Heavy
MLP | 0.5296 | 0.5624 | 0.5197 | 0.5087 | 0.5438 | 0.4532 | 0.3757 | 0.3893 | 0.6914 | 0.6852 | 0.6620 | 0.6378
Planetoid [10] | 0.6650 | 0.5890 | 0.5240 | 0.5180 | 0.6720 | 0.5160 | 0.5140 | 0.4880 | 0.744 | 0.7160 | 0.6770 | 0.6680
Chebyshev [22] | 0.7407 | 0.7116 | 0.7006 | 0.6809 | $0.7232^{*}$ | 0.6542 | 0.6276 | 0.5920 | 0.7450 | 0.7358 | 0.6862 | 0.6732
SGC [8] | 0.779 | 0.7800 | 0.7800 | 0.7530 | 0.724 | 0.6780 | $0.6730^{*}$ | 0.6200 | 0.781 | $\mathbf{0.7880}^{*}$ | 0.7560 | 0.6800
APPNP [21] | $0.8132^{*}$ | 0.7913 | 0.7689 | 0.7629 | 0.6862 | 0.6478 | 0.6052 | 0.5903 | 0.7731 | 0.7639 | 0.7369 | 0.6862
GNM-GCN [20] | 0.7594 | 0.7423 | 0.7531 | 0.7196 | 0.6054 | 0.5793 | 0.5717 | 0.5125 | 0.7654 | 0.7552 | 0.7381 | 0.7072
GNM-GAT [20] | 0.7976 | 0.7875 | 0.7638 | 0.7404 | 0.6832 | 0.6524 | 0.6487 | 0.5865 | 0.7666 | 0.7438 | 0.7568 | 0.6891
GCN [2] | 0.7909 | 0.7851 | 0.7775 | 0.7422 | 0.7075 | 0.6786 | 0.5952 | 0.5551 | $0.7845^{*}$ | 0.7673 | 0.7545 | $0.7247^{*}$
GCN-VD | 0.7980 | 0.7951 | 0.7855 | 0.7522 | 0.7122 | 0.6844 | 0.6676 | 0.6408 | 0.7888 | 0.7727 | 0.7729 | 0.7399
GCN-DVD | 0.7951 | 0.7959 | 0.7885 | 0.7555 | 0.7128 | 0.6908 | 0.6769 | 0.6496 | 0.7874 | 0.7741 | 0.7746 | 0.7542
% gain over GCN | 0.53% | 1.38% | 1.41% | 1.79% | 0.75% | 1.8% | 14.2% | 17.0% | 0.37% | 0.89% | 2.67% | 4.07%
GAT [3] | 0.8100 | $0.8067^{*}$ | $0.8019^{*}$ | 0.7578 | 0.7224 | $0.7033^{*}$ | 0.6683 | $0.6475^{*}$ | 0.7714 | 0.7665 | $0.7579^{*}$ | 0.7068
GAT-VD | 0.8133 | 0.8146 | 0.8079 | 0.7708 | 0.7288 | 0.7149 | 0.6833 | 0.6611 | 0.7732 | 0.7783 | 0.7689 | 0.7149
GAT-DVD | 0.8139 | 0.8179 | 0.8119 | 0.7694 | 0.7294 | 0.7172 | 0.6825 | 0.6627 | 0.7735 | 0.7788 | 0.7723 | 0.7210
% gain over GAT | 0.48% | 1.39% | 1.26% | 1.53% | 0.97% | 1.97% | 2.12% | 2.34% | 0.27% | 1.6% | 1.9% | 2.0%

TABLE IV: Performance of NELL. The ‘*’ indicates the best results of the baselines. Best results of all methods are indicated in bold. ‘Improvement’ means the improvement percent of GCN-VD/DVD (selected better results) against GCN. Dataset | MLP | Planetoid | SGC | GCN | GCN-VD | GCN-DVD | Improvement
---|---|---|---|---|---|---|---
NELL-1 | 0.2385 | 0.3901 | 0.4128 | $0.4416^{*}$ | 0.4652 | 0.4734 | 7.2%
NELL-5 | 0.4938 | 0.3519 | 0.6295 | $0.7030^{*}$ | 0.7424 | 0.7361 | 5.6%
NELL-10 | 0.5838 | 0.5149 | 0.6275 | $0.7615^{*}$ | 0.7734 | 0.7727 | 1.6%

###  IV-D Results on Label Selection Bias Datasets

The results are given in Table III, and we have the following observations.
First, the proposed models (i.e., GCN/GAT with VD/DVD terms) always achieve
the best performances in most cases, which well demonstrates that the
effectiveness of our proposed debiased GNN framework. Second, comparing with
base models, our proposed models all achieve up to 17.0% performance
improvements, and gain larger improvements under heavier bias scenarios. Since
the major difference between our model with base models is the VD/DVD
regularizer, we can safely attribute the significant improvements to the
effective decorrelation term and its seamless joint with GNN models. Third,
GCN/GAT-DVD achieves better results than GCN/GAT-VD in most cases. It
validates the importance and effectiveness of differentiating variables’
weights in the semi-supervised setting. Moreover, our model still outperforms
baselines in the unbiased setting. In real applications, it is hard to control
the collection process without any distribution shift from the training set to
the test set [23]. Therefore, the problem we study is ubiquitous in reality
and our method is effective in most scenarios.

###  IV-E Results on Small Sample Selection Bias Datasets

As NELL is a large-scale graph, we cannot run GAT on a single GPU with 16GB
memory. We only perform GCN-VD/DVD and compare with representative methods
which can perform on this dataset. The results are shown in Table IV. First,
GCN-VD/DVD achieves significant improvements over GCN. It indicates that
selection bias could be induced by a small number of labeled nodes and our
proposed method relieve the estimation bias. Moreover, with fewer labeled
nodes, i.e., larger selection bias, our methods achieve larger improvements
over base models. It further validates our method is an effective method
against heavy bias. Moreover, GCN-DVD further improves GCN-VD with a large
margin on NELL-1 dataset. It means that decorrelating all the variable pairs
equally is suboptimal, and our differentiated strategy is effective when
labeled nodes are scarce. With the number of labeled nodes increases, it may
not necessary to differentiate the variable weights, but GCN-DVD achieves
competitive results with GCN-VD and still outperforms the base model with a
clear margin.

(a) Cora

(b) Citeseer

(c) Pubmed

Figure 3: Embedding correlation analysis on unweighted and weighted GCN.

###  IV-F Sample Weight Analysis

Here we analyze the effect of sample weights $\mathbf{w}$. We compute the
amount of correlation in the labeled nodes’ embeddings
$\tilde{\mathbf{H}}^{(K-1)}$ learned by standard GCN and the weighted
embeddings of the same layer learned by GCN-DVD. Note that, the weights are
the last iteration of sample weights of GCN-DVD. Following [24, 25], the
amount of correlation of GCN and GCN-DVD is measured by Frobenius norm of
cross-corvairance matrix $||C||_{F}^{2}$ computed from variables of
$\tilde{\mathbf{H}}^{(K-1)}$ and weighted $\tilde{\mathbf{H}}^{(K-1)}$
respectively, where $C_{ij}$ represents the covariance between pairwise
variable $i$ and $j$ and the main diagonal of $C$ is set as zero vector.
Figure 3 shows the amount of correlation in unweighted and weighted
embeddings, and we observe that the embeddings’ correlations in all datasets
are reduced, indicating that the weights learned by GCN-DVD can reduce the
correlation between embedded variables. Moreover, as it is hard to reduce the
correlation to zero, the necessity of differentiating variables’ weights is
further validated.

(a) Light

(b) Medium

(c) Heavy

Figure 4: Accuracy of GCN-DVD with different $\lambda_{1}$ and $\lambda_{2}$
on different biased Cora datasets.

(a) Light

(b) Medium

(c) Heavy

Figure 5: Accuracy of GCN-DVD with different $\lambda_{1}$ and $\lambda_{2}$
on different biased Citeseer datasets.

(a) Light

(b) Medium

(c) Heavy

Figure 6: Accuracy of GCN-DVD with different $\lambda_{1}$ and $\lambda_{2}$
on different biased Pubmed datasets.

###  IV-G Parameter sensitivity

We study the sensitiveness of parameters and report the results of GCN-DVD on
three citation networks in Fig. 4-6. The experimental results show that GCN-
DVD is relatively stable to $\lambda_{1}$ and $\lambda_{2}$ with wide ranges
in most cases, indicating the robustness of our model.

###  IV-H Training time per epoch

We report the results for the mean training time of GCN and GCN-DVD per epoch
(forward pass, cross-entropy calculation, backward pass) for 200 epochs on
Cora, Citeseer and Pubmed datasets, measured in seconds wall-clock time, in
Table V. These methods are performed on a RTX 3090 GPU Card. As we can see,
the training time of GCN-DVD term has the same order of magnitude with GCN.
More importantly, the training time of the DVD term will not be influenced by
the base model we choose, i.e., when the base model is GAT, the running time
of DVD term will not change.

TABLE V: The training time per epoch. | Cora | Citeseer | Pubmed
---|---|---|---
GCN | $1.29\times 10^{-2}$ | $2.00\times 10^{-2}$ | $1.11\times 10^{-1}$
GCN-DVD | $6.19\times 10^{-2}$ | $8.46\times 10^{-2}$ | $2.42\times 10^{-1}$

##  V Related Works

In the past few years, Graph Neural Networks (GNNs) [1, 2, 3, 21, 26, 27, 28,
29, 30] have become the major technology to capture patterns encoded in the
graph due to its powerful representation capacity. Recently, KPGNN [31]
applies GNNs to the social event detection task by preserving the incremental
knowledge emerging in social data. MRFasGCN [32] integrates GCN with a Markov
Random Fields (MRF) model to deal with the semi-supervised community detection
problem. Not only pursuing the performance of GNNs on clean data, [33]
proposes an exploratory adversarial attack method, called EpoAtk, to test
whether existing GNNs are robust with adversarial perturbations on graphs.
Although the current GNNs have achieved great success, when applied to the
inductive setting, they all assume that training nodes and test nodes follow
the same distribution. However, this assumption does not always hold in real
applications. GNM [20] first pays attention to the label selection problem on
graph learning, and it learns an IPW estimator to estimate the probability of
each node to be selected and uses this probability to reweight the labeled
nodes. However, it heavily relies on the accuracy of the IPW estimator, which
depends on the label assignment distribution of the whole graph, hence it is
more suitable for the transductive setting.

To enhance the stability in unseen varied distributions, some literatures [5,
34] have revealed the connection between correlation and prediction stability
under model misspecification. Moreover, a kind of literatures [35, 36, 37]
have studied the problem of removing the features correlation effect in neural
networks, which brings great benefits for deep neural networks. However, these
methods are built on simple regressions or regular neural networks such as
CNNs, but GNNs have more complex architectures and properties needed to be
considered. We also notice that [6] propose a differentiated variable
decorrelation term for linear regression. However, this decorrelation term
requires multiple environments with different correlations between stable
variables and unstable variables available in the training stage while our
method does not require this prior knowledge.

##  VI Conclusion

In this paper, we investigate a general and practical problem: learning GNNs
with agnostic label selection bias. The selection bias will inevitably cause
the GNNs to learn the biased correlation between aggregation mode and class
label and make the prediction unstable. We propose a novel debiased GNN
framework, which combines the decorrelation technique with GNNs in a unified
framework. Extensive experiments well demonstrate the effectiveness and
flexibility of DGNN.

## Acknowledgment

This work is partially supported by the National Natural Science Foundation of
China (No. U20B2045, 62192784, 62172052, 61772082, 62002029) and the
Fundamental Research Funds for the Central Universities (No. 2021RC28). Kun
Kuang is supported in part by National Natural Science Foundation of China
(No. 62006207), Young Elite Scientists Sponsorship Program by CAST and
Zhejiang Province Natural Science Foundation (No. LQ21F020020). Shaohua Fan is
supported by BUPT Excellent Ph.D. Students Foundation (No. CX2021311) and
China Scholarship Council. We thank Dr. Tianchi Yang for discussion on the
proof of Theorem 3.

##  Appendix A Proof of Theorem 2

$\displaystyle\hat{\mathbf{w}}$
$\displaystyle=\arg\min_{\mathbf{w}}\sum^{p}_{j=1}(\alpha^{\mathrm{T}}\text{abs}(\mathbf{H}_{.j}^{\mathrm{T}}\Lambda_{\mathbf{w}}\mathbf{H}_{.-j}/n$
(14)
$\displaystyle-\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{w}/n\cdot\mathbf{H}_{.-j}^{\mathrm{T}}\mathbf{w}/n))^{2}+\frac{\lambda_{1}}{n}\sum_{i=1}^{n}\mathbf{w}_{i}^{2}+\lambda_{2}(\frac{1}{n}\sum_{i=1}^{n}\mathbf{w}_{i}-1)^{2}$

###### Proof.

For simplicity, we denote
$\mathcal{L}_{1}=\sum^{p}_{j=1}(\alpha^{\mathrm{T}}\cdot\text{abs}(\mathbf{H}_{.j}^{\mathrm{T}}\Lambda_{\mathbf{w}}\mathbf{H}_{.-j}/n-\mathbf{H}_{.j}^{\mathrm{T}}\mathbf{w}/n\cdot\mathbf{H}_{.-j}^{\mathrm{T}}\mathbf{w}/n))^{2}$,
$\mathcal{L}_{2}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{w}_{i}^{2}$,
$\mathcal{L}_{3}=(\frac{1}{n}\sum_{i=1}^{n}\mathbf{w}_{i}-1)^{2}$ and
$\mathcal{F}(\mathbf{w})=\mathcal{L}_{1}+\lambda_{1}\mathcal{L}_{1}+\lambda_{2}\mathcal{L}_{2}$.
We first calculate the Hessian matrix of $\mathcal{F}(\mathbf{w})$, denoted as
$\mathbf{H}_{e}$, to prove the uniqueness of the optimal solution
$\hat{\mathbf{w}}$, as follows:

$\mathbf{H}_{e}=\frac{\partial^{2}\mathcal{L}_{1}}{\partial\mathbf{w}^{2}}+\lambda_{1}\frac{\partial^{2}\mathcal{L}_{2}}{\partial\mathbf{w}^{2}}+\lambda_{2}\frac{\partial^{2}\mathcal{L}_{3}}{\partial\mathbf{w}^{2}}$

For the term $\mathcal{L}_{1}$, we can rewrite it as:

$\displaystyle\mathcal{L}_{1}$ $\displaystyle=\sum_{j\neq
k}\alpha_{i}^{2}\alpha_{k}^{2}(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{H}_{i,k}\mathbf{w}_{i}$
$\displaystyle-(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,k}\mathbf{w}_{i}))^{2}$
$\displaystyle=\sum_{j\neq
k}\alpha_{i}^{2}\alpha_{k}^{2}((\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{H}_{i,k}\mathbf{w}_{i})^{2}$
$\displaystyle-(\frac{2}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{H}_{i,k}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,k}\mathbf{w}_{i})$
$\displaystyle+((\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,k}\mathbf{w}_{i}))^{2})$

And when $|\mathbf{H}_{i,j}|\leq c$, for any variable $j$ and $k$, and
$|\mathbf{w}_{i}|\leq c$, we have
$\frac{\partial^{2}}{\partial\mathbf{w}^{2}}(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{H}_{i,k}\mathbf{w}_{i})^{2}=\mathcal{O}(\frac{1}{n^{2}})$,
$\frac{\partial^{2}}{\partial\mathbf{w}^{2}}(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,k}\mathbf{w}_{i})=\mathcal{O}(\frac{1}{n^{2}})$
and
$\frac{\partial^{2}}{\partial\mathbf{w}^{2}}((\frac{2}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{H}_{i,k}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,k}\mathbf{w}_{i}))=\mathcal{O}(\frac{1}{n^{2}})$.
Then with $|\alpha_{i}|\leq c$, we have
$\alpha_{i}^{2}\alpha_{k}^{2}\frac{\partial^{2}}{\partial\mathbf{w}^{2}}(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{H}_{i,k}\mathbf{w}_{i}-(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,k}\mathbf{w}_{i}))^{2}=\mathcal{O}(\frac{1}{n^{2}})$.
$\mathcal{L}_{1}$ is sum of $p(p-1)$ such terms. Then we have

$\frac{\partial^{2}\mathcal{L}_{1}}{\partial\mathbf{w}^{2}}=\mathcal{O}(\frac{p^{2}}{n^{2}}).$

With some algebras, we can also have

$\frac{\partial^{2}\mathcal{L}_{2}}{\partial\mathbf{w}^{2}}=\frac{1}{n}\mathbf{I},\frac{\partial^{2}\mathcal{L}_{3}}{\partial\mathbf{w}^{2}}=\frac{1}{n^{2}}\mathbf{1}\mathbf{1}^{\mathrm{T}},$

thus,

$\mathbf{H}_{e}=\mathcal{O}(\frac{p^{2}}{n^{2}})+\frac{\lambda_{1}}{n}\mathbf{I}+\frac{\lambda_{2}}{n^{2}}\mathbf{1}\mathbf{1}^{\mathrm{T}}=\frac{\lambda_{1}}{n}\mathbf{I}+\mathcal{O}(\frac{p^{2}+\lambda_{2}}{n^{2}}).$

Therefore, if $\frac{\lambda_{1}}{n}\gg\frac{p^{2}+\lambda_{2}}{n^{2}}$,
equivalent to $\lambda_{1}n\gg p^{2}+\lambda_{2}$, $\mathbf{H}_{e}$ is an
almost diagonal matrix. Hence, $\mathbf{H}_{e}$ is positive definite [38].
Then the function $\mathcal{F}(\mathbf{w})$ is convex on
$\mathcal{C}=\\{\mathbf{w}:|\mathbf{w}_{i}|\leq c\\}$, and has unique optimal
solution $\hat{\mathbf{w}}$.

Moreover, because $\mathcal{L}_{1}$ is our major decorrelation term, we hope
$\mathcal{L}_{1}$ to dominate the terms $\lambda_{1}\mathcal{L}_{2}$ and
$\lambda_{2}\mathcal{L}_{3}$. On $\mathcal{C}$, we have
$\mathcal{L}_{1}=\mathcal{O}(1)$, $\mathcal{L}_{2}=\mathcal{O}(1)$, and
$\alpha_{i}^{2}\alpha_{k}^{2}(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{H}_{i,k}\mathbf{w}_{i}-(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,j}\mathbf{w}_{i})(\frac{1}{n}\sum_{i=1}^{n}\mathbf{H}_{i,k}\mathbf{w}_{i}))^{2}=\mathcal{O}(1)$.
Thus $\mathcal{L}_{1}=\mathcal{O}(p^{2})$. When
$p^{2}\gg\max(\lambda_{1},\lambda_{2})$, $\mathcal{L}_{1}$ will dominate the
regularization terms $\mathcal{L}_{2}$ and $\mathcal{L}_{3}$. ∎

## References

  * [1] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, “The graph neural network model,” _IEEE Transactions on Neural Networks_ , vol. 20, no. 1, pp. 61–80, 2008.
  * [2] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _ICLR_ , 2016.
  * [3] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio, and Y. Bengio, “Graph attention networks,” in _ICLR_ , 2017.
  * [4] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in _NeurIPS_ , 2017, pp. 1024–1034.
  * [5] K. Kuang, R. Xiong, P. Cui, S. Athey, and B. Li, “Stable prediction with model misspecification and agnostic distribution shift,” in _AAAI_ , 2020.
  * [6] Z. Shen, P. Cui, J. Liu, T. Zhang, B. Li, and Z. Chen, “Stable learning via differentiated variable decorrelation,” in _KDD_ , 2020, pp. 2185–2193.
  * [7] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Galligher, and T. Eliassi-Rad, “Collective classification in network data,” _AI magazine_ , vol. 29, no. 3, pp. 93–93, 2008.
  * [8] F. Wu, A. Souza, T. Zhang, C. Fifty, T. Yu, and K. Weinberger, “Simplifying graph convolutional networks,” in _ICML_. PMLR, 2019, pp. 6861–6871.
  * [9] B. Zadrozny, “Learning and evaluating classifiers under sample selection bias,” in _ICML_ , 2004, p. 114.
  * [10] Z. Yang, W. W. Cohen, and R. Salakhutdinov, “Revisiting semi-supervised learning with graph embeddings,” in _ICML_ , 2016.
  * [11] M. Nurhonen and S. Puntanen, “A property of partitioned generalized regression,” _Communications in statistics-theory and methods_ , vol. 21, no. 6, pp. 1579–1583, 1992.
  * [12] J. Hainmueller, “Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies,” _Political Analysis_ , vol. 20, no. 1, pp. 25–46, 2012.
  * [13] K. Kuang, L. Li, Z. Geng, L. Xu, K. Zhang, B. Liao, H. Huang, P. Ding, W. Miao, and Z. Jiang, “Causal inference,” _Engineering_ , vol. 6, no. 3, pp. 253–263, 2020.
  * [14] K. Kuang, P. Cui, S. Athey, R. Xiong, and B. Li, “Stable prediction across unknown environments,” in _SIGKDD_ , 2018, pp. 1617–1626.
  * [15] K. Kuang, P. Cui, H. Zou, B. Li, J. Tao, F. Wu, and S. Yang, “Data-driven variable decomposition for treatment effect estimation,” _TKDE_ , 2020.
  * [16] N. Kreif, R. Grieve, I. Díaz, and D. Harrison, “Evaluation of the effect of a continuous treatment: a machine learning approach with an application to treatment for traumatic brain injury,” _Health economics_ , vol. 24, no. 9, pp. 1213–1228, 2015.
  * [17] C. M. Bishop, _Pattern recognition and machine learning_. springer, 2006.
  * [18] H. Zou, P. Cui, B. Li, Z. Shen, J. Ma, H. Yang, and Y. He, “Counterfactual prediction for bundle treatment,” in _NeurIPS_ , vol. 33, 2020.
  * [19] A. Carlson, J. Betteridge, B. Kisiel, B. Settles, E. R. Hruschka, and T. M. Mitchell, “Toward an architecture for never-ending language learning,” in _Twenty-Fourth AAAI Conference on Artificial Intelligence_ , 2010.
  * [20] F. Zhou, T. Li, H. Zhou, H. Zhu, and J. Ye, “Graph-based semi-supervised learning with non-ignorable non-response,” in _NeurIPS_ , 2021.
  * [21] J. Klicpera, A. Bojchevski, and S. Günnemann, “Predict then propagate: Graph neural networks meet personalized pagerank,” in _ICLR_ , 2019.
  * [22] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _NeurIPS_ , 2016, pp. 3844–3852.
  * [23] J. Huang, A. Gretton, K. Borgwardt, B. Schölkopf, and A. Smola, “Correcting sample selection bias by unlabeled data,” in _NeurIPS_ , 2007, pp. 1–8.
  * [24] M. Cogswell, F. Ahmed, R. Girshick, L. Zitnick, and D. Batra, “Reducing overfitting in deep networks by decorrelating representations,” in _ICLR_ , 2016.
  * [25] X. Wang, S. Fan, K. Kuang, C. Shi, J. Liu, and B. Wang, “Decorrelated clustering with data selection bias,” in _IJCAI_ , 2020.
  * [26] K. Xu, W. Hu, J. Leskovec, and S. Jegelka, “How powerful are graph neural networks?” in _ICLR_ , 2019. [Online]. Available: https://openreview.net/forum?id=ryGs6iA5Km
  * [27] D. Jin, X. Song, Z. Yu, Z. Liu, H. Zhang, Z. Cheng, and J. Han, “Bite-gcn: A new gcn architecture via bidirectional convolution of topology and features on text-rich networks,” in _WSDM_ , 2021, pp. 157–165.
  * [28] S. Fan, J. Zhu, X. Han, C. Shi, L. Hu, B. Ma, and Y. Li, “Metapath-guided heterogeneous graph neural network for intent recommendation,” in _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2019, pp. 2478–2486.
  * [29] S. Fan, X. Wang, C. Shi, E. Lu, K. Lin, and B. Wang, “One2multi graph autoencoder for multi-view graph clustering,” in _Proceedings of The Web Conference 2020_ , 2020, pp. 3070–3076.
  * [30] L. Bai, L. Cui, Y. Jiao, L. Rossi, and E. Hancock, “Learning backtrackless aligned-spatial graph convolutional networks for graph classification,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2020.
  * [31] Y. Cao, H. Peng, J. Wu, Y. Dou, J. Li, and P. S. Yu, “Knowledge-preserving incremental social event detection via heterogeneous gnns,” in _WWW_ , 2021, pp. 3383–3395.
  * [32] D. Jin, Z. Liu, W. Li, D. He, and W. Zhang, “Graph convolutional networks meet markov random fields: Semi-supervised community detection in attribute networks,” in _AAAI_ , vol. 33, no. 01, 2019, pp. 152–159.
  * [33] X. Lin, C. Zhou, H. Yang, J. Wu, H. Wang, Y. Cao, and B. Wang, “Exploratory adversarial attacks on graph neural networks,” in _ICDM_. IEEE, 2020, pp. 1136–1141.
  * [34] Z. Shen, P. Cui, T. Zhang, and K. Kuang, “Stable learning via sample reweighting.” in _AAAI_ , 2020, pp. 5692–5699.
  * [35] Z. Ma, J.-H. Xue, A. Leijon, Z.-H. Tan, Z. Yang, and J. Guo, “Decorrelation of neutral vector variables: Theory and applications,” _IEEE transactions on neural networks and learning systems_ , vol. 29, no. 1, pp. 129–143, 2016.
  * [36] P. Rodríguez, J. Gonzalez, G. Cucurull, J. M. Gonfaus, and X. Roca, “Regularizing cnns with locally constrained decorrelations,” in _ICLR_ , 2016.
  * [37] Z. Zhang, Y. Zhang, and Z. Li, “Removing the feature correlation effect of multiplicative noise,” in _NeurIPS_ , 2018.
  * [38] Y. Nakatsukasa, “Absolute and relative weyl theorems for generalized eigenvalue problems,” _Linear Algebra and its Applications_ , vol. 432, no. 1, pp. 242–248, 2010.

|  Shaohua Fan received the B.E. degree in 2015 from Northeastern University
and M.S. degree in 2018 from Beijing University of Posts and
Telecommunications. He is a fourth-year Ph.D. student in the Department of
Computer Science of Beijing University of Posts and Telecommunications and
currently works as a visiting student at MILA. His main research interests
including graph mining and causal machine learning. He has published several
papers in major international conferences, including KDD, WWW, IJCAI, and CIKM
etc.
---|---
|  Xiao Wang is an Associate Professor in the School of Computer Science,
Beijing University of Posts and Telecommunications. He received his Ph.D.
degree from the School of Computer Science and Technology, Tianjin University,
Tianjin, China, in 2016. He was a postdoctoral researcher in Department of
Computer Science and Technology, Tsinghua University, Beijing, China. His
current research interests include data mining, social network analysis, and
machine learning. Until now, he has published more than 70 papers in refereed
journals and conferences.
---|---
|  Chuan Shi received the B.S. degree from the Jilin University in 2001, the
M.S. degree from the Wuhan University in 2004, and Ph.D. degree from the ICT
of Chinese Academic of Sciences in 2007. He is a professor and deputy director
of Beijing Key Lab of Intelligent Telecommunications Software and Multimedia
at present. His research interests are in data mining, machine learning, and
evolutionary computing. He has published more than 100 papers in refereed
journals and conferences.
---|---
|  Kun Kuang received his Ph.D. degree from Tsinghua University in 2019. He is
now an Associate Professor in the College of Computer Science and Technology,
Zhejiang University. He was a visiting scholar with Prof. Susan Athey’s Group
at Stanford University. His main research interests include Causal Inference,
Artificial Intelligence, and Causally Regularized Machine Learning. He has
published over 40 papers in major international journals and conferences,
including SIGKDD, ICML, ACM MM, AAAI, IJCAI, TKDE, TKDD, Engineering, and
ICDM, etc.
---|---
|  Nian Liu received the B.E. degree in 2020 from Beijing University of Posts
and Telecommunications. He is a second-year M.S. student in the Department of
Computer Science of Beijing University of Posts and Telecommunications. His
main research interests including graph mining and contrastive learning.
---|---
|  Bai Wang received the B.S. degree from the Xian Jiaotong University, Xian,
China and Ph.D. degree from the Beijing University of Posts and
Telecommunications, Beijing, China. And she is currently a professor of
computer science in BUPT. She was the director of Beijing Key Lab of
Intelligent Telecommunications Software and Multimedia.
---|---