TensorOperations.jl

Matrix chain multiplication can be seen a special case of tensor contraction, and the https://github.com/Jutho/TensorOperations.jl package has an implementation which tries to optimize the contraction order (macro @tensoropt).

There are a bunch of packages exploring ideas around labelled dimensions, although we are still far from the vision of caring as little as we do about registers: * https://github.com/JuliaArrays/AxisArrays.jl * https://github.com/davidavdav/NamedArrays.jl * https://github.com/invenia/NamedDims.jl * https://github.com/rafaqz/DimensionalData.jl * https://github.com/ITensor/ITensors.jl While these don't pass dimension labels along in the type, there are also a number of Einsum-like packages which should work on any storage order, and could be taught to interface with such types: * https://github.com/Jutho/TensorOperations.jl * https://github.com/under-Peter/OMEinsum.jl * https://github.com/mcabbott/TensorCast.jl

Julia has some really cool libraries for index notation. I've made extensive use of TensorOperations.jl [1] and TensorCast.jl [2]. I find that the capabilities of these two complement each-other really well since TensorOperations.jl is a bit more restricted in the sort of operations it can do, but it does all sorts of really powerful analysis to find the optimal way of doing a given contraction whereas TensorCast.jl is super flexible, permissive and expressive but doesn't always generate an optimal contraction for big high dimensional contractions. I'm excited to see what comes out of Andreas' summer work on his implementation of Einsum! There's definitely room for more cool ideas in this space. In particular, I think baking in AD functionality is a cool approach. [1] https://github.com/Jutho/TensorOperations.jl [2] https://github.com/mc

The julia DSL for this is closer to what you want: "M_ij = A_ik * B_kj" reads @einsum M[i,j] := A[i,k] * B[k,j] Or faster with @tensor from https://github.com/Jutho/TensorOperations.jl & the same notation. I'm not aware of any attempt to distinguish co-vectors; you can have g = Diagonal([-1,1,1,1]) and insert A[i,k]g[k,k′]B[k′,j] and this is handled efficiently.