
Julia Summer of Einsum - selimthegrim
https://nextjournal.com/under-Peter/julia-summer-of-einsum
======
eigenspace
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](https://github.com/Jutho/TensorOperations.jl)

[2]
[https://github.com/mcabbott/TensorCast.jl](https://github.com/mcabbott/TensorCast.jl)

------
skybrian
It sounds neat, but I'm wondering about how people who work with matrices can
make their code readable. Defining everything in terms of array indices seems
a bit finicky and obscure? Suppose you get the index wrong and get the wrong
thing?

It seems like named columns like in SQL might be better?

~~~
eigenspace
There's also a lot of work going on in that direction right now. I'd check out
NamedDims.jl[1] and AxisArrays.jl[2]. This approach is more similar to an
EinSum type apprach than you seem to realize though.

The main difference between this sort of approach and the approach in an
EinSum esque package is that these packages actually encode the name of the
dimensions in the type signature of the array whereas an Einsum-esque
approach, the names are just temporary labels for the purposes of the macro
call.

[1]
[https://github.com/invenia/NamedDims.jl](https://github.com/invenia/NamedDims.jl)

[2]
[https://github.com/JuliaArrays/AxisArrays.jl](https://github.com/JuliaArrays/AxisArrays.jl)

~~~
oxinabox
(Dev of NamedDims here) Currently having a (mostly in my head) debate about
adding automatic permutation based on names to make operations legal if there
is permutation of names that is legel. which would result in `*` act kinda
like tensor contraction (but only for matrix matrix / matrix vectors)

Also towards adding something like, making einsum work where you just specify
the left-hand side of the einsum and then it uses the encoded names. (Rather
than the temporay labels)

PRs to add that feature are welcome

------
svnpenn
Julia syntax is beautiful, but the REPL on Windows is so slow as to be
unusable

[https://github.com/JuliaLang/julia/issues/32198](https://github.com/JuliaLang/julia/issues/32198)

~~~
mruts
Is it? The asymmetry of function definitions seems really ugly to me:

f(x) = x * x

function(x) x * x end

function f(x) x * x end

f = x -> x * x

If Julia just used curly brackets instead, we wouldn't need all of these
different ways to declare functions.

It really have been nicer if I could just do:

{ x -> x * x }

or for named function:

def f(x) { x * x }

~~~
eigenspace
Options two and four do different things than options 1 and 3 so those aren't
really relevant to the discussion.

Options one and three are both good to have because sometimes you want one-
line function definitions (mostly when you want your code to look like math)
and sometimes you want big standard 'code looking' multi-line function
definitions. If you don't want to choose in the future, just choose one now
and stuck with it. Both can be used for everything if you're willing to put up
with some awkward syntax.

f(x) = begin x * x end

If this is your main complaint about julia, I'd say we're doing pretty good.

Curly brackets are useful notation and it'd be a waste to use them up for
something like this.

~~~
mruts
It seems like the core of the issue has to do with array indices vs function
application. Julia uses curly brackets for type parameters because it can't
use brackets because they are taken for array notation. But what's wrong with
less than or greater than for type parameters? If it's followed by colons for
contravariant/covariant types, I don't see why it can't be used for both.

Julia's great, but the Pascal-like 'begin and 'end' delimiters strike me as
incredibly ugly. There are many alternatives that I can see (such as the one I
mentioned above) that would obviate this situation.

I consider myself a 'mathy' person (I work as a quant), but even so, Julia is
a very odd language. All of these compromises or whatever just to please the
math people who have a hard time programming. Some of these things the more I
think about them are great (array indices starting at one make a lot of sense,
actually) and the ability to just write '2x' instead of '2*x' is growing on
me.

I guess I don't know what my point is. Maybe it's that Julia seems a little
alienating for CS types, and that might hamper adoption.

All in all, Julia's pretty great, despite the minor quibbles I mentioned. I
sure as hell wish I could use it for quant work instead of being stuck with
the vastly inferior Python.

~~~
patrick5415
Well the target user of Julia is the mathy/sciency people, not the cs folks,
as far as I know...

------
Blammar
Hopefully, people are also looking at
[https://github.com/dgasmith/opt_einsum](https://github.com/dgasmith/opt_einsum)
for ideas of what could go into a julia implementation.

~~~
under-Peter
This is pretty much what I'm going for :). At this point in the project we
have an einsum that is slow but works for all cases and supports AD. Next
steps are splitting up contractions into pairwise contractions, finding
optimal orders, evaluating the contractions with the best algorithm (so we
might dispatch matrix multiplication to BLAS) and supporting different
backends (probably just GPU but if more general tensor types could be
supported that'd be great).

------
savant_penguin
"It is slow and uses a lot of memory, but it only requires operations that are
implemented for many array types(...)"

Would

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

be slower than direct matrix multiplication?

~~~
under-Peter
Heya, for standard arrays that are small, the einsum-macro from Einsum.jl can
be faster than BLAS. Einsum.jl translates the code simply into three loops and
if cache-access patterns (presumably what most of BLAS optimisation is
concerned with) isn't an issue, the for-loops are very fast. A quick check
shows that for a simple matrix multiplication, Einsum is faster for 5x5
matrices but slower for 10x10. (This is ignoring the compilation-time that
calling a julia function for the first time incurs)

