Learning Attributed Graphlets:
Predictive Graph Mining by Graphlets with Trainable Attribute

We consider a classification problem in which a graph $G$ is an input. A set of nodes and edges of $G$ are written as $V_{G}$ and $E_{G}$, respectively. Each one of nodes $v\in V_{G}$ has a categorical label $L_{v}$ and a continuous attribute vector $\bm{z}^{G}_{v}\in\mathbb{R}^{d}$, where $d$ is an attribute dimension. In this paper, an attribute indicates a continuous attribute vector. We assume that a label and an attribute vector are for a node, but the discussion in this paper is completely same as for an edge label and attribute. A training dataset is $\{(G_{i},y_{i})\}_{i\in[n]}$, in which $y_{i}\in\{-1,+1\}$ is a binary label and $n$ is the dataset size, where $[n]=\{1,\ldots,n\}$. Although we only focus on the classification problem, our framework is also applicable to the regression problem just by replacing the loss function.

We consider extracting important small attributed graphs, which we call attributed graphlets (AGs), that contributes to the classification boundary. Note that throughout the paper, we only consider a connected graph as an AG for a better interpretability (do not consider an AG by a disconnected graph). Let $\psi(G_{i};H)\in[0,1]$ be a feature representing a degree that an input graph includes an AG $H$. We refer to $\psi(G_{i};H)$ as the AG inclusion score (AGIS). Our proposed LAGRA identifies important AGs by applying a feature selection to a model with this AGIS feature.

Suppose that $L(G)$ is a labeled graph having a categorical label $L_{v}$ for each node, and in $L(G)$, an attribute $\bm{z}^{G}_{v}$ for each node is excluded from $G$. We define AGIS so that it has a non-zero value only when $L(H)$ is included in $L(G_{i})$:

where $L(H)\sqsubseteq L(G_{i})$ means that $L(H)$ is a subgraph of $L(G_{i})$, and $\phi_{H}(G_{i})\in(0,1]$ is a function that provides a continuous inclusion score of $H$ in $G_{i}$. The condition $L(H)\sqsubseteq L(G_{i})$ makes AGIS highly interpretable. For example, in the case of chemical composition data, if $L(H)$ represents O-C (oxygen and carbon are connected) and $\psi(G_{i};H)>0$, then we can guarantee that $G_{i}$ must contain O-C. Figure 3(a) shows examples of $L(G_{i})$ and $L(H)$.

The function $\phi_{H}(G_{i})$ needs to be defined so that it can represent how strongly the attribute vectors in $H$ can be matched to those of $G_{i}$. When $L(H)\sqsubseteq L(G_{i})$, there exists at least one injection $m:V_{H}\rightarrow V_{G_{i}}$ in which $m(v)$ for $v\in V_{H}$ preserves node labels and edges among $v\in V_{H}$. Figure 3(b) shows an example of when there exist two injections. Let $M$ be a set of possible injections $m$. We define a similarity between $H$ and a subgraph of $G_{i}$ matched by $m\in M$ as follows

where $\rho>0$ is a fixed parameter that adjusts the length scale. In $\exp$, the sum of squared distances of attribute vectors between matched nodes are taken. To use this similarity in AGIS (1), we take the maximum among all the matchings $M$:

An intuition behind (2) is that it evaluates inclusion of $H$ in $G_{i}$ based on the best macthing in a sense of $\mathrm{Sim}(H,G_{i};m)$ for $m\in M$. If $L(H)\sqsubseteq L(G_{i})$ and there exists $m$ such that $\bm{z}^{H}_{v}=\bm{z}^{G_{i}}_{m(v)}$ for $\forall v\in V_{H}$, then, $\phi_{H}(G_{i})$ takes the maximum value (i.e., $1$).

Our prediction model linearly combines the feature $\psi(G_{i};H)$ as follows:

where $\beta_{H}$ and $\beta_{0}$ are parameters, ${\mathcal{H}}$ is a set of candidate AGs, and $\bm{\beta}$ and $\bm{\psi}_{i}$ are vectors containing $\beta_{H}$ and $\psi(G_{i};H)$ for $H\in{\mathcal{H}}$, respectively. Let ${\mathcal{L}}=\{L\mid L\subseteq L(G_{i}),i\in[n],|L|\leq\text{maxpat}\}$ be a set of all the labeled subgraphs contained in the training input graphs $\{G_{i}\}_{i\in[n]}$, where $|L|$ is the number of nodes in the labeled graph $L$ and maxpat is the user-specified maximum size of AGs. The number of the candidate AGs $|{\mathcal{H}}|$ is set as the same size as $|{\mathcal{L}}|$. We set ${\mathcal{H}}$ as a set of attributed graphs created by giving trainable attribute vectors $\bm{z}_{v}^{H}(v\in V_{H})$ to each one of elements in ${\mathcal{L}}$. Figure 4 shows a toy example. Our optimization problem for $\beta_{H}$, $\beta_{0}$ and $\bm{z}_{v}^{H}$ is defined as the following regularized loss minimization in which the sparse $L1$ penalty is imposed on $\beta_{H}$:

where ${\mathcal{Z}}_{\mathcal{H}}=\{\bm{z}^{H}_{v}\mid v\in V_{H},H\in{\mathcal{H}}\}$, and $\ell$ is a convex differentiable loss function. Here, we employ the squared hinge loss function [10]:

Since the objective function (4) induces a sparse solution for $\beta_{H}$, we can identify a small number of important AGs as $H$ having the non-zero $\beta_{H}$. However, this optimization problem has an intractably large number of optimization variables (${\mathcal{H}}$ contains all the possible subgraphs in the training dataset and each one of $H\in{\mathcal{H}}$ has the attribute vector $\bm{z}_{v}^{H}\in\mathbb{R}^{d}$ for each one of nodes). We propose an efficient optimization algorithm that mitigates this problem.

Before introducing the subset ${\mathcal{W}}$, we first describe update rules of each variable. First, we apply the proximal gradient update to $\bm{\beta}$. Let

be the derivative of the loss term in (4) with respect to $\beta_{H}$. Then, the update of $\bm{\beta}$ is defined by

where $\eta>0$ is a step length, and

is a proximal operator (Note that the proximal gradient for the $L1$ penalty is often called ISTA [9], for which an accelerated variant called FISTA is also known. We here employ ISTA for simplicity). We select the step length $\eta$ by the standard backtrack search.

The bias term $\beta_{0}$ is update by

which is the optimal solution of the original problem (4) for given other variables $\bm{\beta}$ and ${\mathcal{Z}}_{\mathcal{H}}$. Since the objective function of $\beta_{0}$ is a differential convex function, the update rule of $\beta_{0}$ can be derived from the first order condition as

where ${\mathcal{I}}^{\rm(new)}=\{i\mid 1-y_{i}(\bm{\psi}_{i}^{\top}\bm{\beta}+\beta^{\rm(new)}_{0})>0\}$. This update rule contains $\beta^{\rm(new)}_{0}$ in ${\mathcal{I}}^{\rm(new)}$. However, it is easy to calculate the update (6) without knowing $\beta^{\rm(new)}_{0}$ beforehand. Here, we omit detail because it is a simple one dimensional problem (see supplementary appendix A).

For $\bm{z}^{H}_{v}\in{\mathcal{Z}}_{{\mathcal{H}}}$, we employ the standard gradient descent:

where $\alpha>0$ is a step length to which we apply the standard backtrack search.

In every update of $\bm{\beta}$, we incrementally add required $H$ into ${\mathcal{W}}\subseteq{\mathcal{H}}$. For the complement set $\overline{{\mathcal{W}}}={\mathcal{H}}\setminus{\mathcal{W}}$, which contains AGs that have never been updated, we initialize $\beta_{H}=0$ for $H\in\overline{{\mathcal{W}}}$. For the initialization of a node attribute vector $\bm{z}_{v}^{H}\in{\mathcal{Z}}_{\mathcal{H}}$, we set the same initial vector if the node (categorical) labels are same, i.e., $\bm{z}_{v}^{H}=\bm{z}_{v^{\prime}}^{H^{\prime}}$ if $L_{v}=L_{v^{\prime}}$ for $\forall H,H^{\prime}\in{\mathcal{H}}$ (in practice, we use the average of the attribute vectors within each node label). This constraint is required for our pruning criterion, but it is only for initial values. After the update (7), all $\bm{z}^{H}_{v}$ can have different values.

Since $\beta_{H}=0$ for $H\in\overline{{\mathcal{W}}}$, it is easy to derive the following relation from the proximal update (5):

This indicates that if the conditions in the left side hold, we do not need to update $\beta_{H}$ because it remains $0$. Therefore, we set

and apply the update (5) only to $H\in{\mathcal{W}}$. However, evaluating $|g_{H}(\bm{\beta})|>\lambda$ for all $H\in\overline{{\mathcal{W}}}$ can be computationally intractable because it needs to enumerate all the possible subgraphs. The following theorem can be used to avoid this difficulty:

Figure 5 shows an illustration of the forward and backward (gradient) computations of LAGRA. For the gradient pruning, an efficient algorithm can be constructed by combining the rule (10) and a graph mining algorithm. A well-known efficient graph mining algorithm is gSpan [11], which creates the tree by recursively expanding each graph in the tree node as far as the expanded graph is included in a given set of graphs as a subgraph. An important characteristics of the mining tree is that all graphs must contain any graph of its ancestors as subgraphs. Therefore, during the tree traverse (depth-first search) by gSpan, we can prune the entire subtree (all descendant nodes) if $\overline{g}_{H}(\bm{\beta})\leq\lambda$ holds for the AG $H$ in a tree node (Figure 5(b)). This means that we can update ${\mathcal{W}}$ by (9) without exhaustively investigating all the elements in $\overline{{\mathcal{W}}}$.

gSpan has another advantage for LAGRA. To calculate the feature $\phi_{H}(G_{i})$, defined in (2), LAGRA requires a set of injections $M$ (an example is shown in Figure 3(b)). gSpan keeps the information of $M$ during the tree traverse because it is required to expand a subgraph in each $G_{i}$ (see the authors implementation https://sites.cs.ucsb.edu/~xyan/software/gSpan.htm). Therefore, we can directly use $M$ created by gSpan to calculate (2).

We here describe entire procedure of the optimization of LAGRA. We employ the so-called regularization path following algorithm (e.g., [12]), in which the algorithm starts from a large value of the regularization parameter $\lambda$ and gradually decreases it while solving the problem for each $\lambda$. This strategy can start from highly sparse $\bm{\beta}$, in which usually ${\mathcal{W}}$ also becomes small. Further, at each $\lambda$, the solution obtained in the previous $\lambda$, can be used as the initial value by which faster convergence can be expected (so-called warm start).

Algorithm 1 shows the procedure of the regularization path following. We initialize $\bm{\beta}=\bm{0}$, which is obviously optimal when $\lambda=\infty$. In line 2 of Algorithm 1, we calculate $\lambda_{\max}$ at which $\bm{\beta}$ starts having non-zero values: $\lambda_{\max}=\max_{H\in{\mathcal{H}}}\left|\sum_{i\in[n]}^{n}y_{i}\psi(G_{i};H)(1-y_{i}\beta_{0})\right|$, where $\beta_{0}=\sum_{i\in[n]}y_{i}/n$. See supplementary appendix C for derivation. $\lambda_{\rm max}$ can also be written as $\lambda_{\max}=\max_{H\in{\mathcal{H}}}|g_{H}(\bm{0})|$. To find $\max_{H\in{\mathcal{H}}}$, we can use almost the same gSpan based pruning strategy by using an upper bound of $\overline{g}_{H}(\bm{\beta})$ as shown in Section 2.2.2 (the only difference is to search the $\max$ value only, instead of searching all $H$ satisfying $|g_{H}(\bm{\beta})|>\lambda$), though in Algorithm 1, this process is omitted for brevity. After setting $\lambda_{0}\leftarrow\lambda_{\max}$, the regularization parameter $\lambda$ is decreased by using a pre-defined decreasing factor $R$ as shown in line 6 of Algorithm 1. For each $\lambda_{1}>\cdots>\lambda_{K}$, the parameters $\bm{\beta},\beta_{0}$ and ${\mathcal{Z}}_{\mathcal{H}}$ are alternately updated as described in Section 2.2.1 and 2.2.2. We stop the alternate update by monitoring performance on the validation dataset in line 14 (stop by thresholding the decrease of the objective function is also possible).

The algorithm of the pruning strategy described in Section 2.2.2 is shown in Algorithm 2. This function recursively traverses the graph mining tree. At each tree node, first, $\overline{g}_{H}(\bm{\beta})$ is evaluated to prune the subtree if possible. Then, if $|g_{H}(\bm{\beta})|>\lambda_{k}$, $H$ is included in ${\mathcal{W}}$. The expansion from $H$ (creating children of the graph tree) is performed by gSpan, by which only the subgraphs contained in the training set can be generated (see the original paper [11] for detail of gSpan). The initialization of the trainable attribute $\bm{z}_{v}^{H^{\prime}}$ is performed when $H^{\prime}$ is expanded (line 15).