How to initialize a tensor with permutational symmetry?

General Usage Performance
,
1

I am working on a iterative solver for non linear Coupled Cluster equations that would require a computation of lot of tensor elements via tensor contractions. For now, I am using the TensorOperations.jl for this purpose. I am trying to initialize a tensor g where the values of g[I,J,K,L] are given for only certain [I,J,K,L]. And the rest are supposed to be deduced from permutational symmetries. For now, I am using this make-shift code to achieve the above, where the values g[I,J,K,L] had been read from some data file and only the non redundant elements were provided.

for (I,J,K,L) in non_redundant_indices
        g[K,L,I,J]=g[I,J,K,L]
        g[J,I,L,K]=g[I,J,K,L]
        g[L,K,J,I]=g[I,J,K,L]
        g[J,I,K,L]=g[I,J,K,L]
        g[L,K,I,J]=g[I,J,K,L]
        g[I,J,L,K]=g[I,J,K,L]
        g[K,L,J,I]=g[I,J,K,L]
end

Does anyone have an about how may I do the same but maybe a bit more conveniently or elegantly perhaps? I looked a bit into TensorKit.jl and they seem to deal with Tensors as a new data type called TensorMap which does support some symmetry features, but I cannot figure out how to implement this permutational symmetry in those tensors and how to exploit them during tensor operations like contractions.

2

Perhaps a good solution would be to define a new Array type which has its own getindex doing the necessary transformation to arrive at a non redundant index by canonizing a general index.

This would be a commonly occuring problem, so there might be a way to use an existing package for this. I’ve managed to dig up one such package GitHub - JuliaArrays/GetindexArrays.jl: Lazy arrays with arbitrary user-defined transformations .

As an example, a way to use this to represent a triangular 3x3x3 tensor (a simple version of your symmetry requirements):

using Pkg
Pkg.add(url="https://github.com/JuliaArrays/GetindexArrays.jl/")

using GetindexArrays, Random
B = reshape(shuffle(1:9),(3,3))  # some base matrix
M = GetindexArray((1:3,1:3,1:3)) do r,c,d
    B[minimum((r,c,d)),maximum((r,c,d))]
end

and now, M is a bigger tensor, using values from B (above diagonal) as non redundant elements:

julia> B
3×3 Matrix{Int64}:
 4  2  1
 9  3  6
 5  8  7

julia> M
3×3×3 GetindexArray{Any, 3, var"#13#14", Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}} with indices 1:3×1:3×1:3:
[:, :, 1] =
 4  2  1
 2  2  1
 1  1  1

[:, :, 2] =
 2  2  1
 2  3  6
 1  6  6

[:, :, 3] =
 1  1  1
 1  6  6
 1  6  7
1 Like