Decompositions of symmetric tensors#
The examples below demonstrate common decomposition routines for symmetric tensors.
import yastn
import pytest
config_kwargs = {"backend": "np"}
SVD decompositions and truncation#
The example below demonstrates SVD-based decomposition and truncation of a symmetric tensor.
def test_svd_truncate_lowrank(config_kwargs):
"""Check SVD and lowrank SVD combined with truncation."""
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)
# Decompose the tensor and build a reference tensor from the factors.
U, S, V = yastn.linalg.svd(a, axes=((0, 1), (2, 3)), sU=-1)
# Fix the singular values explicitly so that truncation behavior is
# deterministic and easy to reason about in the test.
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)])
# First check that global truncation reduces the spectrum to the
# requested total bond dimension. Here low-rank SVD driver is used.
opts = {'tol': 0.01, 'D_block': 100, 'D_total': 12}
_, S2, _ = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), sU=-1, **opts, policy='lowrank')
assert S2.get_shape() == (12, 12)
# Then verify that truncation by charge produces the same result when
# the left and right factors are split with different charge conventions.
opts = {'D_block': {(0,): 2, (-1,): 0}, 'policy': 'lowrank'}
U1, S1, V1 = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), nU=True, sU=-1, **opts)
assert S1.get_shape() == (2, 2)
a1 = U1 @ S1 @ V1
# Truncation by charge requires care when assigning charges to the
# new factors, especially for different values of nU and sU.
opts = {'D_block': {(1,): 2}, 'policy': 'lowrank'}
U2, S2, V2 = yastn.linalg.svd_with_truncation(a, axes=((0, 1), (2, 3)), nU=False, sU=-1, **opts)
assert S1.get_shape() == (2, 2)
a2 = U2 @ S2 @ V2
assert yastn.norm(a1 - a2) < tol
QR decompositions#
The example below takes a tensor a with four legs, decomposes it using QR, and contracts the resulting Q and R tensors back into a.
def run_qr_combine(a):
"""Decompose and contract tensor ``a`` using a QR decomposition."""
assert a.ndim == 4
def check_diag_R_nonnegative(R):
"""Check that the diagonal of R is chosen to be non-negative."""
for t in R.get_blocks_charge():
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 the signature of the new leg and its placement.
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-based eigs method available in SciPy. Tensor operations can be passed to other SciPy methods in a similar way, although this is currently 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) # e.g., it could be an MPS tensor
a, _ = yastn.qr(a, axes=((0, 1), 2), sQ=-1) # orthonormalize
# Build the dense transfer matrix from a as 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
# Initialize a random tensor matching the transfer matrix from the left.
# Add an extra third leg carrying charges -1, 0, and 1 to solve the
# eigensystem over these three 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 maps a 1-D vector to a YASTN tensor, applies the
# transfer-matrix action, and returns a 1-D vector.
r1d, meta = yastn.split_data_and_meta(v0.to_dict(level=0), squeeze=True)
def f(x):
t = yastn.Tensor.from_dict(yastn.combine_data_and_meta(x, meta))
t2 = yastn.ncon([t, a, a.conj()], [(1, 3, -3), (1, 2, -1), (3, 2, -2)])
t3, _ = yastn.split_data_and_meta(t2.to_dict(level=0, meta=meta), squeeze=True)
return t3
ff = LinearOperator(shape=(len(r1d), len(r1d)), matvec=f, dtype=np.float64)
# Apply SciPy's sparse eigensolver to a YASTN symmetric tensor.
wa, va1d = eigs(ff, v0=r1d, k=1, which='LM', tol=1e-10)
# Transform the eigenvectors into YASTN tensors.
va = [yastn.Tensor.from_dict(yastn.combine_data_and_meta(x, meta)) for x in va1d.T]
# Remove zero blocks now; some eigenvectors have well-defined charge,
# although a superposition of symmetry sectors may appear in degenerate cases.
va = [x.remove_zero_blocks() for x in va]
# We can also restrict the search directly to eigenvectors with the desired
# charge, here n=0.
legs = [a.get_legs(0).conj(),
a.get_legs(0)]
v0 = yastn.rand(config=a.config, legs=legs, n=0)
r1d, meta = yastn.split_data_and_meta(v0.to_dict(level=0), squeeze=True)
def f(x):
t = yastn.Tensor.from_dict(yastn.combine_data_and_meta(x, meta))
t2 = yastn.ncon([t, a, a.conj()], [(1, 3), (1, 2, -1), (3, 2, -2)])
t3, _ = yastn.split_data_and_meta(t2.to_dict(level=0, meta=meta), squeeze=True)
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.Tensor.from_dict(yastn.combine_data_and_meta(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())