use linalg.eigh instead of linalg.eig

src/pygenstability/constructors.py

@@ -31,9 +31,13 @@ def load_constructor(constructor, graph, **kwargs):
     """Load a constructor from its name, or as a custom Constructor class."""
     if isinstance(constructor, str):
         if graph is None:
-            raise Exception(f"No graph was provided with a generic constructor {constructor}")
+            raise Exception(
+                f"No graph was provided with a generic constructor {constructor}"
+            )
         try:
-            return getattr(sys.modules[__name__], f"constructor_{constructor}")(graph, **kwargs)
+            return getattr(sys.modules[__name__], f"constructor_{constructor}")(
+                graph, **kwargs
+            )
         except AttributeError as exc:
             raise Exception(f"Could not load constructor {constructor}") from exc
     if not isinstance(constructor, Constructor):
@@ -53,8 +57,8 @@ def limit(*args, **kwargs):
 
 
 def _compute_spectral_decomp(matrix):
-    """Solve eigenalue problem."""
-    lambdas, v = la.eig(matrix.toarray())
+    """Solve eigenalue problem for symmetric matrix."""
+    lambdas, v = la.eigh(matrix.toarray())
     vinv = la.inv(v.real)
     return lambdas.real, v.real, vinv
 
@@ -67,7 +71,9 @@ def _check_total_degree(degrees):
 
 def _get_spectral_gap(laplacian):
     """Compute spectral gap."""
-    spectral_gap = np.round(max(np.real(sp.linalg.eigs(laplacian, which="SM", k=2)[0])), 8)
+    spectral_gap = np.round(
+        max(np.real(sp.linalg.eigs(laplacian, which="SM", k=2)[0])), 8
+    )
     L.info("Spectral gap = 10^%s", np.around(np.log10(spectral_gap), 2))
     return spectral_gap
 
@@ -80,7 +86,9 @@ class Constructor:
     to return quality matrix, null model, and possible global shift.
     """
 
-    def __init__(self, graph, with_spectral_gap=False, exp_comp_mode="spectral", **kwargs):
+    def __init__(
+        self, graph, with_spectral_gap=False, exp_comp_mode="spectral", **kwargs
+    ):
         """The constructor calls the prepare method upon initialisation.
 
         Args:
@@ -105,7 +113,9 @@ def __init__(self, graph, with_spectral_gap=False, exp_comp_mode="spectral", **k
     def _get_exp(self, scale):
         """Compute matrix exponential at a given scale."""
         if self.exp_comp_mode == "expm":
-            exp = sp.linalg.expm(-scale * self.partial_quality_matrix.toarray().astype(_DTYPE))
+            exp = sp.linalg.expm(
+                -scale * self.partial_quality_matrix.toarray().astype(_DTYPE)
+            )
         if self.exp_comp_mode == "spectral":
             lambdas, v, vinv = self.spectral_decomp
             exp = v @ np.diag(np.exp(-scale * lambdas)) @ vinv
@@ -181,7 +191,9 @@ class constructor_continuous_combinatorial(Constructor):
     @_limit_numpy
     def prepare(self, **kwargs):
         """Prepare the constructor with non-scale dependent computations."""
-        laplacian, degrees = sp.csgraph.laplacian(self.graph, return_diag=True, normed=False)
+        laplacian, degrees = sp.csgraph.laplacian(
+            self.graph, return_diag=True, normed=False
+        )
         _check_total_degree(degrees)
         laplacian /= degrees.mean()
         pi = np.ones(self.graph.shape[0]) / self.graph.shape[0]
@@ -221,7 +233,9 @@ class constructor_continuous_normalized(Constructor):
     @_limit_numpy
     def prepare(self, **kwargs):
         """Prepare the constructor with non-scale dependent computations."""
-        laplacian, degrees = sp.csgraph.laplacian(self.graph, return_diag=True, normed=False)
+        laplacian, degrees = sp.csgraph.laplacian(
+            self.graph, return_diag=True, normed=False
+        )
         _check_total_degree(degrees)
         normed_laplacian = sp.diags(1.0 / degrees).dot(laplacian)
 
@@ -354,7 +368,9 @@ class constructor_directed(Constructor):
     @_limit_numpy
     def prepare(self, **kwargs):
         """Prepare the constructor with non-scale dependent computations."""
-        assert self.exp_comp_mode == "expm"
+        assert (
+            self.exp_comp_mode == "expm"
+        ), 'exp_comp_mode="expm" is required for "constructor_directed"'
 
         alpha = kwargs.get("alpha", 0.8)
         n_nodes = self.graph.shape[0]
@@ -366,11 +382,17 @@ def prepare(self, **kwargs):
 
         self.partial_quality_matrix = sp.csr_matrix(
             alpha * np.diag(dinv).dot(self.graph.toarray())
-            + ((1 - alpha) * np.diag(np.ones(n_nodes)) + np.diag(alpha * (dinv == 0.0))).dot(ones)
+            + (
+                (1 - alpha) * np.diag(np.ones(n_nodes)) + np.diag(alpha * (dinv == 0.0))
+            ).dot(ones)
             - np.eye(n_nodes)
         )
 
-        pi = abs(sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][:, 0])
+        pi = abs(
+            sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][
+                :, 0
+            ]
+        )
         pi /= pi.sum()
         self.partial_null_model = np.array([pi, pi])
 
@@ -424,9 +446,10 @@ def prepare(self, **kwargs):
 
             self.partial_quality_matrix = sp.csr_matrix(
                 alpha * np.diag(dinv).dot(self.graph.toarray())
-                + ((1 - alpha) * np.diag(np.ones(n_nodes)) + np.diag(alpha * (dinv == 0.0))).dot(
-                    ones
-                )
+                + (
+                    (1 - alpha) * np.diag(np.ones(n_nodes))
+                    + np.diag(alpha * (dinv == 0.0))
+                ).dot(ones)
                 - np.eye(n_nodes)
             )
 
@@ -435,7 +458,11 @@ def prepare(self, **kwargs):
                 sp.diags(dinv).dot(self.graph) - sp.diags(np.ones(n_nodes))
             )
 
-        pi = abs(sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][:, 0])
+        pi = abs(
+            sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][
+                :, 0
+            ]
+        )
         pi /= pi.sum()
         self.partial_null_model = np.array([pi, pi])