Source code for netrd.distance.laplacian_spectral_method


Graph distance based on :

author: Guillaume St-Onge
Submitted as part of the 2019 NetSI Collabathon.

import numpy as np
import networkx as nx
from .base import BaseDistance
from ..utilities.graph import unweighted
from scipy.special import erf
from scipy.integrate import quad
from scipy.linalg import eigvalsh
from scipy.sparse.csgraph import csgraph_from_dense
from scipy.sparse.csgraph import laplacian
from scipy.sparse.linalg import eigsh

[docs]class LaplacianSpectral(BaseDistance): """Flexible distance able to compare the spectrum of the Laplacian in many ways."""
[docs] @unweighted def dist( self, G1, G2, normed=True, kernel='normal', hwhm=0.011775, measure='jensen-shannon', k=None, which='LM', ): """Graph distances using different measure between the Laplacian spectra of the two graphs The spectra of both Laplacian matrices (normalized or not) is computed. Then, the discrete spectra are convolved with a kernel to produce continuous ones. Finally, these distribution are compared using a metric. The results dictionary also stores a 2-tuple of the underlying adjacency matrices in the key `'adjacency_matrices'`, the Laplacian matrices in `'laplacian_matrices'`, the eigenvalues of the Laplacians in `'eigenvalues'`. If the networks being compared are directed, the augmented adjacency matrices are calculated and stored in `'augmented_adjacency_matrices'`. Parameters ---------- G1, G2 (nx.Graph) two networkx graphs to be compared. normed (bool) If True, uses the normalized laplacian matrix, otherwise the raw laplacian matrix is used. kernel (str) kernel to obtain a continuous spectrum. Choices available are 'normal', 'lorentzian', or None. If None is chosen, the discrete spectrum is used instead, and the measure is simply the euclidean distance between the vector of eigenvalues for each graph. hwhm (float) half-width at half-maximum for the kernel. The default value is chosen such that the standard deviation for the normal distribution is :math:`0.01`, as in reference [1]_. This option is relevant only if kernel is not None. measure (str) metric between the two continuous spectra. Choices available are 'jensen-shannon' or 'euclidean'. This option is relevant only if kernel is not None. k (int) number of eigenvalues kept for the (discrete) spectrum, also used to create the continuous spectrum. If None, all the eigenvalues are used. k must be smaller (strictly) than the size of both graphs. which (str) if k is not None, this option specifies the eigenvalues that are kept. See the choices offered by `scipy.sparse.linalg.eigsh`. The largest eigenvalues in magnitude are kept by default. Returns ------- dist (float) the distance between G1 and G2. Notes ----- The methods are usually applied to undirected (unweighted) networks. We however relax this assumption using the same method proposed for the Hamming-Ipsen-Mikhailov. See [2]_. References ---------- .. [1] .. [2] """ adj1 = nx.to_numpy_array(G1) adj2 = nx.to_numpy_array(G2) self.results['adjacency_matrices'] = adj1, adj2 directed = nx.is_directed(G1) or nx.is_directed(G2) if directed: # create augmented adjacency matrices N1 = len(G1) N2 = len(G2) null_mat1 = np.zeros((N1, N1)) null_mat2 = np.zeros((N2, N2)) adj1 = np.block([[null_mat1, adj1.T], [adj1, null_mat1]]) adj2 = np.block([[null_mat2, adj2.T], [adj2, null_mat2]]) self.results['augmented_adjacency_matrices'] = adj1, adj2 # get the laplacian matrices lap1 = laplacian(adj1, normed=normed) lap2 = laplacian(adj2, normed=normed) self.results['laplacian_matrices'] = lap1, lap2 # get the eigenvalues of the laplacian matrices if k is None: ev1 = np.abs(eigvalsh(lap1)) ev2 = np.abs(eigvalsh(lap2)) else: # transform the dense laplacian matrices to sparse representations lap1 = csgraph_from_dense(lap1) lap2 = csgraph_from_dense(lap2) ev1 = np.abs(eigsh(lap1, k=k, which=which)[0]) ev2 = np.abs(eigsh(lap2, k=k, which=which)[0]) self.results['eigenvalues'] = ev1, ev2 if kernel is not None: # define the proper support a = 0 if normed: b = 2 else: b = np.inf # create continuous spectra density1 = _create_continuous_spectrum(ev1, kernel, hwhm, a, b) density2 = _create_continuous_spectrum(ev2, kernel, hwhm, a, b) # compare the spectra dist = _spectra_comparison(density1, density2, a, b, measure) self.results['dist'] = dist else: # euclidean distance between the two discrete spectra dist = np.linalg.norm(ev1 - ev2) self.results['dist'] = dist return dist
def _create_continuous_spectrum(eigenvalues, kernel, hwhm, a, b): """Convert a set of eigenvalues into a normalized density function The discret spectrum (sum of dirac delta) is convolved with a kernel and renormalized. Parameters ---------- eigenvalues (array): list of eigenvalues. kernel (str): kernel to be used for the convolution with the discrete spectrum. hwhm (float): half-width at half-maximum for the kernel. a,b (float): lower and upper bounds of the support for the eigenvalues. Returns ------- density (function): one argument function for the continuous spectral density. """ # define density and repartition function for each eigenvalue if kernel == "normal": std = hwhm / 1.1775 f = lambda x, xp: np.exp(-((x - xp) ** 2) / (2 * std**2)) / np.sqrt( 2 * np.pi * std**2 ) F = lambda x, xp: (1 + erf((x - xp) / (np.sqrt(2) * std))) / 2 elif kernel == "lorentzian": f = lambda x, xp: hwhm / (np.pi * (hwhm**2 + (x - xp) ** 2)) F = lambda x, xp: np.arctan((x - xp) / hwhm) / np.pi + 1 / 2 # compute normalization factor and define density function Z = np.sum(F(b, eigenvalues) - F(a, eigenvalues)) density = lambda x: np.sum(f(x, eigenvalues)) / Z return density def _spectra_comparison(density1, density2, a, b, measure): """Apply a metric to compare the continuous spectra Parameters ---------- density1, density2 (function): one argument functions for the continuous spectral densities. a,b (float): lower and upper bounds of the support for the eigenvalues. measure (str): metric between the two continuous spectra. Returns ------- dist (float): distance between the spectra. """ if measure == "jensen-shannon": M = lambda x: (density1(x) + density2(x)) / 2 jensen_shannon = ( _kullback_leiber(density1, M, a, b) + _kullback_leiber(density2, M, a, b) ) / 2 dist = np.sqrt(jensen_shannon) elif measure == "euclidean": integrand = lambda x: (density1(x) - density2(x)) ** 2 dist = np.sqrt(quad(integrand, a, b)[0]) return dist def _kullback_leiber(f1, f2, a, b): def integrand(x): if f1(x) > 0 and f2(x) > 0: result = f1(x) * np.log(f1(x) / f2(x)) else: result = 0 return result return quad(integrand, a, b)[0]