<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Copyright 2021 Google LLC
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
https
Step1:  Part 1
Review of SVD and EVD
Eigenvalue decomposition
Any Hermitian matrix $M$ can be decomposed as
$$
M = U D U^\dagger
$$
where $U$ is a unitary matrix and $D$ is a diagonal matrix. The diagonal entries of $D$ are the eigenvalues of $M$.
Step2:  Singular value decomposition
Any $n \times m$ matrix $M$ with complex entries can be decomposed as
$$
M = U S V^\dagger
$$
where $U$ is an $n \times n$ unitary matrix, $V$ is an $m \times m$ unitary matrix, and $S$ is an $n \times m$ matrix with non-negative entries on the diagonal and zeros everywhere else.
The diagonal entries in $S$ are called singular values. The singular values of $M$ and coincide with the positive square roots of the eigenvalues of $M^\dagger M$ and $M M^\dagger$, which we can deduce from the following equations
Step5:  MPS
tensorSVD is a helper function used for making an MPS. The idea of it is that it does the reshaping of a tensor into a matrix for SVD, and then reshapes the U and V back into tensors for you.
Step6:  We can compress the parity indicator into a MPS form by a sequence of SVDs, as explained in the lecture. The following code block
Step7:  In this scheme, the four MPS tensors consist of each of the $U$ tensors from the SVDs for all but the final tensor. We can see their shapes here
Step8:  We can interpret the function of these MPS tensors by looking at which entries are nonzero. In particular, for $U2$ we have the following
Step9:  Part 2
In this section we will walk through the creation of simple Tree Tensor Neworks (TTNs) using the density matrix method. These examples are phrased as multipart exercises, so you might try to code up the solution yourself before looking ahead.
Problem 1
(i)
We begin by defining the tensor 
$$
A_{ijkm} = \sqrt{i+2j+3k+4l + 5m}
$$
where $i,j,k,l, m$ are each two-dimensional indices. That is, each one can take on the value $0$ or $1$. Initialize this tensor as a numpy array.
Step10:  Solution (i)
Step11:  (ii)
Using the appropriate density matrices, form isometries $W_L$ and $W_R$ that transform $A$ into the tree tensor network shown (truncating internal indices to dimension $\chi=2$).
Store $W_L$, $W_R$, and $B$ as arrays. These form the truncated TTN representaiton of $A$.
Step12:  Solution (ii)
Step13:  (iii)
Check the accuracy of the truncated TTN by computing the truncation error $\epsilon = \|A-A_\text{recover}\|/\|A\|$, with the recovered tensor $A_\text{recover}$ given by contracting the tree.
Here $\| \cdot \|$ is the Frobenius norm that is implemented in numpy as numpy.linalg.norm.
Step14:  Solution (iii)
Step15:  Problem 2
In this problem we'll train a TTN to approximate a batch of images, where again the approximation comes from truncating to bond dimension $\chi =2$. The image set we're going to consider consists of five images with five pixels each, but the method is general and you should try repeating the exercise for a different set of images.
(i)
Using the density matrix approach approximate the given set of images using the TTN shown below (truncating internal indices to bond dimension $\chi=2$).
Step19:  Solution (i)
Step20:  (ii)
Check the accuracy of your result by contracting the TTN into a single dense tensor and comparing with the dense tensor representation of the original data.
Step21:  Solution (ii)
We can check the accuracy using the Frobenius norm, just like in the previous problem.
In this case, the example was chosen to allow perfect compression with $\chi=2$, so the total error just comes from machine precision.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from numpy import linalg as LA
import itertools
Explanation: Copyright 2021 Google LLC
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
https://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
Instructions and Imports
Got to File -> Save a copy in Drive to get your own copy of this notebook to play with.
You may need to go to Edit -> Notebook settings and switch the Runtime type to Python 3 to get all of the code cells to run properly
End of explanation
A = np.random.rand(3,3)
M = A @ A.T
D, U = LA.eigh(M) # D is returned as a vector
print(M)
print()
M2 = U @ np.diag(D) @ np.conj(U.T) # Conjugation not always needed
print(M2)
Explanation: Part 1
Review of SVD and EVD
Eigenvalue decomposition
Any Hermitian matrix $M$ can be decomposed as
$$
M = U D U^\dagger
$$
where $U$ is a unitary matrix and $D$ is a diagonal matrix. The diagonal entries of $D$ are the eigenvalues of $M$.
End of explanation
M = np.random.rand(2, 3) + 1j * np.random.rand(2, 3)
U, singular_values, Vh = LA.svd(M)
S = np.zeros((2, 3))
np.fill_diagonal(S, singular_values)
print(M)
print()
M2 = U @ S @ Vh
print(M2)
print(f'singular values: {singular_values}')
print(f'squares of singlar values: {singular_values**2}')
print()
Mh = np.conjugate(np.transpose(M))
print('M Mh')
D, U = LA.eigh(M @ Mh)
print(D)
print()
print('Mh M')
D, U = LA.eigh(Mh @ M)
print(D)
Explanation: Singular value decomposition
Any $n \times m$ matrix $M$ with complex entries can be decomposed as
$$
M = U S V^\dagger
$$
where $U$ is an $n \times n$ unitary matrix, $V$ is an $m \times m$ unitary matrix, and $S$ is an $n \times m$ matrix with non-negative entries on the diagonal and zeros everywhere else.
The diagonal entries in $S$ are called singular values. The singular values of $M$ and coincide with the positive square roots of the eigenvalues of $M^\dagger M$ and $M M^\dagger$, which we can deduce from the following equations:
$$
M^\dagger M = V S^2 V^\dagger
$$
$$
MM^\dagger = US^2 U^\dagger
$$
End of explanation
def tensorSVD(T,uinds,svd_threshold=1E-16):
  Computes the  SVD of an N-index tensor
  
  Args:
    T: Tensor to decompose
    uinds: List of indices forming the "left" effective index. These indices
      belong to the U tensor, and the rest belong to Vh.
    svd_threshold: Singular values smaller than this are truncated.
  
  Returns:
    tenU: Left tensor of the SVD.
    tenS: Diagonal matrix of singular values.
    tenV: Right tensor of the SVD.
  
  NT = len(T.shape)
  Nu = len(uinds)
  dest = range(Nu) # array 0,1,2,...
  pT = np.moveaxis(T,uinds,dest)
  udims = [pT.shape[n] for n in range(Nu)]
  vdims = [pT.shape[n] for n in range(Nu,NT)]
  uD = np.prod(udims)
  vD = np.prod(vdims)
  rpT = np.reshape(pT,[uD,vD])
  U,S,V = LA.svd(rpT,full_matrices=False)
  # Determine truncation size:
  n_svd = len(S)
  for ix in range(n_svd):
    if S[ix] < svd_threshold:
      n_svd = ix
      break
  # Perform the truncation:
  truncU = U[:,0:n_svd]
  truncV = V[0:n_svd,:]
  truncS = S[0:n_svd]
  
  # Restore tensor structure to truncated U, S, V:
  udims.append(n_svd)
  vdims.insert(0,n_svd)
  tenU = np.copy(np.reshape(truncU,udims))
  tenV = np.copy(np.reshape(truncV,vdims))
  tenS = np.diag(truncS)
  return tenU,tenS,tenV
def makeParityIndicator(num_bits):
  Creates an indicator tensor for the parity dataset. The indicator is a
  tensor with 2**num_bits entries and shape (2,)*num_bits.
  Example: 
    psi4 = makeParityIndicator(4)
    psi[0,0,0,0] = 1
    psi[0,0,0,1] = 0
  Args:
    num_bits: Length of the input strings for the indicator.
  Returns:
    psi: Indicator tensor.
  
  psi = np.zeros((2,)*num_bits)
  for bits in itertools.product([0, 1], repeat=num_bits):
    if sum(bits) % 2 == 0:
      psi[bits] = 1
  #psi += 1E-14*np.random.randn(*psi.shape)
  return psi
Explanation: MPS
tensorSVD is a helper function used for making an MPS. The idea of it is that it does the reshaping of a tensor into a matrix for SVD, and then reshapes the U and V back into tensors for you.
End of explanation
# Create indicator for four-bit parity.
psi4 = makeParityIndicator(4)
# Create the MPS tensors via a sequence of SVDs.
U1,S1,V1 = tensorSVD(psi4,[0])
print("V1.shape = ", V1.shape)
M2 = np.tensordot(S1, V1,[1,0])
print("M2.shape = ", M2.shape)
U2,S2,V2 = tensorSVD(M2, [0,1])
print("V2.shape = ", V2.shape)
M3 = np.tensordot(S2, V2, [1,0])
print("M3.shape = ", M3.shape)
U3,S3,V3 = tensorSVD(M3, [0,1])
print("V3.shape = ", V3.shape)
M4 = np.tensordot(S3, V3, [1,0])
print("M4.shape = ", M4.shape)
Explanation: We can compress the parity indicator into a MPS form by a sequence of SVDs, as explained in the lecture. The following code block
End of explanation
print("U1.shape = ",U1.shape)
print("U2.shape = ",U2.shape)
print("U3.shape = ",U3.shape)
print("M4.shape = ",M4.shape)
Explanation: In this scheme, the four MPS tensors consist of each of the $U$ tensors from the SVDs for all but the final tensor. We can see their shapes here:
End of explanation
print("These are the non-zero elements")
print(round(np.sqrt(2)*U2[0,0,0]))
print(round(np.sqrt(2)*U2[1,1,0]))
print(round(np.sqrt(2)*U2[1,0,1]))
print(round(np.sqrt(2)*U2[0,1,1]))
print("These elements should be zero")
print(round(np.sqrt(2)*U2[0,0,1]))
print(round(np.sqrt(2)*U2[1,1,1]))
print(round(np.sqrt(2)*U2[1,0,0]))
print(round(np.sqrt(2)*U2[0,1,0]))
Explanation: We can interpret the function of these MPS tensors by looking at which entries are nonzero. In particular, for $U2$ we have the following:
End of explanation
# Your code here
Explanation: Part 2
In this section we will walk through the creation of simple Tree Tensor Neworks (TTNs) using the density matrix method. These examples are phrased as multipart exercises, so you might try to code up the solution yourself before looking ahead.
Problem 1
(i)
We begin by defining the tensor 
$$
A_{ijkm} = \sqrt{i+2j+3k+4l + 5m}
$$
where $i,j,k,l, m$ are each two-dimensional indices. That is, each one can take on the value $0$ or $1$. Initialize this tensor as a numpy array.
End of explanation
d = 2 # index dimensions
A = np.zeros((d,d,d,d,d))
for i in range(d):
    for j in range(d):
        for k in range(d):
            for l in range(d):
                for m in range(d):
                    A[i,j,k,l,m] = np.sqrt(i + 2*j + 3*k + 4*l + 5*m)
print(A.shape)
print(A)
 
Explanation: Solution (i)
End of explanation
# Your code here
Explanation: (ii)
Using the appropriate density matrices, form isometries $W_L$ and $W_R$ that transform $A$ into the tree tensor network shown (truncating internal indices to dimension $\chi=2$).
Store $W_L$, $W_R$, and $B$ as arrays. These form the truncated TTN representaiton of $A$.
End of explanation
chi = 2 # Set bond dimension for truncations.
# (a)form left density matrix and isometry
rho_L = A.reshape(d**2,d**3) @ A.reshape(d**2,d**3).T
D_L,U_L = LA.eigh(rho_L)
W_L = (U_L[:,(d**2-chi):]).reshape(d,d,chi)
# (b) form right density matrix and isometry
rho_R = A.reshape(d**3,d**2).T @ A.reshape(d**3,d**2)
D_R,U_R = LA.eigh(rho_R)
W_R = (U_R[:,(d**2-chi):]).reshape(d,d,chi)
# (c) form B tensor
B = np.einsum('ijklm,ijn,lmp->nkp',A,W_L,W_R)
print(B)
Explanation: Solution (ii)
End of explanation
# Your code here
Explanation: (iii)
Check the accuracy of the truncated TTN by computing the truncation error $\epsilon = \|A-A_\text{recover}\|/\|A\|$, with the recovered tensor $A_\text{recover}$ given by contracting the tree.
Here $\| \cdot \|$ is the Frobenius norm that is implemented in numpy as numpy.linalg.norm.
End of explanation
# check truncation error
A_recover = np.einsum('fkg,ijf,lmg->ijklm',B,W_L,W_R)
err_tot = LA.norm(A-A_recover) / LA.norm(A)
print(err_tot)
Explanation: Solution (iii)
End of explanation
chi = 2 # set bond dimension for truncations
# define image data
image_data = np.array([[0,1,1,0,1],[1,1,0,0,1],[1,0,1,0,1],
                      [1,1,0,1,0],[0,1,1,0,1],[1,0,1,0,1]])
n_samples = image_data.shape[0]
n_pixels = image_data.shape[1]
# Your code here
Explanation: Problem 2
In this problem we'll train a TTN to approximate a batch of images, where again the approximation comes from truncating to bond dimension $\chi =2$. The image set we're going to consider consists of five images with five pixels each, but the method is general and you should try repeating the exercise for a different set of images.
(i)
Using the density matrix approach approximate the given set of images using the TTN shown below (truncating internal indices to bond dimension $\chi=2$).
End of explanation
###############################################################
def images_to_sparse(images):
  Creates a sparse tensor representation of image data (one-hot encoding)
  
  Args:
    images: Array of shape (M,N). M is the number of samples, N is the number
      of pixels
      
  Returns:
    pixel_tensors: list of length N, where pixel_tensors[k] is an array of 
      shape (P,M) representing the one-hot encoding of the kth pixels.
  
  
  pixel_tensors = []
  for pixel_loc in range(images.shape[1]):
      pixel_tensors.append(np.array([(1-images[:,pixel_loc]),images[:,pixel_loc]]))
  
  return pixel_tensors
###############################################################
def compute_scalar_products(pixel_tensors):
  Creates the matrix of scalar products from a list of pixel_tensors.
  
  Args:
    pixel_tensors: list of tensors, where pixel_tensors[k] is an array of 
        shape (P,M) representing the one-hot encoding of the kth pixels.
      
  Returns:
    scalar_matrix: array of shape (M,M), formed from the product over pixels k
      of the outer product of pixel_tensors[k] with itself.
  
  
  M = pixel_tensors[1].shape[1]
  scalar_matrix = np.ones((M,M))
  for pixel in pixel_tensors:
      scalar_matrix = scalar_matrix * (pixel.T @  pixel)
  return scalar_matrix
###############################################################
def sparse_to_dense(pixel_tensors):
  Creates a dense tensor from a list of pixel_tensors.
  
  Args:
    pixel_tensors: list of N tensors, where pixel_tensors[k] is an array of 
        shape (P,M) representing the one-hot encoding of the kth pixels.
      
  Returns:
    dense_tensor: array with N dimensions, where dim of kth index kth is equal 
      to the dim P of the kth tensor.
  
  
  dense_tensor = pixel_tensors[0]
  dims_P = [dense_tensor.shape[0]]
  for pixel in pixel_tensors[1:]:
      dim_temp = pixel.shape[0]
      dense_tensor = np.kron(dense_tensor,np.ones((dim_temp,1))) * np.kron(np.ones((np.prod(dims_P),1)),pixel)
      dims_P.append(dim_temp)
  return np.sum(dense_tensor,1).reshape(dims_P)
###############################################################
# export image data to sparse tensor format
V = images_to_sparse(image_data)
# (a) compute left density matrix
gamma_L = compute_scalar_products(V[2:4])
V_L = np.einsum('im,jm->ijm',V[0],V[1])
rho_L = V_L.reshape((4,n_samples)) @ gamma_L @ (V_L.reshape((4,n_samples))).T
# compute left isometry
D_L, U_L = LA.eigh(rho_L)
W_L = U_L[:,(4-chi):].reshape(2,2,chi)
# (b) compute right density matrix 
gamma_R = compute_scalar_products(V[0:2])
V_R = np.einsum('im,jm->ijm',V[3],V[4])
rho_R = V_R.reshape((4,n_samples)) @ gamma_R @ (V_R.reshape((4,n_samples))).T
# compute right isometry
D_R, U_R = LA.eigh(rho_R)
W_R = U_R[:,(4-chi):].reshape(2,2,chi)
# (c) coarse-grain image data using isometries to form B tensor
V_Ltemp = np.einsum('ijk,ijl->kl',W_L,V_L)
V_Rtemp = np.einsum('ijk,ijl->kl',W_R,V_R)
B = sparse_to_dense([V_Ltemp,V[2],V_Rtemp])
Explanation: Solution (i)
End of explanation
# Your code here
Explanation: (ii)
Check the accuracy of your result by contracting the TTN into a single dense tensor and comparing with the dense tensor representation of the original data.
End of explanation
# check accuracy
A_recover = np.einsum('fkg,ijf,lmg->ijklm',B,W_L,W_R)
A_initial = sparse_to_dense(V)
err_tot = LA.norm(A_recover-A_initial)
print(err_tot)
Explanation: Solution (ii)
We can check the accuracy using the Frobenius norm, just like in the previous problem.
In this case, the example was chosen to allow perfect compression with $\chi=2$, so the total error just comes from machine precision.
End of explanation
<END_TASK>