Improve eigen-decomposition accuracy for "adjacency" Laplacian.

sgw_tools/graph.py

@@ -20,10 +20,11 @@ def compute_fourier_basis(self, recompute=False, spectrum_only=False):
                                 'large matrix ({0} x {0}) may take some '
                                 'time.'.format(self.N))
 
+        e, U = self._eigL(spectrum_only=spectrum_only)
         if spectrum_only:
-            self._e = eigh(self.L.toarray(order='F'), eigvals_only=True, overwrite_a=True, driver='ev')
+            self._e = e
         else:
-            self._e, self._U = eigh(self.L.toarray(order='F'), overwrite_a=True, driver='evr')
+            self._e, self._U = e, U
             self._mu = np.max(np.abs(self._U))
 
         self._e[np.isclose(self._e, 0)] = 0
@@ -43,6 +44,9 @@ def compute_fourier_basis(self, recompute=False, spectrum_only=False):
         assert e_max <= e_bound, "Largest eigenvalue was {} but upper bound is {}".format(e_max, e_bound)
         self._lmax = e_max
 
+    def _eigL(self, spectrum_only):
+        return util.eigh(self.L, spectrum_only=spectrum_only)
+
     def _get_upper_bound(self):
         return np.inf
 
@@ -428,6 +432,16 @@ def dwq(self):
             self._Wq, self._dwq = util.magneticAdjacencyMatrix(self, self.q)
         return self._dwq
 
+    def _eigL(self, spectrum_only):
+        if self.lap_type == "adjacency":
+            # more accurate to decompose the adjacency matrix itself then shift the result
+            e, U = util.eigh(self.W, spectrum_only=spectrum_only)
+            e = np.flip(1 - e/np.max(e))
+            U = np.flip(U, axis=1) if U else None
+            return e, U
+        else:
+            return super()._eigL(spectrum_only=spectrum_only)
+
     @property
     def lmin(self):
         """

sgw_tools/util.py

@@ -64,6 +64,13 @@ def extract_components(G):
     return subgraphs
 
 
+def eigh(M, spectrum_only=False):
+    if spectrum_only:
+        return eigh(M.toarray(order='F'), eigvals_only=True, overwrite_a=True, driver='ev')
+    else:
+        return eigh(M.toarray(order='F'), overwrite_a=True, driver='evr')
+
+
 def _estimate_lmin(G, maxiter):
     N = G.N
     if N > 100: