9. Tensor decomposition

As mention in the Manipulating a UniTensor, the specification of rowrank makes it convenient to apply linear algebra operations on UniTensors.

9.1. Singular value decomposition

Here is an example where a singular value decomposition (SVD) is performed on a UniTensor:

  • In Python:

 1T = cytnx.UniTensor.ones([5,5,5,5,5], rowrank = 3)
 2S, U, Vt = cytnx.linalg.Svd(T)
 3U.set_name('U')
 4S.set_name('S')
 5Vt.set_name('Vt')
 6
 7T.print_diagram()
 8S.print_diagram()
 9U.print_diagram()
10Vt.print_diagram()

Output >>

-----------------------
tensor Name :
tensor Rank : 5
block_form  : False
is_diag     : False
on device   : cytnx device: CPU
          ---------
         /         \
   0 ____| 5     5 |____ 3
         |         |
   1 ____| 5     5 |____ 4
         |         |
   2 ____| 5       |
         \         /
          ---------
-----------------------
tensor Name : S
tensor Rank : 2
block_form  : False
is_diag     : True
on device   : cytnx device: CPU
               -----------
              /           \
   _aux_L ____| 25     25 |____ _aux_R
              \           /
               -----------
-----------------------
tensor Name : U
tensor Rank : 4
block_form  : False
is_diag     : False
on device   : cytnx device: CPU
          ----------
         /          \
   0 ____| 5     25 |____ _aux_L
         |          |
   1 ____| 5        |
         |          |
   2 ____| 5        |
         \          /
          ----------
-----------------------
tensor Name : Vt
tensor Rank : 3
block_form  : False
is_diag     : False
on device   : cytnx device: CPU
               ----------
              /          \
   _aux_R ____| 25     5 |____ 3
              |          |
              |        5 |____ 4
              \          /
               ----------

When calling Svd, the first three legs of T are automatically reshaped into one leg according to rowrank=3. After the SVD, the matrices U and Vt are automatically reshaped back into the corresponding index form of the original tensor. This way, we can reconstruct the original UniTensor T by contracting \(U \cdot S \cdot Vt\):

  • In Python:

1Tsymm_diff=T-cytnx.Contract([U,S,Vt])
2print(Tsymm_diff.Norm()/T.Norm())

Output >>

Total elem: 1
type  : Double (Float64)
cytnx device: CPU
Shape : (1)
[7.87947e-16 ]

If we contract \(U \cdot S \cdot Vt\), we get a tensor of the same shape as T and we can subtract the two tensors. The error \(\frac{|T-U \cdot S \cdot Vt|}{|T|}\) is of the order of machine precision, as expected.

If the singular values have a decaying hierarchy, then the smallest values can be omitted without introducing a large error to the reconstructed tensor. More precisely, the 2-norm of the error between the original tensor and the reconstructed tensor is minimized. This truncated SVD is the basis of many tensor network algorithms.

Here we demonstrate the usage of Svd_truncate(), which implements this truncation. In this example we print the singular values obtained by Svd(), and compare them to the result of Svd_truncate():

  • In Python:

1T = cytnx.UniTensor.ones([5,5,5,5,5], rowrank = 3)
2S1, _, _ = cytnx.linalg.Svd(Tin = T)
3print(S1)
4S2, _, _, err = cytnx.linalg.Svd_truncate(Tin = T, keepdim = 10, err = 1e-12, is_UvT = True, return_err = 1)
5print(S2)
6print(err)

Output >>

-------- start of print ---------
Tensor name:
is_diag    : True
contiguous : True

Total elem: 25
type  : Double (Float64)
cytnx device: CPU
Shape : (25)
[5.59017e+01 4.65490e-15 9.57650e-46 3.07595e-77 8.59824e-108 6.24040e-138 1.64397e-168 2.17184e-199 6.65227e-230 2.54214e-261 8.00572e-292 1.46017e-308 5.43472e-323 5.43472e-323 5.43472e-323 4.94066e-323 4.94066e-323 4.94066e-323 4.94066e-323 4.94066e-323 4.94066e-323 4.94066e-323 4.94066e-323 4.94066e-323 4.94066e-323 ]




-------- start of print ---------
Tensor name:
is_diag    : True
contiguous : True

Total elem: 1
type  : Double (Float64)
cytnx device: CPU
Shape : (1)
[5.59017e+01 ]




-------- start of print ---------
Tensor name:
is_diag    : False
contiguous : True

Total elem: 1
type  : Double (Float64)
cytnx device: CPU
Shape : (1)
[4.65490e-15 ]

We note that the singular values obtained by Svd_truncate() are truncated according to our err requirement, a keepdim argument is also passed to specify a maximum desired dimension. Finally we note that by setting return_err = 1 we can get the largest truncated singular values, it is also possible to obtain all truncated value by passing any int >1.

Note

The Svd() and Svd_truncate() calls xgesdd routines from LAPACKE internally, for other algorithms like xgesvd routines one may consider Gesvd() and Gesvd_truncate() functions.

9.2. Eigenvalue decomposition

  • In Python:

 1# Create a rank-2 UniTensor which represents a square matrix
 2uT = cytnx.UniTensor.arange(4*4).reshape(4,4) \
 3                    .set_rowrank(1) \
 4                    .relabel(["in","out"]) \
 5                    .set_name("initial tensor")
 6# Eigenvalue decomposition
 7uD, uV = cytnx.linalg.Eig(uT)
 8# Create uV, uD, uV_inv = Inv(V) with names and labels
 9uV.relabel_(["in","a"]).set_name("eigenvectors")
10uD.relabel_(["a","b"]).set_name("eigenvalues")
11uV_inv = cytnx.linalg.InvM(uV).relabel_(["b","out"]) \
12                     .set_name("inverted eigenvectors")
13# Compare uT with Uv * uD * uV_inv
14uT_new = cytnx.Contract([uV,uD,uV_inv]) \
15              .permute(uT.labels()) \
16              .set_name("reconstruction from eigenvalues")
17diff = uT_new - uT
18print(diff.Norm()) # 1.53421e-14

9.2.1. Truncated Hermitian eigenvalue decomposition

For a Hermitian matrix, Eigh() returns the eigenvalues in ascending order; on a UniTensor with symmetries it does so block by block, one block per charge sector. Eigh_truncate() is the counterpart of Svd_truncate() for this case: it keeps the keepdim eigenpairs of largest magnitude \(|\lambda|\), ranked across all blocks, and returns the truncated eigenvalue UniTensor, the eigenvectors with a truncated eigen bond, and, with return_err, the discarded magnitudes. The arguments keepdim, err, mindim, min_blockdim and return_err have the same meaning as in Svd_truncate(). A typical use is the truncation of a reduced density matrix to its dominant eigenbasis, sector by sector, in one call.

  • In Python:

 1# A Hermitian matrix with Z2 symmetry: one real symmetric block per charge sector,
 2# with eigenvalues of both signs.
 3bond = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0) >> 3, cytnx.Qs(1) >> 2], [cytnx.Symmetry.Zn(2)])
 4M = cytnx.UniTensor([bond, bond.redirect()], labels=["r", "c"], rowrank=1)
 5M.put_block_(cytnx.from_numpy(np.array([[5., 0., 0.], [0., -3., 0.], [0., 0., 0.5]])), 0)
 6M.put_block_(cytnx.from_numpy(np.array([[4., 0.], [0., 1.]])), 1)
 7# The full spectrum, block by block (ascending in each block)
 8D, V = cytnx.linalg.Eigh(M)
 9print([blk.numpy() for blk in D.get_blocks()])
10# Keep the 3 eigenpairs of largest |eigenvalue| over both sectors: 5, 4, -3
11D3, V3, err = cytnx.linalg.Eigh_truncate(M, keepdim=3, return_err=1)
12print([blk.numpy() for blk in D3.get_blocks()])
13print(D3.bonds()[0].getDegeneracies())  # kept dimension per sector
14print(err.item())  # the largest discarded |eigenvalue|: 1

Output >>

[array([-3. ,  0.5,  5. ]), array([1., 4.])]
[array([-3.,  5.]), array([4.])]
[2, 1]
1.0

The ranking is global: eigenvalue magnitudes are compared across all sectors. Without per-block floors that change this selection, the block-sparse and dense paths keep the same eigenvalue magnitudes. When the retained count is exactly keepdim, the reconstruction \(V D V^\dagger\) is an optimal rank-keepdim approximation in spectral norm (the Hermitian case of the Eckart-Young theorem). The err cutoff can reduce the retained count; per-block floors or an exact magnitude degeneracy at the cut can increase it. Floors constrain the selection and can exclude the unconstrained dense optimum: sectors with spectra [-3, 0.5, 5] and [1, 4], keepdim = 2 and floors [0, 2] retain [1, 4] (spectral error 5), whereas the unconstrained optimum retains [5, 4] (error 3). return_err reports the magnitudes actually discarded, with mode 1 reporting their maximum (the spectral-norm reconstruction error).

For a dense UniTensor, the one-element min_blockdim acts as mindim and prevents cutoff-based removal; it does not override keepdim. For block inputs, floors are reserved before applying the global budget. Consequently keepdim = 1, min_blockdim = [3] retains one eigenpair for a dense input but three for a one-sector block input with at least three eigenpairs. Dense eigenvectors also inherit Eigh()’s flattened row space and rowrank 1; block eigenvectors preserve the row-space bonds and rowrank. For example a dense shape [1, 3, 1, 3] with rowrank 2 returns eigenvectors of shape [3, 2] with rowrank 1 when two eigenpairs are kept.

Dtypes follow Eigh(). For fermionic inputs, materializing pending sign flips can promote a real Float sector to Double; truncation inherits the dtype of the selected sectors.

Note

Eigh_truncate() selects by magnitude, not by algebraic value. To select eigenpairs by value, for instance the lowest or the highest eigenvalues of a Hamiltonian, call Eigh() and keep the wanted part of its ascending spectrum (for a symmetric UniTensor, block by block), or use Lanczos() with which = “SA” (smallest algebraic) or “LA” (largest algebraic), which computes only the wanted eigenpairs.

9.3. QR decomposition

The QR decomposition decomposes a matrix M to the form M = QR, where Q is an orthogonal matrix (Q Q^T = I), and R is a upper-right triangular matrix. One can perform a QR decomposition by using Qr().

Qr(Tin, is_tau)
Parameters:

  • Tin (cytnx.UniTensor) – input tensor

  • is_tau (bool) – If is_tau=True, the function returns an additional one-dimensional tensor tau that contains the scaling factors of the Householder reflectors that generate Q along with R. See [ABB+99] for details. Default: is_tau=False

Here is an example of a QR decomposition:

  • In Python:

 1uT = cytnx.UniTensor.ones([5,5,5,5,5], rowrank = 3, name="uT")
 2Q, R = cytnx.linalg.Qr(uT)
 3Q.set_name("Q")
 4R.set_name("R")
 5
 6Q.print_diagram()
 7R.print_diagram()
 8
 9# Verify the recomposition
10print((cytnx.Contract(Q,R)-uT).Norm())

Output >>

-----------------------
tensor Name : Q
tensor Rank : 4
block_form  : False
is_diag     : False
on device   : cytnx device: CPU
          ----------
         /          \
   0 ____| 5     25 |____ _aux_
         |          |
   1 ____| 5        |
         |          |
   2 ____| 5        |
         \          /
          ----------
-----------------------
tensor Name : R
tensor Rank : 3
block_form  : False
is_diag     : False
on device   : cytnx device: CPU
              ----------
             /          \
   _aux_ ____| 25     5 |____ 3
             |          |
             |        5 |____ 4
             \          /
              ----------

Total elem: 1
type  : Double (Float64)
cytnx device: CPU
Shape : (1)
[3.64214e-14 ]
[ABB+99]

E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users' Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999. ISBN 0-89871-447-8 (paperback). https://www.netlib.org/lapack/lug/node69.html.