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.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.