Einstein Summation in Numpy | Hacker News Reader

Einstein Summation in Numpy

obilaniu6266h16.wordpress.com

81 pointsovi25641 comments

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.

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.

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.

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...

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.

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->"

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

Vide for PyTorch: 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....

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...

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.

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.

>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-...

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.

I agree about the value of the article. I guess I got tilted as I was hoping for another good, more fresh, post of similar quality!

Using the Optimize kewarg is a must for any summation involving a few inputs, which is why I brought it up. It is one of the things which is not discussed in a post from 2016.

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.

Python has statically typed arrays in the standard library: in the array module (e.g. array.array).

It’s not used much for scientific computing because numpy made a better tool.

Also, the list data structure in Python is a contiguous memory array... it’s just designed to store a contiguous array of PyObject pointers, and the use case for list is exactly when the flexibility and abstraction of using PyObject data structures is more important than raw performance. You would use numpy arrays or special jitted functions or custom C-level data structures you write in Cython, etc. etc., when raw performance is the higher priority.

My point exactly. Python was not thought from the beginning as a language for numerical computing, and when you try to force it for that, there is a lot of friction. Just having to "import" a library to have sane support for arrays is silly.

I feel like that's the point of Python - only provide the barebones of what you might need initially, and then make importing the rest fast and easy. For most use cases Pythons default collections work just fine (what does an array have that a list doesn't for everyday use?), but if you need more power, just import and go.

From my point of view, numeric arrays are barebones, and strings and dictionaries are fancy data structures that I never need, but they are nevertheless elevated as features of the language itself. It would seem much more natural if it was the other way round. In my everyday life I solve linear systems in octave by doing "x=A\b", but in python it becomes very cumbersome.

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

> 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.

Semi-inactive Octave dev here.

Even though I use Python every day, and I have used scipy often, I still reach for Octave whenever I need a quick calculation that I don't feel like turning into a big project.

And yeah, it's nice that matrices are first-class objects. Sometimes I miss 1-dimensional arrays or sensible strings, but it's not a big deal, and for quick, short things, the Octave syntax is still a lot easier to deal with than Python.

Most of Octave's users describe themselves as "not a programmer", despite using a programming language, so I think there's still room for this kind of software in the world.

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

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!

The big challenge with leaving core data structures and interfaces to the community is fragmentation. If you look at the rust/go numerical communities the big struggle has been a lack of universal ndarray implementation.

Most numeric libraries can be Described as functions of/on ndarrays that return ndarrays, not having a language wide ndarray is the equivalent of having every library carry it's own unique list interface/implementation.

Well, you could still have a set of "core" packages for basic data structures, or an official library. For example, C++ has STL, and most libraries build on top of that.

Probably explains Fortans continued hold on numerical computing

> Probably explains Fortans continued hold on numerical computing

indeed! a sizeable part of numpy is written in fortran

Only insofar as numpy relies on BLAS and LAPACK for certain functionality, and most (all?) implementations are written in Fortran. The actual numpy source code contains no Fortran; you don't need a Fortran compiler to build numpy. Scipy does contain Fortran, though.

"certain functionality" is a bit of an understatement, here. Those are precisely the core algorithms, that allow you to do things with matrices, solve systems, etc. As you put it, numpy is simply a wrapper for this.

I think that's an unfair characterization of numpy's functionality. Numpy introduces data structures not present in the core python library (another one of your complaints about python, though it's not clear to me why it matters that numpy is a separate library), and certain manipulations do not require BLAS/LAPACK acceleration, eg. multiplying an array by a scalar (yes, I know dscal exists in BLAS) or array slicing. Further, it is possible to use numpy without any BLAS/LAPACK libraries at all, and it continues to work (albeit more slowly).

Note that most accelerated implementations of BLAS are not written in Fortran, but use e.g. SIMD instructions directly.

I'm not sure I understand your point. What's the downside of numpy arrays being "non-native"?

Dunno man. I used MATLAB for most of my thesis, and once I discovered NumPy, I switched all my code to it in a heartbeat. Yes, MATLAB is better for matrix manipulation. For everything else, Python is better. My productivity went up quite a bit.

The secret of compiled languages (like Rust and Go) is that they are all mainly designed to implement compilers. After all, self-hosting is almost always the first and most formative software engineering goal. ;)

Would this be solved by inclusion of numpy into python’s standard library? Is there a chance it might? I know it was behind syntax changed in 3.

You don't have to include all of numpy into python standard library. Unifying python's existing 'array' package with numpy arrays and allowing logical and arithmetic operation and the functions in the standard lib 'math' package to operate on those arrays would get you very far. Then you can still have numpy as a separate package for people who wanted more advanced functions.

What do you mean exactly? Can I do matrix products and slices without importing anything? Yes, that would solve the problem. But I do not think that this would be possible without major changes in the language itself.

This is not a library problem, it is a focus problem of the language. The language supports strings and dictionaries natively, as basic data types, but not numeric arrays.

yes it does. nobody actually uses them, but a non-list ordered sequence of values is actually part of the library.

https://docs.python.org/3/library/array.html

For one-dimensional arrays, python's array can be faster. It's also directly accessible through the Python C API/Cython's wrapping of it, which makes it convenient to use.

> “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!

Reply on news.ycombinator.com