diff --git a/sgw_tools/graph.py b/sgw_tools/graph.py index 12547fd..b195b76 100644 --- a/sgw_tools/graph.py +++ b/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): """ diff --git a/sgw_tools/util.py b/sgw_tools/util.py index f311919..fa4d1e6 100644 --- a/sgw_tools/util.py +++ b/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: