6 ms·

Einstein Summation in Numpy

81 points by ovi256 8y ago

s9w 8y ago
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 8y ago
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 8y ago
  >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-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 https://github.com/dgasmith/opt_einsum#news-opt_einsum-will-...
  - rax 8y ago
    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.
  - dgasmith 8y ago
    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 8y ago
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.
- stochastic_monk 8y ago
  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.
  - enriquto 8y ago
    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.
    - hprotagonist 8y ago
      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 https://docs.python.org/3/library/array.html
      - stochastic_monk 8y ago
        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.
  - dagw 8y ago
mlthoughts2018 8y ago
> “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!
BlackFly 8y ago
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 8y ago
  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 8y ago
    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.
singularity2001 8y ago
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 8y ago
  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 8y ago
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-neural-networks-39c097d2925b https://medium.com/inbrowserai/simple-diagrams-of-convoluted....
- p1esk 8y ago
  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/capsnet_cifar.py#L68 https://github.com/michaelklachko/CapsNet/blob/master/capsne...

81 points by ovi256 8y ago

s9w 8y ago

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 8y ago

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 8y ago
>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-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 https://github.com/dgasmith/opt_einsum#news-opt_einsum-will-...
- rax 8y ago
  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.
- dgasmith 8y ago
  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 8y ago

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.

stochastic_monk 8y ago
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.
- enriquto 8y ago
  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.
  - hprotagonist 8y ago
    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 https://docs.python.org/3/library/array.html
    - stochastic_monk 8y ago
      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.
- dagw 8y ago

mlthoughts2018 8y ago

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

BlackFly 8y ago

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 8y ago
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 8y ago
  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.

singularity2001 8y ago

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 8y ago
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 8y ago

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-neural-networks-39c097d2925b https://medium.com/inbrowserai/simple-diagrams-of-convoluted....

p1esk 8y ago
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/capsnet_cifar.py#L68 https://github.com/michaelklachko/CapsNet/blob/master/capsne...