{'fraction_non_alphanumeric': 0.04545795701883866, 'fraction_numerical': 0.014498011058017062, 'mean_word_length': 4.302666489785089, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 21, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Distance Metric Learning for Graph Structured Data

Graphs are versatile tools for representing structured data. As a result, a variety of machine learning methods have been studied for graph data analysis. Although many such learning methods depend on the measurement of differences between input graphs, defining an appropriate distance metric for graphs remains a controversial issue. Hence, we propose a supervised distance metric learning method for the graph classification problem. Our method, named interpretable graph metric learning (IGML), learns discriminative metrics in a subgraph-based feature space, which has a strong graph representation capability. By introducing a sparsity-inducing penalty on the weight of each subgraph, IGML can identify a small number of important subgraphs that can provide insight into the given classification task. Because our formulation has a large number of optimization variables, an efficient algorithm that uses pruning techniques based on safe screening and working set selection methods is also proposed. An important property of IGML is that solution optimality is guaranteed because the problem is formulated as a convex problem and our pruning strategies only discard unnecessary subgraphs. Furthermore, we show that IGML is also applicable to other structured data such as itemset and sequence data, and that it can incorporate vertex-label similarity by using a transportation-based subgraph feature. We empirically evaluate the computational efficiency and classification performance of IGML on several benchmark datasets and provide some illustrative examples of how IGML identifies important subgraphs from a given graph dataset.


Introduction
Because of the growing diversity of data science applications, machine learning methods must adapt to a variety of complicated structured data, from which it is often difficult to obtain typical numerical vector representations of input objects. A standard approach to modeling structured data is to employ graphs. For example, graph-based representations are prevalent in domains such as chemo-and bio-informatics. In this study, we particularly focus on the case in which a data instance is represented as a pair of a graph and its associated class label.
Although numerous machine learning methods explicitly or implicitly depend on how to measure differences between input objects, defining an appropriate distance metric on graphs remains a controversial issue in the field. A widely accepted approach is the graph kernel (Gärtner et al., 2003;Vishwanathan et al., 2010), which enables to apply machine learning methods to graph data without requiring explicit vector representations. Another popular approach would be to use neural networks (Atwood and Towsley, 2016;Narayanan et al., 2017), from which a suitable representation can be learned while avoiding to explicitly define a metric. However, in these approaches, it is difficult to create a metric that explicitly extracts significant sub-structures, i.e., subgraphs. Identifying discriminative subgraphs in an interpretable manner can be insightful for many graph classification tasks. In particular, graph representation is prevalent in scientific data analysis. For example, chemical compounds are often represented by graphs; thus, finding subgraphs that have a strong effect on a target label (e.g., toxicity) is informative. Other examples of graph representations are protein 3D structures and crystalline substances (e.g., Brinda and Vishveshwara, 2005;Xie and Grossman, 2018), where the automatic identification of important sub-structures is expected to provide an insight behind correlation between structures and target labels. Further details of the previous studies are discussed in Section 2.
We propose a supervised method that obtains a metric for graphs, thereby achieving both high predictive performance and interpretability. Our method, named interpretable graph metric learning (IGML), combines the concept of metric learning (e.g., Weinberger and Saul, 2009;Davis et al., 2007) with a subgraph representation, where each graph is represented by a set of its subgraphs. IGML optimizes a metric that assigns a weight m i(H) ≥ 0 to each subgraph H contained in a given graph G. Let φ H (G) be a feature of the graph G that is monotonically non-decreasing with respect to the frequency of subgraph H of G. Note that we assume that subgraphs are counted without overlapped vertices and edges throughout the study. We consider the following squared distance between two graphs G and G : where G is the set of all connected graphs. Although it is known that the subgraph approach has strong graph representation capability (e.g. Gärtner et al., 2003), naïve calculation is obviously infeasible unless the weight parameters have some special structure. We formulate IGML as a supervised learning problem of the distance function (1) using a pairwise loss function of metric learning (Davis et al., 2007) with a sparse penalty on m i(H) . The resulting optimization problem is computationally infeasible at a glance, because the number of weight parameters is equal to the number of possible subgraphs, which is usually intractable. We overcome this difficulty by introducing safe screening (Ghaoui et al., 2010) and working set selection (Fan et al., 2008) approaches. Both of these approaches can significantly reduce the number of variables, and further, they can be combined with a pruning strategy on the tree traverse of graph mining. These optimization tricks are inspired by two recent studies (Nakagawa et al., 2016) and (Morvan and Vert, 2018), which developed safe screening-and working set-based pruning for a linear prediction model with the LASSO penalty, respectively. By combining these two techniques, we constructed a path-wise optimization method that can obtain a sparse solution of the weight parameter m i(H) without directly enumerating all possible subgraphs.
To the best of our knowledge, no previous studies can provide an interpretable subgraph-based metric learned in a supervised manner. The advantages of IGML can be summarized as follows: • Because IGML is formulated as a convex optimization problem, the global optimal can be found by the standard gradient-based optimization.
• The safe screening-and working set-based optimization algorithms make our problem practically tractable without sacrificing optimality.
• We can identify a small number of important subgraphs that discriminate different classes. This implies that the resulting metric is easy to compute and highly interpretable, making it useful for a variety of subsequent data analyses. For example, applying the nearest neighbor classification or decision tree on the learned space would be effective.
Moreover, we propose three extensions of IGML. First, we show that IGML is directly applicable to other structured data, such as itemset and sequence data. Second, its application to a triplet based loss function is discussed. Third, we extend IGML to allow similarity information of vertex-labels to be incorporated. We empirically verify the superior or comparable prediction performance of IGML to other existing graph classification methods (most of which are not interpretable). We also show some examples of extracted subgraphs and data analyses on the learned metric space. The reminder of this paper is organized as follows. In Section 2, we review previous studies on graph data analysis. In Section 3, we introduce a formulation of our proposed IGML. Section 4 discusses strategies to reduce the size of the IGML optimization problem. The detailed computational procedure of IGML is described in Section 5. Three extensions of IGML are presented in Section 6. Section 7 reports our empirical evaluation of the effectiveness of IGML on several benchmark datasets.
Note that this paper is an extended version of a preliminary conference paper (Yoshida et al., 2019a). The source code of the program used in our experiments is available at https://github.com/takeuchi-lab/Learning-Interpretable-Metric-between-Graphs.

Related Work
Kernel-based approaches have been widely studied for graph data analysis, and they can provide a metric of graph data in a reproducing kernel Hilbert space. In particular, subgraph-based graph kernels are closely related to our study. The graphlet kernel ) creates a kernel through small subgraphs with only about 3-5 vertices, which are called graphlets. The neighborhood subgraph pairwise distance kernel (Costa and Grave, 2010) selects pairs of subgraphs from a graph and counts the number of pairs identical to those in another graph. The subgraph matching kernel (Kriege and Mutzel, 2012) identifies common subgraphs based on cliques in the product graph of two graphs.
The feature space created by these subgraph-based kernels is easy to interpret. However, because the above approaches are unsupervised, it is fundamentally impossible to eliminate subgraphs that are unnecessary for a specific target classification task. Therefore, for example, to create the entire kernel matrix of training data, all the candidate subgraphs in the data must be enumerated once, which becomes intractable even for small-sized subgraphs. In contrast, we consider dynamically "pruning" unnecessary subgraphs through a supervised formulation of metric learning. As we will demonstrate in our later experimental results, this significantly reduces the enumeration cost, allowing our proposed algorithm to deal with the larger size of subgraphs than the simple subgraph based kernels.
There are many other kernels including the shortest path (Borgwardt and Kriegel, 2005)-, random walk (Vishwanathan et al., 2010;Sugiyama and Borgwardt, 2015;Zhang et al., 2018b)-, and spectrum-based (Kondor and Borgwardt, 2008;Kondor et al., 2009;Kondor and Pan, 2016;Verma and Zhang, 2017) approaches. The Weisfeiler-Lehman (WL) kernel Shervashidze et al., 2011), which is based on the graph isomorphism test, is a popular and empirically successful kernel that has been employed in many studies (Yanardag and Vishwanathan, 2015;Niepert et al., 2016;Narayanan et al., 2017;Zhang et al., 2018a). Again, all such approaches are unsupervised, and it is difficult to interpret results from the perspective of sub-structures of a graph. Although several kernels deal with continuous attributes on vertices (Feragen et al., 2013;Orsini et al., 2015;Su et al., 2016;Morris et al., 2016), we only focus on the cases where vertex-labels are discrete due to the associated interpretability.
Because obtaining a good metric is an essential task in data analysis, metric learning has been extensively studied to date, as reviewed in (Li and Tian, 2018). However, due to its computational difficulty, metric learning for graph data has not been widely studied. A few studies have considered the edit distance approaches. For example, Bellet et al. (2012) presented a method for learning a similarity function through an edit distance in a supervised manner. Another approach probabilistically formulates the editing process of the graph and estimates the parameters using labeled data Neuhaus and Bunke (2007). However, these approaches cannot provide any clear interpretation of the resulting metric in term of the subgraphs.
Likewise, the deep neural network (DNN) is a standard approach to graph data analysis. The deep graph kernel (Yanardag and Vishwanathan, 2015) incorporates neural language modeling, where decomposed substructures of a graph are regarded as sentences. The PATCHY-SAN (Niepert et al., 2016) and DGCNN (Zhang et al., 2018a) convert a graph to a tensor by using the WL-Kernel and convolute it. Several other studies also have combined popular convolution techniques with graph data (Tixier et al., 2018;Atwood and Towsley, 2016;Simonovsky and Komodakis, 2017). These approaches are supervised, but the interpretability of these DNNs is obviously relatively low. Attention enhances the interpretability of deep learning, but extracting important subgraphs is difficult because attention algorithms for graphs (Lee et al., 2018) only provides the significance of vertex transition on a graph. Another related DNN approach is representation learning. For example, sub2vec (Adhikari et al., 2018) and graph2vec (Narayanan et al., 2017) can embed graph data into a continuous space, but they are unsupervised, and it is difficult to extract substructures that characterize different classes. There are other fingerprint learning methods for graphs by neural networks (e.g. Duvenaud et al., 2015) where the contribution from each node can be evaluated for each dimension of the fingerprint. Although it is possible to highlight sub-structures for the given input graph, this does not produce important common subgraphs for prediction. Supervised pattern mining (Cheng et al., 2008;Novak et al., 2009;Thoma et al., 2010) can be used for identifying important subgraphs by enumerating patterns with some discriminative score. However, these approaches usually 1) employ a greedy strategy to add a pattern for which global optimality cannot be guaranteed, and 2) do not optimize a metric or representation.
A few other studies (Saigo et al., 2009;Nakagawa et al., 2016) have considered optimizing a linear model on the subgraph features with the LASSO penalty using graph mining. A common idea of these two methods is to traverse a graph mining tree with pruning strategies derived from optimality conditions. Saigo et al. (2009) employed a boosting-based approach, which adds a subgraph that violates the optimality condition most severely at every iteration. It was shown that the maximum violation condition can be efficiently identified by pruning the tree without losing the final solution optimality. Nakagawa et al. (2016) derived a pruning criterion by extending safe screening (Ghaoui et al., 2010), which can safely eliminate unnecessary features before solving the optimization problem. This approach can also avoid enumerating the entire tree while guaranteeing the optimality, and its efficiency compared with the boosting-based approach was demonstrated empirically, mainly because it requires much fewer tree traversals. Further, Morvan and Vert (2018) proposed a similar pruning extension of working set selection for optimizing a higher-order interaction model. Although this paper was not for the graph data, the technique is applicable to the same subgraph-based linear model as in (Saigo et al., 2009) and (Nakagawa et al., 2016). Working set selection is a heuristic feature subset selection strategy that has been widely used in machine learning algorithms, such as support vector machines (e.g., Hsu and Lin, 2002). Unlike safe screening, this heuristic selection may eliminate necessary features in the middle of the optimization, but the optimality of the final solution can be guaranteed by iterating subset selection repeatedly until the solution converges. However, these methods can only optimize a linear prediction model. In this study, we focus on metric learning of graphs. Therefore, unlike the above mentioned pruning based learning methods, our aim is to learn a "distance function". In metric learning, a distance function is typically learned from a loss function defined over a relative relation between samples (usually, pairs or triplets), by which a discriminative feature space that is generally effective for subsequent tasks, such as classification and similarity-based retrieval, is obtained. Inspired by (Nakagawa et al., 2016) and (Morvan and Vert, 2018), we derive screening and pruning rules for this setting , and further, we combine them to develop an efficient algorithm.

Optimization problem
Suppose that the training dataset {(G i , y i )} i∈ [n] consists of n pairs of a graph G i and a class label y i , where [n] := {1, . . . , n}. Let G be the set of all connected subgraphs of {G i } i∈ [n] . In each graph, vertices and edges can be labeled. If H ∈ G is a connected subgraph of G ∈ G, we write H G. Further, let #(H G) be the frequency of the subgraph H in G. Note that we adopt a definition of frequency that does not allow any vertices or edges among the counted subgraphs to overlap. As a representation of a graph G, we consider the following subgraph-based feature representation: where g is some monotonically non-decreasing and non-negative function, such as the identity function g(x) = x or indicator function g(x) = 1 x>0 , which takes the value 1 if x > 0, and 0 otherwise. It is widely known that subgraph-based features can effectively represent graphs. For example, g(x) = x allows all nonisomorphic graphs to be distinguished. A similar idea was shown in (Gärtner et al., 2003) for a frequency that allows overlaps. However, this feature space is practically infeasible because the possible number of subgraphs is prohibitively large. We focus on how to measure the distance between two graphs, which is essential for a variety of machine learning problems. We consider the following weighted squared distance between two graphs: and D i ⊆ [n] be the subsets of indices that are in the same and different classes to x i , respectively. For each of these sets, we select the K most similar inputs to x i by using some default metric, such as the graph kernel (further details are presented in Section 3.2). As a loss function for x i , we consider where L, U ∈ R + are constant parameters satisfying U ≤ L, and t (x) = [t − x] 2 + is the standard squared hinge loss function with threshold t ∈ R. This loss function is a variant of the pairwise loss functions used in metric learning (Davis et al., 2007). The first term in the loss function yields a penalty if x i and x l are closer than L for l ∈ D i , and the second term yields a penalty if x i and x j are more distant than U for j ∈ S i . Let R(m) = m 1 + η 2 m 2 2 = m 1 + η 2 m 2 2 be an elastic-net type sparsity-inducing penalty, where η ≥ 0 is a non-negative parameter. We define our proposed IGML (interpretable graph metric learning) as the following regularized loss minimization problem: where λ > 0 is the regularization parameter. The solution of this problem can provide not only a discriminative metric but also insight into important subgraphs because the sparse penalty is expected to select only a small number of non-zero parameters. Let α ∈ R 2nK + be the vector of dual variables where α il and α ij for i ∈ [n], l ∈ D i , and j ∈ S i are concatenated. The dual problem of (4) is written as follows (see Appendix A for derivation): where . .] ∈ R p×2nK . Then, from the optimality condition, we obtain the following relationship between the primal and dual variables: where t (x) = −2[t − x] + is the derivative of t . When the regularization parameter λ is larger than certain λ max , the optimal solution is m = 0. Then, the optimal dual variables are α il = − L (0) = 2L and α ij = − −U (0) = 0. By substituting these equations into (6), we obtain λ max as 3.2 Selection of S i and D i For K = |S i | = |D i |, in the experiments reported later, we employed the small number K = 10 and used a graph kernel to select samples in S i and D i . Although we simply used a pre-determined kernel, selecting the kernel (or its parameter) through cross-validation beforehand is also possible. Using only a small number of neighbors is a common setting in metric learning. For example, Davis et al. (2007), which is a seminal work on the pairwise approach, only used 20c 2 pairs in total, where c is the number of classes. A small K setting has two aims. First, particularly S i , adding pairs that are too far apart can be avoided under this setting. Even for a pair of samples with the same labels, enforcing such distant pairs to be close may cause over-fitting (e.g., when the sample is an outlier). Second, a small K reduces the computational cost. Because the number of pairs is O(n 2 ), adding all of them into the loss term requires a large computational cost. In fact, these two issues are not only for the pairwise formulation but also for other relative loss functions such as the standard triplet loss, for which there exist O(n 3 ) triplets. One potential difficulty in selecting D i and S i is the discrepancy between the initial and the optimal metric. The loss function is defined through D i and S i , which are selected based on the neighbors in the initial metric, but the optimization of the metric may change the nearest neighbors of each sample. A possible remedy for this problem is to adaptively change D i and S i in accordance with the updated metric (Takeuchi and Sugiyama, 2011), though the resulting optimality of this approach is still not known. To the best of our knowledge, this is still an open problem in metric learning, which we consider beyond the scope of this paper. In the experiments (Section 7), we show that a nearest-neighbor classifier in the learned metric with this heuristics selection of D i and S i shows better or comparable performance to standard graph classification methods, such as a graph neural network.

Creating a Tractable Sub-problem
Because the problems of (4) and (5) are convex, the local solution is equivalent to the global optimal. However, naïvely solving these problems is computationally intractable because of the high dimensionality of m. In this section, we introduce several useful rules for restricting candidate subgraphs while maintaining the optimality of the final solution. Note that the proofs for all the lemmas and theorems are provided in the appendix.
To make the optimization problem tractable, we work with only a small subset of features during the optimization process. Let F ⊆ [p] be a subset of features. By fixing m i = 0 for i / ∈ F, we define sub-problems of the original primal P λ and dual D λ problems as follows: where m F , c ijF , and m λ (α) F are sub-vectors specified by F. If the size of F is moderate, these sub-problems are significantly computationally easier to solve than the original problems. We introduce several criteria that determine whether the feature k should be included in F using the techniques of safe screening (Ghaoui et al., 2010) and working set selection (Fan et al., 2008). A general form of our criteria can be written as C k,: q + r C k,: 2 ≤ T, where q ∈ R 2nK + , r ≥ 0, and T ∈ R are constants that assume different values for each criterion. If this inequality holds for k, we exclude the k-th feature from F. An important property is that although our algorithm only solves these small sub-problems, we can guarantee the optimality of the final solution, as shown later.
However, selecting F itself is computationally expensive because the evaluation of (11) requires O(n) computations for each k. Thus, we exploit a tree structure of graphs for determining F. Figure 1 shows an example of such a tree, which can be constructed by a graph mining algorithm, such as gSpan (Yan and Han, 2002). Suppose that the k-th node corresponds to the k-th dimension of x (note that the node index here is not the order of the visit). If a graph corresponding to the k-th node is a subgraph of the k -th node, the node k is a descendant of k, which is denoted as k ⊇ k. Then, the following monotonic relation is immediately derived from the monotonicity of φ H : ... Because any parent node is a subgraph of its children in the gSpan tree Figure 1, the non-overlapped frequency #(H G) of subgraph H in G is monotonically non-increasing while descending the tree node. Then, the condition (12) is obviously satisfied because for a sequence of H H H · · · in the descending path of the tree, is a feature corresponding to H in G i . Based on this property, the following lemma enables us to prune a node during the tree traversal.
be a pruning criterion. Then, if the inequality holds, for any descendant node k ⊇ k, the following inequality holds: where q ∈ R 2nK + and r ≥ 0 are an arbitrary constant vector and scalar variable, respectively.
This lemma indicates that if the condition (14) is satisfied, we can say that none of the descendant nodes are included in F. Assuming that the indicator function g(x) = 1 x>0 is used in (2), a tighter bound can be obtained through the following lemma.
Lemma 4.2. If g(x) = 1 x>0 is set in (2), the pruning criterion (13) can be replaced with By comparing the first terms of Lemmas 4.1 and 4.2, we see that Lemma 4.2 is tighter when g(x) = 1 x>0 as follows:

Step1. Traverse Graph Mining Tree
Step2. Solve subproblem  A schematic illustration of the optimization algorithm for IGML is shown in Figure 2 (for further details, see Section 5). To generate a subset of features F, we first traverse the graph mining tree during which the safe screening/working set selection procedure and their pruning extensions are performed (Step1). Next, we solve the sub-problem (9) with the generated F using a standard gradient-based algorithm (Step2). Safe screening is also performed during the optimization iteration in Step2, which is referred to as dynamic screening. This further reduces the size of F.
Before moving onto detailed formulations, we summarize our rules to determine F in Table 1. The columns represent the different approaches to evaluating the necessity of features, i.e., safe and working set approaches. For the safe approaches, there are further 'single λ' (described in Section 4.1.2) and 'range of λ' (described in Section 4.1.3) approaches. The single λ approach considers safe rules for a specific λ, while the range of λ approach considers safe rules that can eliminate features for a range of λ (not just a specific value). Both the single and range approaches are based on the bounds of the region in which the optimal solution exists, for which details are given in Section 4.1.1. The rows of Table 1 indicate the variation of rules to remove one specific feature and rules to prune all features in a subtree.

Safe Screening
Safe screening (Ghaoui et al., 2010) was first proposed to identify unnecessary features in LASSO-type problems. Typically, this approach considers a bounded region of dual variables in which the optimal solution must exist. Then, we can eliminate dual inequality constraints that are never violated given that the solution exists in that region. The well-known Karush-Kuhn-Tucker (KKT) conditions show that this is equivalent to the elimination of primal variables that take value 0 at the optimal solution. In Section 4.1.1, we first derive a spherical bound for our optimal solution, and then in Section 4.1.2, a rule for safe screening is shown. Section 4.1.3 extends rules that are specifically useful for the regularization path calculation.

Sphere Bound for Optimal Solution
The following theorem provides a hyper-sphere containing the optimal dual variable α .
Theorem 4.1 (DGB). For any pair of m ≥ 0 and α ≥ 0, the optimal dual variable α must satisfy ). This bound is called the duality gap bound (DGB), and the parameters m and α used to construct the bound are referred to as the reference solution. This inequality reveals that the optimal α should be in the inside of the sphere whose center is the reference solution α and radius is 2 P λ (m) − D λ (α), i.e., twice the square root of the duality gap. Therefore, if the quality of the reference solution m and α is better, a tighter bound can be obtained. When the duality gap is zero, meaning that m and α are optimal, the radius is shrunk to zero.
If the optimal solution for λ 0 is available as a reference solution to construct the bound for λ 1 , the following bound, called regularization path bound (RPB), can be obtained.
Theorem 4.2 (RPB). Let α 0 be the optimal solution for λ 0 and α 1 be the optimal solution for λ 1 . Then, This inequality indicates that the optimal dual solution for λ 1 (α 1 ) should be in the sphere whose center is λ0+λ1 2λ0 α 0 and radius is λ0−λ1 2λ0 α 0 2 . However, RPB requires the exact solution, which is difficult to obtain in practice due to numerical errors. The relaxed RPB (RRPB) extends RPB to incorporate the approximate solution as a reference solution.
Theorem 4.3 (RRPB). Assuming that α 0 satisfies α 0 − α 0 2 ≤ , the optimal solution α 1 for λ 1 must satisfy In Theorem 4.1, the reference α 0 is only assumed to be close to α 0 within the radius instead of assuming that α 0 is available. For example, can be obtained using the DGB (Theorem 4.1). Similar bounds to those derived here were previously considered for the triplet screening of metric learning on usual numerical data (Yoshida et al., 2018(Yoshida et al., , 2019b. Here, we extend a similar idea to derive subgraph screening.

Safe Screening and Safe Pruning Rules
Theorem 4.1 and 4.3 identify the regions where the optimal solution exists using a current feasible solution α. Further, from (6), when C k,: α ≤ λ, we have m k = 0. This indicates that where B is a region containing the optimal solution α , i.e., α ∈ B. As we derived in Section 4.1.1, the sphere-shaped B can be constructed using feasible primal and dual solutions. By solving this maximization problem, we obtain the following safe screening (SS) rule.
Theorem 4.4 (SS Rule). If the optimal solution α exists in the bound B = {α | α − q 2 2 ≤ r 2 }, the following rule holds Theorem 4.4 indicates that we can eliminate unnecessary features by evaluating the condition shown in (16).
Here, the theorem is written in a general form, and in practice, q and r can be defined by the center and a radius of one of the sphere bounds, respectively. An important property of this rule is that it guarantees optimality, meaning that the sub-problems (9) and (10) have the exact same optimal solution to the original problem if F is defined through this rule. However, it is still necessary to evaluate the rule for all p features, which is currently intractable. To avoid this problem, we derive a pruning strategy on the graph mining tree, which we call the safe pruning (SP) rule.

Theorem 4.5 (SP Rule). If the optimal solution α is in the bound
This theorem is a direct consequence of Lemma 4.1. If this condition holds for a node k during the tree traversal, a subtree below that node can be pruned. This means that we can safely eliminate unnecessary subgraphs even without enumerating them. In this theorem, note that B has an additional non-negative constraint q ≥ 0, but this is satisfied by all the bounds in Section 4.1.1 because of the non-negative constraint in the dual problem.

Range-based Safe Screening & Safe Pruning
The SS and SP rules apply to a fixed λ. The range-based extension identifies an interval of λ for which the satisfaction of SS/SP is guaranteed. This is particularly useful for the path-wise optimization or regularization path calculation, where the problem must be solved with a sequence of λ. We assume that the sequence is sorted in descending order, as optimization algorithms typically start from the trivial solution m = 0. Let λ = λ 1 ≤ λ 0 . By combining RRPB with the rule (16), we obtain the following theorem.
Theorem 4.6 (Range-based Safe Screening (RSS)). For any k, the following rule holds where This rule indicates that we can safely ignore m k for λ ∈ [λ a , λ 0 ], while if λ a > λ 0 , the weight m k cannot be removed by this rule. For SP, the range-based rule can also be derived from (17).
Theorem 4.7 (Range-based Safe Pruning (RSP)). For any k ⊇ k, the following pruning rule holds: where This theorem indicates that, while λ ∈ [λ a , λ 0 ], we can safely remove the entire subtree of k. Analogously, if the feature vector is generated from g(x) = 1 x>0 (i.e., binary), the following theorem holds.
Theorem 4.8 (Range-Based Safe Pruning (RSP) for binary feature). Assuming g(x) = 1 x>0 in (2), a and b in theorem 4.7 can be replaced with Because these constants a and b are derived from the tighter bound in Lemma 4.2, the obtained range becomes wider than the range in Theorem 4.7.
Once we calculate λ a and λ a of (18) and (19) for some λ, they are stored at each node of the tree. Subsequently, such λ a and λ a can be used for the next tree traversal with different λ . If the conditions of (18) or (19) are satisfied, the node can be skipped (RSS) or pruned (RSP). Otherwise, we update λ a and λ a by using the current reference solution.

Working Set Method
Safe rules are strong rules in the sense that they can completely remove features; thus, they are sometimes too conservative to fully accelerate the optimization. In contrast, the working set selection is a widely accepted heuristic approach to selecting a subset of features.

Working Set Selection & Working Set Pruning
The working set (WS) method optimizes the problem with respect to only selected working set features. Then, if the optimality condition for the original problem is not satisfied, the working set is reselected and the optimization on the new working set restarts. This process iterates until optimality on the original problem is achieved.
Besides the safe rules, we use the following WS selection criterion, which is obtained directly from the KKT conditions: If this inequality is satisfied, the k-th dimension is predicted as m k = 0. Hence, the working set is defined by Although m i = 0 for i / ∈ W is not guaranteed, the final convergence of the procedure is guaranteed by the following theorem.
Theorem 4.9 (Convergence of WS). Assume that there is a solver for the sub-problem (9) (or equivalently (10)) that returns the optimal solution for given F. The working set method, which iterates optimizating the sub-problem with F = W and updating W alternately, returns the optimal solution of the original problem in finite steps.
However, here again, the inequality (20) needs to be evaluated for all features, which is computationally intractable.
The same pruning strategy as for SS/SP can be incorporated into working set selection. The criterion (20) is also a special case of (11), and Lemma 4.1 indicates that if the following inequality Prune WP (k) := Prune(k|α, 0) ≤ λ, holds, then no k ⊇ k is included in the working set. We refer to this criterion as working set pruning (WP).

Relation with Safe Rules
Note that for the working set method, we may need to update W multiple times, unlike in the safe screening approaches, as shown by Theorem 4.9. Instead, the working set method can usually exclude a larger number of features compared with safe screening approaches. In fact, when the condition of the SS rule (16) is satisfied, the WS criterion (20) must likewise be satisfied. Because all the spheres (DGB, RPB and RRPB) contain the reference solution α, which is usually the current solution, the inequality holds, where B is a sphere created by DGB, RPB or RRPB. This indicates that when the SS rule excludes the k-th feature, the WS also excludes the k-th feature. However, to guarantee convergence, WS needs to be fixed until the sub-problem (9)-(10) is solved (Theorem 4.9). In contrast, the SS rule is applicable anytime during the optimization procedure without affecting the final optimality. This enables us to apply the SS rule even to the sub-problem (9)-(10), where F is defined by WS as shown in Step 2 of Figure 2 (dynamic screening).
For the pruning rules, we first confirm the following two properties: is the center of the sphere, r ≥ 0 is the radius, and C ∈ R is a constant. In the case of DGB, the center of the sphere is the reference solution α itself, i.e., q = α. Then, the following relation holds between the SP criterion Prune(k|q, r) and WP criterion Prune WP (k): This once more indicates that when the SP rule is satisfied, the WP rule must be satisfied as well. When the RPB or RRPB sphere is used, the center of the sphere is q = λ0+λ1 2λ0 α 0 . Assuming that the solution for λ 0 is used as the reference solution, i.e., α = α 0 , we obtain Using this inequality, we obtain indicating that the pruning of WS is always tighter than that of the safe rule. However, in our algorithm presented in Section 5, λ 1 < λ 0 holds because we start from a larger value of λ and gradually decrease it. Then, this inequality does not hold, and Prune(k|q, r) < Prune WP (k) becomes a possibility. When the WS and WP rules are strictly tighter than the SS and SP rules, respectively, using both of WS/WP and SS/SP rules is equivalent to using WS/WP only (except for dynamic screening). Even in this case, the range-based safe approaches (the RSS and RSP rules) can still be effective. When the range-based rules are evaluated, we obtain the range of λ such that the SS or SP rule is satisfied. Thus, as long as λ is in that range, we do not need to evaluate any safe or working set rules.

Training with Path-wise Optimization
We employ path-wise optimization (Friedman et al., 2007), where the optimization starts from λ = λ max , which gradually decreases λ while optimizing m. As can be seen from (8), λ max is defined by the maximum of the inner product C k,: α. This value can also be found by a tree search with pruning. Suppose that we calculate C k,: α while traversing the tree andλ max is the current maximum value during the traversal. Using Lemma 4.1, we can derive the pruning rule Prune(k|α, 0) ≤λ max .
If this condition holds, the descendant nodes of k cannot be maximal, and thus we can identify λ max without calculating C k,: α for all candidate features.
Algorithm 1 shows the outer loop of our path-wise optimization. The TRAVERSE and SOLVE functions in Algorithm 1 are shown in Algorithm 2 and 3, respectively. Algorithm 1 first calculates λ max which is the minimum λ at which the optimal solution is m = 0 (line 3). The outer loop in lines 5-14 is the process of decreasing λ with the decreasing rate R. The TRAVERSE function in line 7 determines the subset of features F by traversing tree with SS and WS. The inner loop (line 9-13) alternately solves the optimization problem with the current F and updates F until the duality gap becomes less than the given threshold eps.
Algorithm 2 shows the TRAVERSE function, which recursively visits tree nodes to determine F. The variable node.pruning contains λ a of RSP, and if the RSP condition (19) is satisfied (line 3), the function Algorithm 1: Path-wise Optimization 1 function PathwiseOptimization(R, T, freq, MaxIter, eps) Initialize root node as root.children = empty, root.screening = ∞, and root.pruning = ∞ returns the current F (the node is pruned). The variable node.screening contains λ a of RSS, and if the RSS condition (18) is satisfied (line 5), this node can be skipped, and the function proceeds to the next node. If these two conditions are not satisfied, the function 1) updates node.pruning and node.screening if update is true, and 2) evaluates the conditions of RSP and WP (line 10), and RSS and WS (line 14). In lines 17-18, gSpan expands the children of the current node, and for each child node, the TRAVERSE function is called recursively. Algorithm 3 shows a solver for the primal problem with the subset of features F. Although we employ a simple projected gradient algorithm, any optimization algorithm can be used in this process. In lines 7-10, the SS rule is evaluated at every after freq iterations. Note that this SS is only for the sub-problems (9) and (10) created by the current F (not for the original problems).

Enumerating Subgraphs for Test Data
To obtain a feature vector for test data, we only need to enumerate subgraphs with m k = 0. When gSpan is used as a mining algorithm, a unique code, called minimum DFS code, is assigned to each node. If a DFS code for a node is (a 1 , a 2 , . . . , a n ), a child node is represented by (a 1 , a 2 , . . . , a n , a n+1 ). This enables us to prune nodes that does not generate subgraphs with m k = 0. Suppose that a subgraph (a 1 , a 2 , a 3 ) = (x, y, z) must be enumerated, and that we are currently at node (a 1 ) = (x). Then, a child with (a 1 , a 2 ) = (x, y) should be traversed, but a child with (a 1 , a 2 ) = (x, w) cannot generate (x, y, z), and consequently we can stop the traversal of this node.

Learning Mahalanobis Distance for Selected Subgraphs
Instead of m, the following Mahalanobis distance can be considered where M is a positive definite matrix. Directly optimizing M requires O(p 2 ) primal variables and semidefinite constraint, making the problem computationally expensive, even for relatively small p. Thus, as optional post-processing, we consider optimizing the Mahalanobis distance (22)  Compute α by (7) Above, R : R h×h → R is a regularization term for M , where a typical setting is R(M ) = trM + η 2 M 2 F with tr representing the trace of a matrix. This metric can be more discriminative, because it is optimized to the training data with a higher degree of freedom.

Vector Representation of a Graph
An explicit vector representation of an input graph can be obtained using optimized m as follows: Unlike the original x i , the new representation x i is computationally tractable because of the sparsity of m, and simultaneously, this space should be highly discriminative. This property is beneficial for further analysis of the graph data. We show an example of applying the decision tree to the learned space later in the paper.
In the case of the general Mahalanobis distance given in Section 5.3.1, we can obtain further transformation. Let M = V ΛV be the eigenvalue decomposition of the learned M . By employing the regularization term R(M ) = trM + η 2 M 2 F , some of the eigenvalues of M can be shrunk to 0 because trM is equal to the sum of the eigenvalues. If M has h < h non-zero eigenvalues, Λ can be written as a h × h diagonal matrix, and V is a h × h matrix such that each column is the eigenvector of a non-zero eigenvalue. Then, a representation of the graph is This can be considered as a supervised dimensionality reduction from h-to h -dimensional space. Although each dimension no longer corresponds to a subgraph in this representation, the interpretation remains clear because each dimension of the transformed vector is simply a linear combination of z i .

Extensions
In this section, we consider three extensions of IGML: applications to other data types, employing a triplet loss function, and introducing vertex-label similarity.

Application to Other Structured Data
In addition to graph data, the proposed method can be applied to itemset/sequence data . For an itemset, the Jaccard index, defined as the size of the intersection of two sets divided by the size of the union, is the most popular similarity measure. Although a few studies have considered kernels for an itemset (Zhang et al., 2007), to the best of our knowledge, it remains difficult to adapt a metric on a given labeled dataset in an interpretable manner. In contrast, there are many kernel approaches for sequence data. The spectrum kernel (Leslie et al., 2001) creates a kernel matrix by enumerating all k-length subsequences in the given sequence. The mismatch kernel  enumerates subsequences allowing m discrepancies in a pattern of length k. The gappy kernel (Leslie and Kuang, 2004;Kuksa et al., 2008) counts the number of k-mers (subsequences) with a certain number of gaps g that appear in the sequence. The above kernels require the value of hyperparameter k, although various lengths may in fact be related. The motif kernel (Zhang and Zaki, 2006;Pissis et al., 2013;Pissis, 2014) counts the number of "motifs" appearing in the input sequences, the "motif" must be decided by the user. Because these approaches are based on the idea of the 'kernel', they are unsupervised, unlike our approach. By employing a similar approach to the graph input, we can construct a feature representation φ H (X i ) for both itemset and sequence data. For the itemset data, the i-th input is a set of items X i ⊆ I, where I is a set of all items, e.g., X 1 = {a, b}, X 2 = {b, c, e}, . . . with the candidate items I = {a, b, c, d, e}. The feature φ H (X i ) is defined by 1 H⊆Xi for ∀H ⊆ I. This feature φ H (X i ) also has monotonicity φ H (X i ) ≤ φ H (X i ) for H ⊇ H. In sequence data, the i-th input X i is a sequence of items. Thus, the feature φ H (X i ) is defined from the frequency of a sub-sequence H in the given X i . For example, if we have X i = b, b, a, b, a, c, d and H = b, a , then H occurs twice in X i . For sequence data, the monotonicity property is again guaranteed where H is a sub-sequence of H . Because of these monotonicity properties, we can apply the same pruning procedures to both of itemset and sequence data. Figure 3 shows examples of trees that can be constructed by itemset and sequence mining algorithms (Agrawal et al., 1994;Pei et al., 2001).

Triplet Loss
We formulate the loss function of IGML as the pair-wise loss (3). Triplet loss functions are also widely used in metric learning (e.g., Weinberger and Saul, 2009): where T is an index set of triplets consisting of (i, j, l) satisfying y i = y j , y i = y l . This loss incurs a penalty when the distance between samples in the same class is larger than the distance between samples in different classes. Because the loss is defined by a 'triplet' of samples, this approach can be more time-consuming than the pairwise approach. In contrast, the relative evaluation such as d m (x i , x j ) < d m (x i , x l ) (the j-th sample must be closer to the i-th sample than the l-th sample) can capture the higher-order relations between input objects rather than penalizing the pair-wise distance. A pruning rule can be derived even for the case of triplet loss. By defining c ijl := c il − c ij , the loss function can be written as (i,j,l)∈T 1 (m c ijl ).
Because this has the same form as pairwise loss with L = 1 and U = 0, the optimization problem is reduced to the same form as the pairwise case. We require a slight modification of Lemma 4.1 because of the change of the constant coefficients (i.e., from c ij to c ijl ). The equation (13) is changed to This is easily proven using which is an immediate consequence of the monotonicity inequality (12).

Considering Vertex-Label Similarity
Because IGML is based on the exact matching of subgraphs to create the feature φ H (G), it is difficult to provide a prediction for a graph that does not exactly match many of the selected subgraphs. Typically, this  (26) and (27). The relation L P (u 2 , 3) L G (v 2 , 3) is satisfied between u 2 and v 2 at Step3.
happens when the test dataset has a different distribution of vertex-labels. For example, in the case of the prediction on a chemical compound group whose atomic compositions are largely different from those in the training dataset, the exact match may not be expected as in the case of the training dataset. Therefore, we consider incorporating similarity/dissimilarity information of graph vertex-labels to relax this exact matching constraint. A toy example of vertex-label dissimilarity is shown in Figure 4. In this case, the 'red' vertex is similar to the 'green' vertex, while it is dissimilar to the 'yellow' vertex. For example, we can create this type of table by using prior domain knowledge (e.g., chemical properties of atoms). Even when no prior information is available, a similarity matrix can be inferred using any embedding method (e.g., Huang et al., 2017). Because it is difficult to directly incorporate similarity information into our subgraph isomorphism-based feature φ H (G), we first introduce a relaxed evaluation of inclusion of a graph P in a given graph G. We assume that P is obtained from the gSpan tree of the training data. Our approach is based on the idea of 're-labeling' graph vertex-labels in the Weisfeiler-Lehman (WL) kernel (Shervashidze et al., 2011), which is a well-known graph kernel with an approximate graph isomorphism test. Figure 5 (a) shows an example of the re-labeling procedure, which is performed in a fixed number of recursive steps. The number of steps is denoted as T (T = 3 in the figure) and is assumed to be pre-specified. n ]), which originate from the two vertices u and v in graphs P and G, respectively. We say that L P (v, h) is included in L G (u, h) and denote it by when the following recursive condition is satisfied: , which is assumed to be already obtained on the level h − 1 computation. Because S When the maximum number of matchings is equal to m, this means that there exists an injection σ(i) σ(i) . It is well known that the maximum bipartite matching can be reduced to the maximum flow problem, which can be solved in the polynomial time (Goldberg and Tarjan, 1988). An example of the inclusion relationship is shown in Figure 5 (b).
Let |P | and |G| be the numbers of vertices in P and G, respectively. Then, multisets of the level T hierarchical labels of all the vertices in P and G are written as [L P (u i , T )] i∈[|P |] := [L P (u 1 , T ), L P (u 2 , T ), . . . , L P (u |P | , T )] and [L G (v i , T )] i∈[|G|] := [L G (v 1 , T ), L G (v 2 , T ), . . . , L G (v |G| , T )], respectively. For a feature of a given input graph G, we define the approximate subgraph isomorphism feature (ASIF) as follows: This feature approximately evaluates the existence of a subgraph P in G using the level T hierarchical labels. ASIF satisfies the monotone decreasing property (12), i.e., x P G ≤ x P G if P P , because the number of conditions in (27) only increases when P grows.
To incorporate label dissimilarity information (as shown in Figure 4) into ASIF, we first extend the label inclusion relation (26) by using the concept of optimal transportation cost. As a label similarity-based relaxed evaluation of L P (v, h) L G (u, h), we define an asymmetric cost between L P (u, h) and L G (v, h) as follows where the second term of (29b) is which we refer to as the label transportation cost (LTC) representing the optimal transportation from the multiset S (h) to another multiset S (h) among the set of all injections I := {∀σ : [m] → [n] | σ(i) = σ(j) for i = j}. The equation (29) has a recursive structure similar to that of (26). The first case (29a) occurs when S (h) = S (h) = [ ], which is at the first level h = 1. In this case, cost 1 is defined by dissimilarity(F (1) , F (1) ), which is directly obtained as a dissimilarity between original labels since F (1) and F (1) stem from the original vertex-labels. In the other case (29b), the cost is recursively defined as the sum of the cost from F (h) to F (h) and the optimal-transport cost from S (h) to S (h) . Although this definition is recursive, as in the case of ASIF, the evaluation can be performed by computing sequentially from h = 1 to h = T . Because , which can be solved in polynomial time (Goldberg and Tarjan, 1988).
We define an asymmetric transport cost for two graphs P and G, which we call the graph transportation cost (GTC), as LTC from all level T hierarchical labels of P to those of G: Then, as a feature of the input graph G, we define the following sim-ASIF : where ρ > 0 is a hyperparameter. This sim-ASIF can be regarded as a generalization of (28) based on the vertex-label similarity. When dissimilarity(F (1) , F (1) ) := ∞ × 1 F (1) =F (1) , the feature (31) is equivalent to (28). Similarly to ASIF, GTC(P → G) satisfies the monotonicity property GTC(P → G) ≤ GTC(P → G) for P P because the number of vertices to transport increases as P grows. Therefore, sim-ASIF (31) satisfies the monotonicity property, i.e., x P →G ≤ x P →G if P P . From the definition (31), sim-ASIF always has a positive value x P →G > 0 except when GTC(P → G) = ∞, which may not be suitable for identifying a small number of important subgraphs. Further, in sim-ASIF, the bipartite graph in the minimum-cost-flow calculation LTC(S → S , cost h−1 ) is always a complete bipartite graph, where all vertices in S are connected to all vertices in S . Because the efficiency of most of standard minimum-cost-flow algorithms depends on the number of edges, this may entail a large computational cost. As an extension to mitigate these issues, a threshold can be introduced into sim-ASIF as follows: where t > 0 is a threshold parameter. In this definition, x = 0 when exp{−ρ GTC(P → G)} ≤ t, i.e., GTC(P → G) ≥ −(log t)/ρ. This indicates that if a cost is larger than −(log t)/ρ, we can regard the cost as ∞. Therefore, at any h, if the cost between S is not necessary. Then, the number of matching pairs can be less than m in LTC(·) because of the lack of edges, and in this case, the cost is regarded as ∞. Furthermore, if cost h (F (h) → F (h) ) is larger than −(log t)/ρ in (29b), the computation of LTC(S (h) → S (h) , cost h−1 ) is not required because x = 0 is determined.
Note that transportation-based graph metrics have been studied (e.g., Titouan et al., 2019), but the purpose of such studies was to evaluate the similarity between two graphs (not inclusion). Our (sim-)ASIF provides a feature with the monotonicity property as a natural relaxation of subgraph isomorphism, by which the optimality of our pruning strategy can be guaranteed. In contrast, there have been many studies on inexact graph matching  such as eigenvector- (Leordeanu et al., 2012;Kang et al., 2013), edit distance- (Gao et al., 2010), and random walk-based (Gori et al., 2005;Cho et al., 2010) methods. Some of these methods provide a score for the matching, which can be seen as a similarity score between a searched graph pattern and a matched graph. However, they do not guarantee the monotonicity of the similarity score for pattern growth. If the similarity score satisfies monotonicity, it can be combined with IGML. Although we only consider vertex-labels, edge-labels can also be incorporated into (sim-)ASIF. A simple approach is to transform a labeled-edge into a labeled-node with two unlabeled edges, such that (sim-)ASIF is directly applicable.

Experiments
We evaluate the performance of IGML using the benchmark datasets shown in Table 3. These datasets are available from Kersting et al. (2016). We did not use edge labels because the implementations of compared methods cannot deal with them, and the maximum connected graph is used if the graph is not connected. Note that IGML currently cannot directly deal with continuous attributes, so we did not use them. A possible approach would be to perform discretization or quantization before the optimization, such as taking grid points or applying clustering in the attribute space. Building a more elaborated approach, such as dynamically determining discretization, is a possible future directions. The #maxvertices column in the table indicates the size (number of vertices) of the maximum subgraph considered in IGML. To fully identify important subgraphs, a large value of #maxvertices is preferred, but this can cause a correspondingly large memory requirement to store the gSpan tree. For each dataset, we set the largest value for which IGML could finish with a tractable amount memory. The sets S i and D i were selected as the ten nearest neighborhoods of x i (K = |S i | = |D i | = 10) by using the WL-Kernel. A sequence of the regularization coefficients was created by equally spacing 100 grid points on a logarithmic scale between λ max and 0.01λ max . We set the minimum support in gSpan as 0, meaning that all the subgraphs in a given dataset were enumerated for as far as the graph satisfies the #maxvertices constraint. The gSpan tree is mainly traversed when the beginning of each λ as shown in Algorithm 1 (in the case of WS-based approaches, the tree is also traversed at every working set update). Note that the tree is dynamically constructed during this traversal without constructing the entire tree beforehand. In the working-set method, after convergence, it is necessary to traverse the tree again in order to confirm the overall optimality. If the termination condition is not satisfied, optimization with a new working set must be performed. The termination condition for the optimization is that the relative duality gap is less than 10 −6 . In the experiment, we used g(x) = 1 x>0 in φ H (G) with Lemma 4.2 unless otherwise noted. The dataset was randomly divided in such a way that the ratio of partitioning was train : validation : test = 0.6 : 0.2 : 0.2, and our experimental results were averaged over 10 runs.

Evaluating Computational Efficiency
In this section, we confirm the effect of the proposed pruning methods. We evaluated four settings: Safe Screening and Pruning: "SSP", Range based Safe Screening and Pruning: "RSSP", Working set Selection and Pruning: "WSP", and the combination of WSP and RSSP: "WSP+RSSP". Each method performed dynamic screening with DGB at every update of m. We here used the AIDS dataset, where #maxvertices=30. In this dataset, when we fully traversed the gSpan tree without safe screening/working set selection, the number of tree nodes was more than 9.126 × 10 7 , at which point our implementation with gSpan stopped because we ran out of memory. Figure 6 (a) shows the size of F after the first traverse at each λ, and the number of non-zero m k after the optimization is also shown as a baseline. We first observe that both approaches significantly reduced the number of features. Even for the largest case, where approximately 200 of features were finally selected by m k , only less than 1000 features remained. We observe that WSP exhibited significantly smaller values than SSP. Instead, WSP may need to perform the entire tree search again because it cannot guarantee the sufficiency of the current F, while SSP does not need to search the tree again because it guarantees that F must contain all m k = 0.
The number of visited nodes in the first traversal at each λ is shown in Figure 6 (b). Here, we added RSSP and WSP+RSSP, which are not shown in Figure 6  for SSP and RSSP, and for WSP and WSP+RSSP. Because RSSP is derived from SSP, it does not change the number of screened features. As we discussed in Section 4.2.2, WSP removes more features than RSSP, though it is not safe. We observed that the #visited nodes of SSP was the largest, but it was less than approximately 27000 (27000/9.126 × 10 7 ≈ 0.0003). Comparing SSP and WSP, we see that WSP pruned a larger number of nodes. In contrast, the #visited nodes of RSSP was less than 6000. The difference between SSP and RSSP indicates that a larger number of nodes can be skipped by the range-based method. Therefore, by combining the node skip by RSSP with the stronger pruning of WSP, the #visited nodes was further reduced. RSSP and WSP+RSSP had larger values at λ 0 than the subsequent λ i . This is because of the effect of range-based screening and pruning. At λ 0 , every visited node in the tree calculates the ranges in which the screening and pruning rules are satisfied (i.e., RSS and RSP rules), and as a result, some nodes can be skipped during that λ i is in those ranges. At every λ i for i > 0, the ranges are updated only in the (non-skipped) visited nodes, and thus, the range-based rules take the effect except for λ 0 .
The total time of the path-wise optimization is shown in Table 2. RSSP and WSP+RSSP were fast with regard to the traversing time, and WSP and WSP+RSSP were fast with regard to the solving time. Note that because the tree is dynamically constructed during the traverse, the 'Traverse' time includes the time spent on the tree construction. In total, WSP+RSSP was the fastest. These results indicate that our method only took approximately 1 minute to solve the optimization problem with more than 9.126 × 10 7 variables. We also show the computational cost evaluation for other datasets in the Appendix I.
Although we have confirmed that IGML works efficiently on several benchmark datasets, completely elucidating general complexity of IGML is remains as future work. The practical complexity at least depends on the graph size in the training data, #maxvertices, #samples, and the pruning rate. In terms of the graph size, traversing a large graph dataset using gSpan can be intractable because it requires all matched subgraphs to be maintained at each tree node. Therefore, applying IGML to large graphs, e.g., graphs with more than thousands of nodes, would be difficult. Meanwhile, the scalability of IGML depends not only on the sizes of graphs but also strongly on the performance of the pruning. However, we still do not have any general analytic complexity evaluation for the rate of the pruning that avoids exponential worst-case computations. In fact, we observed that there exist datasets in which efficiency of the pruning is not sufficient. For example, on the IMDB-BINARY and IMDB-MULTI datasets, which are also from , a large number of small subgraphs are shared across all the different classes and instances (i.e., x i,k = 1 for ∀i). Our upper bound in the pruning is based on the fact that x i,k ≤ x i,k for descendant node k in the mining tree. This bound becomes tighter when x i,k = 0 for many i because 0 is the lower bound of x i,k . In contrast, when many instances have x i,k = 1, the bound can be loose, making traversal intractable. This is an important open problems common in predictive mining methods (Nakagawa et al., 2016;Morvan and Vert, 2018).

Predictive Accuracy Comparison
In this section, we compare the prediction accuracy of IGML with those of the Graphlet-Kernel (GK) , Shortest-Path Kernel (SPK) (Borgwardt and Kriegel, 2005), Random-Walk Kernel (RW) (Vishwanathan et al., 2010), Weisfeiler-Lehman Kernel (WL) (Shervashidze et al., 2011), and Deep Graph Convolutional Neural Network (DGCNN) (Zhang et al., 2018a). We used the implementations available at the URLs in the footnote 1 . Note that we mainly compared methods for obtaining a metric between graphs. The graph kernel approach is one of most important existing approaches to defining a metric space of non-vector structured data. Although kernel functions are constructed in an un-supervised manner, their high prediction performance has been widely shown. In particular, the WL kernel is known for its comparable classification performance to recent graph neural networks (e.g., Niepert et al., 2016;Morris et al., 2019). Meanwhile, DGCNN can provide a vector representation of an input graph by using the outputs of some middle layer, which can be interpreted that a metric space is obtained through a supervised learning. We did not compare with (Saigo et al., 2009;Nakagawa et al., 2016;Morvan and Vert, 2018) as they only focused on specific linear prediction models rather than building a general discriminative space. We employed the k-nearest neighbors (k-nn) classifier to directly evaluate the discriminative ability of feature spaces constructed by IGML and each kernel function. We here employed the k-nn classifier for directly evaluating discriminative ability of feature spaces constructed by IGML and each kernel function. A graph kernel can be seen as an inner-product k(G i , G j ) = ϕ(G j ), ϕ(G j ) , where ϕ is a projection from a graph to reproducing kernel Hilbert space. Then, the distance can be written as The values of k for the k-nn were k = 1, 3, 5, 7, ..., 49 and hyperparameters of each method were selected using the validation data, and the prediction accuracy was evaluated on the test data. The graphlet size for GK was set up to 6. The parameter λ RW for RW was set to the recommended λ RW = max i∈Z:10 i <1/d 2 10 i , where d denotes the maximum degree. The loop parameter h of WL was selected from 0, 1, 2, ..., 10 by using the validation data. For DGCNN, the number of hidden units and their sort-pooling were also selected using the validation data, each ranging from 64, 128, 256 and from 40%, 60%, 80%, respectively. The micro-F1 score for each dataset is shown in Table 3. "IGML (Diag)" indicates IGML with the weighted squared distance (1), and "IGML (Diag→Full)" indicates that with post-processing using the Mahalanobis distance (22). "IGML (Diag)" yielded the best or comparable to the best score on seven out of nine datasets. This result is impressive because IGML uses a much simpler metric than the other methods. Among the seven datasets, "IGML (Diag→Full)" slightly improved the mean accuracy on four datasets, but the difference was not significant. This may suggest that the diagonal weighting can have enough performance in many practical settings. WL kernel also exhibited superior performance, showing the best or comparable to the best accuracy on six datasets. DGCNN showed high accuracy with on the DBLP v1 dataset, which has a large number of samples, while its accuracy was low for the other datasets. Figure 7 shows an illustrative example of IGML on the Mutagenicity dataset, where mutagenicity was predicted from a graph representation of molecules. Figure 7 (a) is a graphical representation of subgraphs, each of which has a weight shown in (b). For example, we can clearly see that subgraph #2 is estimated as an important sub-structure to discriminate different classes. Figure 7 (c) shows a heatmap of the transformation Table 3: Comparison of micro-F1 score. OOM means out of memory. ">1week" indicates that the algorithm ran for more than a week. The numbers after the "±" are the standard deviations. Every dataset had two classes. The bold font indicates the best average value and the '*' symbol indicates that it is comparable to the best value in terms of one-sided t-test (significance level 0.05). matrix √ ΛV optimized for the thirteen features, containing three non-zero eigenvalues. For example, we see that the subgraphs of #10 and #12 have similar columns in the heatmap. This indicates that these two similar subgraphs (#10 contains #12) are shrunk to almost same representation by the regularization term R(M ).

Illustrative Examples of Selected Subgraphs
As another example of graph data analysis on the learned representation, we applied the decision tree algorithm to the obtained feature (23) on the Mutagenicity dataset. Although there has been a study constructing a decision tree directly for graph data (Nguyen et al., 2006), it requires a severe restriction on the patterns to be considered for computational feasibility. In contrast, because (23) is a simple vector representation with a reasonable dimension, it is quite easy to apply the decision tree algorithm. We selected two paths from the obtained decision tree as shown in Figure 8. For example, in the path (a), if a given graph contains "O = N ", and does not contain "H−O−C−C = C−C−H", and contains "N−C = C−C = C< C C ", the given graph is predicted as y = 0 with probability 140/146. Both rules clearly separate the two classes, which is highly insightful as we can trace the process of the decision based on the subgraphs.

Experiments for Three Extensions
In this section, we evaluate the performance of the three extensions of IGML described in Section 6.
First, we evaluated the performance of IGML on itemset and sequence data using the benchmark datasets shown in the first two rows of Table 4 and 5. These datasets can be obtained from (Dua and Graff, 2017) and (Chang and Lin, 2011), respectively. We set the maximum-pattern size considered by IGML as 30. Table 4 lists the micro-F1 scores on the itemset datasets. We used k-nn with the Jaccard similarity as a baseline, where k was selected using the validation set, as described in Section 7.2. The scores of both of IGML (Diag) and (Diag→Full) were superior to those of the Jaccard similarity on all datasets. Table 5 lists the micro-F1 scores on the sequence dataset. Although IGML (Diag) did not outperform the mismatch kernel  for the promoters dataset, IGML (Diag→Full) achieved a higher F1-score than the kernel on all datasets. Figure 9 shows an illustrative example of sequences identified by IGML on the promoters dataset, where the task was to predict whether an input DNA sequence stems from a promoter region. Figure 9 (a) is a graphical representation of the sequence, and the corresponding weights are shown in (b). For example, the sub-sequence #1 in (a) can be considered as an important sub-sequence to discriminate different classes.
Second, we show the results of the triplet formulation described in Section 6.2. To create the triplet set T , we followed the approach in Shen et al. (2014), where k neighborhoods in the same class x j and k neighborhoods in different classes x l were sampled for each x i (k = 4). Here, IGML with the pairwise loss is referred to as 'IGML (Pairwise)', and IGML with the triplet loss is referred to as 'IGML (Triplet)'. Table 6 compares the micro-F1 scores of IGML (Pairwise) and IGML (Triplet). IGML (Triplet) showed higher F1-scores than IGML (Pairwise) on three of nine datasets, but it was not computable on the two       datasets due to running out of memory (OOM). This is because the pruning rule in the triplet case (25) was looser than in the pair-wise case. Finally, we evaluated the sim-ASIF (32). We set the scaling factor of the exponential function as ρ = 1, the threshold of the feature as t = 0.7, and the number of re-labeling steps as T = 3. We employed a simple heuristic approach to create a dissimilarity matrix among vertex-labels using labeled graphs in the given dataset. Suppose that the set of possible vertex-labels is L, and f ( , ) is the frequency that ∈ L and ∈ L are adjacent in all graphs of the dataset. By concatenating f ( , ) for all ∈ L, we obtained a vector representation of a label . We normalized this vector representation such that the vector had the unit L2 norm. By calculating the Euclidean distance of this normalized representations, we obtained the dissimilarity matrix of vertex-labels. We are particularly interested in the case where the distribution of the vertex-label frequency is largely different between the training and test datasets, because in this case the exact matching of IGML may not be suitable to provide a prediction. We synthetically emulated this setting by splitting the training and test datasets using a clustering algorithm. Each input graph was transformed into a vector created by the frequencies of each vertex-label ∈ L contained in that graph. Subsequently, we applied the k-means clustering to split the dataset into two clusters, for which C 1 and C 2 denote sets of assigned data points, respectively. We used C 1 for the training and validation datasets and C 2 is used as the test dataset, where |C 1 | ≥ |C 2 |. Following the same partitioning policy as in the above experiments, the size of the validation data was set as the same size of C 2 , resulting from which the size of the training set was |C 1 | − |C 2 |. Table 7 lists the comparison of the micro-F1 scores on the AIDS, Mutagenicity, and NCI1 datasets. We did not consider other datasets as their training set sizes created from the above procedure were too small. We fixed the #maxvertices of sim-ASIF to 8, which was less than the value in our original IGML evaluation Table 3, because sim-ASIF takes more time than the feature without vertex-label similarity. For the original IGML, we show the result for the setting in Table 3 and the results with #maxvertices 8. IGML with sim-ASIF was superior to the original IGML for the both #maxvertices settings on the AIDS and NCI1 datasets, although it has smaller #maxvertices settings, as shown in Table 7. On the Mutagenicity dataset, sim-ASIF was inferior to the original IGML reported in Table 3, but in the comparison under the same #maxvertices value, their performances were comparable. These results suggest that when the exact matching of the subgraph is not appropriate, sim-ASIF can improve the prediction performance of IGML.

Performance on Frequency Feature
In this section, we evaluate IGML with g(x) = log(1 + x) instead of g(x) = 1 x>0 . Note that because computing the frequency without overlapping #(H G) is NP-complete (Schreiber and Schwöbbermeyer, 2005), in addition to the exact count, we evaluated the feature defined by an upper bound of #(H G) (see Appendix J for details). We employed log because the scale of the frequency x is highly diversified. Based on the results in Section 7.1, we used WSP+RSSP in this section. The #maxvertices for each dataset followed those in Table 3.
The comparison of micro-F1 scores for the exact #(H G) and approximation of #(H G) is shown in Table 8. The exact #(H G) did not complete five datasets mainly due to the computational difficulty  of the frequency counting. In contrast, the approximate #(H G) completed on all datasets. Overall, for both the exact and approximate frequency features, the micro-F1 scores were comparable with the case of g(x) = 1 x>0 shown in Table 3. Table 9 lists the total times for the path-wise optimization for the exact #(H G) and the approximation of #(H G). On the AIDS dataset, the exact #(H G) did not complete within a day, while the traversal time using approximate #(H G) was only 8.6 sec. On the BZR dataset, the traversal time using the exact #(H G) was seven times that using the approximate #(H G). The solving time for the approximation was lower because |F| after traversing of the approximation was significantly less than that of the exact #(H G) in this case. Because the approximate #(H G) is an upper bound of the exact #(H G), the variation of the values of the exact #(H G) was smaller than the approximate #(H G). This resulted in higher correlations among features created by the exact #(H G). It is known that the elastic-net regularization tends to select correlated features simultaneously (Zou and Hastie, 2005), and therefore, |F| in the case of the exact #(H G) becomes larger than in the approximate case. Figure 10 shows the number of visited nodes, size of the feature subset |F| after traversal, and the number of selected features on the AIDS dataset with the approximate #(H G). This indicates that IGML keeps the number of subgraphs tractable even if g(x) = log(1 + x) is used as the feature. The #visited nodes was less than 3500, and |F| after traversal was sufficiently close to |{k |m k > 0}|. We see that #visited nodes at λ 0 is larger than many subsequent λ i s, and this is the effect of range-based rules, as shown in the case of Figure 6

Conclusions
In this paper, we proposed an interpretable metric learning method for graph data, named interpretable graph metric learning (IGML). To avoid computational difficulty, we built an optimization algorithm that combines safe screening, working set selection, and their pruning extensions. We also discussed the three extensions of IGML: (a) applications to other structured data, (b) triplet loss-based formulation, and (c) incorporating vertex-label similarity into the feature. We empirically evaluated the performance of IGML compared with existing graph classification methods. Although IGML was the only method with clear interpretability, it showed superior or comparable prediction performance compared to other state-of-the-art methods. The practicality of IGML was further demonstrated through some illustrative examples of identified subgraphs. Although IGML optimized the metric within tractable time in the experiments, the subgraphs were restricted to moderate sizes (up to 30), and a current major bottleneck for extracting larger-sized subgraphs is the memory requirement of the gSpan tree. Therefore, mitigating this memory consumption is an important future directions to apply IGML to a wider class of problems.

A Dual Problem
The primal problem (4) By the definition of the dual function in (33), to minimize L with respect to m, by partially differentiating L, we obtain The convex conjugate function of t is * which can be written as * t (x * ) = 1 4 x 2 * + tx * , (x * ≤ 0).

B Proof of Lemma 4.1
From (12), the value of (x i,k − x j,k ) 2 is bounded as follows: Using this inequality, the inner product C k ,: q is likewise bounded: Similarly, the norm C k ,: 2 is bounded: [ l∈Di max{x i,k , x l,k } 4 + j∈Si max{x i,k , x j,k } 4 ].

C Proof of Lemma 4.2
First, we consider the first term of C k ,: q + r C k ,: 2 : Now, x i,k ∈ {0, 1} is assumed. Then, if x i,k = 0, we obtain diff = l∈Di q il x l,k − j∈Si q ij x j,k ≤ l∈Di q il x l,k .
Meanwhile, if x i,k = 1, we have x i,k = 1 from the monotonicity, and subsequently diff = By using "max", we can unify these two upper bounds into Employing a similar concept, the norm of C k ,: can also be bounded by Thus, we obtain Proof. Because f is γ-strong convex from the assumption, the following inequality holds: At step t, the problem can be written using only the active constraint at the optimal solution x t as follows: From the definition of W t , the working set W t+1 must contain all active constraints {j ∈ W t | h j (x t ) = 0} at the step t and can contain other constraints that are not included in W t . This means that x t+1 must be in the feasible region of the optimization problem at the step t (46): Therefore, from the optimality condition of the optimization problem (46), From the inequalities (45) and (47), we obtain If x t is not optimal, there exists at least one violated constraint h j (x t ) > 0 for some j because otherwise x t is optimal. Then, we see x t+1 = x t because x t+1 should satisfy the constraint h j (x t+1 ) ≤ 0. If x t = x t+1 , from x t+1 − x t 2 > 0, Thus, the objective function always strictly increases (f (x t ) < f (x t+1 )). This indicates that the algorithm never encounters the same working set W t as the set of other iterations t = t. For any step t, the optimal value f (x t ) with a subset of the original constraints W t must be smaller than or equal to the optimal value f (x ) of original problem (44) with all constraints. Therefore, f (x t ) ≤ f (x ) is satisfied, and we obtain f (x T ) = f (x ) at some finite step T . Table 10 lists the computational times on the BZR, DD, and FRANKENSTEIN datasets. We first note that RSSP was approximately 2-4 times faster in terms of the traversal time compared with SSP. Next, comparing RSSP and WSP, we see that RSSP was faster for Traverse, and WSP was faster for Solve, as we observed in Table 2. Thus, the combination of WS&SP and RSSP was the fastest for all three datasets in total.

J Approximating Frequency Without Overlap
Let F G (H) be "frequency without overlap" that is the frequency of a subgraph of a given graph where any shared vertices and edges are disallowed for counting. This F G (H) is non-increasing with respect to the growth of H, but computing it is computationally complicated. Assuming that we know where all the subgraphs H appear in graph G, calculating F G (H) is equivalent to the problem of finding the maximum independent set, which is NP-complete (Schreiber and Schwöbbermeyer, 2005). In this section, using information obtained in the process of generating the gSpan tree, we approximate the frequency without overlap by its upper bound. This upper bound is also a lower bound of the frequency with overlap. Figure 11 shows the process of generating the gSpan tree and frequency. In the figure, we consider the frequency of the subgraph H ( A ○-A ○-B ○) contained in the graph G. The graph H is obtained as a pattern extension of graph A ○-A ○ (green frame) by A ○-B ○ (red frame). gSpan stores the number of these pattern extensions at each traverse node. We define the count by this extension as F max G (H) (e.g., F max G (H) = 5 for A ○-A ○-B ○). Note that F max G (H) is the frequency of H allowing overlap and duplicately counting matches that are equivalent except for the index of nodes (e.g., F max G (H) for A ○-A ○-B ○-A ○-A ○ is two in the figure). Suppose that H currently has e(> 1) edges (for example, e = 2 in A ○-A ○-B ○). We recursively go back the traverse tree (a tree in the right of Figure 11) until we reach e = 1, i.e., the starting edge that generates H (in the case of A ○-A ○-B ○, the starting edge is A ○-A ○). We use the number of unique matches of this starting edges (the number of green frames), which we define as F approx G (H), as an approximation of F G (H). Obviously, F approx