Right-canonical, syn utils, DMRG, starting variational algorithm

__init__.py

@@ -0,0 +1 @@
+from syngular.tensor.utils import mul

tensor/matrix_product_operator.py

@@ -4,6 +4,8 @@
 import gc
 import itertools
 from typing import Tuple, List, Union
+from cv2 import eigen
+from django.forms import PasswordInput
 from numpy.core.numeric import outer
 
 from opt_einsum import contract
@@ -18,6 +20,8 @@ class MatrixProductOperator:
     COMPRESS = True
     DECOMPOSE = False
 
+    MATMUL_MODE = "standard"
+
     def __init__(self, tensor=None, bond_shape: Union[Union[Tuple, List], np.ndarray]=(), verbose=0) -> None:
 
         self.parameters_number = 0
@@ -171,32 +175,118 @@ def __mul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
     def __matmul__(self, mp: Union[MatrixProductOperator, MatrixProductState]):
         min_bond = min(min(self.bond_shape), min(mp.bond_shape))
         max_bond = max(max(self.bond_shape), max(mp.bond_shape))
-        # print("min", min_bond)f
+        # print("min", min_bond)
         sites = []
 
         if isinstance(mp, MatrixProductState):
-            for idx in range(self.sites_number):
-                site = contract(mp.sites[idx], [1,2,3], self.sites[idx], [4,2,6,5], [1,4,6,3,5])
+            if MatrixProductOperator.MATMUL_MODE == "standard":
+                for idx in range(self.sites_number):
+                    site = contract(mp.sites[idx], [1,2,3], self.sites[idx], [4,2,6,5], [1,4,6,3,5])
+
+                    site = site.reshape((
+                        self.shape[idx][0]*mp.shape[idx][0],
+                        self.shape[idx][2],
+                        self.shape[idx][3]*mp.shape[idx][2]
+                    ))
+                    sites.append(site)
+                return MatrixProductState.from_sites(sites) >> min_bond
+            elif MatrixProductOperator.MATMUL_MODE == "opti":
+                blocks = [None for _ in range(self.sites_number)]
+                basis = [None for _ in range(self.sites_number)]
+                computes = [None for _ in range(self.sites_number)]
+
+                for idx in range(self.sites_number-1):
+                    state_site = mp.sites[idx]
+                    op_site = self.sites[idx]
+
+                    struct = [
+                        state_site, [1, 2, 3],
+                        op_site, [4, 2, 5, 6],
+                        op_site, [7, 8, 9, 10],
+                        state_site, [11, 9, 12]
+                    ]
+                    if idx != 0: struct += [blocks[idx-1], [1, 4, 7, 11]]
+                    struct += [[3, 6, 10, 12]]
+
+                    block = contract(*struct)
+                    blocks[idx] = block
+
+                state_site = mp.sites[-1]
+                op_site = self.sites[-1]
+
+                reduced_density_matrix = contract(
+                    blocks[-1], [1,4,7,11],
+
+                    state_site, [1, 2, 3],
+                    op_site, [4, 2, 5, 6],
+                    op_site, [7, 8, 9, 10],
+                    state_site, [11, 9, 12],
 
-                site = site.reshape((
-                    self.shape[idx][0]*mp.shape[idx][0],
-                    self.shape[idx][2],
-                    self.shape[idx][3]*mp.shape[idx][2]
-                ))
-                sites.append(site)
-            return MatrixProductState.from_sites(sites) >> min_bond
-        elif isinstance(mp, MatrixProductOperator):
-            for idx in range(self.sites_number):
-                site = contract(self.sites[idx], [1,2,3,4], mp.sites[idx], [5,3,6,7], [1,5,2,6,4,7])
+                    [5, 8]
+                )
 
-                site = site.reshape((
-                    self.shape[idx][0]*mp.shape[idx][0],
-                    self.shape[idx][1],
-                    self.shape[idx][2],
-                    self.shape[idx][3]*mp.shape[idx][3]
-                ))
-                sites.append(site)
-            return MatrixProductOperator.from_sites(sites) >> min_bond
+                eigenvalues, eigenvectors = np.linalg.eigh(reduced_density_matrix)
+                eigenvalues = eigenvalues[:min_bond]
+                eigenvectors = eigenvectors[:, :min_bond]
+
+                basis[-1] = eigenvectors
+                computes[-1] = contract(
+                    state_site, [1, 2, 3],
+                    op_site, [4, 2, 5, 6],
+                    np.conjugate(basis[-1]), [5, 7],
+                    [1, 4, 7]
+                )
+
+                penult_state_site = mp.sites[-2]
+                penult_op_site = self.sites[-2]
+
+                penult_reduced_density_matrix = contract(
+                    blocks[-2], [1,4,7,11],
+
+                    penult_state_site, [1, 2, 3],
+                    penult_op_site, [4, 2, 5, 6],
+                    penult_op_site, [7, 8, 9, 10],
+                    penult_state_site, [11, 9, 12],
+
+                    # state_site, [3, 13, 14],
+                    # op_site, [6, 13, 15, 16],
+                    # op_site, [10, 17, 18, 19],
+                    # state_site, [12, 18, 20],
+
+                    computes[-1], [3, 6, 15],
+                    np.conjugate(computes[-1]), [12, 10, 17],
+
+                    [5, 17, 10, 8]
+                )
+
+            elif MatrixProductOperator.MATMUL_MODE == "variational":
+
+                left_blocks = [None for _ in range(self.sites_number)]
+                right_blocks = [None for _ in range(self.sites_number)]
+
+                new_state = MatrixProductState.copy(mp)
+
+                for idx in range(self.sites_number-1):
+                    state_site = mp.sites[idx]
+                    op_site = self.sites[idx]
+
+                    # ...
+
+            elif MatrixProductOperator.MATMUL_MODE == "fit-up":
+                pass
+
+        elif isinstance(mp, MatrixProductOperator):
+                for idx in range(self.sites_number):
+                    site = contract(self.sites[idx], [1,2,3,4], mp.sites[idx], [5,3,6,7], [1,5,2,6,4,7])
+
+                    site = site.reshape((
+                        self.shape[idx][0]*mp.shape[idx][0],
+                        self.shape[idx][1],
+                        self.shape[idx][2],
+                        self.shape[idx][3]*mp.shape[idx][3]
+                    ))
+                    sites.append(site)
+                return MatrixProductOperator.from_sites(sites) >> min_bond
 
     def __mod__(self, mode):
         if mode == 'L' or mode == 'left' or mode == 'l' or mode == 'left-canonical':
@@ -349,6 +439,8 @@ def decompose(self, mode="left"):
                     Q = Q[:,:rank_right]
                     R = R[:rank_right, :]
 
+                    print(Q.shape, R.shape)
+
                     self.sites[k] = self.tensoricization(Q, k)
                     T = R
 
@@ -562,6 +654,12 @@ def apply(self, operator, indices):
     def split():
         pass
 
+    def transpose(self):
+        sites = []
+        for idx in range(self.sites_number):
+            sites.append(self.sites[idx].transpose(self.sites[idx], (0,2,1,3)))
+        return MatrixProductOperator.from_sites(sites)
+
     def retrieve(self, input_indices, output_indices):
         einsum_structure = []
 

tensor/matrix_product_state.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import gc
 import itertools
 from typing import Tuple, List, Union
@@ -7,6 +9,8 @@
 
 from opt_einsum import contract
 
+# from syngular.tensor.matrix_product_operator import MatrixProductOperator
+
 class MatrixProductState:
 
     """__init__ method of MatrixProductState
@@ -41,6 +45,8 @@ def __init__(self, tensor=None, bond_shape: Union[Union[Tuple, List], np.ndarray
             self.sites_number = len(self.bond_shape)+1
             self.sites = [None] * self.sites_number
 
+            self.norm = None
+
             if self.order != self.sites_number:
                 raise Exception("dimensions of bond indices do not match order - 1")
 
@@ -228,14 +234,26 @@ def zeros(input_shape, bond_shape):
             sites.append(np.zeros(shape=shape[idx]))
         return MatrixProductState.from_sites(sites)
 
-    """norm method of MatrixProductState
+    """dot method of MatrixProductState
     Compute the norm of the MatrixProductState
 
     @rtype: float
     @returns: the MatrixProductState norm
     """
-    def norm(self):
-        return self | self
+    def dot(self):
+        return np.sqrt(self | self)
+
+    """normalize method of MatrixProductState
+    Normalize the MatrixProductState
+
+    @rtype: MatrixProductState
+    @returns: self
+    """
+    def normalize(self):
+        norm = np.linalg.norm(self.sites[-1].reshape((self.sites[-1].shape[0], -1)))
+        self.sites[-1] /= norm
+
+        return self
 
     """to_tensor method of MatrixProductState
     Retrieve the based tensor of the MatrixProductState
@@ -255,9 +273,26 @@ def to_tensor(self):
         return tensor
 
 
+    """copy method of MatrixProductState
+    Copy the MatrixProductState
+
+    @rtype: MatrixProductState
+    @returns: a copy of the MPS
+    """
     def copy(self):
         return MatrixProductState.from_sites(self.sites)
 
+
+    """decompose method of MatrixProductState
+    Decompose (left/right normalization)
+
+
+    @type mode: "left" | "right"
+    @param mode: mode of orthonormalization (left/rigth)
+
+    @rtype: MatrixProductState
+    @returns: self
+    """
     def decompose(self, mode="left"):
         if not self.decomposed:
             if mode == "left":
@@ -269,9 +304,8 @@ def decompose(self, mode="left"):
                 for k in range(self.sites_number-1):
 
                     L = self.left_matricization(T, index=k)
-
                     Q, R = np.linalg.qr(L, mode="complete")
-                    # print(Q.shape, R.shape)
+
                     rank = self.shape[k][2]
                     Q = Q[:,:rank]
                     R = R[:rank, :]
@@ -288,15 +322,43 @@ def decompose(self, mode="left"):
                 gc.collect()
 
             elif mode == "right":
-                current_matrix = self.tensor.T
-                current_shape = self.shape[0][0]
-                current_rank = 1
 
-                self.decompose(mode="left")
+                n = self.sites_number-1
+                T = self.tensor
+
+                for k in range(self.sites_number-1,0,-1):
+
+                    L = self.right_matricization(T, index=k).T
+                    Q, R = np.linalg.qr(L, mode="complete")
+
+                    rank = self.shape[k][0]
+                    Q = Q[:,:rank].T
+                    R = R[:rank, :].T
+
+                    self.sites[k] = self.tensoricization(Q, k)
+                    T = R
+
+                    self.parameters_number += np.prod(self.sites[k].shape)
+
+                self.sites[0] = self.tensoricization(T, 0)
+                self.parameters_number += np.prod(self.sites[0].shape)
+
+                del self.tensor
+                gc.collect()
 
             self.decomposed = True
         return self
 
+    """compress method of MatrixProductState
+    Compression of the MatrixProductStates QR algorithm
+
+
+    @type dim: int
+    @param dim: dimension of the compression < current bond dimension
+
+    @rtype: MatrixProductState
+    @returns: compressed MatrixProductState
+    """
     def compress(self, dim: int, mode="left", strict=False):
         n = self.sites_number-1
         parameters_number = 0
@@ -405,7 +467,24 @@ def compress(self, dim: int, mode="left", strict=False):
 
             return that
 
-    def apply(self, operator, index, strict=True):
+    """apply method of MatrixProductState
+    Apply an operator (for now MatrixProductOperator) to the MatrixProductState
+
+
+    @type operator: MatrixProductOperator
+    @param operator: the MPO operator to apply
+    
+    @type index: int
+    @param index: the MPO operator to apply
+
+    @type mode: "compress" | "normal" | "density"
+    @param mode: Mode of application of the operator by compressing the bond dimensions or not 
+                 and using an optimized algorithm
+
+    @rtype: MatrixProductState
+    @returns: the transformation of the MatrixProductState by the MatrixProductOperator
+    """
+    def apply(self, operator: MatrixProductOperator, index, strict=True, mode="compress"):
         if strict: that = self.copy()
         # print("> Apply")
         m = len(operator.shape) // 2
@@ -454,6 +533,16 @@ def apply(self, operator, index, strict=True):
 
         return that
 
+    """retrieve method of MatrixProductState
+    Fixe physical indices to retrieve the real/complex value of the tensor  
+
+
+    @type indices: List[int]
+    @param indices: physicial indices list to be fixed
+
+    @rtype: float
+    @returns: the selected coefficient of the MatrixProductState 
+    """
     def retrieve(self, indices):
         einsum_structure = []
         for idx in range(self.sites_number):