From 861a78b529ede29fc9346fe960464cec3981edc7 Mon Sep 17 00:00:00 2001 From: meandmytram Date: Mon, 30 Oct 2023 01:05:07 -0400 Subject: [PATCH] return residual sprectra --- mdopt/utils/utils.py | 35 ++++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/mdopt/utils/utils.py b/mdopt/utils/utils.py index 549f3c98..227094a6 100644 --- a/mdopt/utils/utils.py +++ b/mdopt/utils/utils.py @@ -1,6 +1,6 @@ """This module contains miscellaneous utilities.""" -from typing import Tuple, List, cast +from typing import Tuple, Optional, List, cast from itertools import chain import numpy as np import scipy @@ -12,7 +12,8 @@ def svd( cut: np.float32 = np.float32(1e-12), chi_max: int = int(1e4), renormalise: bool = False, -) -> Tuple[np.ndarray, list, np.ndarray]: + return_residual_spectrum: bool = False, +) -> Tuple[np.ndarray, list, np.ndarray, Optional[List[float]]]: """ Performs Singular Value Decomposition with different features. @@ -26,6 +27,9 @@ def svd( Maximum number of singular values to keep. renormalise : bool Whether to renormalise the singular value spectrum after the cut. + return_residual_spectrum : bool + Whether to return the residual spectrum, i.e., + the singular values we have thrown away. Returns ------- @@ -35,11 +39,8 @@ def svd( The singular values, sorted in non-increasing order. v_r : np.ndarray Unitary matrix having right singular vectors as rows. - - Raises - ------ - ValueError - If the ``np.ndarray`` provided is not two-dimensional. + residual_spectrum : Optional[List[float]] + The discarded singular values. Only returned if `return_residual_spectrum` is True. """ if len(mat.shape) != 2: @@ -64,12 +65,22 @@ def svd( max_num = min(chi_max, np.sum(singular_values > cut)) ind = np.argsort(singular_values)[::-1][:max_num] + residual_spectrum = singular_values[singular_values <= cut].tolist() + u_l, singular_values, v_r = u_l[:, ind], singular_values[ind], v_r[ind, :] if renormalise: singular_values /= np.linalg.norm(singular_values) - return np.asarray(u_l), cast(list, singular_values), np.asarray(v_r) + if return_residual_spectrum: + return ( + np.asarray(u_l), + cast(list, singular_values), + np.asarray(v_r), + residual_spectrum, + ) + + return np.asarray(u_l), cast(list, singular_values), np.asarray(v_r), None def kron_tensors( @@ -209,7 +220,7 @@ def split_two_site_tensor( tensor = tensor.reshape((chi_v_l * phys_l, phys_r * chi_v_r)) if strategy == "svd": - u_l, singular_values, v_r = svd( + u_l, singular_values, v_r, _ = svd( tensor, cut=cut, chi_max=chi_max, renormalise=renormalise ) chi_v_cut = len(singular_values) @@ -486,7 +497,7 @@ def mpo_from_matrix( mpo = [] mps_dim = phys_dim**2 matrix = matrix.reshape((-1, mps_dim)) - matrix, singular_values, v_r = svd(matrix, chi_max=chi_max, renormalise=False) + matrix, singular_values, v_r, _ = svd(matrix, chi_max=chi_max, renormalise=False) if not orthogonalise: v_r = np.matmul(np.diag(np.sqrt(singular_values)), v_r) mpo.append(np.expand_dims(v_r, -1)) @@ -496,7 +507,9 @@ def mpo_from_matrix( matrix = np.matmul(matrix, np.diag(singular_values)) bond_dim = matrix.shape[-1] matrix = matrix.reshape((-1, mps_dim * bond_dim)) - matrix, singular_values, v_r = svd(matrix, chi_max=chi_max, renormalise=False) + matrix, singular_values, v_r, _ = svd( + matrix, chi_max=chi_max, renormalise=False + ) if not orthogonalise: v_r = np.matmul(np.diag(np.sqrt(singular_values)), v_r) v_r = v_r.reshape((-1, mps_dim, bond_dim))