introvoyz041/TensorNetworks

1import numpy as np2 3import decomposition4 5from numpy.typing import NDArray6complex_arr_t = NDArray[np.complex128]7 8 9class MixedCanonical:10    def __init__(self, state: complex_arr_t, j: int, trunc_bond_d=None):11        self.j = j12        self.rank = int(np.log2(state.size))13 14        if self.rank - j < 1:15            raise Exception("Orthogonality center index is out of bounds")16 17        # we extend shape of our tensor for easier calc (we add two 1d legs,one on the left and one on the right)18        extension = np.ones(1)19        state = np.tensordot(extension, state, axes=0)20        state = np.tensordot(state, extension, axes=0)21 22        left_mps_list = []23        right_mps_list = []24 25        # svd decomposition with integrated bond dimension truncation26        def SVD(state0):27            u, s, vh = np.linalg.svd(state0, full_matrices=False)28            d = s.size29 30            # bond dimension truncation31            if trunc_bond_d is not None and d > trunc_bond_d:32                d = trunc_bond_d33                u = u[:, :d]34                vh = vh[:d, :]35                s = s[:d]36 37            return u, s, vh38 39        # calculation of the left canonical part of the MPS40        for _ in range(j):41            shape = state.shape42            state = np.reshape(state, (int(np.prod(shape[0:2])), int(np.prod(shape[2:]))))43 44            u, s, vh = SVD(state)45            bond_d = s.size46 47            u = np.reshape(u, (shape[0], 2, bond_d))48            state = np.diag(s) @ vh49            state = np.reshape(state, np.insert(shape[2:], 0, bond_d))50 51            left_mps_list.append(u)52 53        # calculation of the right canonical part of the MPS54        for _ in range(self.rank - j - 1):55            shape = state.shape56            state = np.reshape(state, (int(np.prod(shape[:-2])), int(np.prod(shape[-2:]))))57 58            u, s, vh = SVD(state)59            bond_d = s.size60 61            vh = np.reshape(vh, (bond_d, 2, shape[-1]))62            state = u @ np.diag(s)63            state = np.reshape(state, np.append(shape[:-2], bond_d))64 65            right_mps_list.insert(0, vh)66 67        self.left_part = left_mps_list68        self.right_part = right_mps_list69        self.ortho_center = state70 71    def is_left_canonical(self):72        return self.j == self.rank - 173 74    def is_right_canonical(self):75        return self.j == 076 77    def norm(self):78        return np.real(np.tensordot(self.ortho_center, np.conj(self.ortho_center), axes=((0, 1, 2), (0, 1, 2))))79 80    # method returns expectation value of provided operator at site j81    def ev_1site(self, operator: complex_arr_t) -> float:82        # from orthogonality center and its conjugate we create one 2-leg tensor83        contraction = np.tensordot(np.conj(self.ortho_center), self.ortho_center, axes=((0, 2), (0, 2)))84        # performing contraction of operator and this tensor85        contraction = np.tensordot(contraction, operator, axes=((0, 1), (0, 1)))86        return np.real(contraction)87 88    # method returns expectation value of provided operator at site j and j+189    def ev_2site(self, operator: complex_arr_t) -> float:90        if self.is_left_canonical():91            raise Exception("j is a last element of the MPS")92 93        next_tensor = self.right_part[0]94 95        # from mps tensors and its conjugate we create "tensor ring"96        # with 4 legs around operator tensor and then perform contraction97        left_part = np.tensordot(self.ortho_center, np.conj(self.ortho_center), axes=((0,), (0,)))98        right_part = np.tensordot(next_tensor, np.conj(next_tensor), axes=((2,), (2,)))99        ring = np.tensordot(left_part, right_part, axes=((1, 3), (0, 2)))100        ring = ring.transpose((0, 2, 1, 3))101 102        return np.real(np.tensordot(operator, ring, axes=((0, 1, 2, 3), (0, 1, 2, 3))))103 104    # method returns expectation value of the operator represented as MPO in our MPS state105    def ev_mpo(self, MPO) -> float:106        L, R = decomposition.mps_mpo_contraction(self, MPO)107        return np.real(np.tensordot(L[0], R[-1], axes=((0, 1, 2), (0, 1, 2))))108 109    def __getitem__(self, i):110        if 0 > i >= self.rank:111            Exception("Index out of bounds")112 113        if i < self.j:114            return self.left_part[i]115        if i == self.j:116            return self.ortho_center117        if i > self.j:118            return self.right_part[i - self.j - 1]119 120    def __setitem__(self, i, value):121        if 0 > i >= self.rank:122            Exception("Index out of bounds")123 124        if i < self.j:125            self.left_part[i] = value126        if i == self.j:127            self.ortho_center = value128        if i > self.j:129            self.right_part[i - self.j - 1] = value130 131    # with this method we transfer our orthogonality center to a new site using gauge transformations132    def change_j(self, j: int):133        if j == self.j:134            return135 136        old_j = self.j137        self.j = j138 139        if j > old_j:140            for _ in range(j - old_j):141                shape = self.ortho_center.shape142                ortho_center = np.reshape(self.ortho_center, (shape[0] * shape[1], shape[2]))143                q, r = np.linalg.qr(ortho_center)144 145                self.left_part.append(np.reshape(q, shape))146                self.ortho_center = np.tensordot(r, self.right_part.pop(0), axes=((1,), (0,)))147        else:148            for _ in range(old_j - j):149                shape = self.ortho_center.shape150                ortho_center = np.reshape(self.ortho_center, (shape[0], shape[1] * shape[2]))151                q, r = np.linalg.qr(ortho_center.transpose())152 153                self.right_part.insert(0, np.reshape(q.transpose(), shape))154                self.ortho_center = np.tensordot(self.left_part.pop(-1), r, axes=((2,), (1,)))155