"""Laplacian matrix of graphs. """ # Copyright (C) 2004-2018 by # Aric Hagberg # Dan Schult # Pieter Swart # All rights reserved. # BSD license. import networkx as nx from networkx.utils import not_implemented_for __author__ = "\n".join(['Aric Hagberg ', 'Pieter Swart (swart@lanl.gov)', 'Dan Schult (dschult@colgate.edu)', 'Alejandro Weinstein ']) __all__ = ['laplacian_matrix', 'normalized_laplacian_matrix', 'directed_laplacian_matrix'] @not_implemented_for('directed') def laplacian_matrix(G, nodelist=None, weight='weight'): """Return the Laplacian matrix of G. The graph Laplacian is the matrix L = D - A, where A is the adjacency matrix and D is the diagonal matrix of node degrees. Parameters ---------- G : graph A NetworkX graph nodelist : list, optional The rows and columns are ordered according to the nodes in nodelist. If nodelist is None, then the ordering is produced by G.nodes(). weight : string or None, optional (default='weight') The edge data key used to compute each value in the matrix. If None, then each edge has weight 1. Returns ------- L : SciPy sparse matrix The Laplacian matrix of G. Notes ----- For MultiGraph/MultiDiGraph, the edges weights are summed. See Also -------- to_numpy_matrix normalized_laplacian_matrix """ import scipy.sparse if nodelist is None: nodelist = list(G) A = nx.to_scipy_sparse_matrix(G, nodelist=nodelist, weight=weight, format='csr') n, m = A.shape diags = A.sum(axis=1) D = scipy.sparse.spdiags(diags.flatten(), [0], m, n, format='csr') return D - A @not_implemented_for('directed') def normalized_laplacian_matrix(G, nodelist=None, weight='weight'): r"""Return the normalized Laplacian matrix of G. The normalized graph Laplacian is the matrix .. math:: N = D^{-1/2} L D^{-1/2} where `L` is the graph Laplacian and `D` is the diagonal matrix of node degrees. Parameters ---------- G : graph A NetworkX graph nodelist : list, optional The rows and columns are ordered according to the nodes in nodelist. If nodelist is None, then the ordering is produced by G.nodes(). weight : string or None, optional (default='weight') The edge data key used to compute each value in the matrix. If None, then each edge has weight 1. Returns ------- N : NumPy matrix The normalized Laplacian matrix of G. Notes ----- For MultiGraph/MultiDiGraph, the edges weights are summed. See to_numpy_matrix for other options. If the Graph contains selfloops, D is defined as diag(sum(A,1)), where A is the adjacency matrix [2]_. See Also -------- laplacian_matrix References ---------- .. [1] Fan Chung-Graham, Spectral Graph Theory, CBMS Regional Conference Series in Mathematics, Number 92, 1997. .. [2] Steve Butler, Interlacing For Weighted Graphs Using The Normalized Laplacian, Electronic Journal of Linear Algebra, Volume 16, pp. 90-98, March 2007. """ import scipy import scipy.sparse if nodelist is None: nodelist = list(G) A = nx.to_scipy_sparse_matrix(G, nodelist=nodelist, weight=weight, format='csr') n, m = A.shape diags = A.sum(axis=1).flatten() D = scipy.sparse.spdiags(diags, [0], m, n, format='csr') L = D - A with scipy.errstate(divide='ignore'): diags_sqrt = 1.0 / scipy.sqrt(diags) diags_sqrt[scipy.isinf(diags_sqrt)] = 0 DH = scipy.sparse.spdiags(diags_sqrt, [0], m, n, format='csr') return DH.dot(L.dot(DH)) ############################################################################### # Code based on # https://bitbucket.org/bedwards/networkx-community/src/370bd69fc02f/networkx/algorithms/community/ @not_implemented_for('undirected') @not_implemented_for('multigraph') def directed_laplacian_matrix(G, nodelist=None, weight='weight', walk_type=None, alpha=0.95): r"""Return the directed Laplacian matrix of G. The graph directed Laplacian is the matrix .. math:: L = I - (\Phi^{1/2} P \Phi^{-1/2} + \Phi^{-1/2} P^T \Phi^{1/2} ) / 2 where `I` is the identity matrix, `P` is the transition matrix of the graph, and `\Phi` a matrix with the Perron vector of `P` in the diagonal and zeros elsewhere. Depending on the value of walk_type, `P` can be the transition matrix induced by a random walk, a lazy random walk, or a random walk with teleportation (PageRank). Parameters ---------- G : DiGraph A NetworkX graph nodelist : list, optional The rows and columns are ordered according to the nodes in nodelist. If nodelist is None, then the ordering is produced by G.nodes(). weight : string or None, optional (default='weight') The edge data key used to compute each value in the matrix. If None, then each edge has weight 1. walk_type : string or None, optional (default=None) If None, `P` is selected depending on the properties of the graph. Otherwise is one of 'random', 'lazy', or 'pagerank' alpha : real (1 - alpha) is the teleportation probability used with pagerank Returns ------- L : NumPy array Normalized Laplacian of G. Raises ------ NetworkXError If NumPy cannot be imported NetworkXNotImplemnted If G is not a DiGraph Notes ----- Only implemented for DiGraphs See Also -------- laplacian_matrix References ---------- .. [1] Fan Chung (2005). Laplacians and the Cheeger inequality for directed graphs. Annals of Combinatorics, 9(1), 2005 """ import scipy as sp from scipy.sparse import identity, spdiags, linalg if walk_type is None: if nx.is_strongly_connected(G): if nx.is_aperiodic(G): walk_type = "random" else: walk_type = "lazy" else: walk_type = "pagerank" M = nx.to_scipy_sparse_matrix(G, nodelist=nodelist, weight=weight, dtype=float) n, m = M.shape if walk_type in ["random", "lazy"]: DI = spdiags(1.0 / sp.array(M.sum(axis=1).flat), [0], n, n) if walk_type == "random": P = DI * M else: I = identity(n) P = (I + DI * M) / 2.0 elif walk_type == "pagerank": if not (0 < alpha < 1): raise nx.NetworkXError('alpha must be between 0 and 1') # this is using a dense representation M = M.todense() # add constant to dangling nodes' row dangling = sp.where(M.sum(axis=1) == 0) for d in dangling[0]: M[d] = 1.0 / n # normalize M = M / M.sum(axis=1) P = alpha * M + (1 - alpha) / n else: raise nx.NetworkXError("walk_type must be random, lazy, or pagerank") evals, evecs = linalg.eigs(P.T, k=1) v = evecs.flatten().real p = v / v.sum() sqrtp = sp.sqrt(p) Q = spdiags(sqrtp, [0], n, n) * P * spdiags(1.0 / sqrtp, [0], n, n) I = sp.identity(len(G)) return I - (Q + Q.T) / 2.0 # fixture for nose tests def setup_module(module): from nose import SkipTest try: import numpy except: raise SkipTest("NumPy not available")