Almost no part of the algorithm as specified in the blog post is GPU friendly, but it can be improved.
The "symbolic" phase is very hard to port to the GPU in the first place and almost always happens on the CPU side. But even then it's hard (but not impossible) even to parallelize it. This is generally OK because it is usually cheaper than the numeric phase anyways. But it can be trouble on workstations where the CPU exists primarily as a device for shoveling data into many GPUs, you don't want to stall their kernel pipelines with your symbolic analysis.
For the numeric phase you need to introduce the concept of either a "front" or a "supernode." This is a technique where multiple columns get batched together and you get a new elimination tree in terms of those batches of columns rather than individual columns. This turns a lot of irregular memory access (gather/scatter style updates) into densely addressed memory patterns, often just calls into level 3 BLAS or LAPACK.
There are some sparsity patterns for which the above supernode/multifrontal approach does not work very well, but for many practical cases like PDE simulations it does. The sparsity patterns which defeat it usually spread nonzeros out all over the matrix rather than containing them within a specific band. What this causes in practice is when you materialize a supernode's dense matrix, most of its numeric entries will just be 0s, so you create a lot of extra work that serves no productive purpose.
I have no submission for this but I joined the hype in my own way by optimizing the training loop. These tiny models are not really well suited to frameworks like pytorch, and with highly patient AI agents we can now just inline the whole thing into C++ just to see what happens, which I do below:
Basically I just don't want to hear about "the state of SIMD in Rust" unless it is about dramatic improvement in autovectorization in the rust compiler.
80%-90% or so of real life vectorization can be achieved in C or C++ just by writing code in a way that it can be autovectorized. Intrinsics get you the rest of the way on harder code. Autovectorization is essentially a solved problem for the vast majority of floating point code.
Not so with rust, because of a dogmatic approach to floating point arithmetic that assumes bitwise reproducibility is the "right" answer for everyone (actually, it's the right answer to almost nobody) to the point of not even allowing a user to even flag on these optimizations. and once you get to the point of writing intrinsics you have to handwrite code for every new architecture when autovectorizers could have gotten you 80%-90% of the way there with a single source and often this is just enough.
the contention with the above is that if a user needs SIMD they can just use some SIMD API and make their intention more clear. this is essentially an argument that we should handwrite intrinsics. well guess what. I'm a programmer and I use compilers because they _do this for me_ and indeed are able to do so very easily in C or C++ when I instruct it that I'm ok with with reordering operations and other "accuracy impacting" optimizations.
The huge joke on us is that these optimizations generally have the effect of _improving_ accuracy because it will reduce the number of rounding steps either by simply reducing the number of operations or by using fused multiply adds which round only once.
There are two sides to numerical linear algebra. The first is the "linear algebra" part, which is very mathematically sophisticated and the language you choose to represent these concepts is not so important as your understanding. pencil and paper is the ideal place to prove out understanding of this.
The "numerical" part is a minefield because it will take all your math and demolish it. just about every theoretical result you proved out above will hold not true and require extra-special-handholding in code to retain _some_ utility.
As such I think a language which enables you to go as fast as possible from an idea to seeing if it crosses the numerical minefield unscathed is the one to use, and these days that is python. It is just so fast to test a concept out, get immediate feedback in the form of plotting or just plain dumb logging if you like, and you can nearly instantly share this with someone even if you're on ARM +linux & they are Intel+windows
The most problematic issue with python&numpy, as it relates to learning _numerical_ side of linear algebra, is making sure you haven't unintentionally promoted a floating point precision somewhere (for example: if your key claim is that an algorithm works entirely in a working precision of 32 bits but python silently promoted your key accumulator to 64 bits, you might get a misleaing idea of how effective the algorithm was) but these promotions don't happen in a vacuum and if you understand how the language works they won't happen.
edit: & I have worked professionally with fortran for a long time, having known some of the committee members and BLAS working group. so I have no particular bias against the language
Thanks! The chunks trick was a fairly straightforward translation of what I would do in C++ if the compiler wouldn't vectorize the reduction for some reason. These days most compilers will do it if you pass enough flags, a fact I really took for granted when doing this because Rust is more conservative.
I've tried using mul_add, but at the moment performance isn't much better. But I also noticed someone else on my machine running a big parallel build, so I'll wait a little later and run the full sweep over the problem sizes with mul_add.
So really the existence of FMA didn't have a performance implication it seems except to confirm that Rust wasn't passing "fast math" to LLVM where Clang was. It just so happens that "fast math" will also allow vectorization of reductions.
Hi I wrote the blog post linked - and I feel a little silly that I didn't check that _both_ loops vectorized. So I fixed the Rust implementation to keep a running vector of partial sums which I finish up at the end - this one did vectorize. The result was a 2X performance bump, which I'm about to include in the blog post as an update.
If it's OK I'll link to this comment as the inspiration.
On the iterators versus loop: for some reason when I use the raw loop _nothing_ vectorizes, not even the obvious loop. What I read online was that bounds checking happens inside the loop body because Rust doesn't know where those indices are coming from. Using iterators instead is supposed to fix this, and it did seem to in my experiments.
I wasn't clear enough. I was responding to the flow of comments of the form: "Riemann hypothesis is hard, this is unlikely to be true." Sure it's true, but doing a little more research could inform that opinion well past the zeroth-order approximation of "it's a hard problem."
I didn't actually take my own advice, I just wait for Terence Tao to write a post then I know it's true :)
Maybe a better term should be "blind" rather than black-box. I think the goal is simply to hold optimization to the same level of reproducibility that is expected of most scientific fields today, and if a researcher is allowed to introduce a hundred tunable parameters that makes their algorithm converge on all the standard test cases then they haven't created a reproducible optimizer - they have created a benchmark solver.
How is this for "performance portability?" I use this solution in C when the function is very expensive, therefore a little extra indirection really won't make any difference - but can potentially improve the reusability a lot. Is inlining calls to function via a pointer a very basic optimization that any self respecting compiler should be able to do, or is it a very advanced optimization that I can't count on working across platforms? Given the possibility of dynamic libraries I don't see how it could be inlined in all cases, therefore at least some kind of analysis must be done before trying it.
The idea of using DSLs in this field is attractive. The amount of math you have to hold in your head to really start talking about large scale simulations is staggering, and it gets more and more complicated as you scale up. There are thousands of moving parts that must be accounted for, any one of them can burn you.
I wonder however, with full respect to other groups who have made good progress in making the idea a reality, why they never attempt collaboration with computer science - a field that has been perfecting language design for so long? I remember speaking to someone doing a similar project, and they didn't even know that parsing was a "thing," yet they were implementing a DSL.
Since I'm more on the math side and not on the computer science side I can't comment on how much duplicated work is actually happening, but it seems like a lot.
Very handy formula. I was recently in the situation of being given a bunch of polygons, each represented in a list-of-edges format and I needed to compute the outward-pointing normals for the edges. It's relatively straightforward to translate list-of-edges into list-of-vertices, but not so straightforward (to me) to ensure that the vertices are in counterclockwise order after this translation.
If you take the shoelace formula and remove the absolute value then you still get the correct area up to floating point errors, just there may be a minus sign in front of it. If there is a minus sign, reverse the order of your vertices and boom you have counterclockwise order. Computing the normals after this is routine.
Symbolic math is not the same thing as formalized math, and doesn't carry with it the same guarantees of correctness. That's the price one pays for having black box functions spit out answers to hard problems in a reasonable time, they won't spit out proofs of correctness alongside their answer.
Any symbolic system should be treated with care, but they can of course be extremely useful. My typical use case is to have it compute very challenging symbolic expressions for me. I then treat it as a "very plausible hypothesis" which I then prove.
Could someone comment on the following C++14 feature on the list here:
"Avoiding/fusing memory allocations"
It lists Clang but nothing else as supporting it.
This sounds like a compiler optimization and not a language feature, and since it's talking about C++14 (not C++11) I know this isn't referring to move semantics/rvalue references.
The lack of price transparency and emergence of 3rd party system was direct result of heavy regulations on health care dating back to WW2, and was slowly solidified into law through tax incentives and then recently became "official" through the ACA.
I don't know how useful it will be in the end. The idea of targeting a powerful (proper) subset of turing completeness to produce fast executables makes a lot of sense to me, and furthermore for reasons that you have already brought up this is the only way it is possible if we truly desire a mechanized solution to this problem.
My concern on the practicality point is that the superoptimizer here really isn't "superoptimizing" so much as super-simplifying. There's a reason modern compilers don't aggressively inline every single function called if they have access to the source code, and for the same reasons they don't aggressively unroll every loop they come across if the loop bounds are known compile time.The tradeoffs of such optimizations are very complicated and hard to reason about because processors are very complicated and hard to reason about.
It's a little bit of a misrepresentation on the part of the presenter. BLAS is a specification for linear algebra software, of which implementers are free to use the most appropriate representation of numbers as they wish - the fastest being floating point because of native hardware support.
Although the naming schemes for BLAS routines do explicitly reference "single precision" and "double precision" so now I'm confused too.
Did the author claim Morte was Turing-complete? I don't think I saw that claim, but I'm a newcomer to all this computer science theory stuff. I think these Wikipedia links answer your question:
The "symbolic" phase is very hard to port to the GPU in the first place and almost always happens on the CPU side. But even then it's hard (but not impossible) even to parallelize it. This is generally OK because it is usually cheaper than the numeric phase anyways. But it can be trouble on workstations where the CPU exists primarily as a device for shoveling data into many GPUs, you don't want to stall their kernel pipelines with your symbolic analysis.
For the numeric phase you need to introduce the concept of either a "front" or a "supernode." This is a technique where multiple columns get batched together and you get a new elimination tree in terms of those batches of columns rather than individual columns. This turns a lot of irregular memory access (gather/scatter style updates) into densely addressed memory patterns, often just calls into level 3 BLAS or LAPACK.
There are some sparsity patterns for which the above supernode/multifrontal approach does not work very well, but for many practical cases like PDE simulations it does. The sparsity patterns which defeat it usually spread nonzeros out all over the matrix rather than containing them within a specific band. What this causes in practice is when you materialize a supernode's dense matrix, most of its numeric entries will just be 0s, so you create a lot of extra work that serves no productive purpose.