7

Coming from this question, I wonder if a more generalized einsum was possible. Let us assume, I had the problem

using PyCall
@pyimport numpy as np

a = rand(10,10,10)
b = rand(10,10)
c = rand(10,10,10)

Q = np.einsum("imk,ml,lkj->ij", a,b,c)

Or something similar, how were I to solve this problem without looping through the sums?

with best regards

Community
  • 1
  • 1
varantir
  • 6,624
  • 6
  • 36
  • 57
  • 1
    Perhaps [this tensor operations library](https://github.com/Jutho/TensorOperations.jl#index-notation) could help. Do mind that it looks abandoned (latest commit March 2015 and currently failing build) – Felipe Lema Aug 19 '15 at 12:50
  • You left off the `c` from the equation. The `python` implementation depends heavily on the `nditer` iterator. – hpaulj Aug 19 '15 at 14:55
  • 1
    Just a note that `einsum` probably hasn't been directly ported to Julia simply because it's just as fast to write the loops out yourself. If I were to write a port of it, I'd probably do so as a macro to decode the subscript string at parse time and directly expand to a bunch of for loops (similar to how `@printf` works). – mbauman Aug 20 '15 at 13:58

1 Answers1

6

Edit/Update: This is now a registered package, so you can Pkg.add("Einsum") and you should be good to go (see the example below to get started).

Original Answer: I just created some very preliminary code to do this. It follows exactly what Matt B. described in his comment. Hope it helps, let me know if there are problems with it.

https://github.com/ahwillia/Einsum.jl

This is how you would implement your example:

using Einsum

a = rand(10,10,10)
b = rand(10,10)
c = rand(10,10,10)
Q = zeros(10,10)

@einsum Q[i,j] = a[i,m,k]*b[m,l]*c[l,k,j]

Under the hood the macro builds the following series of nested for loops and inserts them into your code before compile time. (Note this is not the exact code inserted, it also checks to make sure the dimensions of the inputs agree, using macroexpand to see the full code):

for j = 1:size(Q,2)
    for i = 1:size(Q,1)
        s = 0
        for l = 1:size(b,2)
            for k = 1:size(a,3)
                for m = 1:size(a,2)
                    s += a[i,m,k] * b[m,l] * c[l,k,j]
                end
            end
        end
        Q[i,j] = s
    end
end
digbyterrell
  • 3,449
  • 2
  • 24
  • 24
  • @varantir if this answered your question, please mark it as answered. Thanks! – Chris Rackauckas Aug 10 '16 at 05:53
  • 1
    This seems to implement the summation as a macro, not a function. I'm new to Julia; How would one use this to create a new array with Einstein summation syntax inside an expression. E.g. how would one use this to express what one might write as `z = exp(einsum("j,ijl->il",x,y))` in numpy? – MRule Aug 28 '20 at 15:02