{'fraction_non_alphanumeric': 0.038785986072223416, 'fraction_numerical': 0.010086869121975733, 'mean_word_length': 4.504000790279561, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Kronecker product kernel methods via generalized vec trick

Kronecker product kernel provides the standard approach in the kernel methods literature for learning from graph data, where edges are labeled and both start and end vertices have their own feature representations. The methods allow generalization to such new edges, whose start and end vertices do not appear in the training data, a setting known as zero-shot or zero-data learning. Such a setting occurs in numerous applications, including drug-target interaction prediction, collaborative filtering and information retrieval. Efficient training algorithms based on the so-called vec trick, that makes use of the special structure of the Kronecker product, are known for the case where the training data is a complete bipartite graph. In this work we generalize these results to non-complete training graphs. This allows us to derive a general framework for training Kronecker product kernel methods, as specific examples we implement Kronecker ridge regression and support vector machine algorithms. Experimental results demonstrate that the proposed approach leads to accurate models, while allowing order of magnitude improvements in training and prediction time.


INTRODUCTION
This work concerns the problem of learning supervised machine learning models from labeled bipartite graphs. Given a training set (d i , t j , y h ) n h=1 of edges (d i , t j ), where d i is the start vertex, and t j the end vertex and y h the label, the goal is to learn to predict labels for new edges with unknown labels. We assume that both the start and end vertices have their own feature representations. Further, we assume that the same vertices tend to appear as start or end vertices in multiple edges (for example, (d t , t v ), (d u , t w ), (d t , t w ) and (d u , t v ) might all belong to the same training set). This latter property is what differentiates this learning setting from the standard supervised learning setting, as the data violates the very basic i.i.d. assumption. This overlapping nature of the data has major implications both on the correct way to estimate the generalization performance of learned predictors, and on how to design efficient algorithms for learning from such data.
As discussed by [1], [2], [3], we may divide this type of learning problem into a number of distinct settings, depending whether the training and test graphs are vertex disjoint or not. In this work, we assume that the training and test graphs are vertex disjoint, that is, neither the start vertex d nor the end vertex t of a test edge is part of any edge in the training set. This is in contrast to, for example, the less challenging setting usually considered in recommender systems literature, in which the training and test graphs are only edge disjoint. There the task is to impute missing labels to customers (start vertex) × products (end vertex) matrix, where each row and column already contains at least some known labels (see Section 1.1 and [4] for more discussion). Following works such as [5], [6], we use the term zero-shot learning to describe our setting, as it requires generalization to new graphs that are vertex disjoint with the training graph. This setting is illustrated in Figure 1, where zero-shot learning corresponds to generalizing from a partially observed start vertices × end vertices label matrix to predicting unknown labels at the bottom right for edges connecting new start and end vertices that have no edges connected to them in the training set.
We consider the setting, where both the feature representations of the start and end vertices are defined implicitly via kernel functions [7]. Further, following works such as [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [6], we assume that the feature representation of an edge is defined as the product of the start and end vertex kernels, resulting in the so-called Kronecker product kernel. Kronecker product kernel methods can be trained by plugging this edge kernel to any existing kernel machine solver. However, for problems where the number of edges is large this is not feasible in practice, as typically the runtime and in many cases also the memory use of kernel solvers grows at least quadratically with respect to number of edges.
In this work, we propose the first general Kronecker product kernel based learning algorithm for large-scale problems. This is realized by generalizing the computational short-cut implied by Roth's column lemma [18] (often known as the vec-trick). As an example, we demonstrate in the experiments how the approach outperforms existing state-of-the-art SVM solver by several orders of magnitude when using the Gaussian kernel, the most well-known special case of the Kronecker  product kernel. Further, we show that the proposed methods perform well in comparison to other types of graph prediction methods. This work extends our previous work [19], that specifically considered the case of efficient ridge regression with Kronecker product kernel.

Related work
The bipartite graph learning problem considered in this work appears under various names in the machine learning literature, including link prediction [10], relational learning [15], pair-input prediction [1], and dyadic prediction [4]. The setting provides a unifying framework for prominent machine learning application areas in biology and chemistry dealing with pairwise interaction predictions, such as drug-target [3], and protein-RNA interaction prediction. The setting is prominent also in collaborative filtering, where the central task is to predict ratings or rankings for (customer, product) -pairs (see e.g. [8], [20], [4]). Also, in information retrieval the reranking of search engine results has been popularly cast as a bipartite graph ranking problem over (query, document) -pairs [21]. Many social network prediction tasks such as matching people to advertisements fall also naturally within the considered framework. Applications where the entities whose relations are predicted are of same type (i.e. protein-protein interactions, friendship in social network) can also be modeled as bipartite graph learning problems. In this case two vertices are used to represent each entity, one included in the set of start vertices and one in the set of end vertices.
In this work we consider the supervised graph learning problem, where the training edges are labeled. A limited range of graph learning problems may also be solved using unsupervised clustering methods [22], [23], for example by predicting that a relation holds between nodes assigned to the same cluster.
Kernel methods [7] are among the most widely used approaches in machine learning, both due to their strong foundations in statistical learning theory, as well as their practical applicability for a wide range of learning problems. The methods represent the data using positive semidefinite kernel functions, that allow incorporating prior knowledge about the problem domain, provide implicit non-linear mappings, and may be learned directly from data (see e.g. [24]). Recently, multi-task learning has been a topic of substantial research in the area of kernel methods (see e.g. [25], [26], [27]).
The use of Kronecker product kernel for graph learning problems was originally proposed at the same time in several independent works for collaborative filtering [8] protein-protein interaction prediction [9] and entity resolution [28]. Since then the approach has become especially popular in predicting biological interactions (see e.g. [29], [30], [3]), as well as a standard building block in more theoretical work concerning the development of multi-task learning methods (see e.g. [25], [31], [27]). Concerning specifically the zero-shot learning setting, Kronecker product kernel methods have been shown to outperform a variety of baseline methods in areas such as recommender systems [11], drug-target prediction [3] and image categorization [6].
Theoretical properties of the Kronecker product kernel were recently analyzed by [15]. They showed that if the start and end vertex kernels have the so-called universality property (e.g. Gaussian kernel) [32], then the resulting Kronecker edge kernel is also universal, resulting in universal consistency when training kernel-based learning algorithms, such as support vector machines and ridge regression using the kernel.
Recently, it has been shown that the property that the same vertices tend to appear many times in the training data can be exploited in order to develop much faster Kronecker product kernel method training algorithms. Several efficient machine learning algorithms have been proposed for the special case, where the training sets consists of a complete bipartite graph, meaning that each possible start-end vertex pair appears exactly once, and a ridge regression loss is minimized. Specifically [13], [14], [16], [17], [6] derive closed form solutions based on Kronecker algebraic optimization (see also [33] for the basic mathematical results underlying these studies). Further, iterative methods based on Kronecker product kernel matrix -vector multiplications, have been proposed (see e.g. [10], [14], [16]).
However, the assumption about the training set graph being complete can be considered a major limitation on the applicability of these methods. As far as we are aware, similar efficient training algorithms as ours for regularized kernel methods have not been previously proposed for the case, where the training graph is sparse in the sense that only a subset of the possible edges appear in the training set.
Much of the work done in multi-task learning can be considered to fall under the bipartite graph prediction, and even zero-shot learning framework [25], [34], [14], [31], [17], [6], [35]. For example, in multi-label data one might have as training data images and labels describing the properties of the images (e.g. food, person, outside). In a traditional multi-label learning setting only the images have features describing them, and we assume that the set of possible labels is fixed in advance, and that examples of each label must appear in the training set. However, if the labels have also features describing then, we may relax these assumptions by casting the task as a graph prediction problem over (data point, label) -edges (e.g. predict 1 if label t i should be assigned to data point d j and -1 if not, or rank a set of candidate labels from best matching to worst). This allows predictions even for new data points and labels both not part of the training set.
There exist two related yet distinct graph learning problems, that have been studied in-depth in the machine learning literature. Let D and T denote, respectively, the sets of start and end vertices connected to the training edges, and (d, t) a new edge, for which the correct label needs to be predicted. First, if d ∈ D and t ∈ T , the problem becomes that of matrix completion (predicting the the unknown labels within the center rectangle in Figure 1). Recently this setting has been considered especially in the recommender systems literature, where matrix factorization methods have become the de-facto standard approach for solving these problems (see e.g. [20]). The latent feature representations learned by the factorization methods, however, only generalize to edges that have both vertices connected to at least one edge in the training set. The second setting corresponds to the situation where d ∈ D and t / ∈ T , or alternatively d / ∈ D and t ∈ T (predicting the block on the right or below in Figure 1). Here predictions are needed only for new rows or columns of the label matrix, but not both simultaneously. The problem may be solved by training a standard supervised learning method for each row or column separately, or by using multi-target learning methods that also aim to utilize correlations between similar rows or columns [36].
While simultaneous use of both start and end vertex features can and has been integrated to learning also in these two related settings, they are not essential for generalization, and can even be harmful (so called negative transfer, see e.g. [17], [3]). The learning algorithms considered in this work are applicable in these related settings, provided that both start and end vertex features are available. However, simpler methods that do not need both of these information sources often provide a competitive alternative. Thus we do not consider further these two settings in this work, as these have been already quite thoroughly explored in previous literature.

GENERALIZED VEC-TRICK
In this section, we introduce the novel generalized Vec-trick algorithm (Algorithm 1) for computing such Kronecker products, where a submatrix of a Kronecker product matrix (M ⊗ N) is multiplied to a vector. The algorithm forms the basic building block for developing computationally efficient training methods for Kronecker product kernel methods.
First, we introduce some notation. By [n], where n ∈ N, we denote the index set {1, ..., n}. By A ∈ R a×b we denote an a × b matrix, and by A i,j the i, j:th element of this matrix. By vec(A) we denote the vectorization of A, which is the ab × 1 column vector obtained by stacking all the columns of A in order starting from the first column. A ⊗ C denotes the Kronecker product of A and C. By a ∈ R c we denote a column vector of size c × 1, and by a i its i:th element.
There exist several studies in the machine learning literature in which the systems of linear equations involving Kronecker products have been accelerated with the so-called "vec-trick". This is characterized by the following result known as the Roth's column lemma in the Kronecker product algebra: Lemma 1 (Roth's column lemma; "Vec trick" [18]). Let M ∈ R a×b , Q ∈ R b×c , and N ∈ R c×d be matrices. Then, It is obvious that the right hand side of (1) is considerably faster to compute than the left hand side, because it avoids the direct computation of the large Kronecker product. We next shift our consideration to matrices that correspond to a submatrix of a Kronecker product matrix. First, we introduce some further notation. To index certain of rows or columns of a matrix, the following construction may be used: T ← 0 ∈ R d×a 3: for h = 1, . . . , e do 4: i, j ← r h , t h
where e i is the ith standard basis vector of R a , is a row indexing matrix for M determined by s. The column indexing matrices are defined analogously.
The above type of ordinary indexing matrices can of course be used with Kronecker product matrices as well. However, the forthcoming computational complexity considerations require a more detailed construction in which we refer to the indices of the factor matrices rather than to those of the product matrix. This is characterized by the following lemma which follows directly from the properties of the Kronecker product: Lemma 2 (Kronecker product index matrix). Let M ∈ R a×b , N ∈ R c×d and let S be an index matrix for the Kronecker product M ⊗ N. Then, S can be expressed as where p = (p 1 , . . . , p f ) T ∈ [a] f and q = (q 1 , . . . , q f ) T ∈ [c] f are sequences of row indices of M and N, respectively. The entries of p and q are given as where i ∈ [ac] denotes a row index of M ⊗ N, and p ∈ [a] and q ∈ [c] map to the corresponding rows in M and N, respectively. The column indexing matrices are defined analogously.
Proof. The claim is an immediate consequence of the definition of the Kronecker product, where the relationship between the row indices of M and N with those of their Kronecker product is exactly the one given in (2).

Theorem 1.
Let M ∈ R a×b , N ∈ R c×d , and let R ∈ {0, 1} f ×ac and C ∈ {0, 1} e×bd to be, respectively, a row index matrix and a column index matrix of M ⊗ N, such that R is determined by the sequences p = (p 1 , . . . , n number of labeled edges in training set m number of unique start vertices connected to training edges q number of unique end vertices connected to training edges d number of start vertex features r number of end vertex features and C by r = (r 1 , . . . , r e ) T ∈ [b] e and t = (t 1 , . . . , t e ) T ∈ [d] e . Further, we assume that the sequences p, q, r and t when considered as mappings, are all surjective. Further, let v ∈ R e . The product can be computed in O(min(ae + df, ce + bf )) time using the generalized Vec trick (Algorithm 1).
Proof. Since the sequences p, q, r and t, determining R and C are all surjective on their co-domains, we imply that Then, according to Lemma 1,(3) can be rewritten as which does not involve a Kronecker product. Calculating the right hand side can be started by first computing either S = NV or T = VM T requiring O(ce) and O(ae) time, respectively, due to V having only at most e nonzero entries. Then, each of the f elements of the product can be computed either by calculating the inner product between a row of S and a row of M or between a row of N and a column of T, depending whether S or T was computed. The former inner products require b and the latter d flops, and hence the overall complexity is O(min(ae + df, ce + bf )). The pseudocode for the algorithm following the ideas of the proof is presented in Algorithm 1. We note that the index matrices R and C are not explicitly given as input to the algorithm, because the elements of these sparse matrices are encoded in the vectors p, q, r and t.
Remark 1 (Vec trick as special case). If R = C = I, the generalized vec trick reduces to the standard vec trick, and the complexity of Algorithm 1 is O(min(abd + acd, bcd + abc)), the same as for the standard vec trick in Lemma 1.

GENERAL FRAMEWORK FOR KRONECKER PRODUCT KERNEL METHODS
In this section, we connect the above abstract considerations with practical learning problems with bipartite graph data. We start by defining some notation which will be used through the forthcoming considerations. As observed from Figure 1, the relevant data dimensions are the sizes of the two data matrices and the number of edges with known labels. The notations are summarized in Table 1.
The training set consists of a labeled sequence of edges (d r h , t s h , y h ) n h=1 ∈ (D × T × Y) n of a labeled bipartite graph, each edge consisting of a start vertex d r h ∈ D, end vertex t s h ∈ T and an associated label y h , where D and T are the spaces of start and end vertices, respectively, and Y is the label space (Y = R for regression and Y = {−1, 1} for binary classification). Moreover, let {d i } m i=1 ⊂ D and {t i } q i=1 ⊂ T denote the sets of such start vertices and end vertices, respectively, that are encountered in the training set as part of at least one edge. Further, let r = (r 1 , . . . , r n ) T ∈ [q] n and s = (s 1 , . . . , s n ) T ∈ [m] n be index sequences that map the training edges to their corresponding start and end vertices, respectively. Looking at Figure 1, one may think of r i and s i as row and column indices for the known entries in the start vertices times end vertices matrix, that correspond to the training data.
Let us assume a randomly drawn new sequence of edges , and the labels are unknown. Our aim is to learn from the training set a prediction function f : D × T → Y, such that can correctly predict the labels of such edges. We further assume that the graphs corresponding to the train and test data are completely disconnected, that is, the sets of their start vertices are mutually disjoint and so are the sets of their end vertices. This is the main difference of our setting to the types of problems solved typically for example within the recommender systems literature, rather our setting corresponds to the so-called zero-shot, or cold-start problem where test vertices do not appear in the training set.
Let k : D × D → [0, ∞) and g : T × T → [0, ∞) denote positive semi-definite kernel functions [7] defined for the start and end vertices respectively. Further, the Kronecker product kernel k ⊗ (d, t, d , t ) = k(d, d )g(t, t ) for the edges is defined as the product of these two base kernels. Let K i,j = k(d i , d j ) and G i,j = g(t i , t j ) denote the start and end vertex kernel matrices for the training data, respectively. Further, let R ∈ {0, 1} n×mq be the Kronecker product index matrix determined by the start and end vertex index sequences r and s, as defined in Lemma 2. Then R(G ⊗ K)R T is the Kronecker product kernel matrix containing the kernel evaluations corresponding to the training edges. Note that the methods to be introduced will never explicitly form this prohibitively large Kronecker product matrix, rather kernel matrix multiplications are implemented with Algorithm 1, that requires as input G and K, as well as the index sequences r and s that implicitly define R. 6 We consider the regularized risk minimization problem with where p ∈ R n and y ∈ R n are the predicted and correct outputs for the training set, L a convex nonnegative loss function and λ > 0 a regularization parameter. Choosing as H the reproducing kernel Hilbert space defined by k ⊗ , the generalized representer theorem [37] guarantees that the optimal solutions are of the form where we call a ∈ R n the vector of dual coefficients. Now let us consider the special case, where the start and end vertex spaces are real vector spaces, that is, D = R d and T = R r , and hence both the start and end vertices have a finite dimensional feature representation. Further, if the start and end vertex kernels are linear, the Kronecker product kernel can be written open as the inner product that is the edges have an explicit feature representation with respect to this kernel as the Kronecker product of the two feature vectors. Let D ∈ R m×d and T ∈ R q×r , respectively, contain the feature representations of the training set start and end vertices. Now the joint Kronecker feature representation for the training data can be expressed as X = R(T ⊗ D). In this case we can equivalently define the regularized risk minimizer as where we call w ∈ R dr the vector of primal coefficients.

Efficient prediction
Next, we consider how predictions can be efficiently computed for graph data both with dual and primal predictors. This operation is necessary both during training when training set predictions are needed to evaluate the quality of the current predictor, as well as when applying the final trained predictor on new edges. Let us assume a new sequence of edges where 1 < i < u, 1 < j < v, and d i and t j may or may not overlap with the training set start and end vertices. Let K ∈ R u×m and G ∈ R v×q contain the kernel evaluations between the training start and end vertices, and the new start and end vertices, respectively. Further, the Kronecker product index matrix R ∈ {0, 1} t×uv encodes those new edges, for which the predictions are needed for.
For the dual model, the predictions can be computed as resulting in a computational complexity of O(min(vn + mt, un + qt)). Further, we note that for sparse models where a large portion of the a i coefficients have value 0 (most notably, support vector machine predictors), we can further substantially speed up prediction by removing these entries from a and R and correspondingly replace the term n in the complexity with the number of non-zero coefficients. In this case, the prediction complexity will be where a 0 is the zero-norm measuring the number of non-zero elements in a. This is in contrast to the explicit computation by forming the full test kernel matrix, which would result in complexity.
In the primal case the predictions can be computed as resulting in a computational complexity of O(min(vdr + dt, udr + rt)). For high-dimensional data (dr >> n) one will save substantial computation by using the dual model instead of the primal. 2 Loss functions and their (sub)gradients and (generalized) Hessians. For binary classification losses marked * we assume that y i ∈ {−1, 1}, while for the rest y i ∈ R.

Efficient learning
The regularized risk minimization problem provides a convex minimization problem, whose optimum can be located with (sub)gradient information. Next, we show how to efficiently compute the gradient of the objective function, and for twice differentiable loss functions Hessian-vector products, when using the Kronecker product kernel. For non-smooth losses or losses with non-smooth derivatives, we may instead consider subgradients and the generalized Hessian matrix (see e.g. [40], [43], [44]). In this section, we will make use of the denominator-layout notation.
First, let us consider making predictions on training data in the dual case. Based on previous considerations, we can compute The regularizer can be computed as λ 2 f 2 H = λ 2 a T R(G ⊗ K)R T a. We use the following notation to denote the gradient and the Hessian (or a subgradient and/or generalized Hessian) of the loss function, with respect to p: The exact form of g and H depends on the loss function, Table 2 contains some common choices in machine learning. While maintaining the full n × n Hessian would typically not be computationally feasible, for univariate losses the matrix is diagonal. Further, for many multivariate losses efficient algorithms for computing Hessian-vector products are known (see e.g. [42] for examples of efficiently decomposable ranking losses). The gradient of the empirical loss can be decomposed via chain rule as The gradient of the regularizer can be computed as λR(G ⊗ K)R T a. Thus, the gradient of the objective function becomes The Hessian of L with respect to a can be defined as Hessian for the regularizer is defined as λR(G ⊗ K)R T . Thus the Hessian for the objective function becomes For commonly used univariate losses, assuming p has been computed for the current solution, g, and H can be computed in O(n) time. The overall cost of the loss, gradient and Hessian-vector product computations will then be dominated by the efficiency of the Kronecker-kernel vector product algorithm, resulting in O(qn + mn) time.
Conversely. in the primal case we can compute the predictions as p = R(T ⊗ D)w, loss gradient as (T T ⊗ D T )R T g and Hessian for the loss as (T T ⊗ D T )R T HR(T ⊗ D), (both w.r.t. w). For the regularizer the corresponding values are λ 2 w T w, λw and λI. The overall complexity of the loss, gradient and Hessian-vector product computations is O(min(qdr + dn, mdr + rn)). Require: R ∈ {0, 1} n×mq , K ∈ R m×m , G ∈ R q×q , y ∈ R n , λ > 0 1: a ← 0 2: repeat 3: p ← R(G ⊗ K)R T a 4: COMPUTE H, g 5: SOLVE (HR(G ⊗ K)R T + λI)x = g + λa 6: a ← a − δx (δ constant or found by line search) 7: until convergence Algorithm 3 Primal Truncated Newton optimization for regularized risk minimization 6: w ← w − δx (δ constant or found by line search) 7: until convergence

Optimization framework
As a simple approach that can be used for training regularized risk minimization methods with access to gradient and (generalized) Hessian-vector product operations, we consider a truncated Newton optimization scheme (Algorithms 2 and 3; these implement approaches similar to [40], [45]). This is a second order optimization method that on each iteration computes the Newton step direction ∂ 2 J(f )x = ∂J(f ) approximately up to a pre-defined number of steps for the linear system solver. The approach is well suited for Kronecker product kernel method optimization, since while computing the Hessians explicitly is often not computationally feasible, computing Hessian-vector products can be done efficiently using the generalized vec trick algorithm.
Alternative optimization schemes can certainly be used, such as the Limited-memory BFGS algorithm [46] or trustregion Newton optimization [43]. Further, for non-smooth losses such as the hinge loss, methods tailored for non-smooth optimization such as the bundle method [44] should be used rather than the truncated Newton method. Such works are orthogonal to our work, as the proposed algorithm can be used to speed up computations for any optimization method that is based on (sub)gradient, and possibly (generalized) Hessian-vector product computations. Optimization methods that process either the edges or the model coefficients individually or in small batches (e.g. stochastic gradient descent [47], coordinate descent [48], SVM decomposition methods [49]) may however not be a good match for the proposed algorithm, as the largest speed-up is gained when doing the computations for all of the training data in one single batch.
Each iteration of the Truncated Newton algorithm starts by computing the vector of training set predictions p for the current solution, after which the gradient g and Hessian H can be computed (see Table 2 for examples on how to compute these for several widely used loss functions). After this, the algorithm solves approximately the linear system to find the next direction of descent (for the primal case we solve analogously ∂ 2 J ∂w 2 x = ∂J ∂w ). Here, we may use methods developed for solving linear systems, such as the Quasi-Minimal Residual iteration method [50] used in our experiments. For the dual case, we may simplify this system of equations by removing the common term R(G ⊗ K)R T from both sides, resulting in the system (HR(G ⊗ K)R T + λI)x = g + λa .
While the optimization approaches can require a large number of iterations in order to converge, in practice often good predictive accuracy can be obtained using early stopping both in solving the system of linear equations on line 5 of Algorithm 2 and line 5 of Algorithm 3 and the outer truncated Newton optimization loops, as there is no need to continue optimization once the error of the prediction function stops decreasing on a separate validation set (see e.g. [51], [52], [53]).

Complexity comparison
The operations that dominate the computational cost for Algorithms 2 and 3 are matrix-vector products of the form R(G ⊗ K)R T v and (T T ⊗ D T )R T v. A natural question to ask is, how much does the proposed fast generalized vec Trick algorithm (Algorithm 1) speed these operations up, compared to explicitly forming the Kronecker product matrices (here denoted the 'Baseline' method). 9 TABLE 3 Dual case: complexity comparison of the proposed method and the baseline approach that constructs the Kronecker product kernel matrix explicitly. We consider three different settings: • Independent: all edges in the training set are mutually disjoint, indicating that n = m = q.
• Dependent: training edges may share start or end vertices or both with each other, that is max(m, q) < n < mq.
• Complete: the training set is a complete bipartite graph, that is n = mq.
'Dependent' is the default setting assumed in this paper, while 'Independent' and 'Complete' are extreme cases of the setting.
In Table 3 we provide the comparison for the dual case. In the 'Independent' -case the complexity of the proposed algorithm reduces to that of the baseline method, since there are vertices shared between the edges, that the algorithm could make use of. However, beyond that the complexity of the baseline approach grows quadratically with respect to the number of edges, while the complexity of the proposed method grows linearly with respect to the numbers of start vertices, end vertices, and edges. In Table 4 we provide the comparison for the primal case. Again, the complexity is the same for the 'Independent' case, but the proposed method is much more efficient for the other settings assuming m << n and q << n. Finally, we note that for both the primal and the dual settings, for the 'Complete' case where R = I, the multiplication could as efficiently be computed using the 'VecTrick' -method implied by Roth's column lemma (1). To conclude, we observe that the proposed algorithm considerably outperforms previously known methods, if max(m, q) << n < mq.

CASE STUDIES: RIDGE REGRESSION AND SUPPORT VECTOR MACHINES
Next, as an example on how our framework may be applied we consider two commonly used losses, deriving fast ridge regression and support vector machine training algorithms for the Kronecker product kernel.

Kronecker ridge regression
As our first example learning algorithm that uses the above defined concepts, let us consider the well-known ridge regression method (aka regularized least-squares, aka least-squares SVM) [38], [54]. For this specific case, the resulting algorithm is simpler than the general optimization framework of Algorithms 2 and 3, as the linear system appearing in the algorithms needs to be solved only once.
In this case, L = 1 2 (p − y) T (p − y), H = I and g = p − y (see Table 2). For the dual case, based on (8) it can be observed that the full gradient for the ridge regression is Writing p open (7) and setting the gradient to zero, we recover a linear system A solution to this system can be recovered by canceling out R(G ⊗ K)R T from both sides, resulting in the linear system (R(G ⊗ K)R T + λI)a = y This system can be solved directly via standard iterative solvers for systems of linear equations. Combined with the generalized vec trick algorithm, this results in O(mn + qn) cost for each iteration of the method.
For the primal case, the gradient can be expressed as Again, writing p open and setting the gradient to zero, we recover a linear system Solving this system with a linear system solver together with the generalized Kronecker product algorithm, each iteration will require O(min(mdr + nr, drq + dn)) cost. The naive approach of training ridge regression with the explicitly formed Kronecker product kernel or data matrix using standard solvers of systems of linear equations would result in O(n 2 ) and O(ndr) cost per iteration for the dual and primal cases, respectively.

Kronecker support vector machine
Support vector machine is one of the most popular classification methods in machine learning. Two of the most widely used variants of the method are the so-called L1-SVM and L2-SVM (see Table 2). Following works such as [40], [45], [43] we consider the L2-SVM variant, since unlike L1-SVM it has an objective function that is differentiable and yields a non-zero generalized Hessian matrix, making it compatible with the presented Truncated Newton optimization framework.
The loss can be defined as } denote the subset of training set that have non-zero loss for given prediction function. Further, let S S denote the index matrix corresponding to this index set. Now we can re-write the loss in a least-squares form as L = 1 2 (S S p − S S y) T (S S p − S S y), its gradient as g = S T S (S S p−S S y), and the Hessian as H = S T S S S , which is a diagonal matrix with entry 1 for all members of S, and 0 otherwise.
By inserting g to (8) and writing p open (7) we recover the gradient of the L2-SVM objective function, with respect to a, as Inserting g and H to (9), we see that ∂J ∂a∂a x = ∂J ∂a can be solved from: The gradient and Hessian-vector products can be computed again at O(qn + mn) time. However, it should be noted that the support vector machine algorithm encourages sparsity, meaning that as the optimization progresses a may come to contain many zero coefficients, and the set S may shrink (at the optimal solution these two sets coincide, denoting the so-called support vectors). Thus given a new solution a, we may compute p = R(G ⊗ K)R T a and thus also the right-hand size of (4.2) in O(min(q a 0 + mn, m a 0 + qn) time, by removing the zero-coefficients from a and the corresponding columns from R T . Further, when solving (4.2) with Truncated Newton optimization, we can in each inner iteration of the method compute matrix-vector multiplication In the primal case the gradient can be written as

generalized Hessian matrix can be defined as
The most efficient existing SVM solvers can be expected at best to scale quadratically with respect to training set size when solving the dual problem, and linearly with respect to training set size and number of features when solving the primal case [55], [56]. Thus using the full Kronecker product kernel or data matrices, such solvers would have O(n 2 ) and O(ndr) scaling for the dual and primal cases, respectively.

EXPERIMENTS
Computational costs for iterative Kronecker product kernel method training depend on two factors: the cost of a single iteration, and the number of iterations needed to reach a good solution. The cost of a single iteration is dominated by gradient computations and Hessian-vector products. These can be efficiently performed with Algorithm 1. Further, as discussed at the end of Section 3.3, the number of iterations needed can be limited via early stopping once a predictor that works well on independent validation data has been reached. For prediction times, the dominating costs are the data matrix or kernel matrix multiplications with the primal or (possibly sparse) dual coefficient vectors, both operations that can be speeded up with Algorithm 1.
In the experiments, we demonstrate the following: 1) Parameter selection: we show how the hyperparameters of the methods can be tuned using validation data.
2) Fast training: the combination of the proposed short-cuts and early stopping allows orders of magnitude faster training than with regular kernel method solvers. 3) Fast prediction: proposed short-cuts allow orders of magnitude faster prediction for new data, than with regular kernel predictors. 4) Accurate predictions: the predictive accuracy of Kronecker models is competitive compared to alternative types of scalable graph learning methods.
We implement the Kronecker methods in Python, they are made freely available under open source license as part of the RLScore machine learning library 1 [57]. For comparison, we consider the LibSVM software [58], which implements a highly efficient support vector machine training algorithm [49]. Further, we compare our approach to stochastic gradient descent and K-nearest neighbor based graph prediction methods. As a typical bipartite graph learning problem, we consider the problem of predicting drug-target interactions. From a data base of known drugs, targets and their binary interactions, the goal is to learn a model that can for new previously unseen drugs and targets predict whether they interact. We consider four drug-target interaction data sets, the GPCR, IC, E data sets [59], and the Ki data [60]. The Ki-data is a naturally sparse data set where the class information is available only for a subset of the edges. For the other four drug-target data sets we use for the experiments a subset of all the possible drug-target interactions. We use exactly the same preprocessing of the data as [3]. The characteristics of the data sets are described in Table 5. The edges in the data are (d i , t j , y h ) triplets, consisting of a pair of vectors encoding the features of the drug and the target (see [3] for details of the features), and a label having value 1 if the drug and target interact, and −1 if they do not. Mapping the problem to the terminology used in this paper, we may consider the drugs as start vertices, the targets as end vertices, and the drug-target pairs with known interaction or non-interaction as edges. Further, we generated two data sets using a variant of the Checkerboard simulation, a standard non-linear problem for benchmarking large scale SVM solvers (see e.g. [61]). In our simulation both start and end vertices have a single feature describing them, drawn from continuous uniform distribution in range (0, 100). The output assigned to an edge (d, t) is +1 whenever both d and t are either odd or even, and -1 when one of them is odd and the other even. When plotting the output against d and t, this results in a highly non-linear checkerboard-type of pattern. Finally, random noise is introduced to the data by flipping with 0.2 probability the class of each edge. The numbers of start and end vertices are same, and labels are assigned for 25% of all the possible edges (i.e. m = q and n = 0.25m 2 ).

Data and general setup
In the zero-shot learning setting the aim is to generalize to such (d i , t j )-edges that are vertex disjoint with the training graph. Therefore, cross-validation with graph data is more complicated than for standard i.i.d. data, since this aim must be reflected in the train-test split (see e.g. [1], [3]). The split is illustrated in Figure 2. To ensure that the training and test graphs are vertex disjoint, the edges are divided into training and test edges as follows. First, both the start vertex-indices [1, ..., m] and end vertex-indices [1, ..., q] are randomly divided into a training and test part. Then, an edge (d i , t j , y h ) is assigned to training set if i belongs to the training start vertex indices and j belongs to the training end vertex indices. It is assigned to the test set if i belongs to the test start vertex indices and j belongs to test end vertex indices. Finally, the rest of the edges are discarded, that is, they belong neither to the training nor test part (the greyed out blocks in Figure 2). To tune the hyper-parameters without the risk of overfitting, one can split the data into training, validation and test parts in an analogous way to the above described train-test split. Combining this with cross-validation is in the literature known as nested cross-validation (for a detailed description of this approach for graph learning, see [3]).
Ridge regression is trained with the minimum residual iteration algorithm [62] implemented in the scipy.sparse.linalg.minres package, while the inner optimization loop of the SVM training algorithm uses the scipy.sparse.linalg.qmr implementing of quasi-minimal residual iteration [50] (SciPy version 0.14.1), with δ = 1. Regularization parameters from the grid [2 −20 , ..., 2 20 ] were tested. We restrict our plots to values [2 −10 , 2 −5 , 2 0 , 2 5 , 2 10 ], as these allow representing all the main trends in the experiments. The classification performance on test data was measured with area under ROC curve (AUC). The experiments were run on a desktop computer with Intel Core i7-3770 CPU (3.40GHz) running Ubuntu Linux 15.04.
Early stopping experiments were run with the linear kernel, in order to allow comparing dual and primal optimization. The LibSVM comparison was done using the Gaussian kernel. LibSVM does not directly support the Kronecker product kernel, this issue was resolved as follows. If both start vertex and end vertex kernels are Gaussian with width γ, then , that is, the Kronecker product kernel is equal to using the Gaussian kernel with concatenated features of the start and end vertex. Fig. 2. Ninefold cross-validation. The matrix rows represent start vertices and columns end vertices, the edges with known labels corresponds to a subset of the elements of this matrix. The row-and column indices are both divided into three sets. The Cartesian product sets of the row and column sets index nine blocks of possible edges. On each round, the test fold is formed from all the edges that belong to one of the blocks, and further have also known label. The training folds are formed from the four blocks, that share no common rows or columns with the test fold. Four blocks are left unused, as their edges connect vertices belonging to training and test folds.

Choosing hyperparameters
In this section, we set up guidelines for selecting the values of the hyperparameters based on experimental verification. The Kronecker algorithms have the following hyperparameters: the regularization parameter, the number of iterations for ridge regression, and both the number of inner and outer iterations for the SVM. We consider only the dual optimization, the observed behavior was very similar also for the primal case.
We run the optimization up to 100 iterations measuring regularized risk J(f ), and AUC on the test set (in selected experiments with 500 iterations we found little improvements). Ridge regression results are presented in Figure 3, while SVM results with inner optimization loop terminated after 10 and 100 iterations are presented in Figures 4 and 5.
In all the experiments the optimal test set AUC is obtained within tens of iterations. Beyond this point reduction in regularized risk no longer translates into better predictions. Further, for SVMs increasing the number of inner iterations from 10 to 100 allows achieving much faster decrease in regularized risk. However, this comes at the cost of having to perform ten times more computation each iteration, and does not lead into faster increase in test set AUC. Thus it can be observed, that rather than having to solve exactly the large optimization problems corresponding to training Kronecker product kernel methods, often only a handful of iterations need to be performed in order to obtain maximal predictive accuracy.
Several observations can be made based on the experiments. The regularized risk decreases quite quickly even if the SOLVE operation in Algorithms 2 and 3 is terminated after a small number of iterations. Moreover, early termination of SOLVE provides us more fine-grained control of the degree of fitting (contrast Figures 4 and 5, where for 100 inner iterations test performance in some cases starts immediately decreasing).
To conclude, the results suggest that on large data sets a good strategy is to start by setting the number of iterations to a small constant (e.g. 10 inner and outer iterations), and increase these parameters only if the predictive accuracy keeps increasing beyond the first iterations on independent validation set. Furthermore, we note that one could sidestep the selection of regularization parameter λ by setting it to a small constant and regularizing only with early stopping, however tuning also λ on separate validation data can sometimes yield even better performance.

Training time
In order to demonstrate the improvements in training speed that can be realized using the sparse Kronecker product algorithm, we compare our Kronecker SVM algorithm to the LibSVM solver on the Ki-data set. Based on preliminary tests, we set γ = 10 −5 , as this value produces informative (not too close to identity matrix, or to matrix full of ones) kernel matrices for both the start and end vertices. For the Kronecker SVM, we perform 10 inner and 10 outer iterations with λ = 2 −5 . For LibSVM, we present the values for the grid [2 −7 , 2 −5 , 2 −3 , 2 −1 ], as results for it vary more based on the choice of regularization parameter. We perform 9-fold cross-validation on the K i data, for various training set sizes. In Figure 6 (left) we present the running times for training the methods for different numbers of training edges, as well as the corresponding cross-validated AUCs. The KronSVM results, while similar as before, are not directly comparable to those in the earlier experiments due to different kernel function being used.
As can be expected based on the computational complexity considerations, the KronSVM algorithm has superior scalability compared to the regular SVM implemented in the LibSVM package. In the experiment on 42000 edges the 2 · 10 6 4 · 10 6 6 · 10 6 8 · 10 6 10 7 Training edges difference is already 25 seconds versus 15 minutes. The LibSVM runtimes could certainly be improved for example by using earlier termination of the optimization. Still, this would not solve the basic issue that without using computational shortcuts such as the generalized Kronecker product algorithm (Algorithm 1) proposed in our work, one will need to construct a significant part of the kernel matrix for the edges. For LibSVM, the scaling is roughly quadratic in the number of edges, while for KronSVM it is linear. Considering the cross-validated AUCs (Figure 6 (middle)), it can be observed that the AUC of the KronSVM with early stopping is very much competitive with that of LibSVM. To conclude, the proposed Kronecker product algorithm allows substantially faster training for graph data than if one uses current state-of-the-art solvers, that are however not able to make use of the shared structure.

Prediction time
Next, we show how the sparse Kronecker product algorithm can be used to accelerate predictions made for new data. As in previous experiment, we train LibSVM on the Ki-data for varying data set sizes, with the training and test split done as in the 9-fold cross-validation experiment. We use λ = −5, and the same kernel and kernel parameters as before. We use the predictor learned by the SVM to make predictions for the 10000 drug-target pairs in the test set. We compare the running times for two approaches for doing the predictions. 'Baseline' refers to the standard decision function implemented in LibSVM. 'Kronecker', refers to implementation that computes the predictions using the sparse Kronecker product algorithm. The 'Kronecker' method is implemented by combining the LibSVM code with additional code that after training LibSVM, reads in the dual coefficients of the resulting predictor and generates a new predictor that uses the shortcuts proposed in this paper. Both predictors are equivalent in the sense that they produce (within numerical accuracy) exactly the same predictions.
The results are plotted in Figure 6 (right). For both methods the prediction time increases linearly with the training set size (see equations (5) and (6)). However, the Kronecker method is more than 1000 times faster than Baseline, demonstrating that significant speedups can be realized by using the sparse Kronecker product algorithm for computing predictions.

Large scale experiments with simulated data
In order to demonstrate the scaling of the proposed Kronecker product kernel methods to larger problem sizes than those encountered in the drug-target data sets, we implemented scalability experiments on various sized subsets of the checkerboard data. The training set is formed as follows. First, we generate the same number of start vertices m and end vertices q, then generate The independent test set is generated in the same way. We use the Gaussian kernel, based on preliminary tests we set λ = 2 −7 and γ = 1, as parameters around this range allow learning the simulated function. As before, KronSVM uses 10 inner and 10 outer iterations. We train both KronSVM and LibSVM for varying data set sizes, up until they reach the point where training takes more than 24 hours to complete. We also measure how long computing predictions takes for a test set of the same size as the training set. Finally, we also measure test set AUC in order to show that the learners can really learn to solve the simulated problem (note that due to random flipping of classes, even the optimal predictor would have only 0.8 AUC). The results are presented in Figure 7. Again, KronSVM outperforms LibSVM by several orders of magnitude. KronSVM can be trained in 24 hours on approximately 10 million edges (correspondingly, with 6400 start and end vertices). The LibSVM experiments were discontinued after training on 64000 edges (correspondingly 1600 start and end vertices) took more than 27 hours. For the same training set size, KronSVM can be trained in 23 minutes. Regarding prediction times, for LibSVM model trained with 64000 edges it took 93 minutes to compute predictions for a test set of same size, whereas with the generalized Kronecker product shortcuts the same computations can be done in 7 seconds. KronSVM can with a model trained on 10 million edges, make predictions for a test set of also 10-million edges in 7 minutes. When trained on 10 million edges, the KronSVM implementation used roughly 1.5 Gigabytes of memory.

Comparison of graph learning methods
Finally, we compare the proposed approach to alternative scalable graph learning approaches, such that can generalize to making predictions for edges that are vertex disjoint with the training graph. We consider the following baseline methods, that use as feature representation the concatenation [d, t] of the start and end vertex features: • Linear model, stochastic gradient descent (SGD) [47]: We fit a linear model f (d, t) = w, [d, t] to the data using stochastic gradient descent over the edges. The approach is extremely scalable, and it is not necessary to load all data into memory at once. Previously, [4] has used sgd with logistic regression for cold start learning with recommender systems. We consider both the logistic and the hinge loss.
• K-nearest neighbors (KNN): KNN methods have enjoyed substantial popularity in graph prediction applications such as biological interaction prediction [63] and recommender systems [64]. The method can model highly nonlinear functions and can scale well especially to low-dimensional problems by using efficient data structures for speeding up the neighborhood-search.
For the baseline methods, we use the implementations from the scikit-learn package [65]. The sgd regularization parameter and KNN number of neighbors parameters are selected with internal 3-fold cross-validation. The number of stochastic gradient descent updates is set to 10 6 (or at minimum one full pass through data). For KronSVM and KronRidge we set λ = 0.0001, as the methods were not very sensitive to amount of regularization when optimization was terminated early. KronSVM uses 10 inner and 10 outer iterations on each data set, while KronRidge uses 100 iterations. As before, for Checker data sets the Kronecker methods use Gaussian kernel with γ = 1 for both vertex kernels, for other data sets we use linear vertex kernels. For the drug-target interaction data sets we use 3x3 fold cross-validation as before, for the checkerboard data sets we generate separate test set with 6250000 edges.
In Table 6 we present the AUCs for the compared methods. Overall, KronSVM performs the best, yielding the best performance on Ki, E, Checker and Checker+ data sets. KronRidge performs slightly worse than KronSVM on most data sets, possibly due to the fact that the squared loss is not as well suited for classification as the squared hinge loss. While the linear SGD methods provide a surprisingly competitive baseline, they do not quite reach the performance of the best methods, and it is impossible for them to outperform random guessing for the Checker data sets, due to the non-linearity of the task. The KNN method performs reliably over all the data sets, but does not yield the best performance on any of them.
Regarding the runtimes presented in Table 7, KronSVM and KronRidge strike a good balance between accurate predictions and the ability to scale to all the data sets. The linear SGD methods provide overall the best scalability, but at the cost of not being able to model nonlinearities in the data. The scalability of the KNN depends on the dimensionality of the data; on Checker and Checker+ the method excels because there are only 2 features, whereas on Ki, IC, E, and GPCR, the method is not competitive.

CONCLUSION
In this work, we have proposed a generalized Kronecker product algorithm. A simple optimization framework is described in order to show how the proposed algorithm can be used to develop efficient training algorithms for pairwise kernel methods. Both computational complexity analysis and experiments show that the resulting algorithms can provide order of magnitude improvements in computational efficiency both for training and making predictions, compared to existing kernel method solvers. Further, we show that the approach compares favorably to other types of graph learning methods. The implementations for the generalized Kronecker product algorithm, as well as for Kronecker product kernel learners are made freely available under open source license.