Merge pull request #44 from pmacg/iss27-visualise-spectrum

Issue #27 - add option to plot eigenvalues in spectrum methods

docs/my_graph.edgelist

@@ -0,0 +1,3 @@
+0 1
+1 2
+2 0

docs/source/changelog.rst

@@ -7,6 +7,7 @@ Unreleased
 **Added**
 
 * `Issue #32 <https://github.com/pmacg/py-sgtl/issues/32>`_ - add tensor product method for combining graphs
+* `Issue #27 <https://github.com/pmacg/py-sgtl/issues/27>`_ - add option to plot spectrum of graph
 
 0.4.4 - 2022-01-07
 ------------------

docs/source/generated/sgtl.graph.Graph.rst

@@ -26,6 +26,7 @@ sgtl.graph.Graph
       ~Graph.inverse_degree_matrix
       ~Graph.inverse_sqrt_degree_matrix
       ~Graph.laplacian_matrix
+      ~Graph.normalised_adjacency_matrix
       ~Graph.normalised_laplacian_matrix
       ~Graph.number_of_edges
       ~Graph.number_of_vertices

docs/source/generated/sgtl.spectrum.normalised_adjacency_spectrum.rst

@@ -0,0 +1,6 @@
+sgtl.spectrum.normalised\_adjacency\_spectrum
+=============================================
+
+.. currentmodule:: sgtl.spectrum
+
+.. autofunction:: normalised_adjacency_spectrum

docs/source/generated/sgtl.spectrum.rst

@@ -16,6 +16,7 @@
 
       adjacency_spectrum
       laplacian_spectrum
+      normalised_adjacency_spectrum
       normalised_laplacian_spectrum
 
 

docs/source/getting-started.rst

@@ -42,29 +42,27 @@ For example:
            [0., 0., 1., 0., 1.],
            [1., 0., 0., 1., 0.]])
 
-You can also generate some graphs with standard shapes.
-For example:
+Or, you can load a more complex graph from an edgelist file like this.
 
     >>> import sgtl.graph
-    >>> graph = sgtl.graph.cycle_graph(5)
-    >>> graph.adjacency_matrix().toarray()
-    array([[0., 1., 0., 0., 1.],
-           [1., 0., 1., 0., 0.],
-           [0., 1., 0., 1., 0.],
-           [0., 0., 1., 0., 1.],
-           [1., 0., 0., 1., 0.]])
+    >>> graph = sgtl.graph.from_edgelist("my_graph.edgelist")
 
-You can also generate some graphs with standard shapes.
-For example:
+See the documentation of the :any:`sgtl.graph.from_edgelist` method for more information on the
+required format of the edgelist file.
+
+Viewing the spectrum of a graph
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Investigating the spectrum of a graph is very simple. For example, you could visualise
+the spectrum of the tensor product of two graphs like this.
+
+.. code-block:: python
 
     >>> import sgtl.graph
-    >>> graph = sgtl.graph.cycle_graph(5)
-    >>> graph.adjacency_matrix().toarray()
-    array([[0., 1., 0., 0., 1.],
-           [1., 0., 1., 0., 0.],
-           [0., 1., 0., 1., 0.],
-           [0., 0., 1., 0., 1.],
-           [1., 0., 0., 1., 0.]])
+    >>> import sgtl.spectrum
+    >>> g1 = sgtl.graph.cycle_graph(5)
+    >>> g2 = sgtl.graph.complete_graph(5)
+    >>> g3 = g1.tensor_product(g2)
+    >>> spectrum = sgtl.spectrum.normalised_adjacency_spectrum(g3, plot=True)
 
 The stochastic block model
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~

sgtl/graph.py

@@ -116,6 +116,17 @@ def adjacency_matrix(self):
         """
         return self.adj_mat
 
+    def normalised_adjacency_matrix(self):
+        """
+        Return the normalised adjacency matrix of the graph. The normalised adjacency matrix is defined to be
+
+        .. math::
+            \\mathcal{A} = D^{-1/2} A D^{-1/2}
+
+        where :math:`A` is the unnormalised adjacency matrix and :math:`D` is the diagonal degree matrix of the graph.
+        """
+        return self.inverse_sqrt_degree_matrix() @ self.adjacency_matrix() @ self.inverse_sqrt_degree_matrix()
+
     def laplacian_matrix(self):
         """
         Construct the Laplacian matrix of the graph. The Laplacian matrix is defined to be
@@ -498,7 +509,7 @@ def rbf_graph(data, variance=1, threshold=0.1) -> Graph:
     return Graph(adj_mat)
 
 
-def from_edgelist(filename: str, directed=False, num_vertices=None, **kwargs) -> Graph:
+def from_edgelist(filename, directed=False, num_vertices=None, **kwargs) -> Graph:
     """
     Construct an ``sgtl.Graph`` object from an edgelist file.
 

sgtl/spectrum.py

@@ -5,19 +5,36 @@
 import scipy
 import scipy.linalg
 import scipy.sparse.linalg
+import matplotlib.pyplot as plt
 import sgtl
 
 
-def adjacency_spectrum(graph: sgtl.Graph, num_eigenvalues=None) -> List[float]:
+def _plot_spectrum(spectrum):
+    """
+    Plot the given list of floats as a matplotlib point plot. If the eigenvalues given are complex, we
+    ignore the imaginary part!
+    """
+    # If the given eigenvalues are complex, ignore the imaginary part!
+    if len(spectrum) > 0 and isinstance(spectrum[0], complex):
+        spectrum = [x.real for x in spectrum]
+
+    plt.plot(sorted(spectrum), '.')
+    plt.show()
+
+
+def adjacency_spectrum(graph: sgtl.Graph, num_eigenvalues=None, plot=False) -> List[complex]:
     """
     Return a list of the eigenvalues of the adjacency matrix of the given graph.
 
     Computing all of the eigenvalues of the adjacency matrix is an expensive operation. You should set num_eigenvalues
     to control the number of eigenvalues you would like to compute.
 
+    If the plot argument is true, then the spectrum is plotted using matplotlib before being returned.
+
     :param graph: The SGTL graph object on which to operate.
     :param num_eigenvalues: How many eigenvalues to return. If this value is not none, the method will return the
                             eigenvalues with the maximum absolute values.
+    :param plot: Whether to plot the spectrum before returning it
     """
     if num_eigenvalues is None:
         num_eigenvalues = graph.number_of_vertices()
@@ -26,10 +43,44 @@ def adjacency_spectrum(graph: sgtl.Graph, num_eigenvalues=None) -> List[float]:
         eigvals, _ = scipy.linalg.eig(graph.adjacency_matrix().toarray())
     else:
         eigvals, _ = scipy.sparse.linalg.eigs(graph.adjacency_matrix(), k=num_eigenvalues)
+
+    # Plot if required
+    if plot:
+        _plot_spectrum(eigvals)
+
+    return list(eigvals)
+
+
+def normalised_adjacency_spectrum(graph: sgtl.Graph, num_eigenvalues=None, plot=False) -> List[complex]:
+    """
+    Return a list of the eigenvalues of the normalised adjacency matrix of the given graph.
+
+    Computing all of the eigenvalues of the adjacency matrix is an expensive operation. You should set num_eigenvalues
+    to control the number of eigenvalues you would like to compute.
+
+    If the plot argument is true, then the spectrum is plotted using matplotlib before being returned.
+
+    :param graph: The SGTL graph object on which to operate.
+    :param num_eigenvalues: How many eigenvalues to return. If this value is not none, the method will return the
+                            eigenvalues with the maximum absolute values.
+    :param plot: Whether to plot the spectrum before returning it
+    """
+    if num_eigenvalues is None:
+        num_eigenvalues = graph.number_of_vertices()
+
+    if num_eigenvalues >= graph.number_of_vertices() - 1:
+        eigvals, _ = scipy.linalg.eig(graph.normalised_adjacency_matrix().toarray())
+    else:
+        eigvals, _ = scipy.sparse.linalg.eigs(graph.normalised_adjacency_matrix(), k=num_eigenvalues)
+
+    # Plot if required
+    if plot:
+        _plot_spectrum(eigvals)
+
     return list(eigvals)
 
 
-def laplacian_spectrum(graph: sgtl.Graph, num_eigenvalues=None, magnitude='smallest') -> List[float]:
+def laplacian_spectrum(graph: sgtl.Graph, num_eigenvalues=None, magnitude='smallest', plot=False) -> List[complex]:
     """
     Return a list of the eigenvalues of the laplacian matrix of the given graph.
 
@@ -39,10 +90,13 @@ def laplacian_spectrum(graph: sgtl.Graph, num_eigenvalues=None, magnitude='small
     When computing only a subset of the eigenvalues, the 'magnitude' parameter controls whether the eigenvalues with
     largest or smallest magnitude will be returned. This defaults to 'smallest'.
 
+    If the plot argument is true, then the spectrum is plotted using matplotlib before being returned.
+
     :param graph: The SGTL graph object on which to operate.
     :param num_eigenvalues: How many eigenvalues to return.
     :param magnitude: Should be 'smallest' or 'largest' - whether to return the eigenvalues with smallest or largest
                       magnitude.
+    :param plot: Whether to plot the spectrum before returning it
     """
     if num_eigenvalues is None:
         num_eigenvalues = graph.number_of_vertices()
@@ -56,10 +110,18 @@ def laplacian_spectrum(graph: sgtl.Graph, num_eigenvalues=None, magnitude='small
         eigvals, _ = scipy.linalg.eig(graph.laplacian_matrix().toarray())
     else:
         eigvals, _ = scipy.sparse.linalg.eigs(graph.laplacian_matrix(), k=num_eigenvalues, which=mag_sp)
+
+    # Plot if required
+    if plot:
+        _plot_spectrum(eigvals)
+
     return list(eigvals)
 
 
-def normalised_laplacian_spectrum(graph: sgtl.Graph, num_eigenvalues=None, magnitude='smallest') -> List[float]:
+def normalised_laplacian_spectrum(graph: sgtl.Graph,
+                                  num_eigenvalues=None,
+                                  magnitude='smallest',
+                                  plot=False) -> List[complex]:
     """
     Return a list of the eigenvalues of the normalised laplacian matrix of the given graph.
 
@@ -69,10 +131,13 @@ def normalised_laplacian_spectrum(graph: sgtl.Graph, num_eigenvalues=None, magni
     When computing only a subset of the eigenvalues, the 'magnitude' parameter controls whether the eigenvalues with
     largest or smallest magnitude will be returned. This defaults to 'smallest'.
 
+    If the plot argument is true, then the spectrum is plotted using matplotlib before being returned.
+
     :param graph: The SGTL graph object on which to operate.
     :param num_eigenvalues: How many eigenvalues to return.
     :param magnitude: Should be 'smallest' or 'largest' - whether to return the eigenvalues with smallest or largest
                       magnitude.
+    :param plot: Whether to plot the spectrum before returning it
     """
     if num_eigenvalues is None:
         num_eigenvalues = graph.number_of_vertices()
@@ -86,4 +151,9 @@ def normalised_laplacian_spectrum(graph: sgtl.Graph, num_eigenvalues=None, magni
         eigvals, _ = scipy.linalg.eig(graph.normalised_laplacian_matrix().toarray())
     else:
         eigvals, _ = scipy.sparse.linalg.eigs(graph.normalised_laplacian_matrix(), k=num_eigenvalues, which=mag_sp)
+
+    # Plot if required
+    if plot:
+        _plot_spectrum(eigvals)
+
     return list(eigvals)