From 9a79ca01e42ba0bb4f01120fbcd39db77ecd51ad Mon Sep 17 00:00:00 2001 From: Antoine Debouchage Date: Sat, 12 Mar 2022 17:00:35 +0100 Subject: [PATCH] Right-canonical, syn utils, DMRG, starting variational algorithm --- __init__.py | 1 + tensor/matrix_product_operator.py | 140 +++++++++++++++++++++++++----- tensor/matrix_product_state.py | 109 ++++++++++++++++++++--- tensor/utils.py | 100 +++++++++++++++++++++ test/test_dmrg.py | 88 ++++++++++++++++--- test/test_mpo.py | 16 +++- test/test_mps.py | 65 +++++++++++++- test/test_syn.py | 30 +++++++ variational/dmrg.py | 20 ++++- variational/optimizer.py | 11 ++- 10 files changed, 530 insertions(+), 50 deletions(-) create mode 100644 __init__.py create mode 100644 tensor/utils.py create mode 100644 test/test_syn.py diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..4609bd7 --- /dev/null +++ b/__init__.py @@ -0,0 +1 @@ +from syngular.tensor.utils import mul \ No newline at end of file diff --git a/tensor/matrix_product_operator.py b/tensor/matrix_product_operator.py index 47d90b2..731b81d 100644 --- a/tensor/matrix_product_operator.py +++ b/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 = [] diff --git a/tensor/matrix_product_state.py b/tensor/matrix_product_state.py index 6a4943b..194099a 100644 --- a/tensor/matrix_product_state.py +++ b/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): diff --git a/tensor/utils.py b/tensor/utils.py new file mode 100644 index 0000000..8c59849 --- /dev/null +++ b/tensor/utils.py @@ -0,0 +1,100 @@ +from syngular.tensor import MatrixProductOperator +from syngular.tensor import MatrixProductState + +from typing import Union + +from opt_einsum import contract + +def mul(op1: Union[MatrixProductOperator, MatrixProductState], op2: Union[MatrixProductOperator, MatrixProductState], mode="standard"): + n = op1.sites_number + m = op2.sites_number + + min_bond = min(min(op1.bond_shape), min(op2.bond_shape)) + max_bond = max(max(op1.bond_shape), max(op2.bond_shape)) + + if n != m: + raise Exception("both operator do not have the same number of sites") + + ''' + Standard mode: + - basic contractions of physicial legs + - compression algorithm afterwards + + > MPO / MPS + > MPS / MPS + > MPO / MPO + ''' + if mode == "standard": + sites = [] + + if isinstance(op1, MatrixProductState) and isinstance(op2, MatrixProductOperator): + if not op1.decomposed or not op2.decomposed: + raise Exception("Operators and States must be decomposed") + + for idx in range(n): + site = contract(op1.sites[idx], [1,2,3], op2.sites[idx], [4,2,5,6], [1,4,5,3,6]) + + site = site.reshape(( + op2.shape[idx][0]*op1.shape[idx][0], + op2.shape[idx][2], + op2.shape[idx][3]*op1.shape[idx][2] + )) + sites.append(site) + return MatrixProductState.from_sites(sites) >> min_bond + elif isinstance(op1, MatrixProductOperator) and isinstance(op2, MatrixProductState): + if not op1.decomposed or not op2.decomposed: + raise Exception("Operators and States must be decomposed") + + for idx in range(n): + site = contract(op2.sites[idx], [1,2,3], op1.sites[idx], [4,2,5,6], [1,4,5,3,6]) + + site = site.reshape(( + op1.shape[idx][0]*op2.shape[idx][0], + op1.shape[idx][2], + op1.shape[idx][3]*op2.shape[idx][2] + )) + sites.append(site) + return MatrixProductState.from_sites(sites) >> min_bond + elif isinstance(op1, MatrixProductState) and isinstance(op2, MatrixProductState): + return op1 | op2 + elif isinstance(op1, MatrixProductOperator) and isinstance(op2, MatrixProductOperator): + for idx in range(n): + site = contract(op2.sites[idx], [1,2,3,4], op1.sites[idx], [5,3,6,7], [1,5,2,6,4,7]) + + site = site.reshape(( + op1.shape[idx][0]*op2.shape[idx][0], + op1.shape[idx][1], + op2.shape[idx][2], + op1.shape[idx][3]*op2.shape[idx][3] + )) + sites.append(site) + return MatrixProductOperator.from_sites(sites) >> min_bond + else: + raise Exception("`syn.mul` should be provided MatrixProductState or MatrixProductOperator objects only") + + + ''' + Variational mode: + - compute left and right blocks + - update random state + - sweep + ''' + if mode == "variational": + pass + + ''' + Standard mode: + - basic contractions of physicial legs + - compression algorithm afterwards + ''' + if mode == "optimized": + pass + + ''' + Standard mode: + - basic contractions of physicial legs + - compression algorithm afterwards + ''' + if mode == "fitup": + pass + diff --git a/test/test_dmrg.py b/test/test_dmrg.py index c050df8..e3f1107 100644 --- a/test/test_dmrg.py +++ b/test/test_dmrg.py @@ -1,42 +1,106 @@ +from syngular.tensor.matrix_product_state import MatrixProductState from syngular.variational import DMRG from syngular.tensor import MatrixProductOperator from syngular.variational import Lanczos import numpy as np -if __name__ == "__main__": - +def dmrg_lanczos_1(): + input_shape = (3,3,3,3) #,3,3) output_shape = (3,3,3,3) #,3,3) - bond_shape = (2,2,2) #,2,2) + bond_shape = (9,9,9) #,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) + # 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()# + mpo = MatrixProductOperator.random(input_shape=input_shape, output_shape=output_shape, bond_shape=bond_shape) # MatrixProductOperator(B, bond_shape=(6,6)).decompose()# - print('C', np.around(mpo.to_tensor(), decimals=2)) + diff = np.abs(B - mpo.to_tensor()) + print('Min', np.min(diff)) + print('Max', np.max(diff)) + # 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("Eigenstate", eigenstate.shape) + print("Eigentensor", eigentensor.shape) - print(mpo.to_tensor().shape) + print("MPO Shape", mpo.to_tensor().shape) eigenvalues, eigenvectors = np.linalg.eig(mpo.to_tensor().reshape((3**4, 3**4))) - eigenvalues2, eigenvectors2 = np.linalg.eig(A) + # eigenvalues2, eigenvectors2 = np.linalg.eig(A) eigenvalue = eigenvectors[0] eigenvalue /= np.linalg.norm(eigenvalue) + print(eigenvalues) - print(eigenvalues2) + # print(eigenvalues2) print(eigenvalue - eigentensor.reshape((81,))) - # print(eigenvalues[0] - eigenstate.to_tensor()) \ No newline at end of file + print(np.abs(eigenvalues[0] - eigenstate.to_tensor())) + +def dmrg_lanczos_random(): + + k = 3 + input_shape = (k,k,k,k) + output_shape = (k,k,k,k) + bond_shape = (k**2,k**2,k**2) + + print((*input_shape,*output_shape,)) + + tensor = np.random.rand(*(*input_shape,*output_shape,)) + tensor = (tensor+tensor.T)/2 + # n = k**4 + # A = np.zeros(n**2).reshape((n,n)) + # np.fill_diagonal(A, list(range(n))) + # tensor = A.reshape((*input_shape,*output_shape,)) + operator = MatrixProductOperator.random(input_shape=input_shape, output_shape=output_shape, bond_shape=bond_shape).decompose()#(tensor, bond_shape=bond_shape).decompose() + + diff = np.abs(tensor - operator.to_tensor()) + print('Min', np.min(diff)) + print('Max', np.max(diff)) + print("Param", operator.parameters_number, "True param", operator.real_parameters_number) + + + eigenstate = DMRG.solve(operator, optimizer=Lanczos) + eigenstate.left_orthonormalization() + eigenstate.normalize() + + eigentensor = eigenstate.to_tensor() + eigentensor_reshaped = eigentensor.reshape(81) + + print("Norm eigentensor", np.linalg.norm(eigentensor), eigenstate.dot()) + + eig = np.linalg.eig(tensor.reshape((3**4, 3**4))) + ground_eigenstate = eig[1][0] + # print(eig[0]) + ground_eigenstate /= np.linalg.norm(ground_eigenstate) + # print(np.linalg.norm(ground_eigenstate)) + + print("Eigenstate DMRG Shape", eigentensor.shape) + print("Eigenstate Ground Shape", ground_eigenstate.shape) + print("Norm", np.linalg.norm(ground_eigenstate), np.linalg.norm(eigentensor_reshaped)) + + print(ground_eigenstate) + print("\n") + print(eigentensor_reshaped) + print("\n") + + for gst in eig[1]: + gst /= np.linalg.norm(gst) + print(np.sum(np.abs(gst - eigentensor_reshaped))) + + print((operator @ eigenstate).dot(), eigenstate.dot()) + + +if __name__ == "__main__": + dmrg_lanczos_random() + diff --git a/test/test_mpo.py b/test/test_mpo.py index 7033cee..50e1bdb 100644 --- a/test/test_mpo.py +++ b/test/test_mpo.py @@ -294,15 +294,27 @@ def test_to_tensor(): print("-------------------------------------------------------") def test_compare(): - y = np.ones(2**4).reshape((4,4)) + y = np.zeros(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() + + # y = np.array([ + # [1,0,0,0], + # [0,2,0,0], + # [0,0,3,0], + # [0,0,0,2], + # ]).reshape((2,2,2,2)) + Y = MatrixProductOperator(y, bond_shape=(4,)).decompose().right_orthonormalization() print(y) print(np.around(Y.to_tensor(), decimals=2)) + print('-----') + for site in Y.sites: + print(site) + print("") + print("---------- Compare with to_tensor() ----------") start = time.time() diff --git a/test/test_mps.py b/test/test_mps.py index d37095c..7b0994e 100644 --- a/test/test_mps.py +++ b/test/test_mps.py @@ -97,10 +97,73 @@ def test_random(): print(Z) print(Z.to_tensor()) +def test_dot(): + x = np.arange(4**3).reshape((4,4,4)) + X = MatrixProductState(x, bond_shape=(4,4)).decompose() + + print("---------- MatrixProductOperator.dot() ----------") + start = time.time() + print('Norm', np.linalg.norm(x)**2) + print('MPS Norm', X.dot()) + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + +def test_normalize(): + x = np.arange(4**3).reshape((4,4,4)) + X = MatrixProductState(x, bond_shape=(4,4)).decompose() + + print("---------- MatrixProductOperator.norm() ----------") + start = time.time() + print('Norm (before normalization)', X.dot()) + X.normalize() + print('Norm (after normalization)', X.dot()) + + print('Norm (after tensorization)', np.linalg.norm(X.to_tensor())) + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + +def test_left_canonical(): + x = np.arange(4**3).reshape((4,4,4)) + X = MatrixProductState(x, bond_shape=(2,2)) + + print("---------- MatrixProductOperator.decompose(mode='left') ----------") + start = time.time() + + X.decompose(mode="left") + + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + print(X.to_tensor()) + +def test_right_canonical(): + x = np.arange(4**3).reshape((4,4,4)) + X = MatrixProductState(x, bond_shape=(2,2)) + + print("---------- MatrixProductOperator.decompose(mode='right') ----------") + start = time.time() + + X.decompose(mode="right") + + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + print(X.to_tensor()) # test_add() # test_compress() # test_random() # test_zeros() # test_augment() -test_apply() \ No newline at end of file +# test_apply() +# test_dot() + +# test_normalize() +test_left_canonical() +test_right_canonical() diff --git a/test/test_syn.py b/test/test_syn.py new file mode 100644 index 0000000..92895eb --- /dev/null +++ b/test/test_syn.py @@ -0,0 +1,30 @@ +from colorlog import WARNING +import syngular as syn + +from syngular.tensor import MatrixProductOperator +from syngular.tensor import MatrixProductState + +import numpy as np + + +def test_mul(): + + x = np.arange(2**2) + X = MatrixProductState(x.reshape((2,2)), bond_shape=(2,)).decompose() + + w = np.arange(2**4).reshape((4,4)) + W = MatrixProductOperator(w.reshape((2,2,2,2)), bond_shape=(2,)).decompose() + + print('------') + print(np.matmul(x, w)) + print(syn.mul(W, X, mode="standard").to_tensor()) + print('------') + print(np.matmul(x, w)) + print(syn.mul(X, W, mode="standard").to_tensor()) + print('------') + print(np.matmul(x, x)) + print(syn.mul(X, X, mode="standard")) + print('------') + + +test_mul() \ No newline at end of file diff --git a/variational/dmrg.py b/variational/dmrg.py index 23d5ed8..e549d86 100644 --- a/variational/dmrg.py +++ b/variational/dmrg.py @@ -18,10 +18,24 @@ def solve(operator: MatrixProductOperator, optimizer: Type[Optimizer]=Lanczos): n = operator.sites_number input_shape = operator.input_shape - bond_shape = operator.bond_shape + bond_shape = operator.bond_shape if input_shape[0] >= operator.bond_shape[0] else input_shape[:-1] + + print(input_shape, bond_shape) + + tensor = np.random.random(size=input_shape) + tensor /= np.linalg.norm(tensor) + + print("Shape", tensor.shape) + print(np.linalg.norm(tensor)) - state = MatrixProductState.random(input_shape=input_shape, bond_shape=bond_shape) + state = MatrixProductState(tensor, bond_shape=bond_shape).decompose() #.random(input_shape=input_shape, bond_shape=bond_shape)# + state.normalize()#MatrixProductState.random(input_shape=input_shape, bond_shape=bond_shape) state.right_orthonormalization() + print("Random state norm", np.linalg.norm(state.to_tensor())) + + # state = MatrixProductState.random(input_shape=input_shape, bond_shape=bond_shape) + # state.right_orthonormalization() + print('[DMRG] Initialization') @@ -153,9 +167,11 @@ def __right_sweep(operator: MatrixProductOperator, state: MatrixProductState, le _, nsite = optimizer.fit(site, init_vec=init_vec) init_vec = nsite[0] + print(nsite) 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]) diff --git a/variational/optimizer.py b/variational/optimizer.py index 36c4778..eac1251 100644 --- a/variational/optimizer.py +++ b/variational/optimizer.py @@ -10,7 +10,7 @@ def fit(tensor): class Lanczos(Optimizer): @staticmethod - def fit(tensor, init_vec = None, iteration=100, ): + def fit(tensor, init_vec = None, iteration=100) -> tuple[np.ndarray]: 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) @@ -38,7 +38,7 @@ def fit(tensor, init_vec = None, iteration=100, ): T[0,0] = alpha - for j in range(1, m-1): + for j in range(1, m): beta = np.sqrt(np.dot(w,w)) vj = w/beta @@ -58,3 +58,10 @@ def fit(tensor, init_vec = None, iteration=100, ): T[j,j-1] = beta return T, V + + +class Davidson: + + @staticmethod + def fit(tensor, init_vec = None, iteration=100) -> tuple[np.ndarray]: + pass \ No newline at end of file