Decompositions of symmetric tensors

import yastn
import pytest
config_kwargs = {"backend": "np"}

SVD decompositions and truncation

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)  # e.g., it could be an 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.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)
    # 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.Tensor.from_dict(yastn.combine_data_and_meta(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 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())