[ITensors] [BUG] Issue with transpose

[kmp5VT]

Description of bug

A transposed ITensor does not return the transposed data when converted

Minimal runnable code

julia> m = random_itensor(Index.((5,6)))
julia> transpose(array(m)) == array(transpose(m))
false

[mtfishman]

transpose of ITensors is a bit subtle and works differently from tranpose of Arrays. Right now it is defined as swapping primed and unprimed indices (though that can be customized with keyword arguments). Additionally, it works on n-dimensional arrays, not just matrices and vectors, where the interpretation is that we are transposing a linear map between two sets of indices. So it is better to think about it as a generalization of matrix transpose to n-dimensional arrays and where the transposition is handled lazily by renaming the indices.

Maybe this demonstration helps:

julia> i = Index(2)
(dim=2|id=609)

julia> a = random_itensor(i', i)
ITensor ord=2 (dim=2|id=609)' (dim=2|id=609)
NDTensors.Dense{Float64, Vector{Float64}}

julia> transpose(a)
ITensor ord=2 (dim=2|id=609) (dim=2|id=609)'
NDTensors.Dense{Float64, Vector{Float64}}

julia> transpose(array(a))
2×2 transpose(::Matrix{Float64}) with eltype Float64:
 0.925524  -0.10343
 0.985219  -0.503981

julia> array(transpose(a))
2×2 Matrix{Float64}:
  0.925524   0.985219
 -0.10343   -0.503981

julia> a + transpose(a)
ITensor ord=2 (dim=2|id=609)' (dim=2|id=609)
NDTensors.Dense{Float64, Vector{Float64}}

julia> array(a + transpose(a))
2×2 Matrix{Float64}:
 1.85105    0.881789
 0.881789  -1.00796

julia> j = Index(2)
(dim=2|id=115)

julia> a = random_itensor(i', j', i, j)
ITensor ord=4 (dim=2|id=609)' (dim=2|id=115)' (dim=2|id=609) (dim=2|id=115)
NDTensors.Dense{Float64, Vector{Float64}}

julia> array(a)
2×2×2×2 Array{Float64, 4}:
[:, :, 1, 1] =
 0.200765   -0.902874
 0.0431985  -0.844064

[:, :, 2, 1] =
  1.09159   -0.0412995
 -0.699256   0.766109

[:, :, 1, 2] =
 1.43898  0.249222
 2.07979  1.4075

[:, :, 2, 2] =
 -0.740071  -1.01018
  0.349539  -0.202951

julia> transpose(a)
ITensor ord=4 (dim=2|id=609) (dim=2|id=115) (dim=2|id=609)' (dim=2|id=115)'
NDTensors.Dense{Float64, Vector{Float64}}

julia> array(transpose(a))
2×2×2×2 Array{Float64, 4}:
[:, :, 1, 1] =
 0.200765   -0.902874
 0.0431985  -0.844064

[:, :, 2, 1] =
  1.09159   -0.0412995
 -0.699256   0.766109

[:, :, 1, 2] =
 1.43898  0.249222
 2.07979  1.4075

[:, :, 2, 2] =
 -0.740071  -1.01018
  0.349539  -0.202951

julia> a + transpose(a)
ITensor ord=4 (dim=2|id=609)' (dim=2|id=115)' (dim=2|id=609) (dim=2|id=115)
NDTensors.Dense{Float64, Vector{Float64}}

julia> array(a + transpose(a))
2×2×2×2 Array{Float64, 4}:
[:, :, 1, 1] =
 0.401531   0.536106
 1.13479   -1.58414

[:, :, 2, 1] =
  1.13479  2.0385
 -1.39851  1.11565

[:, :, 1, 2] =
 0.536106  0.498443
 2.0385    0.397318

[:, :, 2, 2] =
 -1.58414   0.397318
  1.11565  -0.405902

Note that with the NamedDimsArrays rewrite we will probably change the interface a little bit through more explicit use of an ITensorOperator type, which would represent an ITensor that is a linear map between two sets of its indices, and transposing would simply be swapping the domain and codomain indices.