
Grassmann.jl A\b 3x faster than Julia's StaticArrays.jl - DreamScatter
https://discourse.julialang.org/t/grassmann-jl-a-b-3x-faster-than-julias-staticarrays-jl/41451?u=chakravala
======
ChrisFoster
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
This algorithm is _not_ Cramer's rule, it is analogous to it, but based on the
original exterior algebra of Hermann Grassmann.

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

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

------
dklend122
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
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
The allocations for the StaticArrays solve must come from somewhere

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

~~~
marmaduke
So the real question is why a library with the word "static" in its title does
30 thousand allocations..

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

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

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

~~~
nextos
It has a lot of potential. You can often write code in Julia that is really
performant by just following a few simple rules. Most notably, keeping type
stability. Libraries are now catching up. Some are great, some are unique,
others are still missing.

I have generated SIMD code that outperforms C and C++, even after some hand
optimizations. As a matter of fact, some people are re-writing a BLAS drop-in
replacement in pure Julia, and the performance is not going to be too far off
from BLAS [1]. If you work in scientific computing, you can probably
appreciate how amazing that is.

The combination of multiple dispatch, a simple type system and an LLVM backend
is wonderful. In the long run, Julia will probably grab lots of Python's and
R's marketshare because having a single language instead of two languages
(Python + Cython/Pythran/... or R + C++) is a huge advantage.

[1] [https://discourse.julialang.org/t/we-can-write-an-
optimized-...](https://discourse.julialang.org/t/we-can-write-an-optimized-
blas-library-in-pure-julia-please-skip-op-and-jump-to-post-4/11634/4)

~~~
krapht
Maybe Julia will become more popular, but I'd bet on the 800 pound gorilla
Python getting faster through broader use of Dask, Cython, and Numba. I'll
never bet against the network effects of language and library development.

~~~
ta17711771
That's why Go will supersede Python. All the new cool and secure small shit is
being done in Go.

~~~
dunefox
Please no.

~~~
ta17711771
Elaborate?

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

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

