_laplacian.py 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128
  1. """
  2. Laplacian of a compressed-sparse graph
  3. """
  4. # Authors: Aric Hagberg <hagberg@lanl.gov>
  5. # Gael Varoquaux <gael.varoquaux@normalesup.org>
  6. # Jake Vanderplas <vanderplas@astro.washington.edu>
  7. # License: BSD
  8. from __future__ import division, print_function, absolute_import
  9. import numpy as np
  10. from scipy.sparse import isspmatrix
  11. ###############################################################################
  12. # Graph laplacian
  13. def laplacian(csgraph, normed=False, return_diag=False, use_out_degree=False):
  14. """
  15. Return the Laplacian matrix of a directed graph.
  16. Parameters
  17. ----------
  18. csgraph : array_like or sparse matrix, 2 dimensions
  19. compressed-sparse graph, with shape (N, N).
  20. normed : bool, optional
  21. If True, then compute normalized Laplacian.
  22. return_diag : bool, optional
  23. If True, then also return an array related to vertex degrees.
  24. use_out_degree : bool, optional
  25. If True, then use out-degree instead of in-degree.
  26. This distinction matters only if the graph is asymmetric.
  27. Default: False.
  28. Returns
  29. -------
  30. lap : ndarray or sparse matrix
  31. The N x N laplacian matrix of csgraph. It will be a numpy array (dense)
  32. if the input was dense, or a sparse matrix otherwise.
  33. diag : ndarray, optional
  34. The length-N diagonal of the Laplacian matrix.
  35. For the normalized Laplacian, this is the array of square roots
  36. of vertex degrees or 1 if the degree is zero.
  37. Notes
  38. -----
  39. The Laplacian matrix of a graph is sometimes referred to as the
  40. "Kirchoff matrix" or the "admittance matrix", and is useful in many
  41. parts of spectral graph theory. In particular, the eigen-decomposition
  42. of the laplacian matrix can give insight into many properties of the graph.
  43. Examples
  44. --------
  45. >>> from scipy.sparse import csgraph
  46. >>> G = np.arange(5) * np.arange(5)[:, np.newaxis]
  47. >>> G
  48. array([[ 0, 0, 0, 0, 0],
  49. [ 0, 1, 2, 3, 4],
  50. [ 0, 2, 4, 6, 8],
  51. [ 0, 3, 6, 9, 12],
  52. [ 0, 4, 8, 12, 16]])
  53. >>> csgraph.laplacian(G, normed=False)
  54. array([[ 0, 0, 0, 0, 0],
  55. [ 0, 9, -2, -3, -4],
  56. [ 0, -2, 16, -6, -8],
  57. [ 0, -3, -6, 21, -12],
  58. [ 0, -4, -8, -12, 24]])
  59. """
  60. if csgraph.ndim != 2 or csgraph.shape[0] != csgraph.shape[1]:
  61. raise ValueError('csgraph must be a square matrix or array')
  62. if normed and (np.issubdtype(csgraph.dtype, np.signedinteger)
  63. or np.issubdtype(csgraph.dtype, np.uint)):
  64. csgraph = csgraph.astype(float)
  65. create_lap = _laplacian_sparse if isspmatrix(csgraph) else _laplacian_dense
  66. degree_axis = 1 if use_out_degree else 0
  67. lap, d = create_lap(csgraph, normed=normed, axis=degree_axis)
  68. if return_diag:
  69. return lap, d
  70. return lap
  71. def _setdiag_dense(A, d):
  72. A.flat[::len(d)+1] = d
  73. def _laplacian_sparse(graph, normed=False, axis=0):
  74. if graph.format in ('lil', 'dok'):
  75. m = graph.tocoo()
  76. needs_copy = False
  77. else:
  78. m = graph
  79. needs_copy = True
  80. w = m.sum(axis=axis).getA1() - m.diagonal()
  81. if normed:
  82. m = m.tocoo(copy=needs_copy)
  83. isolated_node_mask = (w == 0)
  84. w = np.where(isolated_node_mask, 1, np.sqrt(w))
  85. m.data /= w[m.row]
  86. m.data /= w[m.col]
  87. m.data *= -1
  88. m.setdiag(1 - isolated_node_mask)
  89. else:
  90. if m.format == 'dia':
  91. m = m.copy()
  92. else:
  93. m = m.tocoo(copy=needs_copy)
  94. m.data *= -1
  95. m.setdiag(w)
  96. return m, w
  97. def _laplacian_dense(graph, normed=False, axis=0):
  98. m = np.array(graph)
  99. np.fill_diagonal(m, 0)
  100. w = m.sum(axis=axis)
  101. if normed:
  102. isolated_node_mask = (w == 0)
  103. w = np.where(isolated_node_mask, 1, np.sqrt(w))
  104. m /= w
  105. m /= w[:, np.newaxis]
  106. m *= -1
  107. _setdiag_dense(m, 1 - isolated_node_mask)
  108. else:
  109. m *= -1
  110. _setdiag_dense(m, w)
  111. return m, w