Decompositions of symmetric tensors#
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
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) < tol # == 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) < tol # == 0.0
assert Q2.is_consistent()
assert R2.is_consistent()
check_diag_R_nonnegative(R2.fuse_legs(axes=((0, 1), 2), mode='hard'))
SVD decompositions and truncation#
def test_svd_truncate():
#
# start with random tensor with 4 legs
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
# creat 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() < tol 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() < tol and S2.size == 0
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.
@pytest.mark.skipif(not config_U1.backend.BACKEND_ID=="numpy", reason="uses scipy for raw data")
def test_eigs_simple():
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=tol) == 1.0 for x in (w_ref, wa, wb))
print("va -> ", va.pop())
print("vb -> ", vb.pop())