Vectorization: Introduction
cvw.cac.cornell.edu
cvw.cac.cornell.edu
I guess this is a scientific computing course or something but I feel like even so it’s important to point out that most processors have many vector instructions that operate on integers, bitfields, characters, and the like. The fundamental premise of “do a thing on a bunch of data at once” isn’t limited to just floating point operations.
> Array programming, a style of computer programming where operations are applied to whole arrays instead of individual elements
> Automatic vectorization, a compiler optimization that transforms loops to vector operations
> Image tracing, the creation of vector from raster graphics
> Word embedding, mapping words to vectors, in natural language processing
> Vectorization (mathematics), a linear transformation which converts a matrix into a column vector
Vector (disambiguation) https://en.wikipedia.org/wiki/Vector
> Vector (mathematics and physics):
> Row and column vectors, single row or column matrices
> Vector space
> Vector field, a vector for each point
And then there are a number of CS usages of the word vector for 1D arrays.
Compute kernel: https://en.wikipedia.org/wiki/Compute_kernel
GPGPU > Vectorization, Stream Processing > Compute kernels: https://en.wikipedia.org/wiki/General-purpose_computing_on_g...
sympy.utilities.lambdify.lambdify() https://github.com/sympy/sympy/blob/a76b02fcd3a8b7f79b3a88df... :
> """Convert a SymPy expression into a function that allows for fast numeric evaluation [with the CPython math module, mpmath, NumPy, SciPy, CuPy, JAX, TensorFlow, SymPt, numexpr,]
pyorch lambdify PR, sympytorch: https://github.com/sympy/sympy/pull/20516#issuecomment-78428...
sympytorch: https://github.com/patrick-kidger/sympytorch :
> Turn SymPy expressions into PyTorch Modules.
> SymPy floats (optionally) become trainable parameters. SymPy symbols are inputs to the Module.
sympy2jax https://github.com/MilesCranmer/sympy2jax :
> Turn SymPy expressions into parametrized, differentiable, vectorizable, JAX functions.
> All SymPy floats become trainable input parameters. SymPy symbols become columns of a passed matrix.
I'd say auto-vectorization should still be learned by modern high-performance programmers, because its very low-hanging fruit. You barely have to do anything and suddenly your for-loops are AVX512 optimized, though maybe not to the fullest extent possible.
Still, I suggest that programmers also learn how to properly make SIMD code. Maybe intrinsics are too hard in practice, but ISPC, CUDA, and other SIMD-programming environments make things far easier to learn than you might expect.
------------
ISPC in particular is Intel's SIMD programming language, much akin to CUDA except it compiles into AVX512. So for AVX512-like code execution environments, using the ISPC language/compiler is exceptionally useful.
Its harder to learn a new language than to learn a few compiler-options to enable auto-vectorization however. So in practice, auto-vectorization will continue to be used. But for tasks that specifically would benefit from SIMD-thinking, the dedicated ISPC language should be beneficial.
Somewhat related: 55 GiB/s FizzBuzz https://codegolf.stackexchange.com/questions/215216/high-thr... (discussion: https://news.ycombinator.com/item?id=29031488)
If one wants to sum two arrays instead of for looping through the elements of two arrays one might instead use iterate by chunks so that the compiler can easily tie them all together as a single operation which can then easily vectorize.
The trick with coding for auto-vectorization is to keep your loops small and free of clutter.
I don't have the documentation handy but I think you only need to follow a couple rules:
- loop must have a defined size (for-loop instead of while-loop)
- don't muck with pointers inside the loop (simple pointer increment is okay)
- don't modify other variables (only the array should be modified)
Microsoft, GCC and Clang all do this too, though with different compiler flags and messages.
I'd say that the whole point of this document listed here is to build up the programmer to understanding these error messages and specifically know how to fix the errors that causes a autovectorization-fail.
E.g. when you do summation using a for loop, and the compiler is trying to preserve the rounding-correctness your float operations. Or because you created a dependency graph with a for loop where each iteration requires the previous iteration's result.
TLDR: SSE1 and SSE2 are parts of the base AMD64 ISA, they are always available while compiling for 64 bits. SSE3, SSE4.1 are old, market penetration is very close to 100%, always available in practice. For the newer extensions like AVX and FMA, the optimal tradeoff depends on the product and the market.
If they were compiling for least-common denominator like SSE2, recompiling with -march=native might make a difference. The exact amount of that difference depends heavily on the code.
Automatic vectorizers are limited, they only generating good vectorized code for pure vertical operations which handle long dense arrays of FP32 or FP64 values.
OTOH, some widely-used libraries like Eigen https://eigen.tuxfamily.org/index.php?title=Main_Page are using manually-written intrinsics under the hood, selecting an implementation with preprocessor macros. If your R library uses Eigen, compiling for AVX2 or AVX512 will probably cause a non-trivial performance win.
If your library uses runtime dispatch to select best implementation in runtime, and is compiled recently (i.e. you don’t have new instruction sets in your CPU which weren’t available back then), you’re unlikely to get much profit.
BTW, I’ve noticed C standard library implemented by VC++ uses runtime dispatch, functions like memcpy() are using AVX vectors internally, when supported.
github.com/google/highway supports compiling your single source code once into several versions, and then dispatching to the best one.
Disclosure: I am the main author.
One of the common interview questions I ask is computing the sum of an array of doubles. Barely 1% of applicants can even get the basics correctly such as having multiple partial sums.
Once you do the work such that your code actually allows execution to happen in parallel, then leveraging superscalar and simd happens automatically.
Also they get asked multiple exercises, not doing great at one is not eliminatory.
Being able to quickly identify whether an individual is exceptionally good is very important to ensure the employer can fast-track them and be responsive quickly enough that he doesn't accept another offer.
if you're working in those contexts then it's all fine, but it's important to acknowledge that most software development is bottlenecked by concerns far above this level of concern
Ehhh... I've got multiple hats from my college classes.
Yes, optimizations say that multiple partial sums is fastest. But... my floating-point analysis classes has my brain yelling "cancellation error". The only way to minimize cancellation error is to sort the numbers from biggest to smallest, and then separate them from positives-and-negatives. You add the numbers from smallest to largest in magnitude, only matching signs, and then finally combine them at the end.
The only two partial sums you should be keeping are "positive numbers" and "negative numbers", and everything else needs to be folded together to minimize rounding and cancellation errors.
--------
You know. For stability and minimization of errors. There's other, faster, methods to possibly be employed. But optimization of partial sums is probably the last thing I'm thinking if I see "float" or "double".
If things are "int" or "long", I'll have a higher chance of actually thinking (and caring) about optimization, because stability/cancellation error isn't a thing with int or long.
Yes, floating-point isn't associative, so the results aren't numerically equivalent, and that is why you need to do that optimization yourself since the compiler isn't allowed to do it for you.
As per the exercise, there is no information about the distribution of the data, so there is no indication that any order would be any better than another for numerical stability.
If you care about stability you'd implement a kahan sum, which isn't really affected by changing the order of the sum so that you can parallelize it.
The only way you fix cancellation error is summing the positives in one variable, and the negatives in another, and do a big cancellation at the end.
Anyway, my overall point is not to criticize Kahan summing or other techniques. My point is that these techniques are incompatible with optimizations. The faster way is often times a less stable operation.
What is a "proper" SIMD tool?
As such, it's easier to write high performance code in practice. Intel's ISPC is the closest tool that replicates this effect for AVX512.
The one about avoiding function calls is a no-go, since things like + and * are function calls too. But is it still required for all function calls to be inline-able?
Indirect indexing isn't mentioned in Julia's case, but is it one of the things that would disable auto-vectorization and require an explicit @simd annotation?
[1] https://viralinstruction.com/posts/hardware/#74a3ddb4-8af1-1... [2] https://cvw.cac.cornell.edu/vector/coding_vectorizable
`@simd ivdep` may be needed when you have indirect indexing. Do not forget the `ivdep`, but only if you know it is actually safe. This can lead to hard to track down bugs if it is not.
RISC-V has, or at least had, some classic vector extension proposal. Not sure if it actually got implemented.
(I'm sure you know this; this is just for the benefit of anyone else perusing the comment section.)
> "Perhaps the most interesting part of the open RISC-V instruction set architecture (ISA) is the vector extension (RISC-V "V"). In contrast to the average single-instruction multipe-data (SIMD) instruction set, RISC-V vector instructions are vector length agnostic (VLA). Thus, a RISC-V "V" CPU is flexible in choosing a vector register size while RISC-V "V" binary code is portable between different CPU implementations."
RISC-V also has the P extension, which follows this pattern. It will be very cool if RISC-V proves flexible enough to really enable a resurgence of classic vectors though.
I think the RISC-V V extension, as well as ARM SVE, both have something similar as well. What RVV and SVE have, that Cray vector ISA didn't, is the ability to support different HW vector register lengths with the same ISA.
After all most tensors are implemented as 1D arrays with strides for each dimension with the innermost dimension always contiguous.
Sure! Let me break it down for you:
When we talk about "tensor ops," we're referring to operations performed on tensors, which are multidimensional arrays of numbers commonly used in mathematics and computer science.
Now, these tensor operations are typically implemented using a technique called "vectorizing code." In simple terms, vectorizing code means performing operations on entire arrays of data instead of looping through each element one by one. It's like doing multiple calculations at once, which can be more efficient and faster.
Tensors are usually represented as 1D arrays, meaning all the elements are arranged in a single line. Each dimension of the tensor has a concept called "stride," which represents how many elements we need to skip to move to the next element in that dimension. This helps us efficiently access and manipulate the data in the tensor.
Additionally, the innermost dimension of a tensor is always "contiguous," which means the elements are stored sequentially without any gaps. This arrangement also aids in efficient processing of the tensor data.
So, in summary, tensor ops involve performing operations on multidimensional arrays, and we use vectorized code to do these operations more efficiently. Tensors are represented as 1D arrays with strides for each dimension, and the innermost dimension is always stored sequentially without any gaps.
If the OP cares about implementation details about how an API like PyTorch is made, I think the MiniTorch 'book' is a pretty good intro:
I think accelerators mostly do computations over small matrices and not vectors, so it is "matricized" code, but I am not an expert in this area.
Tensors are really different mathematical objects with a far more rich structure than those used in the Deep Learning context.
https://en.wikipedia.org/wiki/Vector_(mathematics_and_physic...
In linear algebra you should learn this as soon as you hear about it. Is the english Wikipedia page correct in that you always use nabla for the gradient? https://en.wikipedia.org/wiki/Gradient#Generalizations
Because `grad f = ∇f` holds in scalar fields f only (not in vector or tensor fields (tensors that aren't scalars)).
In the ML context it generally doesn't matter because you don't do these transformations (and more generally there is no metric for your space). So you can just treat a gradient as an array of numbers. But in a physical context the distinction starts to matter, at least on a manifold with curvature.
The one place it does matter in ML (that I'm aware of) is information geometry, where you try to do a transformation from the normal coordinate space to a more "natural" coordinate system that is based on the Fisher information, and this introduces curvature into model.
Would indeed be cool to have true tensor processing units :-)
PKD level of cool as you could argue those chips directly process spacetime chunks
Just calling it an array might be underselling it, at this point. Perhaps a tensor in the context of CompSci is just a particular type of array with certain expected properties.
I like the NDArray terminology used by numpy. I think it gets the point across more clearly than “tensor”, and conceptually they’re pretty much equivalent
This is overly simplified though. Things get different when we start talking about {S,M}I{S,M}D (page addresses SIMD) or GPUs. Parallel computing is a whole other game, and CUDA takes it to the next level (lots of people can write CUDA kernels, not a lot of people can write GOOD CUDA kernels (I'm definitely not a pro, but know some)). Parallelism adds another level of complexity due to being able to access different memory locations simultaneously. If you think concurrent optimization is magic, parallel optimization is wizardry (I still believe this is true)
A good basic thing to know is a tensor is a single dimensional array under the hood. At least as far as I know in my mental model - I have not yet looked at the code :-)
Therefore "converting" a 2 * 5 * 10 tensor into a 2 * 50 one is an immediate operation, it doesn't have to scan. Just change some metadata.
A 2-dimensional array A will, in the object program, be stored sequentially in
the order Ai,i, A2,i Am l , A| i2 , A2f2 Am, 2> , Am,„. Thus
it is stored “columnwise”, with the first of its subscripts varying most
rapidly, and the last varying least rapidly. The same is true of 3-dimensional arrays.
1 -dimensional arrays are of course simply stored sequentially. All arrays are stored backwards in storage; i.e. the above sequence is in the order of decreasing absolute location.
I can understand almost everything in the explanations of how ChatGPT, neural networks and fine-tuning works
But one thing isn’t being explained…
How do you embed the words as vectors????
Do you just make up your own scheme?
I understand about vector databases being used for search, to retrieve snippets and stuff them into a prompt. I don’t mean that.
I mean how did people embed the words WHEN THEY USE THE EMBEDDINGS API?
I think this whole vector database and pinecone stuff will be phased out once the LLM windows get larger, like with the new Claude from Anthropic. And the ability to fine tune the model on your data. Then the LLM can make its own lower-rank model of your entire corpus of data, rather than having to do vector lookups and stuffing it into a prompt.
it is, hence the downvotes. This one is more about CPU operations, AVX/SIMD
https://jvns.ca/blog/2023/02/08/why-does-0-1-plus-0-2-equal-...
And there's nothing in vectorization that requires floating-point numbers; it's just that the algorithms that tend to demand the most out of a computer tend to be large numerical algorithms, which use floating-point.
But 0.999... is 1.
Edit: Nevermind, I see you meant "writing down 1/3 as 0.333" literally.
one person's "wrong" is another person's approximation. by reducing accuracy we can often get very large performance improvements. good engineering makes an appropriate tradeoff.
there's even more scope for this kind of thing inside some numerical algorithms where having a "pretty good" guess gives a massive speed boost, but the algorithm can be iterated to find quite precise answers even if the "pretty good" guess is an approximation that's off by 20%
there's yet more scope for this if we consider if we're trying to mathematically model reality -- even if we perform all the math with arbitrary precision calculations, the theoretical model itself is still an approximation of reality, there's some error there. there's limited value in solving a model exactly, as there's always some approximation error between the model and reality.
it's amazing we have such good support for high-performance scientific quality ieee floating point math everywhere in commodity hardware. it's a great tool to have in the toolbelt