Replace OMEinsum with TensorOperations

[mofeing]

After some benchmarking, I've found out that OMEinsum is terribly slow and would like to see how TensorOperations performs compared to it. This PR replaces OMEinsum with TensorOperations. If it improves runtime, I will replace our backend for it.

@lkdvos @Jutho I'm running into some edge cases that are managed by OMEinsum but not TensorOperations. Would you mind lending me a hand?

[lkdvos]

I'm obviously also curious about the benchmarks, and would gladly help out a bit. This being said, I think one of the main differences between OMEinsumeinsum and Tensoroperations implementation is that TensorOperations really requires strict Einstein summation, every index must appear exactly twice: either once in the input and once in the output, or twice in the inputs such that it is contracted. From a quick first glance, these are the errors you are having?

I am not sure if this is easily resolved, both the contraction order determination and the elementary operations make use of this assumption. For instance, I think the hyper indices would require a major rewrite as we really have that contractions happen pairwise.

I'm not familiar enough with OMEinsum, but maybe there is a way to only dispatch the implementation of certain pairwise contractions to TensorOperations?

[mofeing]

The edge cases that are currently failing are the following:

1. A[j, k] = C[i, j, k] Sum over dim (i)
2. A[i, j] = C[i, j, i] Diagonal selection (i)
3. C[i, k, l] = A[i, j, k] * B[k, l, j] Contraction over j & Hadamard product over k
4. C[i, x] = A[i, x] * B[x] Hadamard product over x

Cases (1) and (2) are managed by tensortrace! and cases (3) and (4) are managed by tensorcontract!. Case (3) can be complicated to fix, but I believe cases (1), (2) and (4) could be potentially managed by TensorOperations.

I am not sure if this is easily resolved, both the contraction order determination and the elementary operations make use of this assumption. For instance, I think the hyper indices would require a major rewrite as we really have that contractions happen pairwise.

We look for the contraction path with EinExprs and then call pairwise contractions with tensorcontract! or single contractions (i.e. sum over dim(s) or diagonal selection) with tensortrace!. If case (4) (and maybe case (3)) is fixed, then hyperindices shouldn't be a problem.

I'm not familiar enough with OMEinsum, but maybe there is a way to only dispatch the implementation of certain pairwise contractions to TensorOperations?

I'm trying to remove OMEinsum completely since I try to keep Tenet as lightweight as posible, OMEinsum has terrible performance and compile times, and loads some unneeded packages like ChainRules (when they could be loaded as package extensions).

[lkdvos]

I don't have anything against supporting that, but note that none of these cases are really true Einstein summation things, hence they don't currently work.
I also have to say that I don't have too much time myself to invest into this right now, and at least for us there is not too many reasons to add this: our ncon and @tensor implementation will not generate calls to this and error earlier. Additionally, within the context of symmetric tensors and TensorKit, none of these cases are very well-defined when considering charge conservation.
This being said, if you would find a nice way to add this, I am not opposed to supporting that for Arrays

[mofeing]

I'm thinking that we should have sth like a "EinsumInterface.jl", similar to what "DifferentiationInterface.jl" does for AD engines.

[lkdvos]

Maybe, but it's hard to cleanly define what you want to include. TensorOperations does have an interface and backend selection mechanisms, so you could plug in whatever implementation you want quite easily, but limits itself to simple einsum expressions. There is definitely a case to be made for hyperindices, summed indices, but at that point you might also want to consider fixed indices, and maybe you don't want to have just contractions with conjugations, but also allow 'abs' or 'max', and have different reduction operations etc. It's not that obvious where to draw the line, and any line seems as arbitrary as any other line, making it hard to establish a good interface...

[mofeing]

@lkdvos I've been thinking how to circunvent these issues and what can we do from Tenet's side. so for each case...:

• case 1: i can just detect and call sum on it
• case 2: i can just detect and call view or view and collect on it
• case 4: i can just detect, reshape and call broadcast(*)

for case 3, it's a bit harder but i believe what we can do is just project k, call tensorcontract! for each projection of k and stack them.

maybe you want to reconsider support for case 3, not for the explanation i gave before, but because it's equivalent to batch tensor contraction, and BLAS implementations are starting to give good support for batch GEMM (but i saw @Jutho's answers to issues asking for batch tensor contract).

[lkdvos]

I think the main issue will be that you cannot intercept the macro without fully rewriting it, so while your functions will work, the actual @tensor calls will not be able to support these edge cases. If you are okay with that, I guess that's not too bad?

[mofeing]

ahh i can't use the @tensor macro because the indices are chosen during runtime, so my contract function directly calls tensortrace! and tensorcontract!

i can detect these edge cases in contract, apply these "solutions" and call tensortrace!/tensorcontract!

[lkdvos]

That indeed sounds plausible, and probably not that hard to implement?