
Einstein Summation in Numpy - ovi256
https://obilaniu6266h16.wordpress.com/2016/02/04/einstein-summation-in-numpy/
======
BlackFly
Where is the slot for the metric so that the summation doesn't blindly treat
my tensor's covector slots like vector slots? Maybe I am being a bit too
cynical, but it is a shame that people have come to know about "Einstein
summation convention" without ever learning that it is a special case of Ricci
calculus and was not invented by Einstein. Every time I see someone do u_i v_i
as an implied sum I know I am going to have a hard time. "The sum over the 0
index has a minus sign," yes because you need to contract the metric there to
convert between a covector and a vector to get a sum.

I made a C++ operator overloading implementation of Ricci calculus once upon a
time. That way you could define your tensors and then write statements like:

    
    
      T = (rho + p)*u("^a")u("^b") + p * g("^a^b");
      traceT = T("^a_a");
    

Anyways, I find the domain specific language of this particular numpy function
to be a bit abhorrent. Why not write "ik,kj->ij" as "M_ij = A_ik * B_kj" and
make it much more obvious what is happening instead of requiring this
understanding. Heck, if you used __kwargs, you could even avoid passing the
same argument multiple times like in the quadratic form example:

    
    
      x = np.ricci_calculus("M_ij v^i v^j", M, v)
    

Here I am leaving out that "x = ..." prefix in the DSL since you can take an
ordering convention and a scalar has no slots anyways.

~~~
mcabbott
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](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.

~~~
BlackFly
Yeah, that's the stuff right there!

Well to distinguish covectors you need to define something like a manifold or
just any canonical mapping between the tangent and cotangent space and then
when you define your first few tensors you define them in the scope of that
mapping.

There is some generalization to gauge fields where a bijective canonical
mapping may not be possible in general, but I never studied that extensively.

~~~
mcabbott
I guess lots of things distinguish rows & columns which, for complex numbers,
you can think of as up & downstairs -- rows are conjugated. But otherwise
everybody seems to live in flat space.

If you had reason to do so, it would be easy to bolt something onto the front
of @tensor etc. which understood say A^[i]_[j] * B^[j,k] . Something similar
would not be so hard to do to np.einsum either, just messing with the
string... it could take exactly your "email Latex notation". In fact I'm a bit
surprised that I can't find someone who's done this, in a few minutes'
googling. I did find that opt-einsum has a notation in which the indices at
least follow the array:

[https://optimized-
einsum.readthedocs.io/en/latest/input_form...](https://optimized-
einsum.readthedocs.io/en/latest/input_format.html#interleaved-input)

~~~
BlackFly
Flat space and cartesian coordinates with a couple of known exceptions is the
norm (we do use some curvilinear coordinates where some differential geometry
methods are not necessary but can be illuminating).

You have honestly piqued my curiosity towards Julia. I will need to find an
opportunity to dig in.

------
singularity2001
For all of the listed applications there is a cleaner function in numpy (and
tensorflow, right?)

    
    
        Vector inner product: "a,a->" (Assumes two vectors of same length)
        Vector element-wise product: "a,a->a" (Assumes two vectors of same length)
        Vector outer product: "a,b->ab" (Vectors not necessarily same length.)
        Matrix transposition: "ab->ba"
        Matrix diagonal: "ii->i"
        Matrix trace: "ii->"
        1-D Sum: "a->"
        2-D Sum: "ab->"
        3-D Sum: "abc->"
        Matrix inner product "ab,ab->" (If you pass twice the same argument, it becomes a matrix L2 norm)
        Left-multiplication Matrix-Vector: "ab,b->a"
        Right-multiplication Vector-Matrix: "a,ab->b"
        Matrix Multiply: "ab,bc->ac"
        Batch Matrix Multiply: "Yab,Ybc->Yac"
        Quadratic form / Mahalanobis Distance: "a,ab,b->"

~~~
blt
Sure. The point of those examples is to use familiar operations to help the
reader understand how the notation works, and show the generality.

------
stared
Vide for PyTorch:
[https://rockt.github.io/2018/04/30/einsum](https://rockt.github.io/2018/04/30/einsum)
("Einsum is All you Need - Einstein Summation in Deep Learning" by Tim
Rocktäschel).

The same convention as in NumPy, but IMHO much nicer visual explanation for
beginners. Also, think about tensor diagrams, as in the Feynman diagram in
"Simple diagrams of convolutional neural networks"
[https://medium.com/inbrowserai/simple-diagrams-of-
convoluted...](https://medium.com/inbrowserai/simple-diagrams-of-convoluted-
neural-networks-39c097d2925b).

~~~
p1esk
Nice! I remember trying to implement Hinton's capsules layer in tensorflow
[1], and stumbling upon Einstein notation syntax. Took me a while to figure it
out.

[1]
[https://github.com/michaelklachko/CapsNet/blob/master/capsne...](https://github.com/michaelklachko/CapsNet/blob/master/capsnet_cifar.py#L68)

------
s9w
This is a very good overview, much better than the official documentation. I
used einsum extensively for research in optical signal transmission and found
it to be insanely fast, but often had problems writing the syntax. This would
have been a godsend.

------
rax
This is a reasonable resource, which I have consulted before, but it's from
2016 so the title should be updated.

OptEinsum is also worth checking out if you use Einsum a lot, it has
additional features to what is in Numpy.

~~~
ovi256
>it's from 2016 so the title should be updated

I will try to update the title. There are no changes since then that diminish
this document's value - the einsum notation is still the same. The biggest
change since then is that the numpy.einsum implementation has been
tremendously optimized and has become the benchmark in capabilities that other
libraries follow.

Edit: Opt_Einsum is amazing and was a big part of the optimization efforts -
it has been merged fully into numpy since v1.12. See
[https://github.com/dgasmith/opt_einsum#news-opt_einsum-
will-...](https://github.com/dgasmith/opt_einsum#news-opt_einsum-will-be-in-
numpy-112-and-blas-features-in-numpy-114-call-opt_einsum-as-npeinsum-
optimizetrue-this-repository-contains-more-advanced-features-such-as-dask-or-
tensorflow-backends-as-well-as-a-testing-ground-for-newer-features-in-this-
ecosystem)

~~~
dgasmith
Note that `opt_einsum` is no longer fully merged with NumPy (and likely will
not be) as we now support a wider range of inputs, backends (TensorFlow, Dask,
etc), and more powerful optimization algorithms for expressions with hundreds
of tensors.

Disclaimer: I am the author of the `opt_einsum` package.

------
enriquto
Sometimes I wake up suddenly at night, sweating, terrified, and crying like a
baby, because python+numpy has become a sort of standard in numerical
computing.

Arrays of numbers should be a native type of a programming language used for
serious numerical computing. Strings and dictionaries are nice, but I can do
without. However, arrays of numbers are fundamental and it is _very_ sad that
they are accessible only via idiosyncratic libraries such as numpy or pytorch.

~~~
amelius
On the contrary, I prefer my languages lightweight, and capable of allowing
arrays (and strings, dictionaries, etc) to be implemented by libraries.

~~~
enriquto
> On the contrary, I prefer my languages lightweight, and capable of allowing
> arrays (and strings, dictionaries, etc) to be implemented by libraries.

I agree with this philosophy! But that is certainly not at all the case with
python, where fancy strings and even fancier dictionaries are readily
available at the language level, instead of being implemented as separate
libraries that you have to import.

If Python did not support strings nor dictionaries as basic data types, you
might have a point. But it does, and it chooses to ignore multidimensional
arrays and their basic operations. This is weird, and even a bit _offensive_
when all you need to do is dealing with such arrays of numbers. You always
feel like an "outsider", unlike when you program in octave or julia, for
example.

~~~
altairiumblue
To give a different perspective:

\- the existence of tools for working with strings and dictionaries doesn't
take away from numpy at all

\- I have never been offended by a programming language or felt like an
outsider while using numpy

~~~
mcabbott
It always seems a bit weird to me that np.einsum needs to parse a string in
its own private syntax. I guess it's a good thing that at least it has tools
for this!

------
mlthoughts2018
> “Despite its generality, it can reduce the number of errors made by computer
> scientists and reduce the time they spend reasoning about linear algebra. It
> does so by being simultaneously clearer, more explicit, more self-
> documenting, more declarative in style and less cognitively burdensome to
> use.”

Einsum will also make you more youthful, you’ll grow back missing hair, be a
big success on the dating scene, get a promotion at work, and feel a new lease
on life!

