How to symmetrize a tensor by adding a tensor to itself after permuting indices ?

[Uranium92238]

I am trying to achieve the ITensor equivalent of the following

@tensor R2[a,b,i,j] := I_vvoo[a,b,i,j] + I_vvoo[b,a,j,i]

in the above snippet, I am using TensorOperations.jl and all the entities are Array{Float64,4}. I am trying to achieve the same when I_vvoo is an ITensor with indices a_1,a_2,i_1,i_2. I came across the permute function which I though was what I needed and wrote

R = I_vvoo + permute(I_vvoo,a_2,a_1,i_2,i_1)

However, these two yeild different results (the later former being what the result should be). A minimal example of the above is the follows

I_vvoo = rand(Float64,2,2,4,4)
begin
    a_1 = Index(2,"a_1")
    a_2 = Index(2,"a_2")
    i_1 = Index(4,"i_1")
    i_2 = Index(4,"i_2")
end
iI_vvoo = ITensor(I_vvoo,a_1,a_2,i_1,i_2)
iR = iI_vvoo + permute(iI_vvoo,a_2,a_1,i_2,i_1)
@tensor R[a,b,i,j] := I_vvoo[a,b,i,j] + I_vvoo[b,a,j,i]


Array(iR,a_1,a_2,i_1,i_2) - R

If everything went as expected, the last expression should have resulted in a zero tensor. But it does not.

[ryanlevy]

I think the nicer way to do this would be with swapinds potentially:

julia> a,b,i,j = Index.((2,2,2,2))
((dim=2|id=288), (dim=2|id=169), (dim=2|id=62), (dim=2|id=523))

julia> A = random_itensor(a,b,i,j);

julia> A[a=>1,b=>1,i=>1,j=>2] - A[a=>1,b=>1,i=>2,j=>1]
-0.8067652838024874

julia> B = (A + swapinds(A,i,j));

julia> B[a=>1,b=>1,i=>1,j=>2] - B[a=>1,b=>1,i=>2,j=>1]
0.0

Permute will just change the memory layout, it won’t change the indices of the tensors (which are “fixed”)

julia> Matrix(C,i,i')
2×2 Matrix{Float64}:
 -0.747951  -1.69175
  0.77861   -0.445552

julia> Matrix(permute(C,i',i),i,i')
2×2 Matrix{Float64}:
 -0.747951  -1.69175
  0.77861   -0.445552

julia> Matrix(swapinds(C,i,i'),i,i')
2×2 Matrix{Float64}:
 -0.747951   0.77861
 -1.69175   -0.445552

[mtfishman]

Also note you should use random_itensor rather than randomITensor, that’s the new notation.

[Uranium92238]

@ryanlevy

I think the nicer way to do this would be with swapinds potentially:

But this only seems to work for a single swap. Doing

iR = iI_vvoo + swapinds(iI_vvoo,a_2,a_1,i_2,i_1)

throws the error

ERROR: MethodError: no method matching swapinds(::NTuple{4, Index{Int64}}, ::Index{Int64}, ::Index{Int64}, ::Index{Int64}, ::Index{Int64})

Closest candidates are:
  swapinds(::ITensor, ::Any...; kwargs...)
   @ ITensors ~/.julia/packages/ITensors/2Xz5q/src/itensor.jl:1211
  swapinds(::ITensors.NDTensors.Tensor, ::Any...; kwargs...)
   @ ITensors ~/.julia/packages/ITensors/2Xz5q/src/itensor.jl:1214
  swapinds(::Function, ::ITensors.NDTensors.Tensor, ::Any...)
   @ ITensors ~/.julia/packages/ITensors/2Xz5q/src/itensor.jl:1203
  ...

Stacktrace:
 [1] swapinds(::ITensors.NDTensors.DenseTensor{…}, ::Index{…}, ::Vararg{…}; kwargs::@Kwargs{})
   @ ITensors ~/.julia/packages/ITensors/2Xz5q/src/itensor.jl:1215
 [2] swapinds(::ITensors.NDTensors.DenseTensor{…}, ::Index{…}, ::Index{…}, ::Index{…}, ::Vararg{…})
   @ ITensors ~/.julia/packages/ITensors/2Xz5q/src/itensor.jl:1214
 [3] swapinds(::ITensor, ::Index{Int64}, ::Vararg{Index{Int64}}; kwargs::@Kwargs{})
   @ ITensors ~/.julia/packages/ITensors/2Xz5q/src/itensor.jl:1211
 [4] top-level scope
   @ ~/Desktop/Andreas/phase_2/julia-code/ccd/itensors/t5.jl:56
Some type information was truncated. Use `show(err)` to see complete types.

[ryanlevy]

You can nest the swaps, passing four things it would ambiguous to tell which swapped with what.

iI_vvoo + swapinds(swapinds(iI_vvoo,i_2,i_1),a_2,a_1)

Or more nicely

iI_vvoo + swapinds(iI_vvoo,(i_1,a_1),(i_2,a_2))

[mtfishman]

You could also use replaceinds:

julia> using ITensors: Index, random_itensor, replaceinds

julia> a, b, i, j = Index.((2, 2, 2, 2))
((dim=2|id=54), (dim=2|id=907), (dim=2|id=731), (dim=2|id=593))

julia> t = random_itensor(a, b, i, j)
ITensor ord=4 (dim=2|id=54) (dim=2|id=907) (dim=2|id=731) (dim=2|id=593)
NDTensors.Dense{Float64, Vector{Float64}}

julia> s = replaceinds(t, a => b, b => a, i => j, j => i)
ITensor ord=4 (dim=2|id=907) (dim=2|id=54) (dim=2|id=593) (dim=2|id=731)
NDTensors.Dense{Float64, Vector{Float64}}

julia> t[a => 1, b => 2, i => 1, j => 2]
1.518138597973099

julia> s[b => 1, a => 2, j => 1, i => 2]
1.518138597973099