Graph Space Embedding

1 Introduction

Learning from interconnected systems can be a particularly difficult task due to the possibly non-linear interaction between the components [Linde et al., 2015; Bereau et al., 2018]. In some cases, these interactions are known and therefore constitute an important source of prior information [Jonschkowski, 2015; Zhou et al., 2018]. Although prior knowledge can be leveraged in a variety of ways [Yu et al., 2010], most of the research involving interactions, is focused on their discovery. One popular approach to deal with feature interactions, is to cast the interaction network as a graph and then use kernel methods based on graph properties, such as walk-lengths or subgraphs [Borgwardt and Kriegel, 2005; Shervashidze et al., 2009; Kriege and Mutzel, 2012] or, more recently, graph deep convolutional methods [Defferrard et al., 2016; Fout et al., 2017; Kipf and Welling, 2017]. In this work however, we focus on the case in which the interactions are feature specific and a universal property of the data instances, which make the pattern search algorithms not suitable for this task. To our knowledge, there is limited research involving this setting, although we suggest many problems can be formulated in the same way (see Figure 1). To address this knowledge gap, we present a novel method: Graph Space Embedding (GSE), an approach related to the 'randomwalk' graph kernel [G¨artner et al., 2003; Kang et al., 2012] with an important difference: it is not limited to the sum of all walks of a given length, but rather compares similar edges in two different graphs, which results in better expressiveness. Our empirical evaluation demonstrates that GSE leads to an improvement in performance compared to other baseline algorithms when plasma protein measurements and their interactions are used to predict ischaemia in patients with Coronary Artery Disease (CAD) [van Nunen et al., 2015; Zimmermann et al., 2015]. Moreover, the kernel can be computed in O(n2 ), where n is the number of features, and its hyperparameters efficiently optimized via maximization of the kernel matrix variation.

1.1 Main Contributions

1. Graph Space Embedding function that efficiently maps input into an "interaction-based" space
2. Novel theoretical result on optimal regime for the GSE, namely feasibility region for its parameters
3. Even Decent Sampling Algorithm: a strategy to gain insight on which interactions are responsible for the certain prediction

2 Approach

A remark on notation: we will use bold capital letters for matrices, bold letters for arrays and lower case letters for scalars/functions/1-d variables (ex. X, x, x).

2.1 Interaction Graphs

Any network can be represented by a graph G = {V,E}, where E is a set of edges, V a set of vertices. Denote by A|V |×|V | (|V | is equal to the number of features N) the adjacency matrix, where Ai,j represents the interaction between feature i and j, and whose value is 0 if there is no interaction.

Let x1×N be an array with measurements of features 1 to N for a given point in the data. In order to construct an instancespecific matrix, one can weigh the interaction between each pair of features with a function of their values' product:

where ϕ(A) is some function of the network interaction matrix A, and the operator ◦ represents the Hadamard product, i.e. (A◦B)i,j = (A)i,j(B)i,j.

2.2 Graph Kernel

Unlike the distance in euclidean geometry, which intuitively represents the length of a line between two points, there is no such tangible metric for graphs. Instead, one has to decide what is a reasonable evaluation for the difference between two graphs in the context of the problem.

A popular approach [G¨artner et al., 2003] is to compare random walks on both graphs. The i, jth entry of the order k power of an adjacency matrix A|V |×|V | : Ak = AA...A ︸ ︷︷ ︸ k times ,

corresponds to the number of walks of length k from i to j. Any function that maps the data into a feature space H: φ : X → H, k(x, y) =< φ(x), φ(y) > is a kernel function. Using the original graph kernel formulation, it is possible to define a kernel that will implicitly map the data into a space where the interactions are incorporated:

where G and G′ correspond to Gx(A) and Gx′(A) (see eq. 1); γi,j is a function that "controls" the mapping φ(·); and n is the maximum allowed "random walks" length. If γ is decomposed into UΛUT , where U is a matrix whose columns are the eigenvectors of γ, and Λ a diagonal matrix with its eigenvalues at each diagonal entry, then equation 2 can be re-factored into:

where φi,k,l(G) = ∑nj=1 [√ΛUT ]i,jGj . Consequently, different forms of the function γ can be chosen, with different interpretations. For the case where γi,j = θi θj , which yields:

the kernel entry can be interpreted as an inner product in a space where there is a feature for every node pair {k, l}, which represents the weighted sum of paths of length 1 to n from k to l (φk,l = ∑ni=1 θiGik,l ) [Tsivtsivadze et al., 2011]. The kernel can then be used with a method that employs the kernel trick, such as support vector machines, kernel PCA or kernel clustering. Another interesting case is when we consider the weighted sum of paths of length 1 to ∞. This can be calculated using:

since eβG = limn→+∞ ∑ni=0 βi i! Gi , where β is a parameter.

2.3 Graph Space Embedding

Since we are dealing with a universal interaction matrix for every data point and the interactions are feature specific, it makes sense to compare the same set of edges for every pair of points. As a consequence, we can also avoid solving timeconsuming graph structure problems. With these two points in mind, we combined the previous graph kernel methods and the radial basis function (RBF) to develop a new kernel which we will henceforth refer to as Graph Space Embedding (GSE). The radial basis function is defined as:

where c = e- ||x||2 σ2 e- ||y||2 σ2 . The GSE uses the distance 〈√γ[G], √γ[G′ ]〉F in the radial basis function:

If we then take the upper term of the fraction in r w to be [2 ∑|E| i=0 γ GiG′i ]n , we can use the multinomial theorem to expand each term of the exponential power series, and the expression for the kernel then becomes:

where Γ is the gamma function, Gi ∈ E is the value of edge i in G and ν = σ2 γ . Here, αn (·) represents a combination of |E| integers: (α1, α2, ..., α|E|), with ∑|E| i αni (·) = n, and the sum in r e is taken over all possible combinations of αn (·). For instance, for n = 3 in a graph with |E| = 5, possible examples of α3 (·) include (0, 1, 1, 1, 0) or (0, 2, 1, 0, 0) (see Figure 2). We begin by noting that since the sum in r e is taken over all combinations (l, k) ∈ V × V of size n, the GSE then represents a mapping from the input space to a space where all combinations of n = 0 → ∞ edges are compared between G and G′ , walks or otherwise (see fig 2). Notice that this is in contrast with the kernel of equation 5, where the comparison is between a sum of all possible walks of length n = 0 → ∞ from one node to another in the two graphs.

Figure 2: The GSE kernel implicitly compares all edge combinations between G and G′ . In this hypothetical graph, we show a sample of four α combinations for n = 3. We denote by r e(α(i)) the value inside the sum r e (see eq. 8) corresponding to the combination α(i). Note that while α(1) is a graph walk and α(2) is not, r e(α(1)) = r e(α(2)). However, due to the repetitions in α(3) and α(4), their value is shrunk in relation to the others. The higher the number of repetitions, the more the value shrinks.

The GSE also allows repeated edges. However, if the data is normalized so that µ(Gi) ≃ 0, σ(Gi) ≃ 1, then both the power in the numerator and the denominator of r e will effectively dampen most combinations with repeated edges, with a higher dampening factor for higher number of repetitions and/or combinations. Even for outlier values, the gamma function will quickly dominate the numerator of r e. The λ factor serves the purpose of shrinking the combinations with higher number of edges for ν > 2. Finally, σ2 now serves a dual purpose: the usual one in RBF to control the influence of points in relation to their distance (see equation 6), while at the same time controlling how much combinations of increasing order are penalized.

2.4 ν Feasibility Region

As discussed in the above section, the hyperparameter ν controls the shrinking of the contribution of higher order edge combinations. Intuitively, not all values of ν will yield a proper kernel matrix since too large of a value will leave out too many edge combinations while one too small will saturate the kernel values. This motivates the search for a ν value feasible operation region, where the kernel incorporates the necessary information for separability. Informally speaking, the kernel entry k(G,G′ ) measures the similarity of G and G′ . In case too few/many edge combinations are considered, the variation of the kernel values will be equal to 1. Therefore, we use the variation of the kernel matrix σ2 (K) as a proxy to detect if ν is within acceptable bounds. We shall refer to the ability of the kernel to map the points in the data into separable images φ(x) as kernel expressiveness.

To determine this region analytically, we find the νmax that yields the largest kernel variation, and then use the loss function around this value to determine in which direction the value ν should take for minimal loss.

Lemma 2.1. maxν σ2 (K(ν)) can be numerically estimated and is guaranteed to converge with a learning rate α ≤

D 2(D-1)dmax , where D is the total number of inter graph combinations and dmax is the largest combination distance.

Proof. The analytical expression for the variance is:

where we used the binomial theorem to expand b, and d = ||G - G′ ||2 . To guarantee the convergence of numerical methods the function derivative must be Lipschitz continuous:

by overloading the notation: K′ (ν) = ∂σ2 (K(ν)) ∂ν to simplify the expression. The left side of equation 10 becomes:

Since 0 ≤ e-β ≤ 1 : ∀β ∈ R, then:

When ϵ = ν - ν ′ → 0 :

and δ tends much faster to 0 then ϵ, since the denominator of δ is the exponential of the sum of ν and ν ′ . Thus, the function k′ (ν) is Lipschitz continuous with constant equal to: L(K′ (ν)) = 2 ( D-1 D ) dmax .

We shall later demonstrate empirically that ν ∗ = maxν σ2 (K(ν)) improves the class separability for our dataset.

2.5 Comparison with Standard Graph Kernels

The original formulation of the graph kernel by Gartner et. al (see eq. 2), multiplies sums of random walks of length i from one edge to another (k → l) by sums of random walks k → l from the other graph being compared of a length not necessarily equal to i:

The infinite length random walk formulation (see eq. 5) behaves in a similar way. Our method though, always compares the same set of edges in the two graphs.

Another important difference is the complexity of our method versus the random-walk graph kernel. For an m×m kernel and n×n graph, the worst-case complexity for a length k′ random walk kernel is O(m2k′n4 ) and O(m2k′n2 ) for dense and sparse graphs, respectively [Vishwanathan et al., 2010]. The GSE, on the other hand, is always O (m2n2 ) since the heaviest operation is the Frobenius inner product in order to compute the distance between G and G′ . Moreover, once this distance is computed, evaluating the kernel for different values of ν is O(1), which combined with the fact that the variance of this kernel is Lipschitz continuous, allows for efficient searching of optimal hyperparameters (see section 2.4).

2.6 Interpretability

How could we better understand what the GSE is doing, when it maps points into an infinite-dimensional space? A successful recent development in explaining black-box models is that of Local Interpretable Model-agnostic Explanations (LIME) [Ribeiro et al., 2016], where a model is interpreted locally by making slight perturbations in the input and building an interpretable model around the new predictions. We too shall monitor our model's response to changes in the input, but instead of making random perturbations, we will perturb the input in the direction of maximum output change.

Given an instance from the dataset x1×N, where N is the number of features, and the function that will incorporate the feature connection network Gx(A) (e.g. Gx(A) = A ◦ x⊤x), we will find the direction to which the model is the most sensitive (positive and negative). Unlike optimization, where the goal is to converge as fast as possible, here we are interested in the intermediate steps of the descent. This is because we shall use the set G = {Gx1 , Gx2 , ..., GxM } and the black-box model's predictions f = {f(x1), f(x2), ..., f(xM)} to fit our interpretable model h(G) ∈ H (where xi is a variation of the original sample x0, and H represents the space of all possible interpretable functions h). This way, we will indirectly unveil the interactions that our model is most sensitive to, and show how these impact the predictions. To penalize complex models over simpler ones, we will introduce a function Ω(h) that measures model complexity. To scale the model complexity term appropriately, we can find a scalar θ so that the expected value of Ω(h) is equal to a fraction ε of the expected value of the loss:

Lastly, for highly non-linear models, the larger the input space the more complex the output explanations are likely to be, so we will weigh the sample deviations the same as the original sample x0 using the model's own similarity measure k(Gxi ,Gx0 ). Putting it all together:

where L(h, f, k(Gxi , Gx0 )) is the loss of h when using Gxi to predict the black-box model output f(xi), weighted by the kernel distance to the original sample k(Gxi ,Gx0 ).

Even Descent Sampling Method

In order to adequately cover the most sensitive regions, we need to take steps with equidistant output values. Thus, we developed a novel adaptive method to sample more in steeper regions and less in flatter ones. The intuition is that we would like to approximate the function values in unexplored regions, so that we choose an appropriate sampling step while considering the uncertainty of the approximation. Due to the stochastic nature of the method, it is able to escape local extremes. Consider the value of function f at a point x0 and its first order Taylor approximation at an arbitrary point x:

The larger the difference δ = x - x0, the less likely it is that the approximation error f(x) - ˆ f(x) is small. Assume we would like to model the random variable F, which takes the value of 1 if the approximation error is small (δ = | ˆ f (x)- f(x)| ≈ 0), and 0 otherwise. We will model the probability density function of F as being:

Consider also the random variable T which takes the value of 1 if the absolute difference in the output for a point x exceeds an arbitrary threshold (|f(x) -f(x0)| > τ ), and 0 otherwise. Assume there is zero probability this event occurs for sufficiently small steps: δ < a(τ), for some value a(τ). Let us further assume that our confidence that |f(x)-f(x0)| > τ increases linearly after the value δ = a(τ), until the maximum confidence level is reached at δ = b. After some value δ = c, we decide not to make any further assumptions about this event, so we attribute zero probability from that point on. This can be modeled as:

where v = 2c - a(τ) - b , u = b - a(τ) and T = 1, if |f(x) - f(x0| > τ and 0 otherwise. The distribution of interest is then pS = p(f = 1 ∩ t = 1|δ). To simplify the calculations, we impose the uncertainty about our approximation (expressed by F) and the likelihood of a sufficiently large output difference (expressed by T) to be independent given δ: p(f = 1 ∩ t = 1|δ) = p(t = 1|δ)p(f = 1|δ), and since the goal is to sample steps from this distribution, we will divide it by the normalization constant: Z = p(f = 1 ∩ t = 1) = ∫ +∞ -∞ p(f = 1 ∩ t = 1|δ)dδ. See Figure 3 for an illustration of the method.

There are a couple of properties that can be manipulated for a successful sampling of the output space:

Controlled Termination

To force the algorithm to terminate after a minimum number of samples Mmin have been sampled, one can decrease the value of a(τ) with each iteration so that it becomes increasingly more likely that a value of δ will be picked such that |f(x) - f(x0)| < τ, terminating the routine. For this purpose, one can compute the estimated threshold value τ0 that will keep the routine running.

where N is the number of features. This is an underdetermined equation, but one possible trivial solution is to set:

where N ′ is the number of non-zero gradient values, then let a(τ) decay with time so that it will reach this limit value after Mmin iterations:

Escaping Local Extrema

To make it more likely to escape local extrema, one possibility is to set the cut-off value c larger when the norm of τ0 (eq. 21) is larger than its expected value, and smaller otherwise:

This formulation allows jumping out from zones where the gradient is locally small, while taking smaller steps where the gradient is larger than expected.

Termination When Too Far from Original Sample

Since we are trying to explain the model locally, the sampling should terminate when the algorithm is exploring too far from the original sample. For that purpose, one can set λ to increase with increasing distance d to the original sample, pushing the probability density towards the left: λ(d) = e- d σ2 .

Putting all of the above design considerations together, you can find the complete routine in algorithm 1.

Algorithm 1 Even Descent algorithm

Input: f, x0,A

Parameter: τ, λ, θa, b, cl,Mmin

Output: X′ , f

• 1: i ← 0, fi ← f(x0), f ← [fi], converged←False
• 2: E[||τ0||] = 0, X′ ← [x0]
• 3: while converged ̸= True do
• 4: i ← i + 1
• 5: ∇f ←ComputePartialDers(x0,A, f)
• 6: τ0 ← τ/|N ′ ∗∇f|
• 7: a, b, c ←UpdatepS(i, θa,Mmin,E[||τ0||], cl)
• 8: E[||τ0||] ← (E[||τ0||](i -1) + ||τ0||)/i
• 9: δ ←EvenSample(λ, a, b, c)
• 10: xi ←xi ± δ ∗∇f
• 11: Append(f, f(xi)), Append(X′ , xi)
• 12: if |fi - fi-1| < τ then
• 13: converged ← True
• 14: end if
• 15: end while
• 16: return X′ , f

3 Experiments

3.1 Materials

For all our analysis, we used plasma protein levels of patients with suspected coronary artery disease who were diagnosed for the presence of ischaemia [Bom et al., 2018]. A total of 332 protein levels were measured using proximity extension arrays [Assarsson et al., 2014], and of the 196 patients, 108 were diagnosed with ischaemia. The protein-protein interactions data is available for download at StringDB [Jensen et al., 2009]. We implemented the GSE and the random walk kernel in python and used sci-kit learn implementation [Pedregosa et al., 2011] for the other algorithms in the comparison.

3.2 Ischaemia Classification Performance

We benchmarked the GSE performance and running time when predicting ischaemia against the random-walk graph kernel, RBF, and random forests. Additionally, in order to test the hypothesis that the protein-interaction information is improving the analysis, we also tested GSE using a constant matrix full of ones as the interaction matrix. For this benchmark, we performed a 10-cycle stratified shuffle crossvalidation split on the normalized protein data and recorded the average ROC area under the curve (AUC). To speed up the analysis, we used a training set of 90 pre-selected proteins using univariate feature selection with the F-statistic [Hira and Gillies, 2015]. The results are shown in table 1. The GSE outperformed all the other compared methods, and the fact that the GSE with a constant matrix (GSE*) had a lower performance increases our confidence that the prior interaction knowledge is beneficial for the analysis. The GSE is also considerably faster than the Random-Walk kernel, as expected. To test how both scale increasing feature size, we compared the running time of both for different pre-selected numbers of proteins. The results are depicted in Figure 4.

3.3 Performance for Different ν Values

Recall from section 2.4 that a feasible operating region for the ν values in the GSE kernel was analytically determined. We wanted to investigate how the loss function performs within this region, and whether it is possible to draw conclusions regarding the GSE kernel behaviour with respect to the interactions. To test this, the ν ∗ = maxν σ2 [k(ν)] was found using a gradient descent (ADAM [Kingma and Ba, 2015]) on the training set over 20 stratified shuffle splits (same preprocessing as in 3.2). We then measured the ROC AUC on the validation set using 12 multiples of ν ∗ . The results can be seen in Figure 5. It is quite interesting that our proxy for measuring kernel expressiveness turns out to be a convex function peaking at ν ∗ .

3.4 Interpretability Test

To test how interpretable our model's predictions are, first we trained the model on a random subset of our data and used the trained model to predict the rest of the data. Then we employed the method described in section 2.6 on a random patient in the test set, using decision trees as the interpretable models h(G) ∈ H, and a linear weighted combination of max depth and min samples per split as the complexity penalization term Ω(h). We then picked the two most important features and made a 3d plot using an interpolation of the prediction space. The result is depicted in Figure 6.

The Even Descent Sampling tests instances which are approximately equidistant in the output values. For this patient, our model 'predicts' its ischaemia risk could be mitigated by lowering protein TIMP metallopeptidase inhibitor 4 ("TIMP4") and the interaction between lipoprotein lipase ("LPL") and renin ("REN").

4 Conclusions

In this paper, we address the problem of analyzing interconnected systems and leveraging the often-known information about how the components interact. To tackle this task, we developed the Graph Space Embedding algorithm and compared it to other established methods using a dataset of proteins and their interactions from a clinical cohort to predict ischaemia. The GSE results outperformed the other algorithms in running time and average AUC. Moreover, we presented an optimal regime for the GSE in terms of a feasibility region for its parameters, which vastly decreases the optimization time. Finally, we developed a new technique for interpreting black-box models' decisions, thus making it possible to inspect which features and/or interactions are the most relevant.