14 ms·

Grassmann.jl A\b 3x faster than Julia's StaticArrays.jl

117 points by DreamScatter 6y ago

ssivark 6y ago
Could somebody motivate the comparison being made? It is not apriori obvious what Grassman numbers have to do with linear solvers. Thanks.
- spacedome 6y ago
  They are solving the linear system Ax=b, written in Julia/Matlab usually as A\b, specifically for a fixed size 5x5 matrix. The Grassmann algebras can give a way of representing the matrix and vectors, and a comparison is being made to the representation (and solver) given by StaticArrays.jl, which just uses typical matrices and arrays but optimized for small sizes. No real motivation to use Grassmann algebras if your only concern is generally solving a linear system, other than this method seeming to be faster, which mostly (imo) points to the StaticArrays.jl method needing to be optimized. Edit: This post also does not mention numerical stability of the Grassmann.jl method, which is a real concern in practice, even for such small systems.
  - deleted 6y ago
    [deleted]
- deltasquared 6y ago
  This is a pretty neat trick, they implemented exterior algebra see https://en.wikipedia.org/wiki/Exterior_algebra https://en.wikipedia.org/wiki/Exterior_algebra You see this kind of math in textbooks focused on computer graphics.
- contravariant 6y ago
  Do you have a Physics background by any chance? Anyway, I believe the module might be named after the mathematician, not the number system called after the same mathematician. Regardless the name Grassman is tied to exterior algebras (and therefore geometric algebras), which are implemented by this library and which can be used to solve matrix equations (although I couldn't recommend using the geometric algebra approach for higher dimensions, unless you can do so symbolically).
  - deltasquared 6y ago
    I do have a bit of physics in my background. I read about Grassmann recently in a game programming book. Foundations of game engine development Volume 1 By Eric Lengyel ends with a nice explanation of Grassmann algebra. The road to reality by Roger Penrose also has a treatment of Grassmann algebra in physics context. I agree the library is likely a reference to the mathematician.
    - contravariant 6y ago
      Interesting. Yeah it's nothing too important it's just that I've only ever seen physicists use Grassman numbers, whereas mathematicians tend to work with the exterior algebra directly (which works in much the same way except the wedge operator is used explicitly, IMHO this makes things clearer).
      - DreamScatter 6y ago
        Grassmann algebra and exterior algebra are the same thing, regardless of whether a mathematician or physicist uses them.
        contravariant 6y ago
        Sure, but the Grassman numbers are a bit of an odd model of Grassman algebra with its own notation. They live in some anonymous vector space that is effectively defined in terms of the Grassman numbers, whereas mathematicians usually have some vector space they're interested in and construct an exterior algebra from its elements. The differences are minor but they're two clearly distinct approaches. Like I said it's nothing terribly important but it's still an interesting difference.
        deltasquared 6y ago
        I agree, the subtle difference between the physics and math mindset is interesting. I wound up on the sidelines as a computer scientist who worked in a physics lab.
- 6gvONxR4sf7o 6y ago
  There's a neat kind of linear algebra (multilinear algebra, technically) that doesn't just deal with vectors (think vectors and planes and volumes), but it's often computationally unweildy. This is a fast library for it that competes with simple matrix code in efficiency.
snicker7 6y ago
The real performance win seems to be reducing the number of allocations (2 vs. 30 thousand).
- doublesCs 6y ago
  That's not uncommon in HPC. I'm surprised if this was inherent to StaticArrays.jl though. Probably some interop edge case, I'm imagining. Shall we summon Chris Rackauckas?
  - ChrisRackauckas 6y ago
    It worked, I answered on Discourse because someone pinged me for an answer, saying "you've been summoned" haha.
- marmaduke 6y ago
  So the real question is why a library with the word "static" in its title does 30 thousand allocations..
  - dnautics 6y ago
    I believe the "static" refers to static, aka, compile-time dimension of the array (yes, julia is a compiling language), versus an array that has a run-time, mutable array dimensionality.
  - ChrisFoster 6y ago
    "Static" is a reference to the size of the arrays being fixed and encoded in the type. There's various allocation strategies, for instance StaticArrays.SizedArray uses a normal heap-allocated Array as the storage. That said, the backing storage for the particular static array types used in this example (SMatrix,SVector) is an immutable Tuple which should generally be stack allocated if the library and optimizer are working together as intended. So the allocations here are a bit of a surprise; they may indicate a lack of inlining in something we expected to be inlined.
  - tkoolen 6y ago
thanatropism 6y ago
When I first heard of Julia in 2016 it was sandwiched between R and Python and thought it had no chance whatsoever of succeeding. I’m increasingly impressed by Julia. It’s a crowded space and they remain relevant. Maybe even best of class for some applications.
- plafl 6y ago
  The same first impression for me. It's great people with more faith than me have pushed forward and actually I'm still waiting to make the jump. It was very similar with python for scientific computing in the beginning. When I wrote my master thesis in python in 2006 python was completely dwarfed by Matlab. I was more enthusiastic on those days trying new things and also had more time. It is yet to be seen if Julia displaces python but at least it has non zero chance of doing it.
  - peatmoss 6y ago
    I’ve wanted Scheme / Racket to make a showing in places where R and Python have dominated the data science and numerical computing spaces (for those of us not writing C/C++/Fortran). But, I think Julia is the closest thing to an honest-to-goodness lisp that has a prayer of mainstream adoption. I quite like Julia and think the community has leveraged the meta programming strengths of the language well enough to convince me to relinquish my love of minimalist lisp syntax. I do wonder if the biggest risk to Julia isn’t that it’s not extreme / superlative in any dimension. If you want raw performance, Rust seems to be capturing the lion’s share of excitement among “new” languages. If you want meta programming, Julia still feels like a compromise compared to Scheme/Racket/Clojure/Common Lisp. If you want libraries, Python is still king. If you want stats libraries, R is still king. Julia is impressively good at all of those things... but it’s not good enough to be the most compelling of any of them. Still, it makes me think that it may be time to dive in head first rather than dabble my toes as I've been doing for the past few years. I guess my main barrier at the moment is compiler / deployment story—it still seems hard to build a deployable binary that fits in, say, a Lambda layer as far as I can tell.
    - 6y ago
dklend122 6y ago
I'd be interested to see a static arrays benchmark run on 1.5 as it can now stack allocate views and other types due to general compiler improvements.
- stabbles 6y ago
  StaticArrays.jl is just linear algebra operations for tuples, and they have always been stack allocated as they are immutable and have a size known at compile-time, so Julia 1.5 will not bring anything new for it. The deal with 1.5 is that the GC has improved in that it does not require structs with references to dynamically allocated objects to be heap allocated themselves.
  - dklend122 6y ago
    The allocations for the StaticArrays solve must come from somewhere
dklend122 6y ago
Turns out the allocations are just from the fact that the computation as written was not type stable. Making it so brings down the allocations to 2 and increases the relative speed to just under 2x. Here's the post from chris: "That’s not type-stable. x = [@SMatrix randn(5,5) for i in 1:10000] is a type stable way to compute an array of random SMatrices, and that nearly doubles the speed of the computation and removes the allocations down to 2, i.e. from: 2.549 ms (29496 allocations: 1.44 MiB) to 1.469 ms (2 allocations: 390.70 KiB) for me. Still a nice performance for Grassmann.jl, but seems to be <2x."
ChrisFoster 6y ago
My main concern with this comparison is that there's no mention of numerical stability / accuracy of this method. Calculating determinants is notoriously prone to cancellation error if done naively, so you've got to be really careful. By extension, the usual text book implementation of Cramer's rule is a rather bad way to solve matrix systems. (I don't know much about Grassman.jl so perhaps they've used some known tricks to make this work better, but you do have to be really careful with this stuff.) Actually in early versions of StaticArrays we were pretty cavalier about numerical accuracy because it was just a library we needed for fast geometric calculations. However as people started using this library for serious numerical work we've had to pay a lot more attention to having good numerical properties. There's still a lot to learn - and no doubt much work to do - in upgrading the implementation to make sure everything is first robust, and second as fast as possible.
- DreamScatter 6y ago
  This algorithm is not Cramer's rule, it is analogous to it, but based on the original exterior algebra of Hermann Grassmann.
  - ChrisFoster 6y ago
    I agree the exterior algebra is beautifully elegant mathematics, but expressive elegance of the formalism is largely irrelevant to numerical robustness. (Cramer's rule itself is a perfect example of this: it's elementary to express, but in simple form is a terrible idea numerically.) To be clear, I'm not claiming there's anything wrong with your algorithm but simply adding a note of caution. Any careful comparison between methods must include both efficiency and numerical accuracy.
canjobear 6y ago
I would like to move from Python to Julia because Julia gives easier access to better abstractions for numerical operations. What keeps me on Python is that there is no real replacement for PyTorch.
- eigenspace 6y ago
  Just a thought, but have you considered simply using PyTorch through PyCall.jl or something like ThArrays.jl? This may not be worth it if the majority of your code is reliant on PyTorch, but it's certainly an option that should be taken seriously imo.

117 points by DreamScatter 6y ago

Could somebody motivate the comparison being made? It is not apriori obvious what Grassman numbers have to do with linear solvers. Thanks.

spacedome 6y ago
They are solving the linear system Ax=b, written in Julia/Matlab usually as A\b, specifically for a fixed size 5x5 matrix. The Grassmann algebras can give a way of representing the matrix and vectors, and a comparison is being made to the representation (and solver) given by StaticArrays.jl, which just uses typical matrices and arrays but optimized for small sizes. No real motivation to use Grassmann algebras if your only concern is generally solving a linear system, other than this method seeming to be faster, which mostly (imo) points to the StaticArrays.jl method needing to be optimized. Edit: This post also does not mention numerical stability of the Grassmann.jl method, which is a real concern in practice, even for such small systems.
- deleted 6y ago
  [deleted]
deltasquared 6y ago
This is a pretty neat trick, they implemented exterior algebra see https://en.wikipedia.org/wiki/Exterior_algebra https://en.wikipedia.org/wiki/Exterior_algebra You see this kind of math in textbooks focused on computer graphics.
contravariant 6y ago
Do you have a Physics background by any chance? Anyway, I believe the module might be named after the mathematician, not the number system called after the same mathematician. Regardless the name Grassman is tied to exterior algebras (and therefore geometric algebras), which are implemented by this library and which can be used to solve matrix equations (although I couldn't recommend using the geometric algebra approach for higher dimensions, unless you can do so symbolically).
- deltasquared 6y ago
  I do have a bit of physics in my background. I read about Grassmann recently in a game programming book. Foundations of game engine development Volume 1 By Eric Lengyel ends with a nice explanation of Grassmann algebra. The road to reality by Roger Penrose also has a treatment of Grassmann algebra in physics context. I agree the library is likely a reference to the mathematician.
  - contravariant 6y ago
    Interesting. Yeah it's nothing too important it's just that I've only ever seen physicists use Grassman numbers, whereas mathematicians tend to work with the exterior algebra directly (which works in much the same way except the wedge operator is used explicitly, IMHO this makes things clearer).
    - DreamScatter 6y ago
      Grassmann algebra and exterior algebra are the same thing, regardless of whether a mathematician or physicist uses them.
      - contravariant 6y ago
        Sure, but the Grassman numbers are a bit of an odd model of Grassman algebra with its own notation. They live in some anonymous vector space that is effectively defined in terms of the Grassman numbers, whereas mathematicians usually have some vector space they're interested in and construct an exterior algebra from its elements. The differences are minor but they're two clearly distinct approaches. Like I said it's nothing terribly important but it's still an interesting difference.
        deltasquared 6y ago
        I agree, the subtle difference between the physics and math mindset is interesting. I wound up on the sidelines as a computer scientist who worked in a physics lab.
6gvONxR4sf7o 6y ago
There's a neat kind of linear algebra (multilinear algebra, technically) that doesn't just deal with vectors (think vectors and planes and volumes), but it's often computationally unweildy. This is a fast library for it that competes with simple matrix code in efficiency.

snicker7 6y ago

The real performance win seems to be reducing the number of allocations (2 vs. 30 thousand).

doublesCs 6y ago
That's not uncommon in HPC. I'm surprised if this was inherent to StaticArrays.jl though. Probably some interop edge case, I'm imagining. Shall we summon Chris Rackauckas?
- ChrisRackauckas 6y ago
  It worked, I answered on Discourse because someone pinged me for an answer, saying "you've been summoned" haha.
marmaduke 6y ago
So the real question is why a library with the word "static" in its title does 30 thousand allocations..
- dnautics 6y ago
  I believe the "static" refers to static, aka, compile-time dimension of the array (yes, julia is a compiling language), versus an array that has a run-time, mutable array dimensionality.
- ChrisFoster 6y ago
  "Static" is a reference to the size of the arrays being fixed and encoded in the type. There's various allocation strategies, for instance StaticArrays.SizedArray uses a normal heap-allocated Array as the storage. That said, the backing storage for the particular static array types used in this example (SMatrix,SVector) is an immutable Tuple which should generally be stack allocated if the library and optimizer are working together as intended. So the allocations here are a bit of a surprise; they may indicate a lack of inlining in something we expected to be inlined.
- tkoolen 6y ago

thanatropism 6y ago

When I first heard of Julia in 2016 it was sandwiched between R and Python and thought it had no chance whatsoever of succeeding. I’m increasingly impressed by Julia. It’s a crowded space and they remain relevant. Maybe even best of class for some applications.

plafl 6y ago
The same first impression for me. It's great people with more faith than me have pushed forward and actually I'm still waiting to make the jump. It was very similar with python for scientific computing in the beginning. When I wrote my master thesis in python in 2006 python was completely dwarfed by Matlab. I was more enthusiastic on those days trying new things and also had more time. It is yet to be seen if Julia displaces python but at least it has non zero chance of doing it.
- peatmoss 6y ago
  I’ve wanted Scheme / Racket to make a showing in places where R and Python have dominated the data science and numerical computing spaces (for those of us not writing C/C++/Fortran). But, I think Julia is the closest thing to an honest-to-goodness lisp that has a prayer of mainstream adoption. I quite like Julia and think the community has leveraged the meta programming strengths of the language well enough to convince me to relinquish my love of minimalist lisp syntax. I do wonder if the biggest risk to Julia isn’t that it’s not extreme / superlative in any dimension. If you want raw performance, Rust seems to be capturing the lion’s share of excitement among “new” languages. If you want meta programming, Julia still feels like a compromise compared to Scheme/Racket/Clojure/Common Lisp. If you want libraries, Python is still king. If you want stats libraries, R is still king. Julia is impressively good at all of those things... but it’s not good enough to be the most compelling of any of them. Still, it makes me think that it may be time to dive in head first rather than dabble my toes as I've been doing for the past few years. I guess my main barrier at the moment is compiler / deployment story—it still seems hard to build a deployable binary that fits in, say, a Lambda layer as far as I can tell.
  - 6y ago

dklend122 6y ago

I'd be interested to see a static arrays benchmark run on 1.5 as it can now stack allocate views and other types due to general compiler improvements.

stabbles 6y ago
StaticArrays.jl is just linear algebra operations for tuples, and they have always been stack allocated as they are immutable and have a size known at compile-time, so Julia 1.5 will not bring anything new for it. The deal with 1.5 is that the GC has improved in that it does not require structs with references to dynamically allocated objects to be heap allocated themselves.
- dklend122 6y ago
  The allocations for the StaticArrays solve must come from somewhere

dklend122 6y ago

Turns out the allocations are just from the fact that the computation as written was not type stable. Making it so brings down the allocations to 2 and increases the relative speed to just under 2x. Here's the post from chris: "That’s not type-stable. x = [@SMatrix randn(5,5) for i in 1:10000] is a type stable way to compute an array of random SMatrices, and that nearly doubles the speed of the computation and removes the allocations down to 2, i.e. from: 2.549 ms (29496 allocations: 1.44 MiB) to 1.469 ms (2 allocations: 390.70 KiB) for me. Still a nice performance for Grassmann.jl, but seems to be <2x."

ChrisFoster 6y ago

My main concern with this comparison is that there's no mention of numerical stability / accuracy of this method. Calculating determinants is notoriously prone to cancellation error if done naively, so you've got to be really careful. By extension, the usual text book implementation of Cramer's rule is a rather bad way to solve matrix systems. (I don't know much about Grassman.jl so perhaps they've used some known tricks to make this work better, but you do have to be really careful with this stuff.) Actually in early versions of StaticArrays we were pretty cavalier about numerical accuracy because it was just a library we needed for fast geometric calculations. However as people started using this library for serious numerical work we've had to pay a lot more attention to having good numerical properties. There's still a lot to learn - and no doubt much work to do - in upgrading the implementation to make sure everything is first robust, and second as fast as possible.

DreamScatter 6y ago
This algorithm is not Cramer's rule, it is analogous to it, but based on the original exterior algebra of Hermann Grassmann.
- ChrisFoster 6y ago
  I agree the exterior algebra is beautifully elegant mathematics, but expressive elegance of the formalism is largely irrelevant to numerical robustness. (Cramer's rule itself is a perfect example of this: it's elementary to express, but in simple form is a terrible idea numerically.) To be clear, I'm not claiming there's anything wrong with your algorithm but simply adding a note of caution. Any careful comparison between methods must include both efficiency and numerical accuracy.

canjobear 6y ago

I would like to move from Python to Julia because Julia gives easier access to better abstractions for numerical operations. What keeps me on Python is that there is no real replacement for PyTorch.

eigenspace 6y ago
Just a thought, but have you considered simply using PyTorch through PyCall.jl or something like ThArrays.jl? This may not be worth it if the majority of your code is reliant on PyTorch, but it's certainly an option that should be taken seriously imo.