Laplacian-Optimized Diffusion For Semi-Supervised Learning

1.2. Graph-based Approaches to SSLFoundations. A commonly-used assumption for input data stored in a meaningful embedding space is that two datapoints are likely to share the same label if they are “close” to each other, or if they lie on a samelow-dimensional sub-manifold (i.e., if points of a same class form low dimensional subspaces locally, which can be seen as tangent spaces of a low-dimensional manifold in the high-dimensionalembedding space); see, e.g., Zhu and Goldberg (2009), Chapelle et al. (2010) and inset. Thesesimple geometric assumptions, that we will respectively refer to as the cluster assumption and the manifold assumptionas various authors have done in the past (see Sindhwani et al. (2010)), are at the root of a large family of “graph-basedmethods”, including Zhu et al. (2003); Zhou et al. (2003); Wang et al. (2006); Karasuyama and Mamitsuka (2017),which are considered state-of-the-art methods in graph-based SSL.Methodology. Graph-based methods begin with the creation of a sparse graph whose vertices are input (labeled andunlabeled) data points. Typically, two data points are connected with an edge if either the distance between them isless than some threshold in the embedding space or if one of them is in the set of the K nearest neighbors of the other.This graph thus connects nearby points, and offers a discrete representation of the “geometry” of the data, whichincludes any sub-manifold around which the data concentrate. Each edge of this graph is then assigned a weightbased on either its length in the embedding space as in Belkin and Niyogi (2004) (sometimes using a data-dependentrescaling of the coordinates, see Karasuyama and Mamitsuka (2017)), or on the distribution of the nearest data pointsaround as in Wang et al. (2006). Equipped with this weighted graph, classification is then performed typically throughpropagation of the known labels to the rest of the graph, using diffusion based on the edge weights or through someform of entropy minimization.Strengths and limitations. One advantage of semi-supervised learning compared to deep learning is that the graph isentirely defined by the data: it therefore dramatically reduces the number of hyperparameters like number of layers,neuron count, connectivity type, etc. While SSL involves a weighted graph over which classification is performed,an obvious limitation is that the edge weights are rarely globally trained (i.e., optimized): they are usually chosento be Gaussian weights à la Belkin and Niyogi (2004) or Locally Linear Embedding (LLE) weights à la Roweis andSaul (2000). Since these “pre-canned” weights strongly influence the behavior of the classifier as they indicate thesimilarity between data points, graph-based approaches to SSL also depend critically on how well the input data fitsthe underlying geometric assumptions in the actual embedding (feature) space. The few methods that try to partiallyoptimize weights based on the input data do rely on nonlinear, nonconvex minimization, which hurts their scalability inpractice. Note that many methods have built over this graph-based diffusion approach, introducing various additionalterms such as iterated Laplacians (Zhou and Belkin (2011)), Tikhonov (as in Belkin et al. (2006)) or 1 (as in Rustamovand Klosowski (2018)) regularization, and diffusion-inspired convolutional filters, see Kipf and Welling (2017).1.3. ContributionsIn this paper, we approach the design of a new Laplacian-optimized diffusion for semi-supervised learning using discrete differential geometry. To exploit the intrinsic structure of the data which is often spatially distributed in clustersand/or on low-dimensional manifolds within the high-dimensional representation space, we propose a biconvex minimization where the edge weights of a data-dependent graph are optimized in alternation with an inference of labelsover the unlabeled data points. This edge weight optimization constructs a Laplacian operator whose quadratic formrepresents a smoothness energy in a metric adapted to the intrinsic distribution of the input in the embedding space.The resulting Laplacian-optimized graph diffusion of known labels exploits the geometry of the input, classifying unlabeled data more accurately than current SSL techniques. The algorithmic simplicity and efficiency of our approachalso make it attractive, despite the seemingly complex task of optimizing all the graph weights.2. Review of Graph-based SSL MethodsBefore introducing our approach, we first delve into the existing graph-based methods to gather some of the keyconcepts we will leverage in our work, but also to point out how usual differential geometry notions (from Laplacianand Dirichlet energy to optimal transport) play a central role throughout this literature.3

2.1. NotationsSemi-supervised learning considers an input dataset D {xi }i 1.n of n points distributed in the D-dimensional Euclidean space D . The first input points are labeled, so that the input set can be written as D DL DU , whereDL {xi }i 1. are the labeled points and DU {xi }i 1.n the unlabeled points, with n. Assuming there areC categories of labels (defined to be the integers {1, 2, .,C} for simplicity), we store label information in an n Cmatrix P, such that 1if xi DL and label(xi ) k,Pik 0otherwise.Note that the Euclidean norm of vectors will be denoted by k.k, while the Frobenius norm of a matrix will be represented as k.kF , with kMk2F Tr(MMt). We will also denote by diag(.) the operator that turns a vector v of size k intoa diagonal k k matrix diag(v) such that [diag(v)]ii [v]i .All graph-based methods proceed quite similarly. They begin by defining a connected graph G (V , E ) with thenodes V corresponding to the n input data points, and where an edge ei j E connecting points xi and x j is creatediff x j is among the K nearest neighbors of xi or vice-versa. With this data representation, they proceed to propagateknown labels through the graph. This label propagation is achieved through various means designed to exploit thegeometric assumptions of the data, as we review next.2.2. Harmonic Gaussian FieldsHarmonic Gaussian fields (HGF, Zhu et al. (2003)), one of the early approaches in graph-based SSL, proposed to usea discrete Laplacian operator over the graph to induce diffusion of the labels. For every graph edge ei j E , they assigna positive Gaussian weight wi j of the form: wi j exp kxi x j k22 /σ 2 ,(1)where σ is a length scale hyperparameter, typically picked based on the average distance between data points connected by an edge. Since these weights are known to give rise to a discrete approximation of the Laplace operator,see Belkin and Niyogi (2004), the inherent smoothness of the data can be leveraged by forming a harmonic interpolation of the labels based on these weight expressions. More precisely, let’s define the matrix Q of size n C, where Qidenotes its i-th row with C values, representing the likelihood scores for node xi to be in each of the C classes. TheHGF approach computes matrix Q as the minimizer of the following discrete Dirichlet energy:(2)arg min wi j kQi Q j k2 s.t. Qi Pi xi DL .Qei j EThis minimization calculates C harmonic functions (one per class, stored in columns of Q) subject to label constraints;the inferred label of a given unlabeled node is then defined as the index of the largest element in Q’s i-th row:label(xi ) arg max Qik .(3)k 1.CThis approach can be interpreted as a Gaussian regression, since it corresponds to the conditional expectation of labelassignment assuming that the latter is a realization of a Gaussian process with a covariance function given by theLaplace operator; it also shares obvious connections to the graph mincut approach of Blum and Chawla (2001). Notethat the optimal Q is easily solved for via a sparse linear system, and by definition of Eq. (2), this harmonic labelpropagation is interpolatory, i.e., known input labels are kept unaltered. Prior knowledge of class sizes can also beaccommodated through post-processing using class mass normalization as in de Sousa (2015).2.3. Learning with Local and Global ConsistencyLearning with Local and Global Consistency method (LLGC, Zhou et al. (2003)) relies on the same Gaussian edgeweights, but modifies the variational problem for Q to mimic spreading activation networks: first, they replace theinterpolation of P with a fitting penalty, which only favors the initial labeling (in lieu of enforcing it exactly) to resultin a globally smoother assignment; second, they propose to rescale the value of Q at a node xi in the Dirichlet energyby dividing by the sum of the edge weights emanating from that node in order to compensate for unavoidable densityvariation in real datasets. Classification is therefore achieved by solving:Q(4)arg min kQ Pk2F η wi j k Qi j k2 .Qei j E4 k wik k w jk

Note that this normalization of the Dirichlet energy corresponds precisely to what is called the “normalized Laplacian”quadratic form (Lnorm ) in Normalized Cuts of Shi and Malik (2000), a well-known approach to spectral clustering. Thetrade-off between the two competing energy terms in Eq. (4) is captured by a positive hyperparameter η. Finally,minimizing the resulting energy turns out to be an implicit integration step of diffusion as in Desbrun et al. (1999);indeed, differentiating the objective function shows that optimization of Eq. (4) amounts to solving:(Id ηLnorm ) Q P.(5)Unlike HGF, LLGC does not guarantee label preservation, and mislabeling happens in practice for most choices ofsize η. Still, the normalization of the Dirichlet energy allows for better results overall compared to previous methods.2.4. Linear Neighborhood PropagationThe Linear Neighborhood Propagation (LNP, Wang and Zhang (2008)) departs from previous methods by rejectingthe use of Gaussian weights since the Laplacian is certainly not the only operator to capture smoothness of labelassignments. Instead, they take a cue from work on dimensionality reduction (namely, Locally Linear Embedding(LLE), Roweis and Saul (2000)), and propose to locally optimize edge weights to exploit the local linearity assumption. For a given node xi , weights of emanating edges ei j are found through constrained quadratic optimization:arg min kxi wi j x j k2 s.t. wi j 0 and{wi j } j N(i) wi j 1,(6)j N(i)j N(i)where N (i) represents the indices of the vertices connected to xi in the graph G . Since this optimization is performedindependently for each node, note that wi j 6 w ji , and we should consider the graph as directed in this case. Theinterpretation of this optimization, also adopted in de Sousa (2015), is clear: a weight wi j measures (between 0 and 1)how similar x j is to xi , independent of the actual distances unlike what is used for Gaussian weights. It is thus morelikely to build strong edges between quite separate manifolds, whereas Gaussian weights would just assume they areonly weakly linked simply based on large distances.With these (directed) edge weights, the labels are inferred like in LLGC, by solving for the matrix Q that minimizes:kQ Pk2F η wi j Qi Q j2.(7)i, jThe resulting implicit diffusion thus uses a locally-optimized “composite” Laplacian, and this labeling was shown toimprove accuracy sometimes.2.5. Measure-based ClassifiersPrevious classifiers are, in essence, using C different scalar functions, the largest of which indicates the predicted labelfor a given node. Subramanya et al. Subramanya and Bilmes (2011) noted that these C values can also be seen asproviding a histogram (or probability distribution) for each data point, measuring the probability for this point to beof class c. As a consequence, they proposed to consider the smoothness of an assignment of probability distributionsover the graph (instead of the smoothness of the numerical labels) to further improve accuracy. They formulated aloss function based on Kullback-Leibler (KL) divergence1 which, after several tweaks, results in a convex programming problem. This problem was then solved via an alternating minimization (AM) approach, akin to the traditionalExpectation-Maximization of Dempster et al. (1977) method, with closed-form updates. Improvements in classification accuracy were demonstrated from this measure-based approach. Further noting that the KL divergence cannotcapture interactions between histogram bins, Solomon et al. (2014) suggested using optimal transport instead, andproposed a generalized Dirichlet energy over distribution-valued nodes of a graph to induce information propagation.However, they neither discussed computational complexity nor demonstrated improvement in classification accuracy.2.6. Learning with weight optimizationFinally, the best results to date in graph-based SSL involve optimizing the edge weights globally, to best adapt themto the data. In fact, the idea dates back to HGF, where the authors noted that the Gaussian weights from Eq. (1) canbe extended to read: D(xi,d x j,d )2(8)wσi j exp σd2d 11 KLdivergence is also called “relative entropy” or “information gain’, depending on the scientific field where it is invoked.5

Figure 2: Swirls. For a 2D datasets with points on 4 spiraling arms with 100 points and one label (marked as a colored cross) per spiral, previousSSL approaches fail to capture the obvious manifold structure of the data (each point is colored based on the inferred labeling). The three leftmostresults (HGF from Zhu et al. (2003), LLGC from Zhou et al. (2003)) and LNP from Wang et al. (2006)) rely on local choices of edge weights, withthe leftmost two overemphasizing label proximity and the third one focusing on local linearity instead. While the AEW approach from Karasuyamaand Mamitsuka (2017) optimizes the weights globally (before using label propagation through HGF or LLGC), it properly discovers only twospirals at best. Our approach distinguishes the four spirals even with only one label per class.Figure 3: Weighted graph. We visualize the different weighted graphs of the Swirls example in Fig. 2 as 2D matrices where element (i, j) is wi j(for LNP, (wi j w ji )/2). Points are arranged by class (first 100 for class 1, ., last 100 for class 4), but shuffled within each class. Since weightsmeasure similarity between vertices of a common edge, one should expect a block diagonal structure, with weights corresponding to interclassedges near zero. Weights below 10 5 are not shown, and values above 0.1 are all displayed in a dark blue color.where now the D scaling hyperparameters σ (σ1 , .σD ) can be optimized based on the data to “learn” the relevant dimensions of the dataset. Unfortunately, Zhu et al. (2003) proposed to minimize label entropy to learn thesehyperparameters, which led to very limited improvements: first, their optimization required regularization to preventthe hyperparameters from going to 0, otherwise it would end up selecting non-zero weights only between labels andtheir spatially closest unlabeled nodes; second, their entropy-based functional is highly nonlinear and nonconvex,hampering its efficient minimization.Recently, Karasuyama and Mamitsuka Karasuyama and Mamitsuka (2017) revisited this idea by optimizing parameters σ (σ1 , .σD ) so as to minimize an LLE-type energy instead: they compute σ through1arg min kxi wσi j x j k2 . k wσik j N(i){σ1 ,.,σD } iThat is, they mixed the LLE weights Roweis and Saul (2000) advocated in LNP with the weight normalization proposed in Normalized Cuts Shi and Malik (2000), and used the above “normalized LLE” energy to find the hyperparameters {σd }d 1.D that make the weighted graph best fit the input data manifolds (i.e., for which every point can be bestrepresented as a linear combination of its neighbors). The authors dubbed their approach Adaptive Edge Weighting(AEW). Unfortunately, the energy they use to optimize the weights is not convex. Minimizing it requires some formof gradient descent, and each gradient evaluation can be quite costly for high-dimensional input spaces since it is ofcomplexity O(nKD2 ). Yet, diffusing using HGF or LLGC the known labels through the graph with AEW-optimizedweights was shown to beat previous methods on many datasets.2.7. DiscussionThe best current diffusion-based SSL approaches rely on graphs with learned edge weights (and their correspondingLaplacian operators) to proceed with label propagation in a way that most effectively leverages oft-present geometricproperties of large data sets. However, they are often hampered by highly nonlinear optimizations of the weights,preventing their scalability and reliability; moreover, the loss function used to learn the weights is totally independent of the subsequent diffusion or even of the labels that are known. Instead, we propose to construct a consistentoptimization of edge weights which, with the corresponding anisotropic diffusion they define, allows to adapt to the6

Figure 4: Weights visualization. For the swirls dataset (see Fig. 2), we display the top 10% (in red) and the bottom 10% (in blue) of edge weightsin terms of their magnitude, for the various types of Laplacian used in SSL. The graph was constructed with each node connected to its 20 nearestneighbors. Optimally, edges connecting nearby points of the same class should have larger weights, while edges connecting different classes shouldhave vanishing weights; our approach fares much better than the other SSL techniques.intrinsic geometric structure of the data in order to capture the cluster and manifold assumptions. Our variationalapproach thus improves on previous state-of-the-art semi-supervised learning methods in a number of ways: Edge weights are globally optimized in order to give rise to a (possibly anisotropic) diffusion of labels that bestadapts to the intrinsic geometry of the data; Our loss function is biconvex in edge weights and label likelihood scores, rendering its optimization via alternating minimization particularly efficient; The resulting label propagation is interpolatory (i.e., is guaranteed to keep the known labels intact) and not biasedby large density variations in the input data.More importantly, our approach outperforms previous methods on typical datasets used to evaluate labeling accuracy.3. Our Laplacian-Optimized ClassifierIn this section, we present our variational approach to design a classifier through optimization of the Laplace operator.3.1. NotationsWe will use the same notations as defined earlier in Sec. 2.1. We also assemble the coordinates of the input data pointsinto a n D matrix X (x1 , .xn ), and we denote by S the label selection matrix, that is, the n n diagonal matrixsuch that Sii 1 if xi DL (i.e., if i ), and Sii 0 otherwise. Based on the data-dependent graph G {V , E } ofa set of n vertices V , edges E and a (fixed, but arbitrary) orientation of its edges, we also assemble the “adjacency”matrix d of size E n containing 1 for the endpoints of a given edge with the sign based on edge orientation; thatis, dev 1 if v refers to the index of the target vertex of the edge with index e, dev 1 if v refers to the index of thesource vertex of the edge with index e, and dev 0 otherwise. This matrix is the classical discrete exterior derivativeon 0-forms (see Desbrun et al. (2008)), i.e., a sparse matrix that is completely determined by the graph topology andthe chosen edge orientations. With a slight notation abuse, we will also denote by d the same matrix but where theabsolute value of each element is used instead; i.e., this matrix contains as many 1’s on the row i as there are neighborsof the corresponding vertex xi in G . Finally, we will use the n C matrix Q already mentioned in Sec. 2.2 to store Cdiscrete functions on graph G that indicate the likelihood score for each node to be of a given class.3.2. Parameterization of the LaplacianA reasonable assumption for SSL is that the data in the embedding space is concentrated around a manifold or aset of manifolds. For example, scans of handwritten digits possess certainly less possible variations in style thanthe dimensionality of the ambient representation space which is the number of pixels the scans have. Exploiting thegeometric structures of input datasets is thus crucial: in effect, it brings down the dimensionality of the problem athand. Since the most successful approaches to graph-based SSL thus far have been based on some form of labeldiffusion to best leverage the structure of data, we adopt a similar approach here.Failures of previous Laplacians. However, selecting a good Laplace operator to induce label propagation is a difficulttask. As we reviewed in Sec. 2, the use of “pre-canned” weights on edges such as the Gaussian weights from Eq. (1)relies too much on pure distances in the embedding space to be useful, strongly biasing diffusion in high densityregions; even variants taking local densities into account such as Zelnik-Manor and Perona (2004) to try to decreasethis high dependence to irregular sampling suffer from this issue. Locally-optimized weights like in Eq. (6) fare sometimes better, but double the number of variables needed and totally ignore distances within kNN neighborhoods. Other7

data-dependent choices of weights such as AEW showed improved results by globally optimizing a more general parameterization of the wei

As machine learning evolves, the importance of geometry and topology of data is increasingly acknowledged, seeSindhwani et al.(2010). . Corresponding author Email addresses: mbudnin@caltech.edu (Max Budninskiy), aabdelaz@caltech.edu (Ameera Abdelaziz), ytong@msu.edu (Yiying Tong), mathieu@caltech.edu (Mathieu Desbrun) Preprint submitted to .