From 0fb4d520d8ba200e6265e347a7c4ed1e5352d48e Mon Sep 17 00:00:00 2001 From: Antoine Debouchage Date: Mon, 20 Dec 2021 15:29:05 +0100 Subject: [PATCH] First qunatum glimpse & optimization of einsum with opt_einsum package --- quantum/__init__.py | 2 ++ quantum/circuit.py | 30 +++++++++++++++++ quantum/gate.py | 27 +++++++++++++++ quantum/qbit.py | 20 +++++++++++ tensor/matrix_product_operator.py | 13 ++++--- tensor/matrix_product_state.py | 56 ++++++++++++++++++------------- test/test_mpo.py | 21 ++++++------ test/test_mps.py | 34 ++++++++++++++++++- test/test_quantum.py | 26 ++++++++++++++ trash/deep.py | 2 +- 10 files changed, 190 insertions(+), 41 deletions(-) create mode 100644 quantum/__init__.py create mode 100644 quantum/circuit.py create mode 100644 quantum/gate.py create mode 100644 quantum/qbit.py create mode 100644 test/test_quantum.py diff --git a/quantum/__init__.py b/quantum/__init__.py new file mode 100644 index 0000000..ca150a7 --- /dev/null +++ b/quantum/__init__.py @@ -0,0 +1,2 @@ +from syngular.quantum.circuit import Circuit +from syngular.quantum.qbit import Qbit \ No newline at end of file diff --git a/quantum/circuit.py b/quantum/circuit.py new file mode 100644 index 0000000..82541ff --- /dev/null +++ b/quantum/circuit.py @@ -0,0 +1,30 @@ + + +from syngular.tensor.matrix_product_state import MatrixProductState + + +class Circuit: + + def __init__(self, size, bond=2, initializer="ground", structure=[]): + self.initializer = initializer + self.size = size + + self.structure = structure + + self.current_step = 0 + self.current_state = None + self.states = [] + + def run(self): + pass + + def reset(self): + if self.initializer == 'ground': + self.current_state = MatrixProductState.zeros( + input_shape=(2,)*self.size, + bond_shape=(2,)*(self.size-1) + ) + self.states.append(self.current_state) + + def step(self): + pass \ No newline at end of file diff --git a/quantum/gate.py b/quantum/gate.py new file mode 100644 index 0000000..22e2b6f --- /dev/null +++ b/quantum/gate.py @@ -0,0 +1,27 @@ +import numpy as np + + +I = np.array([ + [1., 0.], + [0., 1.] +]) + +X = np.array([ + [0., 1.], + [1., 0.] +]) + +Y = np.array([ + [0., -1.j], + [1.j, 0.] +]) + +Z = np.array([ + [1., 0.], + [0., -1.] +]) + +H = 1/np.sqrt(2) * np.array([ + [1., 1.], + [1., -1.] +]) \ No newline at end of file diff --git a/quantum/qbit.py b/quantum/qbit.py new file mode 100644 index 0000000..b818fbe --- /dev/null +++ b/quantum/qbit.py @@ -0,0 +1,20 @@ + + +from syngular.tensor.matrix_product_state import MatrixProductState +from syngular.quantum import gate + +class Qbit: + + def __init__(self, size): + self.size = size + self.dim = 2**self.size + + self.state = MatrixProductState.zeros((2,)*size, (2,)*(size-1)).decompose() + self.state.real_parameters_number = self.dim + + for idx in range(self.state.sites_number-1): + self.state.sites[idx][0] = gate.I + self.state.sites[-1][0][0] = 1. + + def to_tensor(self): + return self.state.to_tensor().reshape(self.dim) \ No newline at end of file diff --git a/tensor/matrix_product_operator.py b/tensor/matrix_product_operator.py index 7fb1056..ae0d130 100644 --- a/tensor/matrix_product_operator.py +++ b/tensor/matrix_product_operator.py @@ -5,6 +5,8 @@ import itertools from typing import Tuple, List, Union +from opt_einsum import contract + import numpy as np from scipy import linalg @@ -160,7 +162,7 @@ def __matmul__(self, mp: Union[MatrixProductOperator, MatrixProductState]): if isinstance(mp, MatrixProductState): for idx in range(self.sites_number): - site = np.einsum(mp.sites[idx], [1,2,3], self.sites[idx], [4,2,6,5], [1,4,6,3,5]) + 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], @@ -171,7 +173,7 @@ def __matmul__(self, mp: Union[MatrixProductOperator, MatrixProductState]): return MatrixProductState.from_sites(sites) elif isinstance(mp, MatrixProductOperator): for idx in range(self.sites_number): - site = np.einsum(self.sites[idx], [1,2,3,4], mp.sites[idx], [5,3,6,7], [1,5,2,6,4,7]) + 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], @@ -288,6 +290,7 @@ def to_tensor(self): for inp in itertools.product(*range_inp): for out in itertools.product(*range_out): + print(inp, out) tensor[(*inp, *out)] = self[inp, out] return tensor @@ -461,7 +464,7 @@ def apply(self, operator, indices): struct += output_shape # print(struct) - T = np.squeeze(np.einsum(*struct), axis=(1, len(output_shape)-2)) + T = np.squeeze(contract(*struct), axis=(1, len(output_shape)-2)) for idx in range(m-1): jdx = indices[idx] @@ -505,7 +508,7 @@ def apply(self, operator, indices): # ] # struct += output_shape # print(struct) - T = np.einsum(*struct) + T = contract(*struct) print(T.shape) else: @@ -522,7 +525,7 @@ def retrieve(self, input_indices, output_indices): einsum_structure.append(self.sites[idx][:, input_indices[idx], output_indices[idx], :]) einsum_structure.append([idx, idx+1]) - return np.einsum(*einsum_structure) + return contract(*einsum_structure) def left_orthonormalization(self, bond_shape=()): diff --git a/tensor/matrix_product_state.py b/tensor/matrix_product_state.py index 3277f66..263b1be 100644 --- a/tensor/matrix_product_state.py +++ b/tensor/matrix_product_state.py @@ -5,7 +5,7 @@ import numpy as np from scipy import linalg - +from opt_einsum import contract class MatrixProductState: @@ -115,7 +115,7 @@ def __or__(self, mp): self.sites[idx], [idx,n+idx+1,idx+1], mp.sites[idx], [2*n+idx+1, n+idx+1, 2*n+idx+2] ] - r = np.einsum(*struct) + r = contract(*struct) return r.reshape(1)[0] else: raise Exception("right-hand site must be a MatrixProductState") @@ -186,7 +186,6 @@ def from_sites(sites: List[np.ndarray], orthogonality=None, real_parameters_numb mp.sites = sites mp.sites_number = len(sites) - mp.decomposed = True mp.decomposed = True mp.orthonormalized = orthogonality @@ -210,6 +209,19 @@ def from_sites(sites: List[np.ndarray], orthogonality=None, real_parameters_numb def empty(): return MatrixProductState() + @staticmethod + def zeros(input_shape, bond_shape): + n = len(input_shape) + + shape = [(1, input_shape[0],bond_shape[0])] + shape += [(bond_shape[i-1], input_shape[i], bond_shape[i]) for i in range(1, n-1)] + shape += [(bond_shape[n-2], input_shape[n-1], 1)] + + print(shape) + sites = [] + for idx in range(n): + sites.append(np.zeros(shape=shape[idx])) + return MatrixProductState.from_sites(sites) """norm method of MatrixProductState Compute the norm of the MatrixProductState @@ -244,32 +256,29 @@ def copy(self): def decompose(self, mode="left"): if not self.decomposed: if mode == "left": - current_matrix = self.tensor - current_shape = self.shape[0][1] - current_rank = 1 + n = self.sites_number-1 + T = self.tensor + + for k in range(self.sites_number-1): - unfolding_matrix = np.reshape(current_matrix, newshape=(current_shape*current_rank, -1)) - rank = self.shape[k][2] + L = self.left_matricization(T, index=k) - Q, R = np.linalg.qr(unfolding_matrix, mode="complete") - + Q, R = np.linalg.qr(L, mode="complete") + print(Q.shape, R.shape) + rank = self.shape[k][2] Q = Q[:,:rank] - Q = np.reshape(Q, newshape=(current_rank, current_shape, -1)) R = R[:rank, :] - if self.verbose: print(f"Core {k} with {current_rank} , {current_shape}") + self.sites[k] = self.tensoricization(Q, k) + T = R - self.sites[k] = Q + self.parameters_number += np.prod(self.sites[k].shape) - current_shape = self.shape[k+1][1] - current_rank = rank - current_matrix = R + self.sites[-1] = self.tensoricization(T, n) + self.parameters_number += np.prod(self.sites[-1].shape) - current_matrix = current_matrix[:, :, np.newaxis] - self.sites[-1] = current_matrix - del self.tensor gc.collect() @@ -352,7 +361,7 @@ def compress(self, dim: int, mode="left", strict=False): else: that = self.copy() - print(that) + # print(that) if mode == 'left': # if self.orthonormalized != 'left': self.left_orthonormalization() @@ -367,8 +376,8 @@ def compress(self, dim: int, mode="left", strict=False): W = S @ R - print(Q.shape, L.shape, S.T.shape, R.shape) - print(W.shape) + # print(Q.shape, L.shape, S.T.shape, R.shape) + # print(W.shape) l1 = list(that.shape[k]) l2 = list(that.shape[k+1]) @@ -392,8 +401,7 @@ def retrieve(self, indices): for idx in range(self.sites_number): einsum_structure.append(self.sites[idx][:, indices[idx], :]) einsum_structure.append([idx, Ellipsis, idx+1]) - - return np.einsum(*einsum_structure) + return contract(*einsum_structure) def left_orthonormalization(self, bond_shape=()): diff --git a/test/test_mpo.py b/test/test_mpo.py index e786fe1..ac35bf7 100644 --- a/test/test_mpo.py +++ b/test/test_mpo.py @@ -246,7 +246,7 @@ def test_random(): print("---------- MatrixProductOperator.random(input_shape, output_shape, bond_shape) ----------") start = time.time() - MatrixProductOperator.random((4,4,4), (2,2,2), (3,3,)) + MatrixProductOperator.random((4,4,4), (2,2,2), (8,8,)) end = time.time() print(f"> Execution time : {end-start:.8f}sec") print("-------------------------------------------------------") @@ -284,6 +284,7 @@ def test_to_tensor(): tensor = np.arange((2*2*2*2*2*2*2*2)*(2*2*2*2*2*2*2*2)).reshape(2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2) mpo = MatrixProductOperator(tensor, bond_shape=(2,2,2,2,2,2,2,)).decompose() print(mpo) + print(mpo.real_parameters_number, mpo.parameters_number) print("---------- MatrixProductOperator.to_tensor() ----------") start = time.time() mpo.to_tensor() @@ -293,14 +294,14 @@ def test_to_tensor(): -# test_decompose() -# test_from_sites() -# test_random() -# test_to_tensor() -# test_add() +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() \ No newline at end of file +test_dot() +test_compress() +test_orthogonality() +test_apply() \ No newline at end of file diff --git a/test/test_mps.py b/test/test_mps.py index 086577a..7471654 100644 --- a/test/test_mps.py +++ b/test/test_mps.py @@ -7,6 +7,35 @@ np.set_printoptions(suppress=True) +def test_augment(): + x = np.arange(8).reshape((2,2,2)) + X = MatrixProductState(x, bond_shape=(2,2)).decompose() + + Y = MatrixProductState.zeros((2,2,2),(2,2,)) + + print("---------- MatrixProductOperator.__add__(a, b) ----------") + start = time.time() + Z = X + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + Y + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + + print(Z) + print(Z.to_tensor()) + print(X.to_tensor()) + print((Z >> 2).to_tensor()) + + +def test_zeros(): + print("---------- MatrixProductOperator.__add__(a, b) ----------") + start = time.time() + Z = MatrixProductState.zeros((2,2,2),(2,2,)) + end = time.time() + print(f"> Execution time : {end-start:.8f}sec") + print("-------------------------------------------------------") + print(Z.to_tensor()) + print(Z.shape) + def test_add(): x = np.arange(8).reshape((2,2,2)) y = np.arange(8).reshape((2,2,2)) @@ -58,4 +87,7 @@ def test_random(): # test_add() # test_compress() -test_random() \ No newline at end of file +# test_random() +# test_zeros() + +test_augment() \ No newline at end of file diff --git a/test/test_quantum.py b/test/test_quantum.py new file mode 100644 index 0000000..2807d80 --- /dev/null +++ b/test/test_quantum.py @@ -0,0 +1,26 @@ +from syngular.quantum import Circuit, Qbit +import syngular.quantum.gate as gate + +circ = Circuit(initializer="ground", size=2, structure=[ + (gate.X, 0) +]) + +size = 30 + +ground = Qbit(size) + +# print(ground.to_tensor()) +print(ground.state[(1,)*size]) +print(ground.state[(0,)*size]) +# print(ground.state.real_parameters_number, ground.state.parameters_number) +# import numpy as np +# from syngular.tensor import MatrixProductState +# ground = np.zeros(shape=(2**size)) +# ground[0] = 1 +# print(ground) +# qbit2 = MatrixProductState(ground.reshape((2,)*size), (3,)*(size-1)).decompose() +# print(qbit2) +# for site in qbit2.sites: +# print('---') +# print(site) +# print(qbit2.to_tensor().reshape(2**size)) \ No newline at end of file diff --git a/trash/deep.py b/trash/deep.py index e045b94..4fe28d5 100644 --- a/trash/deep.py +++ b/trash/deep.py @@ -1,7 +1,7 @@ import itertools from PIL.Image import new from numpy.core.fromnumeric import size -from syngular.tensor.tensor_train import MatrixProductOperator, MatrixProductState +from syngular.trash.tensor_train import MatrixProductOperator, MatrixProductState import numpy as np