Support for permutational index symmetry

[Krzmbrzl]

In #50 (comment) I found

symmetries under index permutations are currently not supported. In general, it's quite hard to take them into account, in particular in a way that leads to a large speedup.

I just wanted to come back to it (given that this statement is from 2021) and ask: Has there been any work in this direction in the meantime?

FWIW I think already support in the sense of figuring out that if $f_{ab} = f_{ba}$, then it is enough to compute and store the elements for $b > a$ (which especially for higher-dimensional tensors can lead to significant memory and also processing time savings).

[Jutho]

Dear @Krzmbrzl , support for this is still lacking, as it is not really of interest to us. What kind of permutation symmetry does your use case require?

In quantum many body physics, one has symmetric or antisymmetric wave functions in first quantisation. We tend to work in second quantisation. However, if the total particle number is enforced, both descriptions should be equivalent, and have an equal number of degrees of freedom. I am wondering if something can be done with this. For example, an antisymmetric tensor of N indices where every index can range over 1:L is equivalent to a tensor with L indices, where every index only ranges over 0:1 with the added constraint that the only nonzero values of the tensor are entries T[i1,i2,...,iL] for which i1 + i2 + ... + iL = N. The latter can be represented by TensorKit. So maybe this is a way to then also support tensors with index symmetries, by writing a wrapper that converts the data in second quantisation structure to first quantisation structure.

[Krzmbrzl]

Thank you very much for your feedback @Jutho 👍

What kind of permutation symmetry does your use case require?

I'm working in the field of quantum chemistry, mostly with Coupled Cluster methods. The resulting tensors are almost always "column-symmetric", which comes about due to the particle-1,2 symmetry (i.e. the indistinguishability of electrons). For instance the tensors representing the "coupled cluster amplitudes" have the symmetry

t^{ab}_{ij} = t^{ba}_{ji}

these things can become rather large in typical applications and can therefore end up becoming a memory bottleneck in practical applications (especially once triple excitations $t^{abc}_{ijk}$ come in (which have 6 equivalent index permutations)).

In terms of where computational effort could be saved, in Coupled Cluster methods we do non-linear optimization to find a set of cluster amplitudes for which the "residual" tensor becomes zero. That means that in every iteration we compute this residual tensor by means of a series of tensor contractions. This residual (in most cases) also has the previously mentioned column-symmetry. A contribution to such a residual could be

O^{ij}_{ab} = g^{ij}_{kl} t^{kl}_{ab}

and thus by restricting the "external indices" a, b, i and j to the non-redundant part of O, we could actually save on necessary computations.

Intermediates in our expressions (appearing because we factorize expressions into a series of binary tensor contractions) can in principle show arbitrary index permutation symmetry (it depends on the exact expressions that contribute to the intermediate). However, I believe that being able to support the basic column-symmetry would be a big gain for us already.

---

if the total particle number is enforced, both descriptions should be equivalent, and have an equal number of degrees of freedom. I am wondering if something can be done with this. For example, an antisymmetric tensor of N indices where every index can range over 1:L is equivalent to a tensor with L indices, where every index only ranges over 0:1 with the added constraint that the only nonzero values of the tensor are entries T[i1,i2,...,iL] for which i1 + i2 + ... + iL = N

That's an interesting suggestion. For the kinds of column symmetry that I'd need, this approach would (to my understanding) require the introduction of composite indices though.
Assuming a,b run over 1-3 and i,j of 1-2 respectively we'd need to introduce

A = a + 3 * (i - 1)
B = b + 3 * (j - 1)

and then express the tensor $O^{ij}{ab}$ as a tensor $O{AB}$ which now has the symmetry

O_{AB} = O_{BA}

which can be represented by a tensor T[i1, i2, i3, i4, i5, i6] where all is can take values 0, 1 or 2. The sum of the indices would then be restricted to be equal to 2.

The problem that I see with this approach is that we can't use the composite index for the expressions that actually compute O as there the indices a, b, i and j may appear on different tensors that are contracted with each other (via different indices).

Finally, wouldn't there be performance issues by introducing tensors that have a huge amount of indices but where each indices runs only over two or three values? I imagine that caching could become problematic for such a data layout 🤔