tensor contractions etc

[Jutho]

Dear Julia developers,

I was thinking of writing some code, or maybe even a package, that can compute generic tensor contractions, traces etc.

A contraction of two tensors would be any operation where you start with arrays A and B with NA and NB dimensions, and

some indices of A will be contracted with equally sized indices in B (like the inner index in matrix multiplication) and others will

be present in the output.

I am still thinking of the best syntax, but something like

C["i,j,k"]=A["k,m,n,i"]B["n,j,m"]

seems to be very popular in similar C++ style libraries.

A naive approach is to use permutedims to reshuffle the dimensions of A and B such that the open indices can be grouped together and the contraction indices can be grouped together, then reshape to matrices, use matrix multiplication, reshape to get the original dimension back, and do a final permutedims to get the dimensions as requested in the output tensor. This requires in general three permutedims and thus three internal objects that require memory allocation and deallocation, so I would like to avoid that by writing a better approach that does not require temporise and is cache friendly,  as inspired by the fallback method general_matmatmul in Base.LinAlg.  I would of course also like to reuse as much of the code in Base as possible, and not have to rewrite almost identical code.

In that respect, I was wondering, however, why generic_matmatmul, which is a fallback for Base.LinAlg.BLAS.gemm!, was chosen not to exactly mimic gemm! so as to allow also to add the result of A*B to C (with some coefficients) rather than just overwriting C. It would be a near-zero cost operation to modify this function. If it is then also called gemm! instead of generic_matmatmul, it would be completely transparant.

Similarly, it would be nice to have a permutedims! which adds rather then overwrites (like a combination of permutedims! and axpy!).

I can of course modify this myself, but want to check first the original motivation behind the current approach.

Best regards,

Jutho

[Miles Gould]

On 21/03/14 20:17, Jutho wrote:
> I am still thinking of the best syntax, but something like
> C["i,j,k"]=A["k,m,n,i"]B["n,j,m"]
> seems to be very popular in similar C++ style libraries.

It's been a while since I worked with tensors, but would it be useful to
distinguish between co- and contravariant indices as in Einstein
summation notation?

Miles

[Tim Holy]

On Friday, March 21, 2014 01:17:11 PM Jutho wrote:
> In that respect, I was wondering, however, why generic_matmatmul, which is
> a fallback for Base.LinAlg.BLAS.gemm!, was chosen not to exactly mimic
> gemm! so as to allow also to add the result of A*B to C (with some
> coefficients) rather than just overwriting C. It would be a near-zero cost
> operation to modify this function. If it is then also called gemm! instead
> of generic_matmatmul, it would be completely transparant.

Order in which various things were written. gemm! was "recently" modified to
mimic axpy!.

[Jiahao Chen]

Two other serious projects that come to mind (albeit very
young) are the Cyclops tensor framework (http://ctf.eecs.berkeley.edu)
and BTAS (http://itensor.org/btas).

I'm not necessarily suggesting to wrap these libraries, but rather to
examine how they are implemented.

One area I have been concerned about is the efficient representation
and manipulation of tensors that are symmetric or antisymmetric in
certain indices. In physics and chemistry applications, such
symmetries are prevalent, and in higher dimensions, the redundancy
factor in storing such objects grows exponentially in the number of
(anti)symmetric indices.

The co- vs contravariant index issue I suspect will turn out to be a
red herring in the end, considering that all you need to convert one
to the other is to contract with g or g^-1 in the requisite index.
Thanks,

Jiahao Chen
Staff Research Scientist
MIT Computer Science and Artificial Intelligence Laboratory

[Jutho]

By coincidence, I had a meeting today with the developer of Cyclops (I am in Berkeley right now) and I know personally one of the developers of BTAS, although I had no idea he was working on this. I will contact him. We are working in a similar area of physics, where tensor contractions are used a lot (but permutation symmetry is typically not present in the tensors, unlike in coupled cluster applications in chemistry).

Op vrijdag 21 maart 2014 17:48:25 UTC-7 schreef Jiahao Chen:

[Jutho]

It would in many physical applications be meaningful to distinguish between co- and contravariant indices (relativity) or incoming and outgoing indices (other contexts), but this just requires a new type that can store the type of the indices. I was in the first place talking about just providing some additional functionality for the existing types.

Op vrijdag 21 maart 2014 14:17:15 UTC-7 schreef Miles Gould: