GNNs have two inputs: first, a feature matrix $X\in {?}_N\subset {ℝ}^{N\times n}$ where $n$ is the number of different node features and $N$ is the total number of nodes, second, a graph $G=(V,E)$ with adjacency matrix $A\in {\mathscr{H}}_N={\{0,1\}}^{N\times N}$. Depending on the particular downstream application of interest, GNNs may have different outputs; we address node level semi-supervised classification problems, whereby the task is to learn a function that predicts node labels in a given set $?$ (e.g. topics of a scientific paper). An example of such type of models, proposed by [8], is the following two hidden layer graph convolutional network (GCN) that computes the class probabilities as

where $w=(W_1,W_2)\in {ℝ}^d$ are the model parameters. The process of learning $f_w$ from a small subset of training nodes $V_{\mathrm{?????}}$ (for which the ground truth label is known) is most often cast to an optimization problem of the form

where $\mathrm{\ell }:?\times ?\to {ℝ}^{+}$ is a point-wise loss function, and $\mathrm{\Omega }$ is a regularizer.

Now, let us suppose that the information about the true adjacency matrix $A$ is missing or incomplete. Since, ultimately, we are interested in finding a model that minimizes the generalization error, we assume the existence of a second subset of instances with known target, $V_{\mathrm{???}}$ (the validation set), from which we can estimate the generalization error. Hence, we propose to find $A\in {\mathscr{H}}_N$ that minimizes the function

where $w_A$ is the minimizer, assumed unique, of $L$ for a fixed adjacency matrix$A$. We can then consider Equations (2) and (3) as the inner and outer objective of a mixed-integer bilevel programming problem where the outer objective aims to find an optimal discrete graph structure and the inner objective the optimal parameters of a GCN given a graph.

The resulting bilevel problem is intractable to solve exactly even for small graphs due to its high dimensional combinatorial nature. To circumvent this difficulty we introduce a generative model and reformulate the bilevel program in terms of the (continuous) parameters of the resulting distribution over discrete graphs. Specifically, we propose to model each edge with a Bernoulli random variable. Let ${\overline{\mathscr{H}}}_N=\mathrm{Conv}({\mathscr{H}}_N)$ be the convex hull of the set of all adjacency matrices for $N$ nodes. By modeling all possible edges as a set of mutually independent Bernoulli random variables with parameter matrix $\theta \in {\overline{\mathscr{H}}}_N$ we can sample graphs as ${\mathscr{H}}_N\ni A\sim \mathrm{Ber}(\theta )\equiv P_{\theta }$. Eqs. (2) and (3) can then be replaced by

where we take the expectation of the objectives over graph structures. In this way both the inner and the outer objectives become continuous and smooth functions of the Bernoulli parameters.

The bilevel problem given by Eqs. (4)-(5) is still challenging to solve efficiently. This is because the solution of the inner problem is not available in closed form for GCNs; and the expectations are intractable to compute exactly⁴⁴4 This is different than e.g. model free reinforcement learning, where the objective function is usually unknown, depending in an unknown way from the action and the environment.. We now discuss a practical algorithm to approach the problem. Regarding the inner problem (5), we can choose a tractable approximate learning dynamics $\mathrm{\Phi }$ such as stochastic⁵⁵5 The stochasticity arises from the graph distribution. gradient descent (SGD) [3],

where ${\gamma }_t$ is a learning rate and $A_t\sim P_{\theta }$ is drawn at each iteration.

Let $w_{\theta ,T}$ be an approximate minimizer of $?\left[L\right]$ (where $T$ may depend on $\theta$). We now need to compute an estimator for the hypergradient:

We can swap gradient and expectation operators since $\mathrm{Ber}(\theta )$ is finite and we assume $F$ bounded. Next, we use the so-called straight-through estimator [2] and set ${\nabla }_{\theta }A:=I$ (which would be $0$ a.e. otherwise) and propagate the gradient through (portions of) the learning dynamics to compute ${\nabla }_{\theta }{w}_{\theta ,T}$[10, 5]. Finally, we take the single sample Monte Carlo estimator of (7) to update the parameters $\theta$ with projected gradient descent on the unit hypercube.

A sketch of the method is presented in Algorithm 1, where inputs and operations in squared brackets are optional.

The expected output of a GCN trained with LDS is

which can be cheaply estimated in an unbiased manner as

where $S>0$ is the number of samples we wish to draw. Given the parametrization of the graph generator with Bernoulli variables, one can sample a new graph in $O(N^2)$. Sampling from a large number of Bernoulli variables, however, is highly efficient, trivially parallelizable, and possible at a rate of millions per second.

Learning discrete distributions has a computational advantage over dense graphs if the sampled graphs are sparse enough: indeed for GCNs computing ${\widehat{f}}_w$ has a cost of $O(SCd)$, opposed to $O(N^{2}d)$ for a fully connected graph, where $C={\sum }_{ij}{\theta }_{ij}$ is the expected number of edges, and $d$ is the dimension of the weights. Another advantage of using a graph-generative model is that we can interpret it probabilistically which is not the case when learning a dense adjacency matrix.

On the incomplete graphs scenarios (Figure 2) LDS achieves competitive results in all settings and shows accuracy gains of up to $7$ percentage points. Notably, LDS improves the generalization accuracy of GCN models also when the given graph is that of the respective dataset (100% of edges retained), by learning a handful of additional, helpful, edges. The accuracy of 84.1% and 75.0% for Cora and Citeseer, respectively, exceed all previous state-of-the-art results at the time of publication. The learned graphs are very sparse: e.g. for Cora, on average, less than $0.2\%$ edges are present.

For semi-supervised classification problems (Table 1), we noted that supervised learning baselines work well on some datasets but fail to provide competitive results on others. The semi-supervised learning baselines LP [17], ManiReg [1] and SemiEmb [15] can only improve the supervised learning baselines on $1$, $3$ and $4$ datasets, respectively. The results for the GCN with different input graphs show that $k$NN-GCN works well and provides competitive results compared to the supervised baselines on all datasets. $k$NN-LDS significantly outperforms $k$NN-GCN on $4$ out of the $7$ datasets. In addition, $k$NN-LDS is among the most competitive methods on all datasets and yields the highest gains on datasets that have an underlying graph.