Extracting Singular Values of an entanglement cut within a Quantum Number Block

[amoghanakru]

Hi, I want to ask a question that is sort of a follow-up/continuation of this question: How to realize momentum Entanglement spectrum in ITensor.jl - #3 by kevinh. The goal of that question (which is my goal as well) is to compute, given an entanglement cut of a QN-conserving wavefunction, the pairs/tuples of Schmidt eigenvalues and the quantum numbers.

If I compute the SVD of a QN-conserving MPS psi like so:

psi = orthogonalize(psi, b)
  U, S, V = svd(psi[b], (linkind(psi, b - 1), siteind(psi, b)))

the object S is an ITensor whose entries (singular values) can be accessed with S[n,n] for n in 1:size(S)[1]. The eigenbasis of the singular value matrix is block diagonal in the quantum number basis, and the code above produces an S that does contain the QN information. For example, if we display S in the terminal, the output is

ITensor ord=2
Dim 1: (dim=181|id=976|"Link,u") <Out>
 1: QN(("A",1),("B",32)) => 6
 2: QN(("A",1),("B",31)) => 13
 ...
 13: QN(("A",-1),("B",-32)) => 6
Dim 2: (dim=181|id=747|"Link,v") <Out>
 1: QN(("A",-1),("B",-32)) => 6
 2: QN(("A",-1),("B",-31)) => 13
 ...
 13: QN(("A",1),("B",32)) => 6

where “A”, “B” are two quantum numbers. Doing println(S) actually prints (along with the above code) the matrix entries corresponding to singular value matrix, in Blocks corresponding to the QN sectors:

NDTensors.DiagBlockSparse{Float64, Vector{Float64}, 2}
 181×181
Block(1, 1)
 [1:6, 1:6]
 0.00028918907365869354  0.0                   0.0                   0.0                    0.0                   0.0
 0.0                     7.634975754941061e-5  0.0                   0.0                    0.0                   0.0
 0.0                     0.0                   6.501313181278113e-5  0.0                    0.0                   0.0
 0.0                     0.0                   0.0                   1.5586653810441255e-5  0.0                   0.0
 0.0                     0.0                   0.0                   0.0                    8.979528300027936e-6  0.0
 0.0                     0.0                   0.0                   0.0                    0.0                   2.2726487218800084e-6
#other blocks to follow

From the print output, I can of course recover the singular values within a QN block of my choice. However, this is messy, and I would much rather have a way to input a given QN sector (like ((“A”,1),(“B”,32))) and recover the corresponding block of the SVD matrix (which in my example is a square matrix only because the left/right bond dimensions are equal). I unfortunately couldn’t see how to do this sort of operation by digging through the ITensors.jl code (and through related questions on this forum), and I would definitely appreciate any direction.

Thanks very much,

Amogh

[mtfishman]

I think this is the best way to do that right now:

using ITensors: Block, Index, QN, commonind, dag, random_itensor, space
using LinearAlgebra: svd

i = Index([QN(0) => 2, QN(1) => 3], "i")
j = Index([QN(0) => 2, QN(1) => 3], "j")
A = random_itensor(i, dag(j))
U, S, V = svd(A, i)
u = commonind(S, U)
v = commonind(S, V)
b1 = findfirst(b -> first(b) == QN(1), space(u))
b2 = findfirst(b -> first(b) == QN(1), space(v))
S[Block(b1, b2)]