Decompositions of symmetric tensors#
import yastn
import pytest
config_kwargs = {"backend": "np"}
SVD decompositions and truncation#
def test_svd_truncate(config_kwargs):
#
# Start with random tensor with 4 legs.
config_U1 = yastn.make_config(sym='U1', **config_kwargs)
legs = [yastn.Leg(config_U1, s=1, t=(0, 1), D=(5, 6)),
yastn.Leg(config_U1, s=1, t=(-1, 0), D=(5, 6)),
yastn.Leg(config_U1, s=-1, t=(-1, 0, 1), D=(2, 3, 4)),
yastn.Leg(config_U1, s=-1, t=(-1, 0, 1), D=(2, 3, 4))]
a = yastn.rand(config=config_U1, n=1, legs=legs)
#
# Fixing singular values for testing.
# Create new tensor *a* that will be used for testing.
U, S, V = yastn.linalg.svd(a, axes=((0, 1), (2, 3)), sU=-1)
S.set_block(ts=(-2, -2), Ds=4, val=[2**(-ii - 6) for ii in range(4)])
S.set_block(ts=(-1, -1), Ds=12, val=[2**(-ii - 2) for ii in range(12)])
S.set_block(ts=(0, 0), Ds=25, val=[2**(-ii - 1) for ii in range(25)])
a = yastn.ncon([U, S, V], [(-1, -2, 1), (1, 2), (2, -3, -4)])
#
# Opts for truncation; truncates below (global) tolerance tol with respect to largest singular value
# up to a total of D_total singular values, with at mot D_block singular values for each charge sector.
opts = {'tol': 0.01, 'D_block': 100, 'D_total': 12}
U1, S1, V1 = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=-1, **opts) # nU=True, total charge of *a* goes with U
assert S1.get_blocks_charge() == ((-2, -2), (-1, -1), (0, 0))
assert S1.s[0] == 1 and U1.n == a.n and V1.n == (0,) and S1.get_shape() == (12, 12)
#
# svd_with_truncation is a shorthand for
U, S, V = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=-1) # nU=True
mask = yastn.linalg.truncation_mask(S, **opts)
U1p, S1p, V1p = mask.apply_mask(U, S, V, axes=(2, 0, 0))
assert all((x - y).norm() < 1e-12 for x, y in ((U1, U1p), (S1, S1p), (V1, V1p)))
#
# Specific charges of S1 depend on signature of a new leg comming from U,
# and where charge of tensor 'a' is attached
U1, S1, V1 = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), **opts) # sU=1, nU=True
assert S1.get_blocks_charge() == ((0, 0), (1, 1), (2, 2))
assert S1.s[0] == -1 and U1.n == a.n and V1.n == (0,) and S1.get_shape() == (12, 12)
U1, S1, V1 = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=-1, nU=False, **opts)
assert S1.get_blocks_charge() == ((-1, -1), (0, 0), (1, 1))
assert S1.s[0] == 1 and U1.n == (0,) and V1.n == a.n and S1.get_shape() == (12, 12)
U1, S1, V1 = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=1, nU=False, **opts)
assert S1.get_blocks_charge() == ((-1, -1), (0, 0), (1, 1))
assert S1.s[0] == -1 and U1.n == (0,) and V1.n == a.n and S1.get_shape() == (12, 12)
#
# Different truncation options.
opts = {'tol': 0.02, 'D_block': 5, 'D_total': 100}
_, S1, _ = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=1, **opts)
assert S1.get_shape() == (11, 11)
#
# tol_block decides on block truncation tolerance with respect to largest singular value of each block
opts = {'tol': 0, 'tol_block': 0.2, 'D_total': 100}
_, S1, _ = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=1, **opts)
assert S1.get_shape() == (9, 9)
#
# We can specify D_block and tol_block for each charge sector independently
opts = {'D_block': {(0,): 2, (-1,): 3}, 'tol_block': {(0,): 0.6, (-1,): 0.1}}
_, S1, _ = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=-1, **opts)
assert S1.get_shape() == (4, 4)
#
# Here use truncate_multiplets option to shift the cut to largest D to retain multiplet
opts = {'D_total': 10, "truncate_multiplets": True}
_, S1, _ = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=-1, **opts)
assert S1.get_shape() == (12, 12)
#
# Empty tensor
opts = {'D_total': 0}
_, S2, _ = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), nU=False, sU=-1, **opts)
assert S2.norm() < 1e-12 and S2.size == 0
QR decompositions#
The function below takes tensor a with 4 legs, decompose it using QR and contracts the resulting Q and R tensors back into a.
def run_qr_combine(a):
""" decompose and contracts tensor ``a`` using qr decomposition """
assert a.ndim == 4
def check_diag_R_nonnegative(R):
""" checks that diagonal of R is selected to be non-negative """
for t in R.struct.t:
assert all(R.config.backend.diag_get(R.real()[t]) >= 0)
assert all(R.config.backend.diag_get(R.imag()[t]) == 0)
Q, R = yastn.linalg.qr(a, axes=((3, 1), (2, 0)))
QR = yastn.tensordot(Q, R, axes=(2, 0))
QR = QR.transpose(axes=(3, 1, 2, 0))
assert yastn.norm(a - QR) < 1e-12 # == 0.0
assert Q.is_consistent()
assert R.is_consistent()
check_diag_R_nonnegative(R.fuse_legs(axes=(0, (1, 2)), mode='hard'))
# change signature of new leg; and position of new leg
Q2, R2 = yastn.qr(a, axes=((3, 1), (2, 0)), sQ=-1, Qaxis=0, Raxis=-1)
QR2 = yastn.tensordot(R2, Q2, axes=(2, 0)).transpose(axes=(1, 3, 0, 2))
assert yastn.norm(a - QR2) < 1e-12 # == 0.0
assert Q2.is_consistent()
assert R2.is_consistent()
check_diag_R_nonnegative(R2.fuse_legs(axes=((0, 1), 2), mode='hard'))
Combining with scipy.sparse.linalg.eigs#
Calculate the dominant eigenvector of a transfer matrix by employing the Krylov-base eigs method available in SciPy. Tensor operations can be similarly passed to other SciPy methods, though this is limited to the NumPy backend.
import numpy as np
from scipy.sparse.linalg import eigs, LinearOperator
@numpy_test
def test_eigs_simple(config_kwargs):
config_U1 = yastn.make_config(sym='U1', **config_kwargs)
legs = [yastn.Leg(config_U1, s=1, t=(-1, 0, 1), D=(2, 3, 2)),
yastn.Leg(config_U1, s=1, t=(0, 1), D=(1, 1)),
yastn.Leg(config_U1, s=-1, t=(-1, 0, 1), D=(2, 3, 2))]
a = yastn.rand(config=config_U1, legs=legs) # could be mps tensor
a, _ = yastn.qr(a, axes=((0, 1), 2), sQ=-1) # orthonormalize
# Dense transfer matrix build from a; reference solution
tm = yastn.ncon([a, a.conj()], [(-1, 1, -3), (-2, 1, -4)])
tm = tm.fuse_legs(axes=((2, 3), (0, 1)), mode='hard')
tmn = tm.to_numpy()
w_ref, v_ref = eigs(tmn, k=1, which='LM') # use scipy.sparse.linalg.eigs
# Initializing random tensor matching tm from left.
# We add an extra 3-rd leg carrying charges -1, 0, 1
# to calculate eigs over those 3 subspaces in one go.
legs = [a.get_legs(0).conj(),
a.get_legs(0),
yastn.Leg(a.config, s=1, t=(-1, 0, 1), D=(1, 1, 1))]
v0 = yastn.rand(config=a.config, legs=legs)
# Define a wrapper that goes r1d -> yastn.tensor -> tm @ yastn.tensor -> r1d
r1d, meta = yastn.compress_to_1d(v0)
def f(x):
t = yastn.decompress_from_1d(x, meta=meta)
t2 = yastn.ncon([t, a, a.conj()], [(1, 3, -3), (1, 2, -1), (3, 2, -2)])
t3, _ = yastn.compress_to_1d(t2, meta=meta)
return t3
ff = LinearOperator(shape=(len(r1d), len(r1d)), matvec=f, dtype=np.float64)
# scipy.sparse.linalg.eigs that goes though yastn symmetric tensor.
wa, va1d = eigs(ff, v0=r1d, k=1, which='LM', tol=1e-10)
# Transform eigenvectors into yastn tensors
va = [yastn.decompress_from_1d(x, meta) for x in va1d.T]
# We can remove zero blocks now, as there are eigenvectors with well defined charge
# (though we might get superposition of symmetry sectors in case of degeneracy).
va = [x.remove_zero_blocks() for x in va]
# we can also limit ourselves directly to eigenvectors with desired charge, here 0.
legs = [a.get_legs(0).conj(),
a.get_legs(0)]
v0 = yastn.rand(config=a.config, legs=legs, n=0)
r1d, meta = yastn.compress_to_1d(v0)
def f(x):
t = yastn.decompress_from_1d(x, meta=meta)
t2 = yastn.ncon([t, a, a.conj()], [(1, 3), (1, 2, -1), (3, 2, -2)])
t3, _ = yastn.compress_to_1d(t2, meta=meta)
return t3
ff = LinearOperator(shape=(len(r1d), len(r1d)), matvec=f, dtype=np.float64)
wb, vb1d = eigs(ff, v0=r1d, k=1, which='LM', tol=1e-10) # scipy.sparse.linalg.eigs
vb = [yastn.decompress_from_1d(x, meta) for x in vb1d.T] # eigenvectors as yastn tensors
# dominant eigenvalue should have amplitude 1 (likely degenerate in our example)
assert all(pytest.approx(abs(x), rel=1e-10) == 1.0 for x in (w_ref, wa, wb))
print("va -> ", va.pop())
print("vb -> ", vb.pop())