A question regarding Matrix and ITensor gates

sourav · November 9, 2023, 12:24am

Dear ITensor Team,

The question may be somewhat long --please bear with me–I appologise in adv.

Before I go to why I am asking this question, let me start with an example. Suppose I have the Hamiltonian H=\sum_{i}X_{i}X_{i+1} + Y_{i}Y_{i+1} + Z_{i}Z_{i+1} where X_i=S_x^{i} and so on for S=1/2 particles.

Now, for time evolution (say via TEBD) I consider a two site gate. To this end, I simply use the following code (say for the first bond):

using ITensors
N = 10
s = siteinds("S=1/2", N)
s1, s2 = s[1], s[2]

hj_1 = 1.0*op("Sx", s1)*op("Sx",s2) + 1.0*op("Sy", s1)*op("Sy",s2) + 1.0*op("Sz", s1)*op("Sz", s2)
print(hj_1)

Now, if I print hj_1 in the above code snippet I simply get

ITensor ord=4
Dim 1: (dim=2|id=826|"S=1/2,Site,n=1")'
Dim 2: (dim=2|id=826|"S=1/2,Site,n=1")
Dim 3: (dim=2|id=246|"S=1/2,Site,n=2")'
Dim 4: (dim=2|id=246|"S=1/2,Site,n=2")
NDTensors.Dense{ComplexF64, Vector{ComplexF64}}
 2×2×2×2
[:, :, 1, 1] =
 0.25 + 0.0im    0.0 + 0.0im
  0.0 + 0.0im  -0.25 + 0.0im

[:, :, 2, 1] =
 0.0 + 0.0im  0.5 + 0.0im
 0.0 + 0.0im  0.0 + 0.0im

[:, :, 1, 2] =
 0.0 + 0.0im  0.0 + 0.0im
 0.5 + 0.0im  0.0 + 0.0im

[:, :, 2, 2] =
 -0.25 + 0.0im   0.0 + 0.0im
   0.0 + 0.0im  0.25 + 0.0im

Now, suppose for some reason (will clarify later), I do not want to use the op(“Sz”, site)*op(“Sz”, site+1) type stuffs. Since at each bond the gate operator is X_{i} \otimes X_{i+1} + Y_{i} \otimes Y_{i+1} + Z_{i} \otimes Z_{i+1} I could have simply done the following

using ITensors
import ITensors: op
N=10

const sz = (1/2.0)*[1. 0.;0. -1.]
const sx = (1/2.0)*[0. 1.; 1. 0.]
const sy = (1/2.0)*[0. -1.0im; +1.0im 0.]


s = siteinds("S=1/2", N)
function op(::OpName"test",::SiteType"S=1/2",s1::Index,s2::Index)
    
    h = kron(sz, sz) + kron(sx, sx) + kron(sy, sy)
    return itensor(Matrix(h),s1',s1,s2',s2)
end

print(op("test", s, 1, 2))

Interestingly, now I get a result of the form

ITensor ord=4
Dim 1: (dim=2|id=51|"S=1/2,Site,n=1")'
Dim 2: (dim=2|id=51|"S=1/2,Site,n=1")
Dim 3: (dim=2|id=183|"S=1/2,Site,n=2")'
Dim 4: (dim=2|id=183|"S=1/2,Site,n=2")
NDTensors.Dense{ComplexF64, Vector{ComplexF64}}
 2×2×2×2
[:, :, 1, 1] =
 0.25 + 0.0im  0.0 + 0.0im
  0.0 + 0.0im  0.0 + 0.0im

[:, :, 2, 1] =
   0.0 - 0.0im  0.5 - 0.0im
 -0.25 + 0.0im  0.0 + 0.0im

[:, :, 1, 2] =
 0.0 - 0.0im  -0.25 + 0.0im
 0.5 - 0.0im    0.0 + 0.0im

[:, :, 2, 2] =
 0.0 + 0.0im   0.0 - 0.0im
 0.0 - 0.0im  0.25 + 0.0im

Clearly the matrix structures in the two cases are different although my naive impression was that I am doing the same thing. Could you please clarify what is going on?

Just to clarify why I’m asking this question: I have a scenario where I am Trotterizing a Liouvillean and for boundary bonds I have the matrix for dissipative part of the Liouvillean which I want to convert in Itensor operator and add to the Hamiltonian part of the Liouvillean which can be written in terms of product of conventional S=1/2 operator available in ITensor.

Best,
Sourav.

mtfishman · November 9, 2023, 1:41pm

It should be itensor(h,s1',s2',s1,s2) or itensor(h,s2',s1',s2,s1) (Matrix isn’t needed, h is already a Matrix).

You can think of the ITensor constructor as doing a reshape of the matrix you are inputting, so you have to make sure it is reshaping properly to match the Index ordering you input, or manually reshape and permute the input to match the corresponding Index ordering you want. Here you also have to think carefully about the memory ordering of the matrix output by kron, the h you make is going to be a matrix from the bra space to the ket space, which is why you need to use an Index pattern like s1',s2',s1,s2.

sourav · November 9, 2023, 3:19pm

Hi Matt,

Thanks for your reply. I used itensor(s1', s2' , s1, s2) type indexing as you suggested. But that does not fully resolve the problem. To clarify :

N = 10
s = siteinds("S=1/2", N)
s1, s2 = s[1], s[2]
hj_1 = 1.0*op("Sz", s1)*op("Sz",s2) 
print(hj_1)

yields an output

ITensor ord=4
Dim 1: (dim=2|id=243|"S=1/2,Site,n=1")'
Dim 2: (dim=2|id=243|"S=1/2,Site,n=1")
Dim 3: (dim=2|id=112|"S=1/2,Site,n=2")'
Dim 4: (dim=2|id=112|"S=1/2,Site,n=2")
NDTensors.Dense{Float64, Vector{Float64}}
 2×2×2×2
[:, :, 1, 1] =
 0.25   0.0
 0.0   -0.25

[:, :, 2, 1] =
 0.0  0.0
 0.0  0.0

[:, :, 1, 2] =
 0.0  0.0
 0.0  0.0

[:, :, 2, 2] =
 -0.25  0.0
  0.0   0.25

On the other hand the code

import ITensors: op
function op(::OpName"test",::SiteType"S=1/2",s1::Index,s2::Index)
   h = kron(sz, sz)
   return itensor(h,s1',s2',s1,s2)
end
print(op("test", s, 1, 2))

yields

ITensor ord=4
Dim 1: (dim=2|id=243|"S=1/2,Site,n=1")'
Dim 2: (dim=2|id=112|"S=1/2,Site,n=2")'
Dim 3: (dim=2|id=243|"S=1/2,Site,n=1")
Dim 4: (dim=2|id=112|"S=1/2,Site,n=2")
NDTensors.Dense{Float64, Vector{Float64}}
 2×2×2×2
[:, :, 1, 1] =
 0.25  0.0
 0.0   0.0

[:, :, 2, 1] =
  0.0    0.0
 -0.25  -0.0

[:, :, 1, 2] =
 0.0  -0.25
 0.0  -0.0

[:, :, 2, 2] =
  0.0  -0.0
 -0.0   0.25

And again the matrices are not the same. However, changing one line in the previous code involving permutationdims and reshaping seems to resolve the issue as follows:

import ITensors: op
function op(::OpName"test",::SiteType"S=1/2",s1::Index,s2::Index)
    h = kron(sz, sz)
    h = permutedims(reshape(h, 2, 2, 2, 2), [1, 3, 2, 4])
    return itensor(h,s1',s2',s1,s2)
end
print(op("test", s, 1, 2))

It gives the output same (in terms of matrix structures) as using op("Sz",s1)*op("Sz",s2). I am not sure why one needs to reshuffle the dims in the [1,3,2,4] pattern. Is there some physical reasoning ? Would be nice if you kindly clarify.

Best,
Sourav

mtfishman · November 9, 2023, 3:48pm

The ITensor interface is agnostic about the internal memory ordering/Index ordering, so from the perspective of ITensor operations those tensors are the same. For example they compare as equal:

using ITensors
N = 10
s = siteinds("S=1/2", N)
s1, s2 = s[1], s[2]

const sz = (1/2.0)*[1. 0.;0. -1.]
const sx = (1/2.0)*[0. 1.; 1. 0.]
const sy = (1/2.0)*[0. -1.0im; +1.0im 0.]

function ITensors.op(::OpName"test",::SiteType"S=1/2",s1::Index,s2::Index)
    h = kron(sz, sz) + kron(sx, sx) + kron(sy, sy)
    return itensor(h,s1',s2',s1,s2)
end

h1 = 1.0*op("Sx", s1)*op("Sx",s2) + 1.0*op("Sy", s1)*op("Sy",s2) + 1.0*op("Sz", s1)*op("Sz", s2)
h2 = op("test", s, 1, 2)
h1 == h2

# If you want them to print the same
@show permute(h1, s1', s2', s1, s2)
@show permute(h2, s1', s2', s1, s2)

sourav · November 9, 2023, 4:07pm

Hi Matt,

Got it. Thanks a lot

Best,
Sourav.

sourav · November 14, 2023, 2:17pm

Hi @mtfishman ,

I had the impression that I had understood the previous answer but it seems that I possibly got something wrong. Let me point it out via an example.

I have defined a new Sitetype with name “mpdo=1/2”. The code and some simple operator goes as follows:

using ITensors
const sz = (1/2.0)*[1. 0.;0. -1.]
const sx = (1/2.0)*[0. 1.; 1. 0.]
const sy = (1/2.0)*[0. -1.0im; +1.0im 0.]
const sp = [0. 1. ; 0. 0.]
const sm = [0. 0. ; 1. 0.]
const id_half = [1. 0.; 0. 1.]

ITensors.space(::SiteType"mpdo=1/2") = 4

function ITensors.op!(Op::ITensor,O::OpName"Sz",::SiteType"mpdo=1/2",s::Index)
    A=kron(sz, id_half)
    for i in 1:4
        for j in 1:4
            Op[s'=>i,s=>j] = A[i,j]
        end
    end
end

function ITensors.op!(Op::ITensor,O::OpName"Sx",::SiteType"mpdo=1/2",s::Index)
    A=kron(sx, id_half)
    for i in 1:4
        for j in 1:4
            Op[s'=>i,s=>j] = A[i,j]
        end
    end
end

function ITensors.op!(Op::ITensor,O::OpName"Sy",::SiteType"mpdo=1/2",s::Index)
    A=kron(sy, id_half)
    for i in 1:4
        for j in 1:4
            Op[s'=>i,s=>j] = A[i,j]
        end
    end
end

As, you see here S_x(mpdo=1/2) = S_{x}(spin=1/2) \otimes I and so on. Now suppose I construct the an operator A by the following method.

L = 6
s = siteinds("mpdo=1/2",L)
s1, s2 = s[1], s[2]
A = op("Sz",s1)*op("Sx",s2)

Alternatively, I could have also constructed (using const sz, sx and id_half) the same operator (call it B this time) as follows:

import ITensors: op
function op(::OpName"Test",::SiteType"mpdo=1/2",s1::Index,s2::Index)
    ss1 = kron(kron(sz, id_half), kron(sx,id_half))
    return itensor(ss1,s1',s2',s1,s2)
end
B = op("Test",s,1,2)

Now, I see that A==B gives me false. However, just by some trick if I do a index permutation during reshaping

import ITensors: op
function op(::OpName"Test",::SiteType"mpdo=1/2",s1::Index,s2::Index)
    ss1 = kron(kron(sz, id_half), kron(sx,id_half))
    ss1 = permutedims(reshape(ss1, 4, 4, 4, 4), [4, 3, 2, 1])
    return itensor(ss1,s1',s2',s1,s2)
end
B = op("Test",s,1,2)

I get A==B true—not sure what’s going on. I would have thought anything without reshuffling is already correct. Could you please point out what is going wrong ?

Best,
Sourav

Topic		Replies	Views
Mixed nearest neighbour and next nearest neighbour TEBD ITensor Julia Questions	6	356	January 21, 2024
Time-evolution of electron system and the difference Julia and C++ ITensor Julia Questions julia , mps	2	363	December 17, 2022
Adding ITensors objects with different index sites ITensor Julia Questions julia , gate_evolution	12	659	October 15, 2022
How to correctly do TEBD with ITensor? ITensor Julia Questions julia , mps	1	734	March 14, 2023
Building a Hubbard Hamiltonian for TEBD DMRG and Numerical Methods julia	7	163	June 20, 2024

A question regarding Matrix and ITensor gates

Related topics