# Toycode a_mps.pyΒΆ

"""Toy code implementing a matrix product state."""
# Copyright 2018-2020 TeNPy Developers, GNU GPLv3

import numpy as np
from scipy.linalg import svd
# if you get an error message "LinAlgError: SVD did not converge",
# uncomment the following line. (This requires TeNPy to be installed.)
#  from tenpy.linalg.svd_robust import svd  # (works like scipy.linalg.svd)

class SimpleMPS:
"""Simple class for a matrix product state.

We index sites with i from 0 to L-1; bond i is left of site i.
We *assume* that the state is in right-canonical form.

Parameters
----------
Bs, Ss, bc:
Same as attributes.

Attributes
----------
Bs : list of np.Array[ndim=3]
The 'matrices' in right-canonical form, one for each physical site
(within the unit-cell for an infinite MPS).
Each B[i] has legs (virtual left, physical, virtual right), in short vL i vR
Ss : list of np.Array[ndim=1]
The Schmidt values at each of the bonds, Ss[i] is left of Bs[i].
bc : 'infinite', 'finite'
Boundary conditions.
L : int
Number of sites (in the unit-cell for an infinite MPS).
nbonds : int
Number of (non-trivial) bonds: L-1 for 'finite' boundary conditions
"""
def __init__(self, Bs, Ss, bc='finite'):
assert bc in ['finite', 'infinite']
self.Bs = Bs
self.Ss = Ss
self.bc = bc
self.L = len(Bs)
self.nbonds = self.L - 1 if self.bc == 'finite' else self.L

def copy(self):
return SimpleMPS([B.copy() for B in self.Bs], [S.copy() for S in self.Ss], self.bc)

def get_theta1(self, i):
"""Calculate effective single-site wave function on sites i in mixed canonical form.

The returned array has legs vL, i, vR (as one of the Bs).
"""
return np.tensordot(np.diag(self.Ss[i]), self.Bs[i], [1, 0])  # vL [vL'], [vL] i vR

def get_theta2(self, i):
"""Calculate effective two-site wave function on sites i,j=(i+1) in mixed canonical form.

The returned array has legs vL, i, j, vR.
"""
j = (i + 1) % self.L
return np.tensordot(self.get_theta1(i), self.Bs[j], [2, 0])  # vL i [vR], [vL] j vR

def get_chi(self):
"""Return bond dimensions."""
return [self.Bs[i].shape[2] for i in range(self.nbonds)]

def site_expectation_value(self, op):
"""Calculate expectation values of a local operator at each site."""
result = []
for i in range(self.L):
theta = self.get_theta1(i)  # vL i vR
op_theta = np.tensordot(op, theta, axes=[1, 1])  # i [i*], vL [i] vR
result.append(np.tensordot(theta.conj(), op_theta, [[0, 1, 2], [1, 0, 2]]))
# [vL*] [i*] [vR*], [i] [vL] [vR]
return np.real_if_close(result)

def bond_expectation_value(self, op):
"""Calculate expectation values of a local operator at each bond."""
result = []
for i in range(self.nbonds):
theta = self.get_theta2(i)  # vL i j vR
op_theta = np.tensordot(op[i], theta, axes=[[2, 3], [1, 2]])
# i j [i*] [j*], vL [i] [j] vR
result.append(np.tensordot(theta.conj(), op_theta, [[0, 1, 2, 3], [2, 0, 1, 3]]))
# [vL*] [i*] [j*] [vR*], [i] [j] [vL] [vR]
return np.real_if_close(result)

def entanglement_entropy(self):
"""Return the (von-Neumann) entanglement entropy for a bipartition at any of the bonds."""
bonds = range(1, self.L) if self.bc == 'finite' else range(0, self.L)
result = []
for i in bonds:
S = self.Ss[i].copy()
S[S < 1.e-20] = 0.  # 0*log(0) should give 0; avoid warning or NaN.
S2 = S * S
assert abs(np.linalg.norm(S) - 1.) < 1.e-14
result.append(-np.sum(S2 * np.log(S2)))
return np.array(result)

def correlation_length(self):
"""Diagonalize transfer matrix to obtain the correlation length."""
import scipy.sparse.linalg.eigen.arpack as arp
assert self.bc == 'infinite'  # works only in the infinite case
B = self.Bs[0]  # vL i vR
chi = B.shape[0]
T = np.tensordot(B, np.conj(B), axes=[1, 1])  # vL [i] vR, vL* [i*] vR*
T = np.transpose(T, [0, 2, 1, 3])  # vL vL* vR vR*
for i in range(1, self.L):
B = self.Bs[i]
T = np.tensordot(T, B, axes=[2, 0])  # vL vL* [vR] vR*, [vL] i vR
T = np.tensordot(T, np.conj(B), axes=[[2, 3], [0, 1]])
# vL vL* [vR*] [i] vR, [vL*] [i*] vR*
T = np.reshape(T, (chi**2, chi**2))
# Obtain the 2nd largest eigenvalue
eta = arp.eigs(T, k=2, which='LM', return_eigenvectors=False, ncv=20)
return -self.L / np.log(np.min(np.abs(eta)))

def init_FM_MPS(L, d, bc='finite'):
"""Return a ferromagnetic MPS (= product state with all spins up)"""
B = np.zeros([1, d, 1], np.float)
B[0, 0, 0] = 1.
S = np.ones([1], np.float)
Bs = [B.copy() for i in range(L)]
Ss = [S.copy() for i in range(L)]
return SimpleMPS(Bs, Ss, bc)

def split_truncate_theta(theta, chi_max, eps):
"""Split and truncate a two-site wave function in mixed canonical form.

Split a two-site wave function as follows::
vL --(theta)-- vR     =>    vL --(A)--diag(S)--(B)-- vR
|   |                       |             |
i   j                       i             j

Afterwards, truncate in the new leg (labeled vC).

Parameters
----------
theta : np.Array[ndim=4]
Two-site wave function in mixed canonical form, with legs vL, i, j, vR.
chi_max : int
Maximum number of singular values to keep
eps : float
Discard any singular values smaller than that.

Returns
-------
A : np.Array[ndim=3]
Left-canonical matrix on site i, with legs vL, i, vC
S : np.Array[ndim=1]
Singular/Schmidt values.
B : np.Array[ndim=3]
Right-canonical matrix on site j, with legs vC, j, vR
"""
chivL, dL, dR, chivR = theta.shape
theta = np.reshape(theta, [chivL * dL, dR * chivR])
X, Y, Z = svd(theta, full_matrices=False)
# truncate
chivC = min(chi_max, np.sum(Y > eps))
piv = np.argsort(Y)[::-1][:chivC]  # keep the largest chivC singular values
X, Y, Z = X[:, piv], Y[piv], Z[piv, :]
# renormalize
S = Y / np.linalg.norm(Y)  # == Y/sqrt(sum(Y**2))
# split legs of X and Z
A = np.reshape(X, [chivL, dL, chivC])
B = np.reshape(Z, [chivC, dR, chivR])
return A, S, B