Svd_truncate/Gesvd_truncate return wrong discarded singular values for symmetric UniTensor

[ianmccul]

A bug spotted by my student Jerry Chen:

Summary

cytnx.linalg.Svd_truncate on symmetric (UTenType.Block) UniTensor appears to return the wrong value in the return_err path.

For block/symmetric tensors, the "discarded singular values" output can overlap with or match the kept singular values. This is inconsistent with the dense tensor path and appears to come from
selecting the kept tail of the sorted singular spectrum instead of the discarded prefix.

The same issue appears in both:

• linalg::Svd_truncate(const UniTensor&, ...)
• linalg::Gesvd_truncate(const UniTensor&, ...)

Versions

This was diagnosed looking at current master. However the bug was spotted using version 1.0.1. The bug has existed since the return_err mechanism was added in commit c8a65dc.

Actual behavior

For symmetric/block UniTensor:

• return_err > 1 returns Sall.get({ac::tilend(smidx)}), which is the kept tail of the sorted singular values, not the discarded values.
• return_err == 1 returns Smin, which is the smallest kept singular value.

This behavior is incorrect in both paths.

Expected behavior

The symmetric/block UniTensor path should match the dense tensor implementation:

• return_err > 1 should return the discarded singular values only.
• return_err == 1 should return the first discarded singular value, not the smallest kept one. This is documented in arxiv:2401.01921 as returning "the largest of the truncated singular values".

Relevant source locations

Block path in Svd_truncate:

• src/linalg/Svd_truncate.cpp:312-316

Block path in Gesvd_truncate:

• src/linalg/Gesvd_truncate.cpp:354-358

Dense/tensor reference behavior:

• src/backend/linalg_internal_cpu/memcpyTruncation.cpp:197-207

Root cause

Both block implementations sort all singular values ascending and compute a truncation threshold Smin with index smidx, then do:

outCyT.push_back(UniTensor(Sall.get({ac::tilend(smidx)})));

That slice starts at smidx and runs to the end, which corresponds to the kept values, not the discarded ones.

Similarly, for return_err == 1, the code returns Smin, which is the truncation threshold / smallest kept value, while the dense path returns the first discarded value.

Additional issue noticed while tracing

There also seems to be a separate mindim inconsistency in the block implementations:

• src/linalg/Svd_truncate.cpp:193-198 ignores mindim in one branch
• src/linalg/Gesvd_truncate.cpp:246-250 ignores mindim in the opposite branch

This looks independent of the discarded-spectrum bug, but likely belongs in the same area.

[ianmccul]

Reproducer of the bug:

#!/usr/bin/env python3
import os, sys, numpy as np
sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '../CytnxTools')))
sys.path.insert(0, "/home/jerrychen/Cytnx/Cytnx_lib")

import cytnx

print("cytnx version =", cytnx.__version__)

def to_nparray(T):
    assert type(T) == cytnx.UniTensor
    if T.is_blockform():
        tmp = cytnx.UniTensor.zeros(T.shape())
        tmp.convert_from(T)
        T = tmp
    return T.get_block().numpy()

def show(tag, M, keepdim=20, cutoff=0.0):
    print("\n" + "=" * 80)
    print(tag)
    print("=" * 80)
    print("input: is_blockform =", M.is_blockform(), "norm^2 =", M.Norm().item()**2)

    S, U, Vt, errs = cytnx.linalg.Svd_truncate(M, keepdim=keepdim, err=cutoff, return_err=2)

    S_np = to_nparray(S)
    e_np = to_nparray(errs)

    sv = np.diag(S_np) if (S_np.ndim == 2 and S_np.shape[0] == S_np.shape[1]) else S_np.ravel()
    ev = e_np.ravel()

    print("S:")
    print(S)
    print("errs:")
    print(errs)
    print("sv =", sv)
    print("ev =", ev)
    print("sum(sv^2) =", np.sum(sv**2))
    print("sum(ev^2) =", np.sum(ev**2))
    print("same sorted? =", len(sv) == len(ev) and np.allclose(np.sort(sv), np.sort(ev)))

def make_dense():
    M = cytnx.UniTensor(cytnx.arange(24).reshape(2,2,6))
    M.relabels_(["a","b","c"])
    M.set_rowrank_(2)
    return M

def make_z2():
    Sym = cytnx.Symmetry.Zn(2)
    Bond, Qs = cytnx.Bond, cytnx.Qs

    bl = Bond(cytnx.BD_IN,  [Qs(0), Qs(1)], [2, 2], [Sym])
    br = Bond(cytnx.BD_OUT, [Qs(0), Qs(1)], [2, 2], [Sym])

    T = cytnx.UniTensor([bl, br], labels=["l", "r"])

    rng = np.random.default_rng(1234)
    nblk = T.get_total_blocks() if hasattr(T, "get_total_blocks") else T.Nblocks()

    for b in range(nblk):
        old_blk = T.get_block_(b) if hasattr(T, "get_block_") else T.get_block(b)
        arr = rng.normal(size=old_blk.shape())
        new_blk = cytnx.from_numpy(arr)

        if hasattr(T, "put_block_"):
            T.put_block_(new_blk, b)
        else:
            T.put_block(new_blk, b)

    return T

show("dense, no truncation expected", make_dense(), keepdim=20, cutoff=0.0)
show("dense, truncation expected",    make_dense(), keepdim=2,  cutoff=0.0)
show("Z2, no truncation expected",    make_z2(),    keepdim=20, cutoff=0.0)
show("Z2, truncation expected",       make_z2(),    keepdim=2,  cutoff=0.0)