Implicit SVD for Graph Representation Learning
Year: 2021

Abstract:
Recent improvements in the performance of state-of-the-art (SOTA) methods for Graph Representational Learning (GRL) have come at the cost of significant computational resource requirements for training, e.g., for calculating gradients via backprop over many data epochs. Meanwhile, Singular Value Decomposition (SVD) can find closed-form solutions to convex problems, using merely a handful of epochs. In this paper, we make GRL more computationally tractable for those with modest hardware. We design a framework that computes SVD of implicitly defined matrices, and apply this framework to several GRL tasks. For each task, we derive linear approximation of a SOTA model, where we design (expensiveto-store) matrix M and train the model, in closed-form, via SVD of M, without calculating entries of M. By converging to a unique point in one step, and without calculating gradients, our models show competitive empirical test performance over various graphs such as article citation and biological interaction networks. More importantly, SVD can initialize a deeper model, that is architected to be nonlinear almost everywhere, though behaves linearly when its parameters reside on a hyperplane, onto which SVD initializes. The deeper model can then be fine-tuned within only a few epochs. Overall, our procedure trains hundreds of times faster than state-of-the-art methods, while competing on empirical test performance. We open-source our implementation at: https://github.com/samihaija/isvd

heading: 1 Introduction
text: Truncated Singular Value Decomposition (SVD) provides solutions to a variety of mathematical problems, including computing a matrix rank, its pseudo-inverse, or mapping its rows and columns onto the orthonormal singular bases for low-rank approximations. Machine Learning (ML) software frameworks (such as TensorFlow) offer efficient SVD implementations, as SVD can estimate solutions for a variaty of tasks, e.g., in computer vision [Turk and Pentland, 1991], weather prediction [Molteni et al., 1996], recommendation [Koren et al., 2009], language [Deerwester et al., 1990, Levy and Goldberg, 2014], and more-relevantly, graph representation learning (GRL) [Qiu et al., 2018].
SVD’s benefits include training models, without calculating gradients, to arrive at globally-unique solutions, optimizing Frobenius-norm objectives (§2.3), without requiring hyperparameters for the learning process, such as the choice of the learning algorithm, step-size, regularization coefficient, etc. Typically, one constructs a design matrix M, such that, its decomposition provides a solution to a task of interest. Unfortunately, existing popular ML frameworks [Abadi et al., 2016, Paszke et al., 2019] cannot calculate the SVD of an arbitrary linear matrix given its computation graph: they
∗Part of this work was done during internship at Intel Labs.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
compute the matrix (entry-wise) then2 its decomposition. This limits the scalability of these libraries in several cases of interest, such as in GRL, when explicit calculation of the matrix is prohibitive due to memory constraints. These limitations render SVD as impractical for achieving state-of-the-art (SOTA) for tasks at hand. This has been circumvented by Qiu et al. [2018] by sampling M entrywise, but this produces sub-optimal estimation error and experimentally degrades the empirical test performance (§7: Experiments).
We design a software library that allows symbolic definition of M, via composition of matrix operations, and we implement an SVD algorithm that can decompose M from said symbolic representation, without need to compute M. This is valuable for many GRL tasks, where the design matrix M is too large, e.g., quadratic in the input size. With our implementation, we show that SVD can perform learning, orders of magnitudes faster than current alternatives.
Currently, SOTA GRL models are generally graph neural networks trained to optimize cross-entropy objectives. Their inter-layer non-linearities place their (many) parameters onto a non-convex objective surface where convergence is rarely verified3. Nonetheless, these models can be convexified (§3) and trained via SVD, if we remove nonlinearities between layers and swap the cross-entropy objective with Frobenius norm minimization. Undoubtedly, such linearization incurs a drop of accuracy on empirical test performance. Nonetheless, we show that the (convexified) model’s parameters learned by SVD can provide initialization to deeper (non-linear) models, which then can be fine-tuned on cross-entropy objectives. The non-linear models are endowed with our novel Split-ReLu layer, which has twice as many parameters as a ReLu fully-connected layer, and behaves as a linear layer when its parameters reside on some hyperplane (§5.2). Training on modest hardware (e.g., laptop) is sufficient for this learning pipeline (convexify → SVD → fine-tune) yet it trains much faster than current approaches, that are commonly trained on expensive hardware. We summarize our contributions as:
1. We open-source a flexible python software library that allows symbolic definition of matrices and computes their SVD without explicitly calculating them.
2. We linearize popular GRL models, and train them via SVD of design matrices.
3. We show that fine-tuning a few parameters on-top of the SVD initialization sets state-of-the-art on many GRL tasks while, overall, training orders-of-magnitudes faster.

heading: 2 Preliminaries & notation
text: We denote a graph with n nodes and m edges with an adjacency matrix A ∈ Rn×n and additionally, if nodes have (d-dimensional) features, with a feature matrix X ∈ Rn×d. If nodes i, j ∈ [n] are connected then Aij is set to their edge weight and otherwise Aij = 0. Further, denote the (row-wise normalized) transition matrix as T = D−1A and denote the symmetrically normalized adjacency with self-connections as Â = (D + I)− 1 2 (A + I)(D + I)− 1 2 where I is identity matrix.
We review model classes: (1) network embedding and (2) message passing that we define as follows. The first inputs a graph (A, X) and outputs node embedding matrix Z ∈ Rn×z with z-dimensions per node. Z is then used for an upstream task, e.g., link prediction. The second class utilizes a function H : Rn×n×Rn×d → Rn×z where the function H(A,X) is usually directly trained on the upstream task, e.g., node classification. In general, the first class is transductive while the second is inductive.

heading: 2.1 Network embedding models based on DeepWalk & node2vec
text: The seminal work of DeepWalk [Perozzi et al., 2014] embeds nodes of a network using a two-step process: (i) simulate random walks on the graph – each walk generating a sequence of node IDs then (ii) pass the walks (node IDs) to a language word embedding algorithm, e.g. word2vec [Mikolov et al., 2013], as-if each walk is a sentence. This work was extended by node2vec [Grover and Leskovec, 2016] among others. It has been shown by Abu-El-Haija et al. [2018] that the learning outcome of
2TensorFlow can caluclate matrix-free SVD if one implements a LinearOperator, as such, our code could be re-implemented as a routine that can convert TensorGraph to LinearOperator.
3Practitioners rarely verify that∇θJ = 0, where J is mean train objective and θ are model parameters.
the two-step process of DeepWalk is equivalent, in expectation, to optimizing a single objective4:
min Z={L,R} ∑ (i,j)∈[n]×[n] [ − E q∼Q [T q] ◦ log σ(LR>)− λ(1−A) ◦ log(1− σ(LR>)) ] ij , (1)
where L,R ∈ Rn× z2 are named by word2vec as the input and output embedding matrices, ◦ is Hadamard product, and the log(.) and the standard logistic σ(.) = (1 + exp(.))−1 are applied element-wise. The objective above is weighted cross-entropy where the (left) positive term weighs the dot-product L>i Rj by the (expected) number of random walks simulated from i and passing through j, and the (right) negative term weighs non-edges (1−A) by scalar λ ∈ R+. The context distribution Q stems from step (ii) of the process. In particular, word2vec accepts hyperparameter context window size C for its stochasatic sampling: when it samples a center token (node ID), it then samples its context tokens that are up-to distance c from the center. The integer c is sampled from a coin flip uniform on the integers [1, 2, . . . , C] – as detailed by Sec.3.1 of [Levy et al., 2015]. Therefore, PQ(q | C) ∝ C−q+1C . Since q has support on [C], then PQ(q | C) = ( 2 (C+1)C ) C−q+1 C .

heading: 2.2 Message passing graph networks for (semi-)supervised node classification
text: We are also interested in a class of (message passing) graph network models taking the general form:
for l = 0, 1, . . . L: H(l+1) = σl ( g(A)H(l)W(l) ) ; H(0) = X; H = H(L); (2)
where L is the number of layers, W(l)’s are trainable parameters, σl’s denote element-wise activations (e.g. logistic or ReLu), and g is some (possibly trainable) transformation of adjacency matrix. GCN [Kipf and Welling, 2017] set g(A) = Â, GAT [Veličković et al., 2018] set g(A) = A ◦ MultiHeadedAttention and GIN [Xu et al., 2019] as g(A) = A + (1 + )I with > 0. For node classification, it is common to set σL = softmax (applied row-wise), specify the size of WL s.t. H ∈ Rn×y where y is number of classes, and optimize cross-entropy objective: . min{Wj}Lj=1 [−Y ◦ logH− (1−Y) ◦ log(1−H)] , where Y is a binary matrix with onehot rows indicating node labels. In semi-supervised settings where not all nodes are labeled, before measuring the objective, subset of rows can be kept in Y and H that correspond to labeled nodes.

heading: 2.3 Truncated Singular Value Decomposition (SVD)
text: SVD is an algorithm that approximates any matrix M ∈ Rr×c as a product of three matrices:
SVDk(M) , arg min U,S,V
||M−USV>||F subject to U>U = V>V = Ik; S = diag(s1, . . . , sk).
The orthonormal matrices U ∈ Rr×k and V ∈ Rc×k, respectively, are known as the left- and right-singular bases. The values along diagonal matrix S ∈ Rk×k are known as the singular values. Due to theorem of Eckart and Young [1936], SVD recovers the best rank-k approximation of input M, as measured by the Frobenius norm ||.||F. Further, if k ≥ rank(M)⇒ ||.||F = 0. Popular SVD implementations follow Random Matrix Theory algorithm of Halko et al. [2009]. The prototype algorithm starts with a random matrix and repeatedly multiplies it by M and by M>, interleaving these multiplications with orthonormalization. Our SVD implementation (in Appendix) also follows the prototype of [Halko et al., 2009], but with two modifications: (i) we replace the recommended orthonormalization step from QR decomposition to Cholesky decomposition, giving us significant computational speedups and (ii) our implementation accepts symbolic representation of M (§4), in lieu of its explicit value (constrast to TensorFlow and PyTorch, requiring explicit M).
In §3, we derive linear first-order approximations of models reviewed in §2.1 & §2.2 and explain how SVD can train them. In §5, we show how they can be used as initializations of non-linear models.
4Derivation is in [Abu-El-Haija et al., 2018]. Unfortunately, matrix in Eq. 1 is dense with O(n2) nonzeros.

heading: 3 Convex first-order approximations of GRL models

heading: 3.1 Convexification of Network Embedding Models
text: We can interpret objective 1 as self-supervised learning, since node labels are absent. Specifically, given a node i ∈ [n], the task is to predict its neighborhood as weighted by the row vector Eq[T q]i, representing the subgraph5 around i. Another interpretation is that Eq. 1 is a decomposition objective: multiplying the tall-and-thin matrices, as LR> ∈ Rn×n, should give a larger value at (LR>)ij = L>j Ri when nodes i and j are well-connected but a lower value when (i, j) is not an edge. We propose a matrix such that its decomposition can incorporate the above interpretations:
M̂(NE) = Eq|C [T q]− λ(1−A) = (
2
(C + 1)C ) C∑ q=1 ( C − q + 1 C ) T q − λ(1−A) (3)
If nodes i, j are nearby, share a lot of connections, and/or in the same community, then entry M̂(NE)ij should be positive. If they are far apart, then M̂(NE)ij = −λ. To embed the nodes onto a low-rank space that approximates this information, one can decompose M̂(NE) into two thin matrices (L,R):
LR> ≈ M̂(NE) ⇐⇒ (LR>)i,j = 〈Li,Rj〉 ≈ M̂(NE)ij for all i, j ∈ [n]. (4)
SVD gives low-rank approximations that minimize the Frobenius norm of error (§2.3). The remaining challenge is computational burden: the right term (1−A), a.k.a, graph compliment, has ≈ n2 nonzero entries and the left term has non-zero at entry (i, j) if nodes i, j are within distance C away, as q has support on [C] – for reference Facebook network has an average distance of 4 [Backstrom et al., 2012] i.e. yielding T 4 with O(n2) nonzero entries – Nonetheless, Section §4 presents a framework for decomposing M̂ from its symbolic representation, without explicitly computing its entries. Before moving forward, we note that one can replace T in Eq. 3 by its symmetrically normalized counterpart Â, recovering a basis where L = R. This symmetric modeling might be emperically preferred for undirected graphs. Learning can be performed via SVD. Specifically, the node at the ith row and the node at the jthth column will be embedded, respectively, in Li and Rj computed as:
U,S,V← SVDk(M̂(NE)); L← US 1 2 ; R← VS 12 (5)
In this k-dim space of rows and columns, Euclidean measures are plausible: Inference of nodes’ similarity at row i and column j can be modeled as f(i, j) = 〈Li,Rj〉 = U>i SVj , 〈Ui,Vj〉S .

heading: 3.2 Convexification of message passing graph networks
text: Removing all σl’s from Eq. 2 and setting g(A) = Â gives outputs of layers 1, 2, and L, respectively,
as: ÂXW(1) , Â2XW(1)W(2) , and ÂLXW(1)W(2) . . .W(L). (6)
Without non-linearities, adjacent parameter matrices can be absorbed into one another. Further, the model output can concatenate all layers, like JKNets [Xu et al., 2018], giving final model output of:
H(NC)linearized = [ X ÂX Â2X . . . ÂLX ] Ŵ , M̂(NC)Ŵ, (7)
where the linearized model implicitly constructs design matrix M̂(NC) ∈ Rn×F and multiplies it with parameter Ŵ ∈ RF×y – here, F = d + dL. Crafting design matrices is a creative process (§5.3). Learning can be performed by minimizing the Frobenius norm: ||H(NC)−Y||F = ||M̂(NC)Ŵ−Y||F. Moore-Penrose Inverse (a.k.a, the psuedoinverse) provides one such minimizer:
Ŵ∗ = argmin Ŵ ∣∣∣∣∣∣M̂Ŵ −Y∣∣∣∣∣∣ F = M̂†Y ≈ VS+U>Y, (8)
with U,S,V ← SVDk(M̂). Notation S+ reciprocates non-zero entries of diagonal S [Golub and Loan, 1996]. Multiplications in the right-most term should, for efficiency, be executed right-to-left.
5Eq[T q]i is a distribution on [n]: entry j equals prob. of walk starting at i ending at j if walk length∼ U [C].
The pseudoinverse M̂† ≈ VS+U> recovers the Ŵ∗ with least norm (§6, Theorem 1). The ≈ becomes = when k ≥ rank(M̂). In semi-supervised settings, one can take rows subset of either (i) Y and U, or of (ii) Y and M, keeping only rows that correspond to labeled nodes. Option (i) is supported by existing frameworks (e.g., tf.gather()) and our symbolic framework (§4) supports (ii) by implicit row (or column) gather – i.e., calculating SVD of submatrix of M without explicitly computing M nor the submatrix. Inference over a (possibly new) graph (A,X) can be calculated by (i) (implicitly) creating the design matrix M̂ corresponding to (A,X) then (ii) multiplying by the explicitly calculated Ŵ∗. As explained in §4, M̂ need not to be explicitly calculated for computing multiplications.

heading: 4 Symbolic matrix representation
text: To compute the SVD of any matrix M using algorithm prototypes presented by Halko et al. [2009], it suffices to provide functions that can multiply arbitrary vectors with M and M>, and explicit calculation of M is not required. Our software framework can symbolically represent M as a directed acyclic graph (DAG) of computations. On this DAG, each node can be one of two kinds:
1. Leaf node (no incoming edges) that explicitly holds a matrix. Multiplications against leaf nodes are directly executed via an underlying math framework (we utilize TensorFlow).
2. Symbolic node that only implicitly represents a matrix as as a function of other DAG nodes. Multiplications are recursively computed, traversing incoming edges, until leaf nodes.
For instance, suppose leaf DAG nodes M1 and M2, respectively, explicitly contain row vector ∈ R1×n and column vector ∈ Rn×1. Then, their (symbolic) product DAG node M = M2@M1 is ∈ Rn×n. Although storing M explicitly requiresO(n2) space, multiplications against M can remain within O(n) space if efficiently implemented as 〈M, .〉 = 〈M2, 〈M1, .〉〉. Figure 1 shows code snippet for composing DAG to represent symbolic node M̂(NE) (Eq. 3), from leaf nodes initialized with in-memory matrices. Appendex lists symbolic nodes and their implementations.

heading: 5 SVD initialization for deeper models fine-tuned via cross-entropy

heading: 5.1 Edge function for network embedding as a (1-dimensional) Gaussian kernel
text: SVD provides decent solutions to link prediction tasks. Computing U,S,V ← SVD(M(NE)) is much faster than training SOTA models for link prediction, yet, simple edge-scoring function f(i, j) = 〈Ui,Vj〉S yields competitive empirical (test) performance. We propose f with θ = {µ, s}:
fµ,s(i, j) = Ex∼N (µ,s)〈Ui,Vj〉Sx = U > i Ex [Sx]Vj = U>i (∫ Ω SxN (x | µ, s) dx ) Vj , (9)
where N is the truncated normal distribution (we truncate to Ω = [0.5, 2]). The integral can be approximated by discretization and applying softmax (see §A.4). The parameters µ ∈ R, σ ∈ R>0 can be optimized on cross-entropy objective for link-prediction:
minµ,s − E(i,j)∈A [log (σ(fµ,s(i, j)))]− k(n)E(i,j)/∈A [log (1− σ(fµ,s(i, j)))] , (10)
where the left- and right-terms, respectively, encourage f to score high for edges, and the low for non-edges. k(n) ∈ N>0 controls the ratio of negatives to positives per batch (we use k(n) = 10). If the optimization sets µ = 1 and s ≈ 0, then f reduces to no-op. In fact, we initialize it as such, and we observe that f converges within one epoch, on graphs we experimented on. If it converges as µ < 1, VS µ > 1, respectively, then f would effectively squash, VS enlarge, the spectral gap.

heading: 5.2 Split-ReLu (deep) graph network for node classification (NC)
text: Ŵ∗ from SVD (Eq.8) can initialize an L-layer graph network with input: H(0) = X, with:
message passing (MP) H(l+1) = [ ÂH(l)W
(l) (p) ] + − [ ÂH(l)W (l) (n) ] + , (11)
output H = l=L∑ l=0 [ H(l)W (l) (op) ] + − [ H(l)W (l) (on) ] +
(12)
initialize MP W(l)(p) ← I;W (l) (n) ← −I; (13)
initialize output W(l)(op) ← Ŵ ∗ [dl : d(l+1)];W (l) (on) ← −Ŵ ∗ [dl : d(l+1)]; (14)
Element-wise [.]+ = max(0, .). Further, W[i : j] denotes rows from (i)th until (j−1)th of W. The deep network layers (Eq. 11&12) use our Split-ReLu layer which we formalize as:
SplitReLu(X;W(p),W(n)) = [XW(p)]+ − [XW(n)]+ , (15)
where the subtraction is calculated entry-wise. The layer has twice as many parameters as standard fully-connected (FC) ReLu layer. In fact, learning algorithms can recover FC ReLu from SplitReLu by assigning W(n) = 0. More importantly, the layer behaves as linear in X when W(p) = −W(n). On this hyperplane, this linear behavior allows us to establish the equivalency: the (non-linear) model H is equivalent to the linear Hlinearized at initialization (Eq. 13&14) due to Theorem 2. Following the initialization, model can be fine-tuned on cross-entropy objective as in §2.2.

heading: 5.3 Creative Add-ons for node classification (NC) models
text: Label re-use (LR): Let M̂(NC)LR , [ M̂ (NC) (Â− (D + I)−1)Y[train] (Â− (D + I)−1)2Y[train] ] . This follows the motivation of Wang and Leskovec [2020], Huang et al. [2021], Wang [2021] and their empirical results on ogbn-arxiv dataset, where Y[train] ∈ Rn×y contains one-hot vectors at rows corresponding to labeled nodes but contain zero vectors for unlabeled (test) nodes. Our scheme is similar to concatenating Y[train] into X, but with care to prevent label leakage from row i of Y to row i of M̂, as we zero-out the diagonal of the adjacency multiplied by Y[train].
Pseudo-Dropout (PD): Dropout [Srivastava et al., 2014] reduces overfitting of models. It can be related to data augmentation, as each example is presented multiple times. At each time, it appears with a different set of dropped-out features – input or latent feature values, chosen at random, get replaced with zeros. As such, we can replicate the design matrix as: M̂> ← [ M̂> PD(M̂)> ] . This
row-wise concatenation maintains the width of M̂ and therefore the number of model parameters.
In the above add-ons, concatenations, as well as PD, can be implicit or explicit (see §A.3).

heading: 6 Analysis & Discussion
text: Theorem 1. (Min. Norm) If system M̂Ŵ = Y is underdetermined6 with rows of M̂ being linearly independent, then solution space Ŵ∗ = { Ŵ ∣∣∣ M̂Ŵ = Y} has infinitely many solutions. Then, for
k ≥ rank(M̂), matrix Ŵ∗, recovered by Eq.8 satisfies: Ŵ∗ = argmin Ŵ∈Ŵ∗ ||Ŵ|| 2 F .
Theorem 1 implies that, even though one can design a wide M̂(NC) (Eq.7), i.e., with many layers, the recovered parameters with least norm should be less prone to overfitting. Recall that this is the goal of L2 regularization. Analysis and proofs are in the Appendix.
Theorem 2. (Non-linear init) The initialization Eq. 13&14 yields H(NC)linearized = H ∣∣ θ←via Eq. 13&14 .
Theorem 2 implies that the deep (nonlinear) model is the same as the linear model, at the initialization of θ (per Eq. 13&14, using Ŵ∗ as Eq. 8). Cross-entropy objective can then fine-tune θ.
This end-to-end process, of (i) computing SVD bases and (ii) training the network fθ on singular values, advances SOTA on competitve benchmarks, with (i) converging (quickly) to a unique solution and (ii) containing merely a few parameters θ – see §7.

heading: 7 Applications & Experiments
text: We download and experiment on 9 datasets summarized in Table 1.
We attempt link prediction (LP) tasks on smaller graph datasets (< 1 million edges) of: Protein-Protein Interactions (PPI) graph from Grover and Leskovec [2016]; as well as ego-Facebook (FB), AstroPh, HepTh from Stanford SNAP [Leskovec and Krevl, 2014]. For these datasets, we use the train-test splits of Abu-El-Haija et al. [2018]. We also attempt semi-supervised node classification (SSC) tasks on smaller graphs of Cora, Citeseer, Pubmed, all obtained from Planetoid [Yang et al., 2016]. For these smaller datasets, we only train and test using the SVD basis (without finetuning)
Further, we attempt on slightly-larger datasets (> 1 million edges) from Stanford’s Open Graph Benchmark [OGB, Hu et al., 2020]. We use the official train-test-validation splits and evaluator of OGB. We attempt LP and SSC, respectively, on Drug Drug Interactions (ogbl-DDI) and ArXiv citation network (ogbn-ArXiv). For these larger datasets, we use the SVD basis as an initialization that we finetune, as described in §5. For time comparisons, we train all models on Tesla K80.

heading: 7.1 Test performance & runtime on smaller datasets from: Planetoid & Stanford SNAP
text: For SSC over Planetoid’s datasets, both A and X are given. Additionally, only a handful of nodes are labeled. The goal is to classify the unlabeled test nodes. Table 2 summarizes the results. For baselines, we download code of GAT [Veličković et al., 2018], MixHop [Abu-El-Haija et al., 2019], GCNII [Chen et al., 2020] and re-ran them with instrumentation to record training time. However, for baselines Planetoid [Yang et al., 2016] and GCN [Kipf and Welling, 2017], we copied numbers
6E.g., if the number of labeled examples i.e. height of M and Y is smaller than the width of M.
from [Kipf and Welling, 2017]. For our models, the row labeled iSVD100(M̂(NC)), we run our implicit SVD twice per graph. The first run incorporates structural information: we (implicitly) construct M̂(NE) with λ = 0.05 and C = 3, then obtain L,R← SVD64(M̂(NE)), per Eq. 5. Then, we concatenate L and R into X. Then, we PCA the resulting matrix to 1000 dimensions, which forms our new X. The second SVD run is to train the classification model parameters Ŵ∗. From the PCA-ed X, we construct the implicit matrix M̂(NC) with L = 15 layers and obtain Ŵ∗ = VS+U>Y[train] with U,S,V ← SVD100(M̂(NC)[train]), per in RHS of Eq. 8. For our second “+ dropout” model variant, we (implicit) augment the data by M̂(NC) > ← [ M̂(NC) > | PD(M̂(NC))> ] , update indices
[train] ← [ train | train ]> then similarly learn as: SVD100(M̂ (NC) [train]) → Ŵ∗ = VS+U>Y[train] Discussion: Our method is competitive yet trains faster than SOTA. In fact, the only method that reliably beats ours on all dataset is GCNII, but its training time is about one-thousand-times longer.
For LP over SNAP and node2vec datasets, the training adjacency A is given but not X. The split includes test positive and negative edges, which are used to measure a ranking metric: ROC-AUC, where the metric increases when test positive edges are ranked above than negatives. Table 3 summarizes the results. For baselines, we download code of WYS [Abu-El-Haija et al., 2018]; we use the efficient implementation of PyTorch-Geometric [Fey and Lenssen, 2019] for node2vec (n2v) and we download code of Qiu et al. [2018] and run it with their two variants, denoting their first variant as NetMF (for exact, explicitly computing the design matrix) as and their second variant as NetM̃F (for approximate, sampling matrix entry-wise) – their code runs SVD after computing either variant. For the first of our models, we compute U,S,V← SVD32(M̂(NE)) and score every test edge as U>i SVj . For the second, we first run SVD256 on half of the training edges, determine the “best” rank ∈ {8, 16, 32, 128, 256} by measuring the AUC on the remaining half of training edges, then using this best rank, recompute the SVD on the entire training set, then finally score the test edges. Discussion: Our method is competitive on SOTA while training much faster. Both NetMF and WYS explicitly calculate a dense matrix before factorization. On the other hand, NetM̃F approximates the matrix entry-wise, trading accuracy for training time. In our case, we have the best
of both worlds: using our symbolic representation and SVD implementation, we can decompose the design matrix while only implicitly representing it, as good as if we had explicitly calculated it.

heading: 7.2 Experiments on Stanford’s OGB datasets
text: We summarize experiments ogbl-DDI and ogbn-ArXiv, respectively, in Tables 4 and 5. For baselines, we copy numbers from the public leaderboard, where the competition is fierce. We then follow links on the leaderboard to download author’s code, that we re-run, to measure the training time. For our models on ogbl-DDI, we (a) first calculate U,S,V← SVD100(M̂(NE)) built only from training edges and score test edge (i, j) using U>i SVj . Then, we (b) then finetune fµ,s (per §5.1, Eq. 9 & 10)) for only a single epoch and score using fµ,s(i, j). Then, we (c) update the SVD basis to include edges from the validation partition and also score using fµ,s(i, j). We report the results for the three steps. For the last step, the rules of OGB allows using the validation set for training, but only after the hyperparameters are finalized. The SVD has no hyperparameters (except the rank, which was already determined by the first step). More importantly, this simulates a realistic situation: it is cheaper to obtain SVD (of an implicit matrix) than back-propagate through a model. For a time-evolving graph, one could run the SVD more often than doing gradient-descent epochs on a model. For our models on ogbn-ArXiv, we (a) compute U,S,V← SVD250(M̂(NC)LR ) where the implicit matrix is defined in §5.3. We (b) repeat this process where we replicate M̂(NC)LR once: in the second replica, we replace the Y[train] matrix with zeros (as-if, we drop-out the label with 50% probability). We (c) repeat the process where we concatenate two replicas of M̂(NC)LR into the design matrix, each with different dropout seed. We (d) fine-tune the last model over 15 epochs using stochastic GTTF [Markowitz et al., 2021]. Discussion: Our method competes or sets SOTA, while training much faster. Time improvements: We replaced the recommended orthonormalization of [Halko et al., 2009] from QR decomposition, to Cholesky decomposition (§A.2). Further, we implemented caching to avoid computing sub-expressions if already calculated (§A.3.4). Speed-ups are shown in Fig. 2.
Table 5: Test classification accuracy over ogbn-arxiv.
Graph dataset: ogbn-ArXiv Baselines: accuracy tr.time
GAT+LR+KD [Ren, 2020] 74.16 ±0.08 (6h) GAT+LR [Niu, 2020] 73.99 ±0.12 (3h)
AGDN [Sun and Wu, 2020] 73.98 ±0.09 (50m) GAT+C&S [Huang et al., 2021] 73.86 ±0.14 (2h)
GCNII [Chen et al., 2020] 72.74 ±0.16 (3h)
Our models: accuracy tr.time
(a) iSVD250(M̂ (NC) LR ) (§3.2, Eq. 8) 68.90 ±0.02 (1s) (b) + dropout(LR) (§5.3) 69.34 ±0.02 (3s) (c) + dropout(M̂(NC)LR ) (§5.3) 71.95 ±0.03 (6s) (d) + finetune H (§5.2, Eq. 11-14) 74.14 ±0.05 (2m)
DDI
FB
AstroPh
PPI
0.00 0.25 0.50 0.75 1.00
(☓; QR) (☓; Cholesky) (✓; Cholesky)
Figure 2: SVD runtime configs of (lazy caching; orthonormalization) as a ratio of SVD’s common default (QR decomposition)

heading: 8 Related work
text: Applications: SVD was used to project rows & columns of matrix M onto an embedding space. M can be the Laplacian of a homogenous graph [Belkin and Niyogi, 2003], Adjacency of user-item bipartite graph [Koren et al., 2009], or stats [Deerwester et al., 1990, Levy and Goldberg, 2014] for word-to-document. We differ: our M is a function of (leaf) matrices – useful when M is expensive to store (e.g., quadratic). While Qiu et al. [2018] circumvents this by entry-wise sampling the (otherwise n2-dense) M, our SVD implementation could decompose exactly M without calculating it.
Symbolic Software Frameworks: including Theano [Al-Rfou et al., 2016], TensorFlow [Abadi et al., 2016] and PyTorch [Paszke et al., 2019], allow chaining operations to compose a computation (directed acyclic) graph (DAG). They can efficiently run the DAG upwards by evaluating (all entries of) matrix M ∈ Rr×c at any DAG node. Our DAG differs: instead of calculating M, it provides product function uM(.) = M× . – The graph is run downwards (reverse direction of edges). Matrix-free SVD: For many matrices of interest, multiplying against M is computationally cheaper than explicitly storing M entry-wise. As such, many researchers implement SVD(uM), e.g. Calvetti et al. [1994], Wu and Simon [2000], Bose et al. [2019]. We differ in the programming flexibility: the earlier methods expect the practitioner to directly implement uM. On the other hand, our framework allows composition of M via operations native to the practitioner (e.g., @, +, concat), and uM is automatically defined.
Fast Graph Learning: We have applied our framework for fast graph learning, but so as samplingbased approaches including [Chen et al., 2018, Chiang et al., 2019, Zeng et al., 2020, Markowitz et al., 2021]. We differ in that ours can be used to obtain an initial closed-form solution (very quickly) and can be fine-tuned afterwards using any of the aforementioned approaches. Additionally, Graphs shows one general application. Our framework might be useful in other areas utilizing SVD.

heading: 9 Conclusion, our limitations & possible negative societal impact
text: We develop a software framework for symbolically representing matrices and compute their SVD. Without computing gradients, this trains convexified models over GRL tasks, showing empirical metrics competitive with SOTA while training significantly faster. Further, convexified model parameters can initialize (deeper) neural networks that can be fine-tuned with cross entropy. Practitioners adopting our framework would now spend more effort in crafting design matrices instead of running experiments or tuning hyperparameters of the learning algorithm. We hope our framework makes high performance graph modeling more accessible by reducing reliance on energy-intensive computation.
Limitations of our work: From a representational prospective, our space of functions is smaller than TensorFlow’s, as our DAG must be a linear transformation of its input matrices, e.g., unable to encode element-wise transformations and hence demanding first-order linear approximation of models. As a result, our gradient-free learning can be performed using SVD. Further, our framework only works when the leaf nodes (e.g., sparse adjacency matrix) fit in memory of one machine.
Possible societal impacts: First, our convexified models directly learn parameters in the feature space i.e. they are more explainable than deeper counterparts. Explainability is a double-edged sword: it gives better ability to interpret the model’s behavior, but also allows for malicious users, e.g., to craft attacks on the system (if they can replicate the model parameters). Further, we apply our method on graphs. It is possible to train our models to detect sensitive attributes of social networks (e.g., ethnicity). However, such ethical concerns exists with any modeling technique over graphs.
Acknowledgments and Disclosure of Funding
This material is based upon work supported by the Defense Advanced Research Projects Agency (DARPA) and the Army Contracting Command- Aberdeen Proving Grounds (ACC-APG) under Contract Number W911NF-18-C-0020.

Review summaries:
The authors suggest to express large matrices as operations over given manageable leaf matrices. They keep track of the resulting abstract syntax tree (AST) and offer a specialised SVD decomposition implementation that can perform the computation on the final expression by traversing the AST and operating on the individual leaf matrices. They then show how using the SVD operation they can perform learning in various inductive and transductive tasks. The advantage of their proposal is that the process is computationally efficient and they can in practice match SOTA results with speedups reaching several orders of magnitude.

The goal of this paper is to speed-up graph representational learning by removing  nonlinearities between layers and replacing them with linear models while also  replacing cross entropy by Frobenius norm minimization. In turn, this allows the  application of SVD, and the authors suggest modifications of randomized SVD  which execute faster and in a matrix-free manner. The above approach is reported  to be much faster than previous state-of-the-art, even though it is required to post-process  the obtained model by feeding it as an initial guess to a more expensive, but also  more accurate, approach. 

In this work, authors design a software framework which allows symbolically representing matrices and computing their SVD implicitly. This paper also derives linear first-order approximations of popular Graph Representation Learning (GRL) models, and they can be trained via SVD. Further, SVD initializations can be used to fine-tune deeper neural networks via cross-entropy on many GRL tasks. Experimental results demonstrate that the proposed model trains orders-of-magnitudes faster than state-of-art methods, while having comparable performance.