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,14 @@ 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)
+ return 1 - e/np.max(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(self, 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: