DMRG & Lanczos implemented (#error : diagonal matrix to mpo unstable)

tensor/matrix_product_operator.py

@@ -573,43 +573,52 @@ def retrieve(self, input_indices, output_indices):
 
 
     def left_orthonormalization(self, bond_shape=()):
-        for k in range(self.sites_number-1):
-            # print(k, k+1)
-            L = self.left_site_matricization(k)
-            R = self.right_site_matricization(k+1)
+        if self.decomposed:
+            for k in range(self.sites_number-1):
+                # print(k, k+1)
+                L = self.left_site_matricization(k)
+                R = self.right_site_matricization(k+1)
 
-            U, B = np.linalg.qr(L)
+                U, B = np.linalg.qr(L)
 
-            W = B @ R
+                W = B @ R
 
-            self.sites[k] = self.tensoricization(U, k)
-            self.sites[k+1] = self.tensoricization(W, k+1)
-
-        # L = self.left_site_matricization(-1)
-        # V, U = np.linalg.qr(L)
-        # self.sites[-1] = self.tensoricization(V, -1)
+                self.sites[k] = self.tensoricization(U, k)
+                self.sites[k+1] = self.tensoricization(W, k+1)
+            
+            # L = self.left_site_matricization(-1)
+            # V, U = np.linalg.qr(L)
+            # self.sites[-1] = self.tensoricization(V, -1)
 
-        self.orthonormalized = 'left'
+            self.orthonormalized = 'left'
+        else:
+            raise Exception("Cannot orthonormalize an undecomposed MatrixProductOperator")
+        return self
 
     def right_orthonormalization(self, bond_shape=()):
-        for k in range(self.sites_number-1, 0, -1):
-
-            R = self.right_site_matricization(k)
-            L = self.left_site_matricization(k-1)
-
-            V, U = np.linalg.qr(R.T)
-            # U, S, V = np.linalg.svd(R, full_matrices=False)
-            W = L @ U.T
-            # W = L @ (U * S)
+        if self.decomposed:
+            for k in range(self.sites_number-1, 0, -1):
+
+                R = self.right_site_matricization(k)
+                L = self.left_site_matricization(k-1)
+
+                V, U = np.linalg.qr(R.T)
+                # U, S, V = np.linalg.svd(R, full_matrices=False)
+                W = L @ U.T
+                # W = L @ (U * S)
 
-            self.sites[k] = self.tensoricization(V.T, k)
-            self.sites[k-1] = self.tensoricization(W, k-1)
+                self.sites[k] = self.tensoricization(V.T, k)
+                self.sites[k-1] = self.tensoricization(W, k-1)
 
-        # R = self.right_site_matricization(0)
-        # V, U = np.linalg.qr(R.T)
-        # self.sites[0] = self.tensoricization(V.T, 0)
+            # R = self.right_site_matricization(0)
+            # V, U = np.linalg.qr(R.T)
+            # self.sites[0] = self.tensoricization(V.T, 0)
 
-        self.orthonormalized = 'right'
+            self.orthonormalized = 'right'
+        else:
+            raise Exception("Cannot orthonormalize an undecomposed MatrixProductOperator")
+
+        return self
 
     def left_site_matricization(self, index):
         return self.left_matricization(self.sites[index], index)

test/test_dmrg.py

@@ -2,12 +2,41 @@
 from syngular.tensor import MatrixProductOperator
 from syngular.variational import Lanczos
 
+import numpy as np
+
 if __name__ == "__main__":
 
-    input_shape = (3,3,3,3,3,3)
-    output_shape = (3,3,3,3,3,3)
-    bond_shape = (2,2,2,2,2)
+    input_shape = (3,3,3,3) #,3,3)
+    output_shape = (3,3,3,3) #,3,3)
+    bond_shape = (2,2,2) #,2,2)
+
+    n = 3*3*3
+    m = 4
+    A = np.zeros(n**2).reshape((n,n))
+    np.fill_diagonal(A, list(range(n)))
+    B = A.reshape((3,3,3,3,3,3))
+    print('B', B)
+
+    mpo = MatrixProductOperator.random(input_shape=input_shape, output_shape=output_shape, bond_shape=bond_shape) # MatrixProductOperator(B, bond_shape=(3,3)).decompose()#
+
+    print('C', np.around(mpo.to_tensor(), decimals=2))
+
+    eigenstate = DMRG.solve(mpo, optimizer=Lanczos)
+    eigentensor = eigenstate.to_tensor()
+    eigentensor /= np.linalg.norm(eigentensor)
+
+    print(eigenstate.shape)
+    print(eigentensor.shape)
+
+    print(mpo.to_tensor().shape)
+    eigenvalues, eigenvectors = np.linalg.eig(mpo.to_tensor().reshape((3**4, 3**4)))
+    eigenvalues2, eigenvectors2 = np.linalg.eig(A)
+
+    eigenvalue = eigenvectors[0]
+    eigenvalue /= np.linalg.norm(eigenvalue)
 
-    mpo = MatrixProductOperator.random(input_shape=input_shape, output_shape=output_shape, bond_shape=bond_shape)
+    print(eigenvalues)
+    print(eigenvalues2)
 
-    DMRG.solve(mpo, optimizer=Lanczos)
+    print(eigenvalue - eigentensor.reshape((81,)))
+    # print(eigenvalues[0] - eigenstate.to_tensor())

test/test_mpo.py

@@ -30,6 +30,7 @@ def test_apply():
     print(W)
     X = Z @ W
     print(X[(1,1),(1,1)])
+    print((z * w)[1,1,1,1])
 
     Z.apply(W, [0,1])
 
@@ -292,16 +293,37 @@ def test_to_tensor():
     print(f"> Execution time : {end-start:.8f}sec")
     print("-------------------------------------------------------")
 
+def test_compare():
+    y = np.ones(2**4).reshape((4,4))
+    np.fill_diagonal(y, list(range(1,5)))
+    y = y.reshape((2,2,2,2))
+    Y = MatrixProductOperator(y, bond_shape=(2,)).decompose().right_orthonormalization()
 
+    print(y)
 
-test_decompose()
-test_from_sites()
-test_random()
-# test_to_tensor() -- very expensive ~1min
-test_add()
-test_mul()
-test_matmul()
-test_dot()
-test_compress()
-test_orthogonality()
-test_apply()
+    print(np.around(Y.to_tensor(), decimals=2))
+
+    print("---------- Compare with to_tensor() ----------")
+    start = time.time()
+
+    diff = np.abs(y - Y.to_tensor())
+    print('Min', np.min(diff))
+    print('Max', np.max(diff))
+
+    end = time.time()
+    print(f"> Execution time : {end-start:.8f}sec")
+    print("-------------------------------------------------------")
+
+
+# test_decompose()
+# test_from_sites()
+# test_random()
+# # test_to_tensor() -- very expensive ~1min
+# test_add()
+# test_mul()
+# test_matmul()
+# test_dot()
+# test_compress()
+# test_orthogonality()
+# test_apply()
+test_compare()

variational/dmrg.py

@@ -1,3 +1,4 @@
+from hashlib import new
 from syngular.tensor import MatrixProductOperator, MatrixProductState
 from syngular.variational.optimizer import Optimizer, Lanczos
 
@@ -111,6 +112,8 @@ def __left_blocks(operator: MatrixProductOperator, state: MatrixProductState):
     '''
     def __right_sweep(operator: MatrixProductOperator, state: MatrixProductState, left_blocks: list, right_blocks: list, optimizer: Type[Optimizer]):
         n = operator.sites_number
+
+        init_vec = None
 
         for k in range(n-1):
             print(f"Step {k+1}/{n-1}")
@@ -146,9 +149,11 @@ def __right_sweep(operator: MatrixProductOperator, state: MatrixProductState, le
                     [11, 7, 9, 10]
                 )
 
-            print(f"Contracted state shape [idx {k}]", merge.shape)
+            print(f"Contracted state shape [idx {k}]", site.shape)
 
-            nsite = optimizer.fit(site)
+            _, nsite = optimizer.fit(site, init_vec=init_vec)
+            init_vec = nsite[0]
+            nsite = np.reshape(nsite, newshape=site.shape)
             state.sites[k], state.sites[k+1] = DMRG.__restore(state=state, site=nsite, index=k)
 
             # print("Shape end step", [site.shape for site in state.sites])
@@ -190,9 +195,11 @@ def __left_sweep(operator: MatrixProductOperator, state: MatrixProductState, lef
                     [11, 7, 9, 10]
                 )
 
-            print(f"Contracted state shape [idx {k}]", merge.shape)
+            print(f"Contracted state shape [idx {k}]", site.shape)
 
-            nsite = optimizer.fit(site)
+            _, nsite = optimizer.fit(site)
+            print(nsite.shape, site.shape)
+            nsite = np.reshape(nsite, newshape=site.shape)
             state.sites[k], state.sites[k+1] = DMRG.__restore(state=state, site=nsite, index=k)
 
 

variational/optimizer.py

@@ -10,19 +10,26 @@ def fit(tensor):
 class Lanczos(Optimizer):
 
     @staticmethod
-    def fit(tensor, iteration=100):
+    def fit(tensor, init_vec = None, iteration=100, ):
         n = tensor.shape[0]*tensor.shape[1] if len(tensor.shape) > 2 else tensor.shape[0]
+        l = tensor.shape[2]*tensor.shape[3] if len(tensor.shape) > 2 else tensor.shape[1]
         m = min(n, iteration)
 
         A = np.reshape(tensor, newshape=(n, -1)) if len(tensor.shape) > 2 else tensor
-
         T = np.zeros((m, m))
-        V = np.zeros((m, n))
+        V = np.zeros((m, l))
+
+        if A.shape[0] != A.shape[1]:
+            A = np.dot(A.conj().T, A)
 
         # First iteration
 
-        v = np.random.rand(n)
-        v /= np.linalg.norm(v)
+
+        if init_vec is None:
+            v = np.random.rand(l)
+            v /= np.linalg.norm(v)
+        else:
+            v = init_vec[:l]
         V[0, :] = v
 
         w = np.dot(A, v)