Algebra#
Canonical form by QR#
There are different algorithms, which can bring MPS/MPO into a canonical form.
The cheapest way is by application of QR decomposition.
# Generate random MPS with no symmetry
psi = mps.random_dense_mps(N=16, D=15, d=2)
# rigth canonical form
#
# --A*-- --
# | | = |Identity
# --A--- --
#
psi.canonize_(to='first')
# left canonical form
#
# --A*-- --
# | | = Identity|
# --A--- --
#
psi.canonize_(to='last')
Check if MPS/MPO is in left/right canonical form by verifying if each tensor forms an isometry after appropriate contraction with its conjugate.
def check_canonize(psi, tol):
""" Canonize mps to left and right, running tests if it is canonical. """
ref_s = (-1, 1, 1) if psi.nr_phys == 1 else (-1, 1, 1, -1)
norm = psi.norm()
for to in ('last', 'first'):
psi.canonize_(to=to, normalize=False)
assert psi.is_canonical(to=to, tol=tol)
assert all(psi[site].s == ref_s for site in psi.sweep())
assert psi.pC is None
assert len(psi.A) == len(psi)
assert abs(psi.factor / norm - 1) < tol
assert abs(mps.vdot(psi, psi) / norm ** 2 - 1) < tol
for to in ('last', 'first'):
phi = psi.shallow_copy()
phi.canonize_(to=to)
assert abs(phi.factor - 1) < tol
assert abs(mps.vdot(phi, phi) - 1) < tol
Canonical form by SVD#
Bringing MPS/MPO into canonical form through SVD decomposition
is computationally more expensive than QR, but allows for truncation.
Truncation is governed by options passed as opts_dict
to yastn.linalg.truncation_mask() that is applied
after each SVD during the sweep through MPS/MPO.
Note
Faithfull SVD truncation requires MPS/MPO to be in the canonical form of the opposite direction to the direction of the truncation sweep.
# There are different options which we can pass, see `truncation_mask`.
# Defaults are assumed for options not explictly specified.
#
opts_svd = {
'D_total': 4, # total number of singular values to keep
'D_block': 2, # maximal number of singular values
# to keep in a single block
'tol': 1e-6, # relative tolerance of singular values
# below which to truncate across all blocks
'tol_blocks': 1e-6 # relative tolerance of singular values
# within individual blocks below which to
}
# Generate random MPS with no symmetry
#
psi = mps.random_dense_mps(N=16, D=15, d=2)
# Bring MPS to canonical form and truncate (here, right canonical form).
# For MPS we usually normalize the state.
#
psi.canonize_(to='last')
psi.truncate_(to='first', opts_svd=opts_svd)
# Generate random MPO with no symmetry
#
H = mps.random_dense_mpo(N=16, D=25, d=2)
# Bring MPO to canonical form and truncate (here, left canonical form).
# For MPO we do not want to change overall scale, thus no normalization.
#
H.canonize_(to='first', normalize=False)
H.truncate_(to='last', opts_svd=opts_svd, normalize=False)
Multiplication#
We can test the multiplication of MPO and MPS using an example of a ground state obtained with DMRG.
def multiplication_example_gs(config=None, tol=1e-12):
"""
Calculate ground state and tests, within eigen-condition,
functions mps.multiply, __mul__, mps.zipper and mps.add
This test presents multiplication as part of the DMRG study.
We use multiplication to get expectation values from a state.
Knowing the exact solution, we will compare it to the value we obtain.
"""
N = 7
Eng = -3.427339492125
t, mu = 1.0, 0.2
#
opts_config = {} if config is None else \
{'backend': config.backend,
'default_device': config.default_device}
# pytest uses config to inject various backends and devices for testing
ops = yastn.operators.SpinlessFermions(sym='U1', **opts_config)
#
# The Hamiltonian is obtained using Hterm.
#
I = mps.product_mpo(ops.I(), N) # identity MPO
terms = []
for i in range(N - 1):
terms.append(mps.Hterm(t, (i, i+1), (ops.cp(), ops.c())))
terms.append(mps.Hterm(t, (i+1, i), (ops.cp(), ops.c())))
for i in range(N):
terms.append(mps.Hterm(mu, [i], [ops.cp() @ ops.c()]))
H = mps.generate_mpo(I, terms)
#
# To standardize this test we fix a seed for random MPS we use.
#
ops.random_seed(seed=0)
#
# In this example we use yastn.Tensor with U1 symmetry.
#
total_charge = 3
psi = mps.random_mps(I, D_total=8, n=total_charge)
#
# We set truncation for DMRG and run the algorithm in '2site' version.
#
opts_svd = {'D_total': 8}
opts_dmrg = {'max_sweeps': 20, 'Schmidt_tol': 1e-14}
out = mps.dmrg_(psi, H, method='2site', **opts_dmrg, opts_svd=opts_svd)
out = mps.dmrg_(psi, H, method='1site', **opts_dmrg)
#
# Test if we obtained the exact solution for energy?:
#
assert abs(out.energy - Eng) < tol
#
# We should have the ground state.
# Now we calculate the variation of energy <H^2>-<H>^2=<(H-Eng)^2>
# to check if DMRG converged properly to tol.
# There are a few ways to do that:
#
# case 1/
# Find approximate Mps representing H @ psi. This is done in two steps.
# First, zipper function gives a good first approximation
# though a series of svd.
#
Hpsi = mps.zipper(H, psi, opts_svd={"D_total": 8}) # here, Hpsi is normalized
#
# Second, compression is optimizing the overlap
# between approximation and a product H @ psi.
#
mps.compression_(Hpsi, (H, psi), method='2site', max_sweeps=5,
Schmidt_tol=1e-6, opts_svd={"D_total": 8},
normalize=False) # norm of Hpsi approximates norm of H @ psi
#
# Use mps.norm() to get variation.
# Norm canonizes a copy of the state and, close to zero,
# is more precise than direct contraction using vdot.
#
p0 = Eng * psi - Hpsi
assert p0.norm() < tol
#
# case 2/
# Here H @ psi is calculated exactly with resulting
# bond dimension being a product of bond dimensions of H and psi.
#
Hpsi = mps.multiply(H, psi)
p0 = Eng * psi - Hpsi
assert p0.norm() < tol
#
# Equivalently we can call.
Hpsi = H @ psi
p0 = Eng * psi - Hpsi
assert p0.norm() < tol