BEB: replace SVD by eigendecomposition

Author: DerWeh

gftool/beb.py

@@ -4,6 +4,7 @@
 
 The implementation is based on a SVD of the `hopping` matrix,
 which is the dimensionless scaling of the hopping of the components. [weh2021]_
+However, we use the unitary eigendecomposition instead of the SVD.
 
 Physical quantities
 -------------------
@@ -75,17 +76,19 @@
 
 """
 # pylint: disable=too-many-locals
+from __future__ import annotations
+
 import logging
 
 from typing import Callable, NamedTuple
 from functools import partial
 
 import numpy as np
 
-from scipy import optimize
 from numpy import newaxis
+from scipy import optimize
 
-from gftool.matrix import decompose_mat
+from gftool.matrix import decompose_her, decompose_mat, UDecomposition
 
 LOGGER = logging.getLogger(__name__)
 
@@ -101,7 +104,7 @@ class SVD(NamedTuple):
     s: np.ndarray
     vh: np.ndarray
 
-    def truncate(self, rcond=None) -> "SVD":
+    def truncate(self, rcond=None) -> SVD:
         """Return the truncated singular values decomposition.
 
         Singular values smaller than `rcond` times the largest singular values
@@ -145,6 +148,55 @@ def partition(self, return_sqrts=False):
         return us, svh
 
 
+class SpecDec(UDecomposition):  # pylint: disable=too-many-ancestors
+    """SVD like spectral decomposition.
+
+    Works only for N×N matrices unlike the `UDecomposition` base class.
+    """
+
+    def truncate(self, rcond=None) -> SpecDec:
+        """Return the truncated spectral decomposition.
+
+        Singular values smaller than `rcond` times the largest singular values
+        are discarded.
+
+        Parameters
+        ----------
+        rcond : float, rcond
+            Cut-off ratio for small singular values.
+
+        Returns
+        -------
+        truncated_svd : SpecDec
+            The truncates the spectral decomposition discarding small singular values.
+
+        """
+        if rcond is None:
+            rcond = np.finfo(self.eig.dtype).eps * max(self.u.shape[-2:])
+        max_eig = np.max(abs(self.eig), axis=-1)
+        significant = abs(self.eig) > max_eig*rcond
+        return self.__class__(rv=self.rv[..., :, significant], eig=self.eig[..., significant],
+                              rv_inv=self.rv_inv[..., significant, :])
+
+    @property
+    def is_trunacted(self) -> bool:
+        """Check if SVD of square matrix is truncated/compact or full."""
+        ushape, uhshape = self.u.shape, self.uh.shape
+        return not ushape[-2] == ushape[-1] == uhshape[-2]
+
+    def partition(self, return_sqrts=False):
+        """Symmetrically partition the spectral decomposition as `u * eig**0.5, eig**0.5 * uh`.
+
+        If `return_sqrts` then `us, np.sqrt(s), suh` is returned,
+        else only `us, suh` is returned (default: False).
+        """
+        sqrt_eig = np.emath.sqrt(self.eig)
+        us, suh = self.u * sqrt_eig[..., newaxis, :], sqrt_eig[..., :, newaxis] * self.uh
+        if return_sqrts:
+            return us, sqrt_eig, suh
+        return us, suh
+
+
 def gf_loc_z(z, self_beb_z, hopping, hilbert_trafo: Callable[[complex], complex],
              diag=True, rcond=None):
     """Calculate average local Green's function matrix in components.
@@ -180,33 +232,33 @@ def gf_loc_z(z, self_beb_z, hopping, hilbert_trafo: Callable[[complex], complex]
     solve_root
 
     """
-    hopping_svd = SVD(*np.linalg.svd(hopping, hermitian=True))
-    LOGGER.info('hopping singular values %s', hopping_svd.s)
-    hopping_svd = hopping_svd.truncate(rcond)
-    LOGGER.info('Keeping %s (out of %s)', hopping_svd.s.shape[-1], hopping_svd.vh.shape[-1])
+    hopping_dec = SpecDec(*decompose_her(hopping))
+    LOGGER.info('hopping singular values %s', hopping_dec.s)
+    hopping_dec = hopping_dec.truncate(rcond)
+    LOGGER.info('Keeping %s (out of %s)', hopping_dec.s.shape[-1], hopping_dec.uh.shape[-1])
     kind = 'diag' if diag else 'full'
 
     eye = np.eye(*hopping.shape)
-    us, sqrt_s, svh = hopping_svd.partition(return_sqrts=True)
+    us, sqrt_s, suh = hopping_dec.partition(return_sqrts=True)
     # [..., newaxis]*eye add matrix axis
     z_m_self = z[..., newaxis, newaxis]*eye - self_beb_z
     z_m_self_inv = np.asfortranarray(np.linalg.inv(z_m_self))
-    dec = decompose_mat(svh @ z_m_self_inv @ us)
+    dec = decompose_mat(suh @ z_m_self_inv @ us)
     diag_inv = 1. / dec.eig
-    if not hopping_svd.is_trunacted:
-        svh_inv = transpose(hopping_svd.vh).conj() / sqrt_s[..., newaxis, :]
-        us_inv = transpose(hopping_svd.u).conj() / sqrt_s[..., :, newaxis]
+    if not hopping_dec.is_trunacted:
+        svh_inv = transpose(hopping_dec.uh).conj() / sqrt_s[..., newaxis, :]
+        us_inv = transpose(hopping_dec.u).conj() / sqrt_s[..., :, newaxis]
         dec.rv = svh_inv @ np.asfortranarray(dec.rv)
         dec.rv_inv = np.asfortranarray(dec.rv_inv) @ us_inv
         return dec.reconstruct(hilbert_trafo(diag_inv), kind=kind)
 
     dec.rv = z_m_self_inv @ us @ np.asfortranarray(dec.rv)
-    dec.rv_inv = np.asfortranarray(dec.rv_inv) @ svh @ z_m_self_inv
+    dec.rv_inv = np.asfortranarray(dec.rv_inv) @ suh @ z_m_self_inv
     correction = dec.reconstruct((diag_inv*hilbert_trafo(diag_inv) - 1) * diag_inv, kind=kind)
     return (diagonal(z_m_self_inv) if diag else z_m_self_inv) + correction
 
 
-def self_root_eq(self_beb_z, z, e_onsite, concentration, hopping_svd: SVD,
+def self_root_eq(self_beb_z, z, e_onsite, concentration, hopping_dec: SpecDec,
                  hilbert_trafo: Callable[[complex], complex]):
     """Root equation r(Σ)=0 for BEB.
 
@@ -220,7 +272,7 @@ def self_root_eq(self_beb_z, z, e_onsite, concentration, hopping_svd: SVD,
         On-site energy of the components.
     concentration : (..., N_cmpt) float array_like
         Concentration of the different components.
-    hopping_svd : SVD
+    hopping_dec : SVD
         Compact SVD decomposition of the (N_cmpt, N_cmpt) hopping matrix in the
         components.
     hilbert_trafo : Callable[[complex], complex]
@@ -240,14 +292,14 @@ def self_root_eq(self_beb_z, z, e_onsite, concentration, hopping_svd: SVD,
     eye = np.eye(e_onsite.shape[-1])  # [..., newaxis]*eye adds matrix axis
     z_m_self = z[..., newaxis, newaxis]*eye - self_beb_z
     # split symmetrically
-    us, svh = hopping_svd.partition()
+    us, suh = hopping_dec.partition()
     # matrix-products are faster if larger arrays are in Fortran order
     z_m_self_inv = np.asfortranarray(np.linalg.inv(z_m_self))
-    dec = decompose_mat(svh @ z_m_self_inv @ us)
+    dec = decompose_mat(suh @ z_m_self_inv @ us)
     dec.rv = us @ np.asfortranarray(dec.rv)
-    dec.rv_inv = np.asfortranarray(dec.rv_inv) @ svh
+    dec.rv_inv = np.asfortranarray(dec.rv_inv) @ suh
     diag_inv = 1. / dec.eig
-    if not hopping_svd.is_trunacted:
+    if not hopping_dec.is_trunacted:
         gf_loc_inv = dec.reconstruct(1./hilbert_trafo(diag_inv), kind='full')
     else:
         gf_loc_inv = z_m_self + dec.reconstruct(1./hilbert_trafo(diag_inv) - diag_inv, kind='full')
@@ -362,12 +414,12 @@ def solve_root(z, e_onsite, concentration, hopping, hilbert_trafo: Callable[[com
     >>> plt.show()
 
     """
-    hopping_svd = SVD(*np.linalg.svd(hopping, hermitian=True))
-    LOGGER.info('hopping singular values %s', hopping_svd.s)
-    hopping_svd = hopping_svd.truncate(rcond)
-    LOGGER.info('Keeping %s (out of %s)', hopping_svd.s.shape[-1], hopping_svd.vh.shape[-1])
+    hopping_dec = SpecDec(*decompose_her(hopping))
+    LOGGER.info('hopping singular values %s', hopping_dec.s)
+    hopping_dec = hopping_dec.truncate(rcond)
+    LOGGER.info('Keeping %s (out of %s)', hopping_dec.s.shape[-1], hopping_dec.uh.shape[-1])
     self_root_part = partial(self_root_eq, z=z, e_onsite=e_onsite, concentration=concentration,
-                             hopping_svd=hopping_svd, hilbert_trafo=hilbert_trafo)
+                             hopping_dec=hopping_dec, hilbert_trafo=hilbert_trafo)
     if self_beb_z0 is None:
         self_beb_z0 = np.zeros(hopping.shape, dtype=complex)
         # experience shows that a single fixed_point is a good starting point
@@ -390,13 +442,13 @@ def solve_root(z, e_onsite, concentration, hopping, hilbert_trafo: Callable[[com
         root_kwds['callback'] = lambda x, f: LOGGER.debug('Residue: %s', np.linalg.norm(f))
 
     sol = optimize.root(root_eq, x0=self_beb_z0, **root_kwds)
-    LOGGER.debug("BEB self-energy root found after %s iterations.", sol.nit)
+    LOGGER.info("BEB self-energy root found after %s iterations.", sol.nit)
 
     if LOGGER.isEnabledFor(logging.INFO):
         # check condition number in matrix diagonalization to make sure it is well defined
-        us, svh = hopping_svd.partition()  # pylint: disable=unbalanced-tuple-unpacking
+        us, suh = hopping_dec.partition()  # pylint: disable=unbalanced-tuple-unpacking
         z_m_self = z[..., newaxis, newaxis]*np.eye(*hopping.shape) - sol.x
-        dec = decompose_mat(svh @ np.linalg.inv(z_m_self) @ us)
+        dec = decompose_mat(suh @ np.linalg.inv(z_m_self) @ us)
         max_cond = np.max(np.linalg.cond(dec.rv))
         LOGGER.info("Maximal coordination number for diagonalization: %s", max_cond)