Why Laplacians?

2. Probabilistic

A justification for Laplacian methods which provides a much stronger justification for the particular form of the operator is that it encodes a particular probabilistic model on the space. In the case of a graph, this is a standard random walk along the edges, with a normalized version of the Laplacian (\(I-D^{-1}L\)) being the Markov transition matrix for a process defined on the vertices and transitioning at random along the edges from that vertex. In the continuum limit, the eigenfunctions of the Laplacian characterize a diffusion process in the ambient space (in particular, they are the eigenfunctions of the associated Fokker-Planck operator), with drift determined by the density of points over the manifold.

The random walk structure helps elucidate a number of properties of the Laplacian in characterizing the structure of the space. First, since a random walk is a local process, over short time scales it moves locally. Over long time scales, over finite graphs, the walk converges to a limit distribution, which describes the global structure. This may fail to be unique if, for example, the graph is not connected. This interpretation also explains the role of the spectrum of the Laplacian. A steady state is an eigenvector of the Markov transition matrix corresponding to the eigenvalue 1, i.e. \((I-D^{-1}L)v=v\), which exists by the Perron-Frobenius theorem. For a Markov chain which does converge to a unique steady state, the speed of this convergence is given by the size of the second eigenvalue: if this is one, the process has multiple steady states, which occurs when there are disconnected components of the graph. If it is near one, convergence occurs very slowly, as would be expected if there are regions of a graph where a walk would remain stuck with high probability because there are few links outside and many links inside. These cases correspond to existence of disconnected or nearly disconnected clusters in a graph, and explain the usefulness of spectral clustering. Quantitatively, the control provided over the structure by this eigenvalue is given by the Cheeger inequality (note that due to sign conventions, the distance of the second largest eigenvalue from 1 of the Markov transition matrix corresponds to the distance of the second smallest eigenvalue from 0 for the Laplacian).

The idea of Laplacian as coming from a random walk also gives, to me, the most intuitive justification for why a second derivative is used: it comes from Itō’s lemma. As the step size goes to zero, a pure random walk on a space converges, by a functional CLT, to a standard Brownian motion, call it \(W_t\). Now consider an arbitrary twice differentiable function on for simplicity, \(\mathbb{R}^n\), \(f(.):\mathbb{R}^n\to\mathbb{R}\). By Itō’s lemma, Taylor expanding \(f(.)\), obtain \(df(W_t)= \frac{1}{2} \Delta f(W_t)dt+\nabla f(W_t)dW_t\), and we see that a random walk induces a drift for any function which is precisely (1/2 of) the Laplacian. Applied to the density \(p(x,t)\), this gives the adjoint Markov operator, or Fokker-Planck equation for the density, \(dp(x,t)=\frac{1}{2} \Delta p(x,t)dt\), the continuous analogue of the transpose of the Markov transition matrix, which instead of mapping state to state maps distribution over states to distributions over states.⁷

So, we can see that one question the Laplacian answers is, given a density over my space, how will it evolve under a random walk. This of course raises the question of why the process needs to be a pure random walk, as opposed to some other process, which might be driven by non-spherical noise, have some drift or jump component, or depend on the location in the space. The answer is that it depends on the application. If nothing about the local structure of the space is known, it makes sense to have a procedure which is locally agnostic. On the other hand, if one is using the induced decomposition of the space to characterize objects arising from a structurally biased process, it may be useful to take that into account. This is especially the case if the methods are being used to describe something which actually does diffuse over a network or manifold-like space, in which case the local distance ought to be in terms of the characteristics of this process, which may be asymmetric or behave differently at boundaries.⁸ Modifications here should be domain specific, and indeed procedures like spectral clustering have been adapted in this way to different characterizations of communities.⁹

As an example, consider the utility of a process consisting of a random walk plus jumps which land at a random point on the space. This will induce a non-local component to the movement, and so induce connections between disconnected or weakly-connected components (which may appear disconnected under sampling). This additional component provides a sort of protection against unmeasured connectedness, and induces the commonly used regularized Laplacian. In particular, this generates the random surfer model underlying Google’s PageRank algorithm, arguably the most successful application of spectral graph methods.

Differential equations

Being a differential operator, probably the most obvious perspective from which to consider the Laplacian is (partial) differential equations, in which second order equations, in which the Laplacian is (part of) the leading term, comprise the core of the subject. The Laplacian shows up in the canonical Poisson equation \(\Delta u(x) = f(x)\), the Heat equation \(\partial_t u(x,t) = \Delta u(x,t)\), the wave equation \(\partial_{tt} u(x,t) = \Delta u(x,t)\), and a large variety of equations with lower order terms, both linear and nonlinear.

In economics (or at least the branches with which I am most familiar), second order PDEs show up most frequently in stochastic optimization problems, in which some noise follows a diffusion process, and an agent chooses a control to maximize a discounted reward over time. With discount rate r, control \(c_t\), noise \(dx_t=\mu (x_t,c_t)dt+\Sigma dW_t\) and reward \(u(x_t,c_t)\), this leads to the HJB equation for the value conditional on any given state as \[rV(x_t,t)=\underset{c_t}{\max} \{u(x_t,c_t) + \partial_t V(x_t,t)+\langle\nabla_x V(x_t,t)), \mu (x_t,c_t)\rangle + \frac{1}{2}Tr(\Sigma^T H_x V(x_t,t)\Sigma)\}\] where the last term, the trace of the Hessian of \(V\), is the weighted Laplacian. One can add boundary conditions (which may also be controlled), constraints, and so on to describe the problem of interest. The simplest case, with a geometric Brownian motion as state variable, no instantaneous payoff, and the only control a binary choice over a terminal boundary condition \(V(x_T,T)=\max\{x_T-p,0\}\) for constant \(p\), gives the well-known Black-Scholes formula for option pricing, though the idea is applicable in more general contexts.

To the extent that solutions inherit features of the space from the second order differential structure, the value function \(V(x,t)\) should be characterizable in terms of the properties of the Laplacian. This idea has attracted attention in the field of reinforcement learning (usually used to describe the discrete time version of this optimization problem), using Laplacian eigenfunctions and derived wavelets as a basis for value function approximation. In cases where the structure of the space is an unknown manifold, this method has the additional advantage of incorporating nonlinear dimension reduction, potentially providing substantial computational advantages over structured domains. As in the probabilistic case, when the driving process is not a pure random walk, it may be advantageous to incorporate structure by using, say, the Fokker-Planck operator and not just the standard Laplacian. However, in the case of substantially nonlinear reward or control, the linear component may provide a poor guide to the behavior of the solution, unlike in the Black-Scholes case where it is exact.

For its uses in more general data analysis, the PDE which is most directly applicable is the Heat equation, \(\partial_t u(x,t) = \Delta u(x,t)\), which, up to a constant, is exactly what I called a Fokker-Planck equation above.¹⁰ Applying standard methods for solution of initial value problems then gives a semigroup which tells us for any t>0 the distribution resulting from any given initial state \(u_0 (x)\), in the form of the Heat Kernel, denoted \(e^{t\Delta}\), which on Euclidean space takes the form \(e^{t\Delta} [u_0](y) = \int_{\mathbb{R}^n}\phi(\frac{|x-y|}{2t})u_0 (x) dx\), where \(\phi\) is the standard Gaussian density on \(\mathbb{R}^n\).

In other words, the solution is convolution with an increasingly dispersed Gaussian, as should be expected from the interpretation of adding Gaussian noise to each point. This gives another way of seeing the connection with Fourier analysis, since by the convolution theorem, any convolution operator acts componentwise on the Fourier frequencies by multiplication. The heat kernel then can be viewed as a form of low-pass filter, which downweights high frequency components by an amount increasing in t. Spectral methods based on the Laplacian, which restrict to the first few eigenfunctions and discard the rest, sharpen this characterization, generating an exact (continuous) low-pass filter. From a statistical perspective, this can be seen as a series estimator of the density.  In contrast, applying the heat kernel can be seen directly as a kernel density estimator: given an initial condition which is noisy (say, an empirical distribution), the heat kernel outputs a regularized version which is smooth: the time parameter plays the role of the bandwidth.

Statistical

The interpretation of Laplacian eigenfunctions and heat kernel smoothing as sieve and kernel methods respectively suggests evaluation of these methods for their statistical properties. Which of these methods is preferable depends on what one believes about the distribution over the graph or manifold, though as is generally the case with sieve versus kernel methods, asymptotic convergence rates are the same for smooth functions, and the decision of which to use usually comes down to other desiderata (e.g. the kernel method outputs a proper density, while the sieve filter is idempotent). More generally, the use of the heat equation picks out a single sieve basis (the Laplacian eigenfunctions) and a single kernel (the Gaussian), which can be compared to other sieves and kernels, which may have more desirable properties depending on the conjectured properties of the objects living on the space. It is well-known, for example, that the Gaussian kernel, as a proper kernel, yields suboptimal convergence rates for highly smooth densities. Since both produce isotropic representations, they also don’t deal with smoothness that may differ by direction (which is especially tricky in the non-Euclidean case since one often lacks canonical directions) or location. Sparse wavelet based methods, as usual, can help here, though it seems to me that proper analysis of optimal wavelet denoising over graphs and manifolds remains understudied.

As for more general alternatives, there are a large number of methods now in the machine learning and statistical literature, many of which rely on Laplacian based methods, via sieve or Reproducing Kernel Hilbert Space methods, the latter including several Bayesian approaches. Estimation over graphs and manifolds remains a hot topic, with a variety of approaches, only some of which rely on Laplacian structure.  As in every other aspect of machine learning these days, neural networks over manifolds provide a promising alternative.  One also has a choice between penalization methods of various kinds for graph-based predictions. The choice between these seems essentially similar to the choice of penalization in standard nonparametric prediction problems, with Laplacian-based methods corresponding to \(L^2\) approaches like splines and support vector machines: a tool in the toolbox and a good first choice for many problems, but not the last.¹¹

Topological

Beyond representing the space on which they live, part of the utility of Laplacian based methods is that important characteristics of the space can easily be “read off” of the Laplacian, especially the spectrum, for some definition of “important.” The use in clustering suggests that this is the degree of connection between the components of the space, a topological feature, or approximately one. In particular, it provides a measure of the 0th homology group of the space. This does not seem to be coincidental. The homology of a topological space (components, holes, voids, etc) can be determined by constructing a simplicial complex from the space, considering not just links between pairs of points, but sets of 3, or 4, and so on, as higher analogues of edges of a graph, which is a first order complex. The homology of each order is determined by the simplices of that order, so graph information is informative about the connected components. By extending the approximating space to higher orders, one can determine also the higher order homologies. Constructing a Laplacian over these higher order connections to obtain the combinatorial Laplacian then in fact gives an operator whose spectrum encodes the higher order homologies in the same way the graph Laplacian describes clustering.

In practice, there doesn’t seem to be a lot in the way of “spectral homology” to give higher order analogues of spectral clustering. Instead, topological data analysis relies mostly on versions of persistent homology, which computes the homology group directly from the complex constructed from connecting nearby points at each length scale. The difference seems to be due to different notions of robustness to noise. In spectral clustering, one claims to have clusters even if the graph is connected so long as there are relatively few links between components, with “few” measurable by a small second eigenvalue of the Laplacian. Persistent homology will instead only show existence of disconnected components if the graph actually is disconnected, for edges connected between sufficiently close points. If I understand this properly (not likely), the analogue for 1st homology group would be to say that there exists an “approximate” hole in the space if there are few paths across a gap relative to paths around it, where few could be measured by the second eigenvalue of a higher order Laplacian. Whether this is useful depends on the model for the space: for sampling exactly from a manifold, one will never get points in the “wrong” place and so persistent homology should describe its holes. If the data is instead drawn from a model with noise (interesting question: is there a simplicial analogue of the stochastic block model?), one will expect outliers and so even if data concentrates on a torus, the persistent homology may not find the bulk of the torus before the hole in the point cloud is completely filled in. To put it visually, the relative merit of these methods for different applications depends on whether you want your algorithm to tell you that a bialy is approximately a bagel, or approximately an onion roll.

Summary

Considering these methods, it gets easier to see why, if you knew not so much about a space, you would first take a look at the Laplacian over it. It describes local and global properties, gives a reasonable definition of smoothness and random processes on a space, and permits translation of the wide variety of second order (Hilbert space) statistical smoothing methods to strange-looking domains. It should be modified if a little more is known about the way things move on the space, and is probably exactly the wrong approach for considering at individual points, but it can tell us a lot if we just want to know how our data are put together.