Skip to content
Snippets Groups Projects
Commit e7a9a080 authored by Giovanni Cataldi's avatar Giovanni Cataldi
Browse files

Upload New File

parent b9387d59
Branches main
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
07_matrix_product_state/MPS.ipynb 0 → 100644
+ 750
− 0
View file @ e7a9a080
%% Cell type:markdown id: tags:
# Matrix product states
Iterativelt applying the Schmidt decomposition we learned in the last lesson we are able to decompose a quantum state in a Matrix Product State (MPS). It is also known as tensor train.
The state of each qubit is represented by a tensor, and the links between the qubits represent the quantum correlation (entanglement) between the qubits. here is the general formula:
$$
\ket{\psi} = \sum_{{s_1,\dots,s_n=0}}^1\, \sum_{{\alpha_1, \dots, \alpha_n = 1}}^\chi \textbf{M}_{1 \alpha_1}^{[1],s_1}\, \textbf{M}_{\alpha_1 \alpha_2}^{[2],s_2} \dots\textbf{M}_{\alpha_{n-2} \alpha_{n-1}}^{[n-1],s_{n-1}}\, \textbf{M}_{\alpha_{n-1} 1}^{[n], s_n}\, \ket{s_1 s_2 \dots s_n}\,
$$
In this hands-on, we will understand how to build and measure local observables on an MPS.
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
yes = "\033[1;32m" + u"\N{check mark}" #+ "\033[1;30m"
no = "\033[1;31m" + u"\N{ballot x}" #+ "\033[1;30m"
```
%% Cell type:code id: tags:
``` python
def tSVD(
tensor,
legs_left,
legs_right,
contract_singvals="N",
max_bond_dimension=128
):
"""Perform a truncated Singular Value Decomposition by
first reshaping the tensor into a legs_left x legs_right
matrix, and permuting the legs of the ouput tensors if needed.
If the contract_singvals = ('L', 'R') it takes care of
renormalizing the output tensors such that the norm of
the MPS remains 1 even after a truncation.
Parameters
----------
tensor : ndarray
Tensor upon which apply the SVD
legs_left : list of int
Legs that will compose the rows of the matrix
legs_right : list of int
Legs that will compose the columns of the matrix
contract_singvals: string, optional
How to contract the singular values.
'N' : no contraction
'L' : to the left tensor
'R' : to the right tensor
Returns
-------
tens_left: ndarray
left tensor after the SVD
tens_right: ndarray
right tensor after the SVD
singvals: ndarray
singular values kept after the SVD
singvals_cutted: ndarray
singular values cutted after the SVD, normalized with the biggest singval
"""
# 1) Reshaping
matrix = np.transpose(tensor, legs_left + legs_right)
shape_left = np.array(tensor.shape)[legs_left]
shape_right = np.array(tensor.shape)[legs_right]
matrix = matrix.reshape(np.prod(shape_left), np.prod(shape_right))
# 2) SVD decomposition
mat_left, singvals_tot, mat_right = np.linalg.svd(
matrix, full_matrices=False
)
# 3) Truncation
lambda1 = singvals_tot[0]
cut = np.nonzero(singvals_tot / lambda1 < 1e-8)[0]
# Confront the cut w.r.t the maximum bond dimension
if len(cut) > 0:
cut = min(max_bond_dimension, cut[0])
else:
cut = max_bond_dimension
cut = min(cut, len(singvals_tot))
singvals_cutted = singvals_tot[cut:]
singvals = singvals_tot[:cut]
mat_left = mat_left[:, :cut]
mat_right = mat_right[:cut, :]
# Contract singular values if requested
if contract_singvals.upper() == "L":
mat_left = np.multiply(mat_left, singvals)
elif contract_singvals.upper() == "R":
mat_right = np.multiply(singvals, mat_right.T).T
elif contract_singvals.upper() != "N":
raise ValueError(
f"Contract_singvals option {contract_singvals} is not "
+ "implemented. Choose between right (R), left (L) or None (N)."
)
# Reshape back to tensors
tens_left = mat_left.reshape(list(shape_left) + [cut])
tens_right = mat_right.reshape([cut] + list(shape_right))
return tens_left, tens_right, singvals, singvals_cutted
def check_isometry(tensor, legs):
id = np.tensordot(tensor, tensor.conj(), (legs, legs))
return np.isclose(np.identity(id.shape[0]), id).all()
```
%% Cell type:markdown id: tags:
## Compressing a general statevector in an MPS
By iteratively using the SVD decomposition decompose the statevector in MPS form.
After this cell there is another to test your code!
<details>
<summary>Solution part 1</summary>
```python
state_tensor = statevector.reshape( [1] + [local_dim] * num_sites + [1] )
```
</details>
<details>
<summary>Solution part 2</summary>
```python
state_tensor, legs[:2], legs[2:], contract_singvals="R",
```
</details>
%% Cell type:code id: tags:
``` python
def from_statevector(
statevector,
local_dim=2,
max_bond_dim=128
):
"""
Initialize the MPS tensors by decomposing a statevector into MPS form.
All the degrees of freedom must have the same local dimension
Parameters
----------
statevector : ndarray of shape( local_dim^num_sites, )
Statevector describing the interested state for initializing the MPS
local_dim : int, optional
Local dimension of the degrees of freedom. Default to 2.
Returns
-------
obj : :py:class:`MPS`
MPS simulator class
Examples
--------
>>> -U1 - U2 - U3 - ... - UN-
>>> | | | |
# For d=2, N=7 and chi=5, the tensor network is as follows:
>>> -U1 -2- U2 -4- U3 -5- U4 -5- U5 -4- U6 -2- U7-
>>> | | | | | | |
# where -x- denotes the bounds' dimension (all the "bottom-facing" indices
# are of dimension d=2). Thus, the shapes
# of the returned tensors are as follows:
>>> U1 U2 U3 U4 U5 U6 U7
>>> [(1, 2, 2), (2, 2, 4), (4, 2, 5), (5, 2, 5), (5, 2, 4), (4, 2, 2), (2, 2, 1)]
"""
if not isinstance(statevector, np.ndarray):
raise TypeError("Statevector must be numpy array")
num_sites = int(np.log(len(statevector)) / np.log(local_dim))
state_tensor = statevector.reshape( ... )
mps = []
for ii in range(num_sites - 1):
legs = list(range(len(state_tensor.shape)))
tens_left, tens_right, singvals, _ = tSVD(
..., ..., ..., ...,
max_bond_dimension=max_bond_dim
)
mps.append( tens_left )
state_tensor = tens_right
mps.append( tens_right )
return mps
```
%% Cell type:code id: tags:
``` python
num_sites = 10
state = np.zeros(2**num_sites)
state[[0, -1]] = 1/np.sqrt(2)
mps = from_statevector(state)
for idx, tt in enumerate(mps[:-1]):
is_iso = check_isometry(tt, [0, 1])
if not is_iso:
print(no + f"Tensor {idx} of the mps is not an isometry" )
break
print(yes + "All tensors are isometries!")
norm = np.tensordot( mps[-1], mps[-1].conj(), ([0, 1, 2], [0, 1, 2]) )
if np.isclose(norm, 1):
print(yes + "Norm of the quantum state is 1!")
else:
print(no + f"Norm of the quantum state is {norm}")
```
%% Cell type:markdown id: tags:
## Moving the isometry in the MPS
The isometry is a very special tensor in the MPS. It can be used to greatly speed up convergence.
It will be better explained at the blackboard!
Here you have to move your isometry from the right to the left!
Don't worry about left-to-right, that will be implemented later on.
Again, check for the cheks afterwards!
<details>
<summary>Solution part 3</summary>
```python
for ii in range(from_idx, to_idx, -1)
```
</details>
<details>
<summary>Solution part 4</summary>
```python
right, left, singvals, _ = tSVD(
tensor, [0], [1, 2], contract_singvals="L"
)
# Update the tensors in the MPS
mps[ii] = left
mps[ii - 1] = np.tensordot(mps[ii - 1], right, ([2], [0]))
```
</details>
%% Cell type:code id: tags:
``` python
def move_iso_to_left(mps, from_idx, to_idx):
"""
Apply a gauge transformation to all bonds between
:py:method:`MPS.num_sites` and `idx`, so that all
sites between the last (rightmost one) and idx
are set to (semi)-unitary tensors.
Parameters
----------
from_idx: int
index of the isometry
to_idx: int
index of the tensor up to which the canonization occurs
"""
for ii in range(from_idx, to_idx, -1):
tensor = mps[ii]
right, left, singvals, _ = tSVD(
..., ..., ..., contract_singvals=...
)
# Update the tensors in the MPS
mps[ii] = left
mps[ii - 1] = np.tensordot(mps[ii - 1], right, (..., ...))
return mps
```
%% Cell type:code id: tags:
``` python
new_idx = 2
old_idx = num_sites-1
mps = move_iso_to_left(mps, old_idx, new_idx)
norm = np.tensordot( mps[new_idx], mps[new_idx].conj(), ([0, 1, 2], [0, 1, 2]) )
if np.isclose(norm, 1):
print(yes + "Norm of the quantum state is 1!")
else:
print(no + f"Norm of the quantum state is {norm}")
```
%% Cell type:markdown id: tags:
## EasyMPS class!
Let's now treat the MPS as a class. It will be much easier to handle stuff!!
In particular, we will store the isometry center, and thus we do not need to remember it.
Furthermore, using the knowledge of HandsOn `03` we will treat the class as a list.
You can access the tensor `idx` of the mps as `mps[idx]`!!
The class is just a container for all the functions we already implemented. The only one
we will not explain from the computational point of view is `meas_tensor_product`.
### Task 1: initialization
Initialize the MPS in the $|00\dots0\rangle$ state. You have to work on the method `__init__`. The cell below will help you in testing it!
<details>
<summary>Solution part 5</summary>
```python
ket0_tensor = np.zeros((1, local_dim, 1))
ket0_tensor[0, 0, 0] = 1
```
</details>
### Task 2: local expectation values
We call local expectation value the expectation value of a local observable. You have to work on the
method `meas_local`.
<details>
<summary>Solution part 6</summary>
```python
first_contraction = np.tensordot(self[idx], operator, ([1], [1]))
first_contraction = first_contraction.transpose(0, 2, 1)
expectation = np.tensordot(first_contraction, self[idx].conj(), ([0, 1, 2], [0, 1, 2]))
```
</details>
%% Cell type:code id: tags:
``` python
class EasyMPS():
def __init__(self, num_sites, local_dim=2, max_bond_dim=128):
self.num_sites = num_sites
self.local_dim = local_dim
self.bd = max_bond_dim
ket0_tensor = ...
...
self.tensors = [ket0_tensor.copy() for _ in range(num_sites)]
self.iso_center = 0
def __len__(self):
"""
Provide number of sites in the MPS
"""
return self.num_sites
def __getitem__(self, key):
"""Overwrite the call for lists, you can access tensors in the MPS using
.. code-block::
MPS[0]
>>> [[ [1], [0] ] ]
Parameters
----------
key : int
index of the MPS tensor you are interested in
Returns
-------
np.ndarray
Tensor at position key in the MPS.tensor array
"""
return self.tensors[key]
def __setitem__(self, key, value):
"""Modify a tensor in the MPS by using a syntax corresponding to lists.
It is the only way to modify a tensor
.. code-block::
tens = np.ones( (1, 2, 1) )
MPS[1] = tens
Parameters
----------
key : int
index of the array
value : np.array
value of the new tensor. Must have the same shape as the old one
"""
self.tensors[key] = value
return None
def __iter__(self):
"""Iterator protocol"""
return iter(self.tensors)
@classmethod
def from_statevector(
cls,
statevector,
local_dim=2,
max_bond_dim=128
):
"""
Initialize the MPS tensors by decomposing a statevector into MPS form.
All the degrees of freedom must have the same local dimension
Parameters
----------
statevector : ndarray of shape( local_dim^num_sites, )
Statevector describing the interested state for initializing the MPS
local_dim : int, optional
Local dimension of the degrees of freedom. Default to 2.
Returns
-------
obj : :py:class:`MPS`
MPS simulator class
Examples
--------
>>> -U1 - U2 - U3 - ... - UN-
>>> | | | |
# For d=2, N=7 and chi=5, the tensor network is as follows:
>>> -U1 -2- U2 -4- U3 -5- U4 -5- U5 -4- U6 -2- U7-
>>> | | | | | | |
# where -x- denotes the bounds' dimension (all the "bottom-facing" indices
# are of dimension d=2). Thus, the shapes
# of the returned tensors are as follows:
>>> U1 U2 U3 U4 U5 U6 U7
>>> [(1, 2, 2), (2, 2, 4), (4, 2, 5), (5, 2, 5), (5, 2, 4), (4, 2, 2), (2, 2, 1)]
"""
if not isinstance(statevector, np.ndarray):
raise TypeError("Statevector must be numpy array")
num_sites = int(np.log(len(statevector)) / np.log(local_dim))
mps = cls(num_sites, local_dim, max_bond_dim)
state_tensor = statevector.reshape([1] + [local_dim] * num_sites + [1])
for ii in range(num_sites - 1):
legs = list(range(len(state_tensor.shape)))
tens_left, tens_right, singvals, _ = tSVD(
state_tensor, legs[:2], legs[2:], contract_singvals="R",
max_bond_dimension=max_bond_dim
)
mps[ii] = tens_left
state_tensor = tens_right
mps[-1] = tens_right
mps.iso_center = num_sites - 1
return mps
def move_iso_to_left(self, to_idx):
"""
Apply a gauge transformation to all bonds between
:py:method:`MPS.num_sites` and `idx`, so that all
sites between the last (rightmost one) and idx
are set to (semi)-unitary tensors.
Parameters
----------
to_idx: int
index of the tensor up to which the canonization occurs
"""
for ii in range(self.iso_center, to_idx, -1):
tensor = mps[ii]
right, left, singvals, _ = tSVD(
tensor, [0], [1, 2], contract_singvals="L"
)
# Update the tensors in the MPS
mps[ii] = left
mps[ii - 1] = np.tensordot(mps[ii - 1], right, ([2], [0]))
self.iso_center = to_idx
return mps
def move_iso_to_right(self, to_idx):
"""
Apply a gauge transformation to all bonds between 0 and `idx`,
so that all sites between the first (òeftmpst one) and idx
are set to (semi)-unitary tensors.
Parameters
----------
to_idx: int
index of the tensor up to which the canonization occurs
"""
for ii in range(self.iso_center, to_idx):
tensor = mps[ii]
left, right, singvals, _ = tSVD(
tensor, [0, 1], [2], contract_singvals="R"
)
# Update the tensors in the MPS
mps[ii] = left
mps[ii + 1] = np.tensordot(mps[ii + 1], right, ([0], [1])).transpose(2, 0, 1)
self.iso_center = to_idx
return mps
def iso_towards(self, new_iso):
"""
Move the isometry center to a new site
Parameters
----------
new_iso : int
New isometry center
"""
self.move_iso_to_left(to_idx=new_iso)
self.move_iso_to_right(to_idx=new_iso)
def meas_local(self, operator, idx):
"""
Measure the local expectation value of the operator `operator`
on the index `idx`
Parameters
----------
operator : np.ndarray
Matrix of the operator
idx : int
Index of the site
Returns
-------
float
Expectation value
"""
self.iso_towards(idx)
first_contraction = np.tensordot(self[idx], operator, ([...], [...]))
first_contraction = first_contraction.transpose(0, 2, 1)
expectation = np.tensordot(..., ..., ([...], [...]))
return np.squeeze(expectation)
def meas_tensor_product(self, ops, idxs):
"""
Measure the tensor products of n operators `ops` acting on the indexes `idxs`
Parameters
----------
ops : list of ndarrays
List of numpy arrays which are one-site operators
idxs : list of int
Indexes where the operators are applied
Returns
-------
measure : float
Result of the measurement
"""
if len(idxs) == 0:
return 1
order = np.argsort(idxs)
idxs = np.array(idxs)[order]
ops = np.array(ops)
ops = ops[order]
self.iso_towards(idxs[0])
transfer_mat = np.eye(self[idxs[0]].shape[0], dtype=np.complex64)
jj = 0
closed = False
for ii in range(idxs[0], self.num_sites):
if closed:
break
# Case of finished tensors
if jj == len(idxs):
# close with transfer matrix of correct size
closing_transfer_mat = np.eye(self[ii].shape[0])
measure = np.tensordot(
transfer_mat, closing_transfer_mat, ([0, 1], [0, 1])
)
closed = True
# Case of operator inside
elif idxs[jj] == ii:
transfer_mat = np.tensordot(transfer_mat, self[ii], ([0], [0]))
transfer_mat = np.tensordot(transfer_mat, ops[jj], ([1], [1]))
transfer_mat = np.transpose(transfer_mat, [0, 2, 1])
transfer_mat = np.tensordot(
transfer_mat, np.conj(self[ii]), ([0, 1], [0, 1])
)
jj += 1
# Case of no operator between the sites
else:
transfer_mat = np.tensordot(transfer_mat, self[ii], ([0], [0]))
transfer_mat = np.tensordot(
transfer_mat, np.conj(self[ii]), ([0, 1], [0, 1])
)
if not closed:
# close with transfer matrix of correct size
closing_transfer_mat = np.eye(self[-1].shape[2])
measure = np.tensordot(transfer_mat, closing_transfer_mat, ([0, 1], [0, 1]))
closed = True
return np.real(measure)
```
%% Cell type:code id: tags:
``` python
zz = np.array([[1, 0], [0, -1]])
xx = np.array([[0, 1], [1, 0]])
```
%% Cell type:code id: tags:
``` python
num_sites = 10
mps = EasyMPS(num_sites)
passed = True
for ii in range(num_sites):
val_z = mps.meas_local(zz, ii)
if not np.isclose(val_z, 1):
print(no + f"Expectation value of Z should be 1, not {val_z}" )
passed = False
break
if passed:
print(yes + f"Expectation value of Z computed correctly" )
passed = True
for ii in range(num_sites):
val_x = mps.meas_local(xx, ii)
if not np.isclose(val_x, 0.0):
print(no + f"Expectation value of X should be 0, not {val_x}" )
passed = False
break
if passed:
print(yes + f"Expectation value of X computed correctly" )
```
%% Cell type:code id: tags:
``` python
num_sites = 7
state = np.zeros(2**num_sites)
state[[0, -1]] = 1/np.sqrt(2)
mps = EasyMPS.from_statevector(state)
print( f"Local expectation value of Z: {mps.meas_local(zz, 0)}" )
print( f"Local expectation value of X: {mps.meas_local(xx, 0)}" )
print( f"Expectation value of the parity: {mps.meas_tensor_product([zz for _ in range(num_sites)], np.arange(num_sites))}")
print( f"Expectation value of the XX...XX: {mps.meas_tensor_product([xx for _ in range(num_sites)], np.arange(num_sites))}")
```
%% Cell type:code id: tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment